brick_mesh.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#include "map_matrix.h"
27#include "brick_mesh.h"
28
29namespace oomph
30{
31 //====================================================================
32 /// Helper namespace for generation of brick from tet mesh
33 //====================================================================
34 namespace BrickFromTetMeshHelper
35 {
36 /// Tolerance for mismatch during setup of boundary coordinates
37 double Face_position_tolerance = 1.0e-12;
38
39 } // namespace BrickFromTetMeshHelper
40
41
42 /// ////////////////////////////////////////////////////////////////////////
43 /// ////////////////////////////////////////////////////////////////////////
44 /// ////////////////////////////////////////////////////////////////////////
45
46
47 //================================================================
48 /// Setup lookup schemes which establish which elements are located
49 /// next to which boundaries (Doc to outfile if it's open).
50 //================================================================
52 {
53 // Martina's fixed and commented version of the code.
54
55 // Define variable doc and set the initial value to False.
56 bool doc = false;
57
58 // If file declared in outfile exists,set doc=true to enable writing to
59 // stream
60 if (outfile) doc = true;
61
62 // Number of boundaries. Gives the value assigned by the function
63 // nboundary()
64 unsigned nbound = nboundary();
65
66 if (doc)
67 {
68 outfile << "The number of boundaries is " << nbound << "\n";
69 }
70
71 // Wipe/allocate storage for arrays
72 // Command clear will reomve all elements of vectors
73 // Command resize will adapt pointer to correct boundary size. Internal data
74 Boundary_element_pt.clear();
76 Boundary_element_pt.resize(nbound);
77 Face_index_at_boundary.resize(nbound);
78
79
80 // Temporary vector of vectors of pointers to elements on the boundaries:
81 // vector_of_boundary_element_pt[i] is the vector of pointers to all
82 // elements that have nodes on boundary i.
83 // This is not a set to ensure UNIQUE ordering on every processor
84 // There are a number of elements with nodes on each boundary.
85 // For each boundary we store the appropriate elements in a vector.
86 // There is therefore a vector for each boundary, which we store in
87 // another vector, i.e. a matrix (vector of vectors).
88 // Then vector_of_boundary_element_pt[0] is a vector of pointers to elements
89 // with nodes on boundary 0 and vector_of_boundary_element_pt[0][0] is the
90 // first element with nodes on boundary 0.
91 Vector<Vector<FiniteElement*>> vector_of_boundary_element_pt;
92 vector_of_boundary_element_pt.resize(nbound);
93
94 // Matrix map for working out the fixed local coord for elements on boundary
95 // Calling pointer to FiniteElement and pointer to Vector<int> defining
96 // the matrix to identify each element on boundary
98
99
100 // Temporary container to store pointers to temporary vectors
101 // so they can be deleted. Creating storage to store these
102 // temporary vectors of vectors of pointers previously defined
103 Vector<Vector<int>*> tmp_vect_pt;
104
105 // Loop over elements
106 //-------------------
107 unsigned nel = nelement();
108 for (unsigned e = 0; e < nel; e++)
109 {
110 // Get pointer to element
111 // and put it in local storage fe_pt.
113
114 // print out values of all elements to doc
115 if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
116
117 // Loop over the element's nodes and find out which boundaries they're on
118 // ----------------------------------------------------------------------
119 // Return number of nodes along one edge of the current element defined by
120 // fe_pt hence, loop over nodes of each element in 3D-dimension
121 unsigned nnode_1d = fe_pt->nnode_1d();
122
123 // Loop over nodes in order
124 for (unsigned i0 = 0; i0 < nnode_1d; i0++)
125 {
126 for (unsigned i1 = 0; i1 < nnode_1d; i1++)
127 {
128 for (unsigned i2 = 0; i2 < nnode_1d; i2++)
129 {
130 // Local node number, which is an assumed ordering defined in our
131 // underlying finite element scheme.
132 unsigned j = i0 + i1 * nnode_1d + i2 * nnode_1d * nnode_1d;
133
134 // Get pointer to the set of boundaries that this node lives on
135 // for each node reset boundaries_pt to 0,
136 // create pointer to each node in each element
137 // and give boundaries_pt a value
138 std::set<unsigned>* boundaries_pt = 0;
139 fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
140
141 // If the node lives on some boundaries....
142 // If not equal to 0, node is on boundary
143 // Loop through values of the current node stored in boundaries_pt
144 if (boundaries_pt != 0)
145 {
146 // Loop over the entries in the set "it" (name we use to refer to
147 // the iterator)
148 for (std::set<unsigned>::iterator it = boundaries_pt->begin();
149 it != boundaries_pt->end();
150 ++it)
151 {
152 // What's the boundary?
153 // Define boundary_id to have the value pointed to by
154 // boundaries_pt
155 unsigned boundary_id = *it;
156
157 // Add pointer to finite element to vector for the appropriate
158 // boundary
159
160 // Does the pointer already exist in the vector
162 std::find(vector_of_boundary_element_pt[*it].begin(),
163 vector_of_boundary_element_pt[*it].end(),
164 fe_pt);
165
166 // Only insert if we have not found it (i.e. got to the end)
167 if (b_el_it == vector_of_boundary_element_pt[*it].end())
168 {
169 vector_of_boundary_element_pt[*it].push_back(fe_pt);
170 }
171
172 // For the current element/boundary combination, create
173 // a vector that stores an indicator which element boundaries
174 // the node is located (boundary_identifier=-/+1 for nodes
175 // on the left/right boundary; boundary_identifier=-/+2 for
176 // nodes on the lower/upper boundary. We determine these indices
177 // for all corner nodes of the element and add them to a vector
178 // to a vector. This allows us to decide which faces of the
179 // element coincide with the boundary since each face of the
180 // element must have all four corner nodes on the boundary.
181 if (boundary_identifier(boundary_id, fe_pt) == 0)
182 {
183 // Here we make our vector of integers and keep track of the
184 // pointer to it The vector stores information about local
185 // node position for each boundary element
186 Vector<int>* tmp_pt = new Vector<int>;
187
188 // Add to the local scheme that will allow us to delete it
189 // later at the very end of the function
190 tmp_vect_pt.push_back(tmp_pt);
191
192 // Add the pointer to the storage scheme defined by
193 // boundary_id and pointer to finite element
194 boundary_identifier(boundary_id, fe_pt) = tmp_pt;
195 }
196
197 // Are we at a corner node? (local for each boundary element)
198 // i0,i1,i2 represent the current 3D coordinates- which local
199 // corner node.
200 if (((i0 == 0) || (i0 == nnode_1d - 1)) &&
201 ((i1 == 0) || (i1 == nnode_1d - 1)) &&
202 ((i2 == 0) || (i2 == nnode_1d - 1)))
203 {
204 // Create index to represent position relative to s_0
205 // s_0,s_1,s_2 are local 3D directions
206 (*boundary_identifier(boundary_id, fe_pt))
207 .push_back(1 * (2 * i0 / (nnode_1d - 1) - 1));
208
209 // Create index to represent position relative to s_1
210 (*boundary_identifier(boundary_id, fe_pt))
211 .push_back(2 * (2 * i1 / (nnode_1d - 1) - 1));
212
213 // Create index to represent position relative to s_2
214 (*boundary_identifier(boundary_id, fe_pt))
215 .push_back(3 * (2 * i2 / (nnode_1d - 1) - 1));
216 }
217 }
218 }
219 }
220 }
221 }
222 }
223
224
225 // Now copy everything across into permanent arrays
226 //-------------------------------------------------
227
228 // Loop over boundaries
229 //---------------------
230 for (unsigned i = 0; i < nbound; i++)
231 {
232 // Loop over elements on given boundary
234 for (IT it = vector_of_boundary_element_pt[i].begin();
235 it != vector_of_boundary_element_pt[i].end();
236 it++)
237 {
238 // Recover pointer to element for each element
239 FiniteElement* fe_pt = (*it);
240
241 // Initialise count for boundary identities (-3,-2,-1,1,2,3)
242 std::map<int, unsigned> count;
243
244 // Loop over coordinates in 3D dimension
245 for (int ii = 0; ii < 3; ii++)
246 {
247 // Loop over upper/lower end of coordinates
248 // count -/+ for each direction separately
249 for (int sign = -1; sign < 3; sign += 2)
250 {
251 count[(ii + 1) * sign] = 0;
252
253 // Initialise map of counts to 0 before loop
254 // count boundary indicator for each element
255 // and its nodes
256 }
257 }
258
259 // Loop over boundary indicators for this element/boundary
260 // count nodes for any one fixed boundary determined by
261 // boundary indicator
262 unsigned n_indicators = (*boundary_identifier(i, fe_pt)).size();
263 for (unsigned j = 0; j < n_indicators; j++)
264 {
265 count[(*boundary_identifier(i, fe_pt))[j]]++;
266 }
267
268 // Determine the correct boundary indicator by checking that it
269 // occurs four times (since four corner nodes of the element's boundary
270 // need to be located on the boundary domain)
271 int indicator = -10;
272
273
274 // Loop over coordinates
275 for (int ii = 0; ii < 3; ii++)
276 {
277 // Loop over upper/lower end of coordinates
278 for (int sign = -1; sign < 3; sign += 2)
279 {
280 // If an index occurs four times then a face is on the boundary
281 // But we can have multiple faces on the same boundary, so add all
282 // the ones that we find!
283 if (count[(ii + 1) * sign] == 4)
284 {
285 indicator = (ii + 1) * sign;
286
287 // Copy into the member data structures
288 Boundary_element_pt[i].push_back(*it);
289 Face_index_at_boundary[i].push_back(indicator);
290 }
291 }
292 }
293 }
294 }
295
296 // Doc?
297 //-----
298 if (doc)
299 {
300 // Loop over boundaries
301 for (unsigned i = 0; i < nbound; i++)
302 {
303 unsigned nel = Boundary_element_pt[i].size();
304 outfile << "Boundary: " << i << " is adjacent to " << nel << " elements"
305 << std::endl;
306
307 // Loop over elements on given boundary
308 for (unsigned e = 0; e < nel; e++)
309 {
311 outfile << "Boundary element:" << fe_pt
312 << " Face index along boundary is "
313 << Face_index_at_boundary[i][e] << std::endl;
314 }
315 }
316 }
317
318
319 // Lookup scheme has now been setup yet
321
322
323 // Cleanup temporary vectors
324 unsigned n = tmp_vect_pt.size();
325 for (unsigned i = 0; i < n; i++)
326 {
327 delete tmp_vect_pt[i];
328 }
329
330 return;
331
332
333 // Doc?
334 /*bool doc=false;
335 if (outfile) doc=true;
336
337 // Number of boundaries
338 unsigned nbound=nboundary();
339
340 // Wipe/allocate storage for arrays
341 Boundary_element_pt.clear();
342 Face_index_at_boundary.clear();
343 Boundary_element_pt.resize(nbound);
344 Face_index_at_boundary.resize(nbound);
345
346
347 // Temporary vector of vectors of pointers to elements on the boundaries:
348 // vector_of_boundary_element_pt[i] is the vector of pointers to all
349 // elements that have nodes on boundary i.
350 // This is not a set to ensure UNIQUE ordering on every processor
351 Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
352 vector_of_boundary_element_pt.resize(nbound);
353
354 // Matrix map for working out the fixed local coord for elements on boundary
355 MapMatrixMixed<unsigned,FiniteElement*,Vector<int>* >
356 boundary_identifier;
357
358 // Tmp container to store pointers to tmp vectors so they can be deleted
359 Vector<Vector<int>*> tmp_vect_pt;
360
361 // Loop over elements
362 //-------------------
363 unsigned nel=nelement();
364 for (unsigned e=0;e<nel;e++)
365 {
366 // Get pointer to element
367 FiniteElement* fe_pt=finite_element_pt(e);
368
369 if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
370
371 // Loop over the element's nodes and find out which boundaries they're on
372 // ----------------------------------------------------------------------
373 unsigned nnode_1d=fe_pt->nnode_1d();
374
375 // Loop over nodes in order
376 for (unsigned i0=0;i0<nnode_1d;i0++)
377 {
378 for (unsigned i1=0;i1<nnode_1d;i1++)
379 {
380 for (unsigned i2=0;i2<nnode_1d;i2++)
381 {
382 // Local node number
383 unsigned j=i0+i1*nnode_1d+i2*nnode_1d*nnode_1d;
384
385 // Get pointer to set of boundaries that this
386 // node lives on
387 std::set<unsigned>* boundaries_pt=0;
388 fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
389
390 // If the node lives on some boundaries....
391 if (boundaries_pt!=0)
392 {
393 for (std::set<unsigned>::iterator it=boundaries_pt->begin();
394 it!=boundaries_pt->end();++it)
395 {
396
397 // What's the boundary?
398 unsigned boundary_id=*it;
399
400 // Add pointer to finite element to vector for the appropriate
401 // boundary
402
403 // Does the pointer already exits in the vector
404 Vector<FiniteElement*>::iterator b_el_it =
405 std::find(vector_of_boundary_element_pt[*it].begin(),
406 vector_of_boundary_element_pt[*it].end(),
407 fe_pt);
408
409 //Only insert if we have not found it (i.e. got to the end)
410 if(b_el_it == vector_of_boundary_element_pt[*it].end())
411 {
412 vector_of_boundary_element_pt[*it].push_back(fe_pt);
413 }
414
415 // For the current element/boundary combination, create
416 // a vector that stores an indicator which element boundaries
417 // the node is located (boundary_identifier=-/+1 for nodes
418 // on the left/right boundary; boundary_identifier=-/+2 for
419 nodes
420 // on the lower/upper boundary. We determine these indices
421 // for all corner nodes of the element and add them to a vector
422 // to a vector. This allows us to decide which face of the
423 element
424 // coincides with the boundary since the (brick!) element must
425 // have exactly four corner nodes on the boundary.
426 if (boundary_identifier(boundary_id,fe_pt)==0)
427 {
428 Vector<int>* tmp_pt=new Vector<int>;
429 tmp_vect_pt.push_back(tmp_pt);
430 boundary_identifier(boundary_id,fe_pt)=tmp_pt;
431 }
432
433 // Are we at a corner node?
434 if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1))
435 &&((i2==0)||(i2==nnode_1d-1)))
436 {
437 // Create index to represent position relative to s_0
438 (*boundary_identifier(boundary_id,fe_pt)).
439 push_back(1*(2*i0/(nnode_1d-1)-1));
440
441 // Create index to represent position relative to s_1
442 (*boundary_identifier(boundary_id,fe_pt)).
443 push_back(2*(2*i1/(nnode_1d-1)-1));
444
445 // Create index to represent position relative to s_2
446 (*boundary_identifier(boundary_id,fe_pt)).
447 push_back(3*(2*i2/(nnode_1d-1)-1));
448 }
449 }
450 }
451 }
452 }
453 }
454 }
455
456
457 // Now copy everything across into permanent arrays
458 //-------------------------------------------------
459
460 // Loop over boundaries
461 //---------------------
462 for (unsigned i=0;i<nbound;i++)
463 {
464 // Loop over elements on given boundary
465 typedef Vector<FiniteElement*>::iterator IT;
466 for (IT it=vector_of_boundary_element_pt[i].begin();
467 it!=vector_of_boundary_element_pt[i].end();
468 it++)
469 {
470 // Recover pointer to element
471 FiniteElement* fe_pt=(*it);
472
473 // Initialise count for boundary identiers (-3,-2,-1,1,2,3)
474 std::map<int,unsigned> count;
475
476 // Loop over coordinates
477 for (int ii=0;ii<3;ii++)
478 {
479 // Loop over upper/lower end of coordinates
480 for (int sign=-1;sign<3;sign+=2)
481 {
482 count[(ii+1)*sign]=0;
483 }
484 }
485
486 // Loop over boundary indicators for this element/boundary
487 unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
488 for (unsigned j=0;j<n_indicators;j++)
489 {
490 count[(*boundary_identifier(i,fe_pt))[j] ]++;
491 }
492
493 // Determine the correct boundary indicator by checking that it
494 // occurs four times (since four corner nodes of the element's boundary
495 // need to be located on the domain boundary
496 int indicator=-10;
497
498 //Check that we're finding exactly one boundary indicator
499 bool found=false;
500
501 // Loop over coordinates
502 for (int ii=0;ii<3;ii++)
503 {
504 // Loop over upper/lower end of coordinates
505 for (int sign=-1;sign<3;sign+=2)
506 {
507 if (count[(ii+1)*sign]==4)
508 {
509 // Check that we haven't found multiple boundaries
510 if (found)
511 {
512 throw OomphLibError(
513 "Trouble: Multiple boundary identifiers!\n",
514 OOMPH_CURRENT_FUNCTION,
515 OOMPH_EXCEPTION_LOCATION);
516 }
517 found=true;
518 indicator=(ii+1)*sign;
519 }
520 }
521 }
522
523 // Check if we've found one boundary
524 if (found)
525 {
526
527 // Store element
528 Boundary_element_pt[i].push_back(*it);
529
530 // Now convert boundary indicator into information required
531 // for FaceElements
532 switch (indicator)
533 {
534 case -3:
535
536 // s_2 is fixed at -1.0:
537 Face_index_at_boundary[i].push_back(-3);
538 break;
539
540 case -2:
541
542 // s_1 is fixed at -1.0:
543 Face_index_at_boundary[i].push_back(-2);
544 break;
545
546 case -1:
547
548 // s_0 is fixed at -1.0:
549 Face_index_at_boundary[i].push_back(-1);
550 break;
551
552
553 case 1:
554
555 // s_0 is fixed at 1.0:
556 Face_index_at_boundary[i].push_back(1);
557 break;
558
559 case 2:
560
561 // s_1 is fixed at 1.0:
562 Face_index_at_boundary[i].push_back(2);
563 break;
564
565 case 3:
566
567 // s_2 is fixed at 1.0:
568 Face_index_at_boundary[i].push_back(3);
569 break;
570
571
572 default:
573
574 throw OomphLibError("Never get here",
575 OOMPH_CURRENT_FUNCTION,
576 OOMPH_EXCEPTION_LOCATION);
577 }
578
579 }
580 }
581 }
582
583 // Doc?
584 //-----
585 if (doc)
586 {
587 // Loop over boundaries
588 for (unsigned i=0;i<nbound;i++)
589 {
590 unsigned nel=Boundary_element_pt[i].size();
591 outfile << "Boundary: " << i
592 << " is adjacent to " << nel
593 << " elements" << std::endl;
594
595 // Loop over elements on given boundary
596 for (unsigned e=0;e<nel;e++)
597 {
598 FiniteElement* fe_pt=Boundary_element_pt[i][e];
599 outfile << "Boundary element:" << fe_pt
600 << " Face index along boundary is "
601 << Face_index_at_boundary[i][e] << std::endl;
602 }
603 }
604 }
605
606
607 // Lookup scheme has now been setup yet
608 Lookup_for_elements_next_boundary_is_setup=true;
609
610
611 // Cleanup temporary vectors
612 unsigned n=tmp_vect_pt.size();
613 for (unsigned i=0;i<n;i++)
614 {
615 delete tmp_vect_pt[i];
616 }*/
617 }
618
619
620} // namespace oomph
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
void setup_boundary_element_info()
Setup lookup schemes which establish whic elements are located next to mesh's boundaries (wrapper to ...
Definition: brick_mesh.h:195
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,...
Definition: mesh.h:91
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
Definition: mesh.h:87
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:95
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1365
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
double Face_position_tolerance
Tolerance for mismatch during setup of boundary coordinates.
Definition: brick_mesh.cc:37
//////////////////////////////////////////////////////////////////// ////////////////////////////////...