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