tet_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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 
31 #include <algorithm>
32 
33 #include "map_matrix.h"
34 #include "tet_mesh.h"
35 
36 
37 namespace oomph
38 {
39 
40 
41 //================================================================
42 /// Global static data that specifies the permitted
43 /// error in the setup of the boundary coordinates
44 //================================================================
46 
47 
48 //================================================================
49 /// Setup lookup schemes which establish which elements are located
50 /// next to which boundaries (Doc to outfile if it's open).
51 //================================================================
52 void TetMeshBase::setup_boundary_element_info(std::ostream &outfile)
53 {
54 
55  //Should we document the output here
56  bool doc = false;
57 
58  if (outfile) doc = true;
59 
60  // Number of boundaries
61  unsigned nbound=nboundary();
62 
63  // Wipe/allocate storage for arrays
64  Boundary_element_pt.clear();
65  Face_index_at_boundary.clear();
66  Boundary_element_pt.resize(nbound);
67  Face_index_at_boundary.resize(nbound);
68 
69  // Temporary vector of vectors of pointers to elements on the boundaries:
70  // This is a vector to ensure UNIQUE ordering in all processors
71  Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
72  vector_of_boundary_element_pt.resize(nbound);
73 
74  // Matrix map for working out the fixed face for elements on boundary
76  face_identifier;
77 
78  // Loop over elements
79  //-------------------
80  unsigned nel=nelement();
81 
82 
83  // Get pointer to vector of boundaries that the
84  // node lives on
85  Vector<std::set<unsigned>*> boundaries_pt(4,0);
86 
87  for (unsigned e=0;e<nel;e++)
88  {
89  // Get pointer to element
91 
92  if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
93 
94  // Only include 3D elements! Some meshes contain interface elements too.
95  if (fe_pt->dim()==3)
96  {
97  // Loop over the element's nodes and find out which boundaries they're on
98  // ----------------------------------------------------------------------
99  //We need only loop over the corner nodes
100  for(unsigned i=0;i<4;i++)
101  {
102  fe_pt->node_pt(i)->get_boundaries_pt(boundaries_pt[i]);
103  }
104 
105  //Find the common boundaries of each face
106  Vector<std::set<unsigned> > face(4);
107 
108  // NOTE: Face indices defined in Telements.h
109 
110  //Face 3 connnects points 0, 1 and 2
111  if(boundaries_pt[0] && boundaries_pt[1] && boundaries_pt[2])
112  {
113  std::set<unsigned> aux;
114 
115  std::set_intersection(boundaries_pt[0]->begin(),boundaries_pt[0]->end(),
116  boundaries_pt[1]->begin(),boundaries_pt[1]->end(),
117  std::insert_iterator<std::set<unsigned> >(
118  aux,aux.begin()));
119 
120  std::set_intersection(aux.begin(),aux.end(),
121  boundaries_pt[2]->begin(),boundaries_pt[2]->end(),
122  std::insert_iterator<std::set<unsigned> >(
123  face[3],face[3].begin()));
124  }
125 
126  //Face 2 connects points 0, 1 and 3
127  if(boundaries_pt[0] && boundaries_pt[1] && boundaries_pt[3])
128  {
129  std::set<unsigned> aux;
130 
131  std::set_intersection(boundaries_pt[0]->begin(),boundaries_pt[0]->end(),
132  boundaries_pt[1]->begin(),boundaries_pt[1]->end(),
133  std::insert_iterator<std::set<unsigned> >(
134  aux,aux.begin()));
135 
136  std::set_intersection(aux.begin(),aux.end(),
137  boundaries_pt[3]->begin(),boundaries_pt[3]->end(),
138  std::insert_iterator<std::set<unsigned> >(
139  face[2],face[2].begin()));
140  }
141 
142  //Face 1 connects points 0, 2 and 3
143  if(boundaries_pt[0] && boundaries_pt[2] && boundaries_pt[3])
144  {
145  std::set<unsigned> aux;
146 
147  std::set_intersection(boundaries_pt[0]->begin(),boundaries_pt[0]->end(),
148  boundaries_pt[2]->begin(),boundaries_pt[2]->end(),
149  std::insert_iterator<std::set<unsigned> >(
150  aux,aux.begin()));
151 
152  std::set_intersection(aux.begin(),aux.end(),
153  boundaries_pt[3]->begin(),boundaries_pt[3]->end(),
154  std::insert_iterator<std::set<unsigned> >(
155  face[1],face[1].begin()));
156  }
157 
158  //Face 0 connects points 1, 2 and 3
159  if(boundaries_pt[1] && boundaries_pt[2] && boundaries_pt[3])
160  {
161  std::set<unsigned> aux;
162 
163  std::set_intersection(boundaries_pt[1]->begin(),boundaries_pt[1]->end(),
164  boundaries_pt[2]->begin(),boundaries_pt[2]->end(),
165  std::insert_iterator<std::set<unsigned> >(
166  aux,aux.begin()));
167 
168  std::set_intersection(aux.begin(),aux.end(),
169  boundaries_pt[3]->begin(),boundaries_pt[3]->end(),
170  std::insert_iterator<std::set<unsigned> >(
171  face[0],face[0].begin()));
172  }
173 
174 
175 
176  //We now know whether any faces lay on the boundaries
177  for(unsigned i=0;i<4;i++)
178  {
179  //How many boundaries are there
180  unsigned count = 0;
181 
182  //The number of the boundary
183  int boundary=-1;
184 
185  //Loop over all the members of the set and add to the count
186  //and set the boundary
187  for(std::set<unsigned>::iterator it=face[i].begin();
188  it!=face[i].end();++it)
189  {
190  ++count;
191  boundary = *it;
192  }
193 
194  //If we're on more than one boundary, this is weird, so die
195  if(count > 1)
196  {
197  std::ostringstream error_stream;
198  error_stream << "Face " << i << " is on " <<
199  count << " boundaries.\n";
200  error_stream << "This is rather strange, so I'm going to die\n";
201  error_stream << "Your mesh may be too coarse or your tetgen mesh\n";
202  error_stream << "may be screwed up...\n";
203  throw OomphLibError(
204  error_stream.str(),
205  OOMPH_CURRENT_FUNCTION,
206  OOMPH_EXCEPTION_LOCATION);
207  }
208 
209  //If we have a boundary then add this to the appropriate set
210  if(boundary >= 0)
211  {
212 
213  // Does the pointer already exits in the vector
215  std::find(vector_of_boundary_element_pt[
216  static_cast<unsigned>(boundary)].begin(),
217  vector_of_boundary_element_pt[
218  static_cast<unsigned>(boundary)].end(),
219  fe_pt);
220 
221  //Only insert if we have not found it (i.e. got to the end)
222  if(b_el_it == vector_of_boundary_element_pt[
223  static_cast<unsigned>(boundary)].end())
224  {
225  vector_of_boundary_element_pt[static_cast<unsigned>(boundary)].
226  push_back(fe_pt);
227  }
228 
229  //Also set the fixed face
230  face_identifier(static_cast<unsigned>(boundary),fe_pt) = i;
231  }
232  }
233 
234  //Now we set the pointers to the boundary sets to zero
235  for(unsigned i=0;i<4;i++) {boundaries_pt[i] = 0;}
236 
237  }
238  }
239 
240  // Now copy everything across into permanent arrays
241  //-------------------------------------------------
242 
243  // Loop over boundaries
244  //---------------------
245  for (unsigned i=0;i<nbound;i++)
246  {
247  // Number of elements on this boundary (currently stored in a set)
248  unsigned nel=vector_of_boundary_element_pt[i].size();
249 
250  // Allocate storage for the coordinate identifiers
251  Face_index_at_boundary[i].resize(nel);
252 
253  unsigned e_count=0;
255  for (IT it=vector_of_boundary_element_pt[i].begin();
256  it!=vector_of_boundary_element_pt[i].end();
257  it++)
258  {
259  // Recover pointer to element
260  FiniteElement* fe_pt=*it;
261 
262  // Add to permanent storage
263  Boundary_element_pt[i].push_back(fe_pt);
264 
265  Face_index_at_boundary[i][e_count] = face_identifier(i,fe_pt);
266 
267  // Increment counter
268  e_count++;
269  }
270  }
271 
272 
273 
274  // Doc?
275  //-----
276  if (doc)
277  {
278  // Loop over boundaries
279  for (unsigned i=0;i<nbound;i++)
280  {
281  unsigned nel=Boundary_element_pt[i].size();
282  outfile << "Boundary: " << i
283  << " is adjacent to " << nel
284  << " elements" << std::endl;
285 
286  // Loop over elements on given boundary
287  for (unsigned e=0;e<nel;e++)
288  {
290  outfile << "Boundary element:" << fe_pt
291  << " Face index of boundary is "
292  << Face_index_at_boundary[i][e] << std::endl;
293  }
294  }
295  }
296 
297 
298  // Lookup scheme has now been setup yet
300 
301 }
302 
303 }
static double Tolerance_for_boundary_finding
Global static data that specifies the permitted error in the setup of the boundary coordinates...
Definition: tet_mesh.h:314
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
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
void setup_boundary_element_info()
Definition: tet_mesh.h:318
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
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 dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806