Public Member Functions | Protected Attributes | List of all members
oomph::TetgenScaffoldMesh Class Reference

#include <tetgen_scaffold_mesh.h>

+ Inheritance diagram for oomph::TetgenScaffoldMesh:

Public Member Functions

 TetgenScaffoldMesh ()
 Empty constructor. More...
 
 TetgenScaffoldMesh (const std::string &node_file_name, const std::string &element_file_name, const std::string &face_file_name)
 Constructor: Pass the filename of the tetrahedra file. More...
 
 TetgenScaffoldMesh (tetgenio &tetgen_data)
 Constructor using direct tetgenio object. More...
 
 ~TetgenScaffoldMesh ()
 Empty destructor. More...
 
unsigned global_node_number (const unsigned &i)
 
unsigned face_boundary (const unsigned &e, const unsigned &i) const
 
unsigned nglobal_edge ()
 Return the number of internal edges. More...
 
bool edge_boundary (const unsigned &i)
 
unsigned edge_index (const unsigned &e, const unsigned &i) const
 
unsigned nglobal_face ()
 Return the number of internal face. More...
 
unsigned face_index (const unsigned &e, const unsigned &i) const
 
double element_attribute (const unsigned &e) const
 
- Public Member Functions inherited from oomph::Mesh
 Mesh ()
 Default constructor. More...
 
 Mesh (const Vector< Mesh * > &sub_mesh_pt)
 
void merge_meshes (const Vector< Mesh * > &sub_mesh_pt)
 
virtual void setup_boundary_element_info ()
 
virtual void setup_boundary_element_info (std::ostream &outfile)
 
virtual void reset_boundary_element_info (Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
 Virtual function to perform the reset boundary elements info rutines. More...
 
template<class BULK_ELEMENT >
void doc_boundary_coordinates (const unsigned &b, std::ofstream &the_file)
 
virtual void scale_mesh (const double &factor)
 
 Mesh (const Mesh &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Mesh &)=delete
 Broken assignment operator. More...
 
virtual ~Mesh ()
 Virtual Destructor to clean up all memory. More...
 
void flush_element_and_node_storage ()
 
void flush_element_storage ()
 
void flush_node_storage ()
 
Node *& node_pt (const unsigned long &n)
 Return pointer to global node n. More...
 
Node * node_pt (const unsigned long &n) const
 Return pointer to global node n (const version) More...
 
GeneralisedElement *& element_pt (const unsigned long &e)
 Return pointer to element e. More...
 
GeneralisedElement * element_pt (const unsigned long &e) const
 Return pointer to element e (const version) More...
 
const Vector< GeneralisedElement * > & element_pt () const
 Return reference to the Vector of elements. More...
 
Vector< GeneralisedElement * > & element_pt ()
 Return reference to the Vector of elements. More...
 
FiniteElement * finite_element_pt (const unsigned &e) const
 
Node *& boundary_node_pt (const unsigned &b, const unsigned &n)
 Return pointer to node n on boundary b. More...
 
Node * boundary_node_pt (const unsigned &b, const unsigned &n) const
 Return pointer to node n on boundary b. More...
 
void set_nboundary (const unsigned &nbound)
 Set the number of boundaries in the mesh. More...
 
void remove_boundary_nodes ()
 Clear all pointers to boundary nodes. More...
 
void remove_boundary_nodes (const unsigned &b)
 
void remove_boundary_node (const unsigned &b, Node *const &node_pt)
 Remove a node from the boundary b. More...
 
void add_boundary_node (const unsigned &b, Node *const &node_pt)
 Add a (pointer to) a node to the b-th boundary. More...
 
void copy_boundary_node_data_from_nodes ()
 
bool boundary_coordinate_exists (const unsigned &i) const
 Indicate whether the i-th boundary has an intrinsic coordinate. More...
 
unsigned long nelement () const
 Return number of elements in the mesh. More...
 
unsigned long nnode () const
 Return number of nodes in the mesh. More...
 
unsigned ndof_types () const
 Return number of dof types in mesh. More...
 
unsigned elemental_dimension () const
 Return number of elemental dimension in mesh. More...
 
unsigned nodal_dimension () const
 Return number of nodal dimension in mesh. More...
 
void add_node_pt (Node *const &node_pt)
 Add a (pointer to a) node to the mesh. More...
 
void add_element_pt (GeneralisedElement *const &element_pt)
 Add a (pointer to) an element to the mesh. More...
 
virtual void node_update (const bool &update_all_solid_nodes=false)
 
virtual void reorder_nodes (const bool &use_old_ordering=true)
 
virtual void get_node_reordering (Vector< Node * > &reordering, const bool &use_old_ordering=true) const
 
template<class BULK_ELEMENT , template< class > class FACE_ELEMENT>
void build_face_mesh (const unsigned &b, Mesh *const &face_mesh_pt)
 
unsigned self_test ()
 Self-test: Check elements and nodes. Return 0 for OK. More...
 
void max_and_min_element_size (double &max_size, double &min_size)
 
double total_size ()
 
void check_inverted_elements (bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
 
void check_inverted_elements (bool &mesh_has_inverted_elements)
 
unsigned check_for_repeated_nodes (const double &epsilon=1.0e-12)
 
Vector< Node * > prune_dead_nodes ()
 
unsigned nboundary () const
 Return number of boundaries. More...
 
unsigned long nboundary_node (const unsigned &ibound) const
 Return number of nodes on a particular boundary. More...
 
FiniteElement * boundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Node * get_some_non_boundary_node () const
 
unsigned nboundary_element (const unsigned &b) const
 Return number of finite elements that are adjacent to boundary b. More...
 
int face_index_at_boundary (const unsigned &b, const unsigned &e) const
 
virtual void dump (std::ofstream &dump_file, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
void dump (const std::string &dump_file_name, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
virtual void read (std::ifstream &restart_file)
 Read solution from restart file. More...
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
void output (std::ostream &outfile)
 Output for all elements. More...
 
void output (std::ostream &outfile, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output (FILE *file_pt)
 Output for all elements (C-style output) More...
 
void output (FILE *file_pt, const unsigned &nplot)
 Output at f(n_plot) points in each element (C-style output) More...
 
void output (const std::string &output_filename)
 Output for all elements. More...
 
void output (const std::string &output_filename, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
 Output a given Vector function at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt)
 
void output_boundaries (std::ostream &outfile)
 Output the nodes on the boundaries (into separate tecplot zones) More...
 
void output_boundaries (const std::string &output_filename)
 
void assign_initial_values_impulsive ()
 Assign initial values for an impulsive start. More...
 
void shift_time_values ()
 
void calculate_predictions ()
 
void set_nodal_and_elemental_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void set_mesh_level_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void set_consistent_pinned_values_for_continuation (ContinuationStorageScheme *const &continuation_stepper_pt)
 Set consistent values for pinned data in continuation. More...
 
bool does_pointer_correspond_to_mesh_data (double *const &parameter_pt)
 Does the double pointer correspond to any mesh data. More...
 
void set_nodal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 Set the timestepper associated with the nodal data in the mesh. More...
 
void set_elemental_internal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void compute_norm (double &norm)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Returns the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
bool is_mesh_distributed () const
 Boolean to indicate if Mesh has been distributed. More...
 
OomphCommunicator * communicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Protected Attributes

unsigned Nglobal_face
 Storage for the number of global faces. More...
 
unsigned Nglobal_edge
 Storage for the number of global edges. More...
 
Vector< unsigned > Global_node
 Storage for global node numbers listed element-by-element. More...
 
std::vector< bool > Edge_boundary
 
Vector< Vector< unsigned > > Face_boundary
 
Vector< Vector< unsigned > > Edge_index
 Vector of vectors containing the global edge index of. More...
 
Vector< Vector< unsigned > > Face_index
 Vector of vectors containing the global edge index of. More...
 
Vector< double > Element_attribute
 
- Protected Attributes inherited from oomph::Mesh
Vector< Vector< Node * > > Boundary_node_pt
 
bool Lookup_for_elements_next_boundary_is_setup
 
Vector< Vector< FiniteElement * > > Boundary_element_pt
 
Vector< Vector< int > > Face_index_at_boundary
 
Vector< Node * > Node_pt
 Vector of pointers to nodes. More...
 
Vector< GeneralisedElement * > Element_pt
 Vector of pointers to generalised elements. More...
 
std::vector< bool > Boundary_coordinate_exists
 

Additional Inherited Members

- Public Types inherited from oomph::Mesh
typedef void(FiniteElement::* SteadyExactSolutionFctPt) (const Vector< double > &x, Vector< double > &soln)
 
typedef void(FiniteElement::* UnsteadyExactSolutionFctPt) (const double &time, const Vector< double > &x, Vector< double > &soln)
 
- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 

Detailed Description

Mesh that is based on input files generated by the tetrahedra mesh generator tetgen.

Constructor & Destructor Documentation

◆ TetgenScaffoldMesh() [1/3]

oomph::TetgenScaffoldMesh::TetgenScaffoldMesh ( )
inline

Empty constructor.

46 {}

◆ TetgenScaffoldMesh() [2/3]

oomph::TetgenScaffoldMesh::TetgenScaffoldMesh ( const std::string &  node_file_name,
const std::string &  element_file_name,
const std::string &  face_file_name 
)

Constructor: Pass the filename of the tetrahedra file.

Constructor: Pass the filename of the tetrahedra file The assumptions are that the nodes have been assigned boundary information which is used in the nodal construction to make sure that BoundaryNodes are constructed. Any additional boundaries are determined from the face boundary information.

42  {
43  // Process the element file
44  // --------------------------
45  std::ifstream element_file(element_file_name.c_str(), std::ios_base::in);
46 
47  // Check that the file actually opened correctly
48 #ifdef PARANOID
49  if (!element_file.is_open())
50  {
51  std::string error_msg("Failed to open element file: ");
52  error_msg += "\"" + element_file_name + "\".";
53  throw OomphLibError(
54  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
55  }
56 #endif
57 
58 
59  // Read in number of elements
60  unsigned n_element;
61  element_file >> n_element;
62 
63  // Read in number of nodes per element
64  unsigned n_local_node;
65  element_file >> n_local_node;
66 
67  // Throw an error if we have anything but linear simplices
68  if (n_local_node != 4)
69  {
70  std::ostringstream error_stream;
71  error_stream
72  << "Tetgen should only be used to generate 4-noded tetrahedra!\n"
73  << "Your tetgen input file, contains data for " << n_local_node
74  << "-noded tetrahedra" << std::endl;
75 
76  throw OomphLibError(
77  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
78  }
79 
80  // Dummy nodal attributes
81  unsigned dummy_attribute;
82 
83  // Element attributes may be used to distinguish internal regions
84  // NOTE: This stores doubles because tetgen forces us to!
85  Element_attribute.resize(n_element, 0.0);
86 
87  // Dummy storage for element numbers
88  unsigned dummy_element_number;
89 
90  // Resize storage for the global node numbers listed element-by-element
91  Global_node.resize(n_element * n_local_node);
92 
93  // Initialise (global) node counter
94  unsigned k = 0;
95  // Are there attributes?
96  unsigned attribute_flag;
97  element_file >> attribute_flag;
98 
99  // If no attributes
100  if (attribute_flag == 0)
101  {
102  // Read in global node numbers
103  for (unsigned i = 0; i < n_element; i++)
104  {
105  element_file >> dummy_element_number;
106  for (unsigned j = 0; j < n_local_node; j++)
107  {
108  element_file >> Global_node[k];
109  k++;
110  }
111  }
112  }
113  // Otherwise read in the attributes as well
114  else
115  {
116  for (unsigned i = 0; i < n_element; i++)
117  {
118  element_file >> dummy_element_number;
119  for (unsigned j = 0; j < n_local_node; j++)
120  {
121  element_file >> Global_node[k];
122  k++;
123  }
124  element_file >> Element_attribute[i];
125  }
126  }
127  element_file.close();
128 
129  // Resize the Element vector
130  Element_pt.resize(n_element);
131 
132  // Process node file
133  //--------------------
134  std::ifstream node_file(node_file_name.c_str(), std::ios_base::in);
135 
136  // Check that the file actually opened correctly
137 #ifdef PARANOID
138  if (!node_file.is_open())
139  {
140  std::string error_msg("Failed to open node file: ");
141  error_msg += "\"" + node_file_name + "\".";
142  throw OomphLibError(
143  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
144  }
145 #endif
146 
147 
148  // Read in the number of nodes
149  unsigned n_node;
150  node_file >> n_node;
151 
152  // Create a vector of boolean so as not to create the same node twice
153  std::vector<bool> done(n_node, false);
154 
155  // Resize the Node vector
156  Node_pt.resize(n_node);
157 
158  // Set the spatial dimension of the nodes
159  unsigned dimension;
160  node_file >> dimension;
161 
162 #ifdef PARANOID
163  if (dimension != 3)
164  {
165  throw OomphLibError("The dimesion of the nodes must be 3\n",
166  OOMPH_CURRENT_FUNCTION,
167  OOMPH_EXCEPTION_LOCATION);
168  }
169 #endif
170 
171  // Flag for attributes
172  node_file >> attribute_flag;
173 
174  // Flag for boundary markers
175  unsigned boundary_markers_flag;
176  node_file >> boundary_markers_flag;
177 
178  // Dummy storage for the node number
179  unsigned dummy_node_number;
180 
181  // Create storage for nodal positions and boundary markers
182  Vector<double> x_node(n_node);
183  Vector<double> y_node(n_node);
184  Vector<double> z_node(n_node);
185  Vector<unsigned> bound(n_node);
186 
187  // Check if there are attributes
188  // If not
189  if (attribute_flag == 0)
190  {
191  // Are there boundary markers
192  if (boundary_markers_flag == 1)
193  {
194  for (unsigned i = 0; i < n_node; i++)
195  {
196  node_file >> dummy_node_number;
197  node_file >> x_node[i];
198  node_file >> y_node[i];
199  node_file >> z_node[i];
200  node_file >> bound[i];
201  }
202  node_file.close();
203  }
204  // Otherwise no boundary markers
205  else
206  {
207  for (unsigned i = 0; i < n_node; i++)
208  {
209  node_file >> dummy_node_number;
210  node_file >> x_node[i];
211  node_file >> y_node[i];
212  node_file >> z_node[i];
213  bound[i] = 0;
214  }
215  node_file.close();
216  }
217  }
218  // Otherwise there are attributes
219  else
220  {
221  if (boundary_markers_flag == 1)
222  {
223  for (unsigned i = 0; i < n_node; i++)
224  {
225  node_file >> dummy_node_number;
226  node_file >> x_node[i];
227  node_file >> y_node[i];
228  node_file >> z_node[i];
229  node_file >> dummy_attribute;
230  node_file >> bound[i];
231  }
232  node_file.close();
233  }
234  else
235  {
236  for (unsigned i = 0; i < n_node; i++)
237  {
238  node_file >> dummy_node_number;
239  node_file >> x_node[i];
240  node_file >> y_node[i];
241  node_file >> z_node[i];
242  node_file >> dummy_attribute;
243  bound[i] = 0;
244  }
245  node_file.close();
246  }
247  }
248 
249  // Determine highest boundary index
250  //------------------------------------
251  unsigned n_bound = 0;
252  if (boundary_markers_flag == 1)
253  {
254  n_bound = bound[0];
255  for (unsigned i = 1; i < n_node; i++)
256  {
257  if (bound[i] > n_bound)
258  {
259  n_bound = bound[i];
260  }
261  }
262  }
263 
264  // Process face file to extract boundary faces
265  //--------------------------------------------
266 
267  // Open face file
268  std::ifstream face_file(face_file_name.c_str(), std::ios_base::in);
269 
270  // Check that the file actually opened correctly
271 #ifdef PARANOID
272  if (!face_file.is_open())
273  {
274  std::string error_msg("Failed to open face file: ");
275  error_msg += "\"" + face_file_name + "\".";
276  throw OomphLibError(
277  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
278  }
279 #endif
280 
281  // Number of faces in face file
282  unsigned n_face;
283  face_file >> n_face;
284 
285  // Boundary markers flag
286  face_file >> boundary_markers_flag;
287 
288  // Storage for the global node numbers (in the tetgen 1-based
289  // numbering scheme!) of the first, second and third node in
290  // each segment
291  Vector<unsigned> first_node(n_face);
292  Vector<unsigned> second_node(n_face);
293  Vector<unsigned> third_node(n_face);
294 
295  // Storage for the boundary marker for each face
296  Vector<unsigned> face_boundary(n_face);
297 
298  // Dummy for global face number
299  unsigned dummy_face_number;
300 
301  // Storage for the (boundary) faces associated with each node.
302  // Nodes are indexed using Tetgen's 1-based scheme, which is why
303  // there is a +1 here
304  Vector<std::set<unsigned>> node_on_faces(n_node + 1);
305 
306  // Extract information for each segment
307  for (unsigned i = 0; i < n_face; i++)
308  {
309  face_file >> dummy_face_number;
310  face_file >> first_node[i];
311  face_file >> second_node[i];
312  face_file >> third_node[i];
313  face_file >> face_boundary[i];
314  if (face_boundary[i] > n_bound)
315  {
316  n_bound = face_boundary[i];
317  }
318  // Add the face index to each node
319  node_on_faces[first_node[i]].insert(i);
320  node_on_faces[second_node[i]].insert(i);
321  node_on_faces[third_node[i]].insert(i);
322  }
323  face_file.close();
324 
325  // Set number of boundaries
326  if (n_bound > 0)
327  {
328  this->set_nboundary(n_bound);
329  }
330 
331  // Create the elements
332  unsigned counter = 0;
333  for (unsigned e = 0; e < n_element; e++)
334  {
335  Element_pt[e] = new TElement<3, 2>;
336  unsigned global_node_number = Global_node[counter];
337  if (done[global_node_number - 1] == false)
338  // ... -1 because node number begins at 1 in tetgen
339  {
340  // If the node is on a boundary, construct a boundary node
341  if ((boundary_markers_flag == 1) && (bound[global_node_number - 1] > 0))
342  {
343  // Construct the boundary ndoe
344  Node_pt[global_node_number - 1] =
345  finite_element_pt(e)->construct_boundary_node(3);
346 
347  // Add to the boundary lookup scheme
348  add_boundary_node(bound[global_node_number - 1] - 1,
349  Node_pt[global_node_number - 1]);
350  }
351  // Otherwise just construct a normal node
352  else
353  {
354  Node_pt[global_node_number - 1] =
355  finite_element_pt(e)->construct_node(3);
356  }
357 
358  done[global_node_number - 1] = true;
359  Node_pt[global_node_number - 1]->x(0) = x_node[global_node_number - 1];
360  Node_pt[global_node_number - 1]->x(1) = y_node[global_node_number - 1];
361  Node_pt[global_node_number - 1]->x(2) = z_node[global_node_number - 1];
362  }
363  // Otherwise just copy the node numbr accross
364  else
365  {
366  finite_element_pt(e)->node_pt(3) = Node_pt[global_node_number - 1];
367  }
368  counter++;
369 
370  // Loop over the other nodes
371  for (unsigned j = 0; j < (n_local_node - 1); j++)
372  {
373  global_node_number = Global_node[counter];
374  if (done[global_node_number - 1] == false)
375  // ... -1 because node number begins at 1 in tetgen
376  {
377  // If we're on a boundary
378  if ((boundary_markers_flag == 1) &&
379  (bound[global_node_number - 1] > 0))
380  {
381  // Construct the boundary ndoe
382  Node_pt[global_node_number - 1] =
383  finite_element_pt(e)->construct_boundary_node(j);
384 
385  // Add to the boundary lookup scheme
386  add_boundary_node(bound[global_node_number - 1] - 1,
387  Node_pt[global_node_number - 1]);
388  }
389  else
390  {
391  Node_pt[global_node_number - 1] =
392  finite_element_pt(e)->construct_node(j);
393  }
394  done[global_node_number - 1] = true;
395  Node_pt[global_node_number - 1]->x(0) =
396  x_node[global_node_number - 1];
397  Node_pt[global_node_number - 1]->x(1) =
398  y_node[global_node_number - 1];
399  Node_pt[global_node_number - 1]->x(2) =
400  z_node[global_node_number - 1];
401  }
402  // Otherwise copy the pointer over
403  else
404  {
405  finite_element_pt(e)->node_pt(j) = Node_pt[global_node_number - 1];
406  }
407  counter++;
408  }
409  }
410 
411 
412  // Resize the "matrix" that stores the boundary id for each
413  // face in each element.
414  Face_boundary.resize(n_element);
415  Face_index.resize(n_element);
416  Edge_index.resize(n_element);
417 
418 
419  // 0-based index scheme used to construct a global lookup for
420  // each face that will be used to uniquely construct interior
421  // face nodes.
422  // The faces that lie on boundaries will have already been allocated so
423  // we can start from there
424  Nglobal_face = n_face;
425 
426  // Storage for the edges associated with each node.
427  // Nodes are indexed using Tetgen's 1-based scheme, which is why
428  // there is a +1 here
429  Vector<std::set<unsigned>> node_on_edges(n_node + 1);
430 
431  // 0-based index scheme used to construct a global lookup for each
432  // edge that will be used to uniquely construct interior edge nodes
433  Nglobal_edge = 0;
434 
435  // Conversion from the edge numbers to the nodes at the end
436  // of each each edge
437  const unsigned first_local_edge_node[6] = {0, 0, 0, 1, 2, 1};
438  const unsigned second_local_edge_node[6] = {1, 2, 3, 2, 3, 3};
439 
440  // Loop over the elements
441  for (unsigned e = 0; e < n_element; e++)
442  {
443  // Each element has four faces
444  Face_boundary[e].resize(4);
445  Face_index[e].resize(4);
446  // By default each face is NOT on a boundary
447  for (unsigned i = 0; i < 4; i++)
448  {
449  Face_boundary[e][i] = 0;
450  }
451 
452  Edge_index[e].resize(6);
453 
454  // Read out the global node numbers of the nodes from
455  // tetgen's 1-based numbering.
456  // The offset is to match the offset used above
457  const unsigned element_offset = e * n_local_node;
458  unsigned glob_num[4] = {0, 0, 0, 0};
459  for (unsigned i = 0; i < 4; ++i)
460  {
461  glob_num[i] = Global_node[element_offset + ((i + 1) % 4)];
462  }
463 
464  // Now we know the global node numbers of the elements' four nodes
465  // in tetgen's 1-based numbering.
466 
467  // Determine whether any of the element faces have already been allocated
468  // an index. This may be because they are located on boundaries, or have
469  // already been visited The global index for the i-th face will be stored
470  // in Face_index[i]
471 
472  // Loop over the local faces in the element
473  for (unsigned i = 0; i < 4; ++i)
474  {
475  // On the i-th face, our numbering convention is such that
476  // it is the (3-i)th node of the element that is omitted
477  const unsigned omitted_node = 3 - i;
478 
479  // Start with the set of global faces associated with the i-th node
480  // which is always on the i-th face
481  std::set<unsigned> input = node_on_faces[glob_num[i]];
482 
483  // If there is no data yet, not point doing intersections
484  // the face has not been set up
485  if (input.size() > 0)
486  {
487  // Loop over the other nodes on the face
488  unsigned local_node = i;
489  for (unsigned i2 = 0; i2 < 3; ++i2)
490  {
491  // Create the next node index (mod 4)
492  local_node = (local_node + 1) % 4;
493  // Only take the intersection of nodes on the face
494  // i.e. not 3-i
495  if (local_node != omitted_node)
496  {
497  // Local storage for the intersection
498  std::set<unsigned> intersection_result;
499  // Calculate the intersection of the current input
500  // and next node
501  std::set_intersection(
502  input.begin(),
503  input.end(),
504  node_on_faces[glob_num[local_node]].begin(),
505  node_on_faces[glob_num[local_node]].end(),
506  std::insert_iterator<std::set<unsigned>>(
507  intersection_result, intersection_result.begin()));
508  // Let the result be the next input
509  input = intersection_result;
510  }
511  }
512  }
513 
514  // If the nodes share more than one global face index, then
515  // we have a problem
516  if (input.size() > 1)
517  {
518  throw OomphLibError(
519  "Nodes in scaffold mesh share more than one global face",
520  OOMPH_CURRENT_FUNCTION,
521  OOMPH_EXCEPTION_LOCATION);
522  }
523 
524  // If the element's face is not already allocated, the intersection
525  // will be empty
526  if (input.size() == 0)
527  {
528  // Allocate the next global index
529  Face_index[e][i] = Nglobal_face;
530  // Associate the new face index with the nodes
531  for (unsigned i2 = 0; i2 < 4; ++i2)
532  {
533  // Don't add the omitted node
534  if (i2 != omitted_node)
535  {
536  node_on_faces[glob_num[i2]].insert(Nglobal_face);
537  }
538  }
539  // Increment the global face index
540  ++Nglobal_face;
541  }
542  // Otherwise we already have a face
543  else if (input.size() == 1)
544  {
545  const unsigned global_face_index = *input.begin();
546  // Set the face index
547  Face_index[e][i] = global_face_index;
548  // Allocate the boundary index, if it's a boundary
549  if (global_face_index < n_face)
550  {
551  Face_boundary[e][i] = face_boundary[global_face_index];
552  // Add the nodes to the boundary look-up scheme in
553  // oomph-lib (0-based) index
554  for (unsigned i2 = 0; i2 < 4; ++i2)
555  {
556  // Don't add the omitted node
557  if (i2 != omitted_node)
558  {
559  add_boundary_node(face_boundary[global_face_index] - 1,
560  Node_pt[glob_num[i2] - 1]);
561  }
562  }
563  }
564  }
565  } // end of loop over the faces
566 
567 
568  // Loop over the element edges and assign global edge numbers
569  for (unsigned i = 0; i < 6; ++i)
570  {
571  // Storage for the result of the intersection
572  std::vector<unsigned> local_edge_index;
573 
574  const unsigned first_global_num = glob_num[first_local_edge_node[i]];
575  const unsigned second_global_num = glob_num[second_local_edge_node[i]];
576 
577  // Find the common global edge index for the nodes on
578  // the i-th element edge
579  std::set_intersection(node_on_edges[first_global_num].begin(),
580  node_on_edges[first_global_num].end(),
581  node_on_edges[second_global_num].begin(),
582  node_on_edges[second_global_num].end(),
583  std::insert_iterator<std::vector<unsigned>>(
584  local_edge_index, local_edge_index.begin()));
585 
586  // If the nodes share more than one global edge index, then
587  // we have a problem
588  if (local_edge_index.size() > 1)
589  {
590  throw OomphLibError(
591  "Nodes in scaffold mesh share more than one global edge",
592  OOMPH_CURRENT_FUNCTION,
593  OOMPH_EXCEPTION_LOCATION);
594  }
595 
596  // If the element's edge is not already allocated, the intersection
597  // will be empty
598  if (local_edge_index.size() == 0)
599  {
600  // Allocate the next global index
601  Edge_index[e][i] = Nglobal_edge;
602  // Associate the new edge index with the nodes
603  node_on_edges[first_global_num].insert(Nglobal_edge);
604  node_on_edges[second_global_num].insert(Nglobal_edge);
605  // Increment the global edge index
606  ++Nglobal_edge;
607  }
608  // Otherwise we already have an edge
609  else if (local_edge_index.size() == 1)
610  {
611  // Set the edge index from the result of the intersection
612  Edge_index[e][i] = local_edge_index[0];
613  }
614  }
615 
616  } // end for e
617 
618 
619  // Now determine whether any edges lie on boundaries by using the
620  // face boundary scheme and
621 
622  // Resize the storage
623  Edge_boundary.resize(Nglobal_edge, false);
624 
625  // Now loop over all the boundary faces and mark that all edges
626  // must also lie on the bounadry
627  for (unsigned i = 0; i < n_face; ++i)
628  {
629  {
630  // Storage for the result of the intersection
631  std::vector<unsigned> local_edge_index;
632 
633  // Find the common global edge index for first edge of the face
634  std::set_intersection(node_on_edges[first_node[i]].begin(),
635  node_on_edges[first_node[i]].end(),
636  node_on_edges[second_node[i]].begin(),
637  node_on_edges[second_node[i]].end(),
638  std::insert_iterator<std::vector<unsigned>>(
639  local_edge_index, local_edge_index.begin()));
640 
641  // If the nodes do not share exactly one global edge index, then
642  // we have a problem
643  if (local_edge_index.size() != 1)
644  {
645  throw OomphLibError(
646  "Nodes in scaffold mesh face do not share exactly one global edge",
647  OOMPH_CURRENT_FUNCTION,
648  OOMPH_EXCEPTION_LOCATION);
649  }
650  else
651  {
652  Edge_boundary[local_edge_index[0]] = true;
653  }
654  }
655 
656  {
657  // Storage for the result of the intersection
658  std::vector<unsigned> local_edge_index;
659 
660  // Find the common global edge index for second edge of the face
661  std::set_intersection(node_on_edges[second_node[i]].begin(),
662  node_on_edges[second_node[i]].end(),
663  node_on_edges[third_node[i]].begin(),
664  node_on_edges[third_node[i]].end(),
665  std::insert_iterator<std::vector<unsigned>>(
666  local_edge_index, local_edge_index.begin()));
667 
668  // If the nodes do not share exactly one global edge index, then
669  // we have a problem
670  if (local_edge_index.size() != 1)
671  {
672  throw OomphLibError(
673  "Nodes in scaffold mesh face do not share exactly one global edge",
674  OOMPH_CURRENT_FUNCTION,
675  OOMPH_EXCEPTION_LOCATION);
676  }
677  else
678  {
679  Edge_boundary[local_edge_index[0]] = true;
680  }
681  }
682 
683  {
684  // Storage for the result of the intersection
685  std::vector<unsigned> local_edge_index;
686 
687  // Find the common global edge index for third edge of the face
688  std::set_intersection(node_on_edges[first_node[i]].begin(),
689  node_on_edges[first_node[i]].end(),
690  node_on_edges[third_node[i]].begin(),
691  node_on_edges[third_node[i]].end(),
692  std::insert_iterator<std::vector<unsigned>>(
693  local_edge_index, local_edge_index.begin()));
694 
695  // If the nodes do not share exactly one global edge index, then
696  // we have a problem
697  if (local_edge_index.size() != 1)
698  {
699  throw OomphLibError(
700  "Nodes in scaffold mesh face do not share exactly one global edge",
701  OOMPH_CURRENT_FUNCTION,
702  OOMPH_EXCEPTION_LOCATION);
703  }
704  else
705  {
706  Edge_boundary[local_edge_index[0]] = true;
707  }
708  }
709  }
710 
711 
712  } // end of constructor
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::FiniteElement::construct_node
virtual Node * construct_node(const unsigned &n)
Definition: elements.h:2509
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::construct_boundary_node
virtual Node * construct_boundary_node(const unsigned &n)
Definition: elements.h:2538
oomph::Mesh::add_boundary_node
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::set_nboundary
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
oomph::TetgenScaffoldMesh::Edge_boundary
std::vector< bool > Edge_boundary
Definition: tetgen_scaffold_mesh.h:130
oomph::TetgenScaffoldMesh::Nglobal_face
unsigned Nglobal_face
Storage for the number of global faces.
Definition: tetgen_scaffold_mesh.h:120
oomph::TetgenScaffoldMesh::face_boundary
unsigned face_boundary(const unsigned &e, const unsigned &i) const
Definition: tetgen_scaffold_mesh.h:71
oomph::TetgenScaffoldMesh::Element_attribute
Vector< double > Element_attribute
Definition: tetgen_scaffold_mesh.h:147
oomph::TetgenScaffoldMesh::Face_index
Vector< Vector< unsigned > > Face_index
Vector of vectors containing the global edge index of.
Definition: tetgen_scaffold_mesh.h:142
oomph::TetgenScaffoldMesh::global_node_number
unsigned global_node_number(const unsigned &i)
Definition: tetgen_scaffold_mesh.h:62
oomph::TetgenScaffoldMesh::Nglobal_edge
unsigned Nglobal_edge
Storage for the number of global edges.
Definition: tetgen_scaffold_mesh.h:123
oomph::TetgenScaffoldMesh::Edge_index
Vector< Vector< unsigned > > Edge_index
Vector of vectors containing the global edge index of.
Definition: tetgen_scaffold_mesh.h:138
oomph::TetgenScaffoldMesh::Global_node
Vector< unsigned > Global_node
Storage for global node numbers listed element-by-element.
Definition: tetgen_scaffold_mesh.h:126
oomph::TetgenScaffoldMesh::Face_boundary
Vector< Vector< unsigned > > Face_boundary
Definition: tetgen_scaffold_mesh.h:134
Eigen::placeholders::end
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
k
char char char int int * k
Definition: level2_impl.h:374
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::Mesh::add_boundary_node(), oomph::FiniteElement::construct_boundary_node(), oomph::FiniteElement::construct_node(), e(), Edge_boundary, Edge_index, Element_attribute, oomph::Mesh::Element_pt, Eigen::placeholders::end, face_boundary(), Face_boundary, Face_index, oomph::Mesh::finite_element_pt(), Global_node, global_node_number(), i, j, k, Nglobal_edge, Nglobal_face, oomph::FiniteElement::node_pt(), oomph::Mesh::Node_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::Mesh::set_nboundary(), and oomph::Global_string_for_annotation::string().

◆ TetgenScaffoldMesh() [3/3]

oomph::TetgenScaffoldMesh::TetgenScaffoldMesh ( tetgenio &  tetgen_data)

Constructor using direct tetgenio object.

Constructor: Pass a tetgenio data structure that represents the input data of the mesh. The assumptions are that the nodes have been assigned boundary information which is used in the nodal construction to make sure that BoundaryNodes are constructed. Any additional boundaries are determined from the face boundary information.

724  {
725  // Find the number of elements
726  unsigned n_element = static_cast<unsigned>(tetgen_data.numberoftetrahedra);
727 
728  // Read in number of nodes per element
729  unsigned n_local_node = static_cast<unsigned>(tetgen_data.numberofcorners);
730 
731  // Throw an error if we have anything but linear simplices
732  if (n_local_node != 4)
733  {
734  std::ostringstream error_stream;
735  error_stream
736  << "Tetgen should only be used to generate 4-noded tetrahedra!\n"
737  << "Your tetgen input data, contains data for " << n_local_node
738  << "-noded tetrahedra" << std::endl;
739 
740  throw OomphLibError(
741  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
742  }
743 
744  // Element attributes may be used to distinguish internal regions
745  // NOTE: This stores doubles because tetgen forces us to!
746  Element_attribute.resize(n_element, 0.0);
747 
748  // Resize storage for the global node numbers listed element-by-element
749  Global_node.resize(n_element * n_local_node);
750 
751  // Initialise (global) node counter
752  unsigned k = 0;
753  // Are there attributes?
754  unsigned attribute_flag =
755  static_cast<unsigned>(tetgen_data.numberoftetrahedronattributes);
756 
757  // Read global node numbers for all elements
758  if (attribute_flag == 0)
759  {
760  for (unsigned i = 0; i < n_element; i++)
761  {
762  for (unsigned j = 0; j < n_local_node; j++)
763  {
764  Global_node[k] =
765  static_cast<unsigned>(tetgen_data.tetrahedronlist[k]);
766  k++;
767  }
768  }
769  }
770  // Otherwise read in the attributes as well
771  else
772  {
773  for (unsigned i = 0; i < n_element; i++)
774  {
775  for (unsigned j = 0; j < n_local_node; j++)
776  {
777  Global_node[k] =
778  static_cast<unsigned>(tetgen_data.tetrahedronlist[k]);
779  k++;
780  }
781  Element_attribute[i] = tetgen_data.tetrahedronattributelist[i];
782  }
783  }
784 
785  // Resize the Element vector
786  Element_pt.resize(n_element);
787 
788  // Read in the number of nodes
789  unsigned n_node = tetgen_data.numberofpoints;
790 
791  // Create a vector of boolean so as not to create the same node twice
792  std::vector<bool> done(n_node, false);
793 
794  // Resize the Node vector
795  Node_pt.resize(n_node);
796 
797  // Flag for boundary markers
798  unsigned boundary_markers_flag = 0;
799  if (tetgen_data.pointmarkerlist != 0)
800  {
801  boundary_markers_flag = 1;
802  }
803 
804  // Create storage for nodal positions and boundary markers
805  Vector<double> x_node(n_node);
806  Vector<double> y_node(n_node);
807  Vector<double> z_node(n_node);
808  Vector<unsigned> bound(n_node);
809 
810  // We shall ignore all point attributes
811  if (boundary_markers_flag == 1)
812  {
813  for (unsigned i = 0; i < n_node; i++)
814  {
815  x_node[i] = tetgen_data.pointlist[3 * i];
816  y_node[i] = tetgen_data.pointlist[3 * i + 1];
817  z_node[i] = tetgen_data.pointlist[3 * i + 2];
818  bound[i] = static_cast<unsigned>(tetgen_data.pointmarkerlist[i]);
819  }
820  }
821  // Otherwise no boundary markers
822  else
823  {
824  for (unsigned i = 0; i < n_node; i++)
825  {
826  x_node[i] = tetgen_data.pointlist[3 * i];
827  y_node[i] = tetgen_data.pointlist[3 * i + 1];
828  z_node[i] = tetgen_data.pointlist[3 * i + 2];
829  bound[i] = 0;
830  }
831  }
832 
833 
834  // Determine highest boundary index
835  //------------------------------------
836  unsigned n_bound = 0;
837  if (boundary_markers_flag == 1)
838  {
839  n_bound = bound[0];
840  for (unsigned i = 1; i < n_node; i++)
841  {
842  if (bound[i] > n_bound)
843  {
844  n_bound = bound[i];
845  }
846  }
847  }
848 
849  // Now extract face information
850  //---------------------------------
851 
852  // Number of faces in face file
853  unsigned n_face = tetgen_data.numberoftrifaces;
854 
855  // Storage for the global node numbers (in the tetgen 1-based
856  // numbering scheme!) of the first, second and third node in
857  // each segment
858  Vector<unsigned> first_node(n_face);
859  Vector<unsigned> second_node(n_face);
860  Vector<unsigned> third_node(n_face);
861 
862  // Storage for the boundary marker for each face
863  Vector<unsigned> face_boundary(n_face);
864 
865  // Storage for the (boundary) faces associated with each node.
866  // Nodes are indexed using Tetgen's 1-based scheme, which is why
867  // there is a +1 here
868  Vector<std::set<unsigned>> node_on_faces(n_node + 1);
869 
870  // Extract information for each segment
871  for (unsigned i = 0; i < n_face; i++)
872  {
873  first_node[i] = static_cast<unsigned>(tetgen_data.trifacelist[3 * i]);
874  second_node[i] =
875  static_cast<unsigned>(tetgen_data.trifacelist[3 * i + 1]);
876  third_node[i] = static_cast<unsigned>(tetgen_data.trifacelist[3 * i + 2]);
877  face_boundary[i] =
878  static_cast<unsigned>(tetgen_data.trifacemarkerlist[i]);
879  if (face_boundary[i] > n_bound)
880  {
881  n_bound = face_boundary[i];
882  }
883  // Add the face index to each node
884  node_on_faces[first_node[i]].insert(i);
885  node_on_faces[second_node[i]].insert(i);
886  node_on_faces[third_node[i]].insert(i);
887  }
888 
889  // Extract hole center information
890  // unsigned nhole=tetgen_data.numberofholes;
891  /*if(nhole!=0)
892  {
893  Hole_centre.resize(nhole);
894 
895  // Coords counter
896  unsigned count_coords=0;
897 
898  // Loop over the holes to get centre coords
899  for(unsigned ihole=0;ihole<nhole;ihole++)
900  {
901  Hole_centre[ihole].resize(3);
902 
903  // Read the centre value
904  Hole_centre[ihole][0]=triangle_data.holelist[count_coords];
905  Hole_centre[ihole][1]=triangle_data.holelist[count_coords+1];
906  Hole_centre[ihole][2]=triangle_data.holelist[count_coords+2];
907 
908  // Increment coords counter
909  count_coords+=3;
910  }
911  }
912  else
913  {
914  Hole_centre.resize(0);
915  }*/
916 
917 
918  // Set number of boundaries
919  if (n_bound > 0)
920  {
921  this->set_nboundary(n_bound);
922  }
923 
924  // Create the elements
925  unsigned counter = 0;
926  for (unsigned e = 0; e < n_element; e++)
927  {
928  Element_pt[e] = new TElement<3, 2>;
929  unsigned global_node_number = Global_node[counter];
930  if (done[global_node_number - 1] == false)
931  // ... -1 because node number begins at 1 in tetgen
932  {
933  // If the node is on a boundary, construct a boundary node
934  if ((boundary_markers_flag == 1) && (bound[global_node_number - 1] > 0))
935  {
936  // Construct the boundary ndoe
937  Node_pt[global_node_number - 1] =
938  finite_element_pt(e)->construct_boundary_node(3);
939 
940  // Add to the boundary lookup scheme
941  add_boundary_node(bound[global_node_number - 1] - 1,
942  Node_pt[global_node_number - 1]);
943  }
944  // Otherwise just construct a normal node
945  else
946  {
947  Node_pt[global_node_number - 1] =
948  finite_element_pt(e)->construct_node(3);
949  }
950 
951  done[global_node_number - 1] = true;
952  Node_pt[global_node_number - 1]->x(0) = x_node[global_node_number - 1];
953  Node_pt[global_node_number - 1]->x(1) = y_node[global_node_number - 1];
954  Node_pt[global_node_number - 1]->x(2) = z_node[global_node_number - 1];
955  }
956  // Otherwise just copy the node numbr accross
957  else
958  {
959  finite_element_pt(e)->node_pt(3) = Node_pt[global_node_number - 1];
960  }
961  counter++;
962 
963  // Loop over the other nodes
964  for (unsigned j = 0; j < (n_local_node - 1); j++)
965  {
966  global_node_number = Global_node[counter];
967  if (done[global_node_number - 1] == false)
968  // ... -1 because node number begins at 1 in tetgen
969  {
970  // If we're on a boundary
971  if ((boundary_markers_flag == 1) &&
972  (bound[global_node_number - 1] > 0))
973  {
974  // Construct the boundary ndoe
975  Node_pt[global_node_number - 1] =
976  finite_element_pt(e)->construct_boundary_node(j);
977 
978  // Add to the boundary lookup scheme
979  add_boundary_node(bound[global_node_number - 1] - 1,
980  Node_pt[global_node_number - 1]);
981  }
982  else
983  {
984  Node_pt[global_node_number - 1] =
985  finite_element_pt(e)->construct_node(j);
986  }
987  done[global_node_number - 1] = true;
988  Node_pt[global_node_number - 1]->x(0) =
989  x_node[global_node_number - 1];
990  Node_pt[global_node_number - 1]->x(1) =
991  y_node[global_node_number - 1];
992  Node_pt[global_node_number - 1]->x(2) =
993  z_node[global_node_number - 1];
994  }
995  // Otherwise copy the pointer over
996  else
997  {
998  finite_element_pt(e)->node_pt(j) = Node_pt[global_node_number - 1];
999  }
1000  counter++;
1001  }
1002  }
1003 
1004 
1005  // Resize the "matrix" that stores the boundary id for each
1006  // face in each element.
1007  Face_boundary.resize(n_element);
1008  Face_index.resize(n_element);
1009  Edge_index.resize(n_element);
1010 
1011 
1012  // 0-based index scheme used to construct a global lookup for
1013  // each face that will be used to uniquely construct interior
1014  // face nodes.
1015  // The faces that lie on boundaries will have already been allocated so
1016  // we can start from there
1017  Nglobal_face = n_face;
1018 
1019  // Storage for the edges associated with each node.
1020  // Nodes are indexed using Tetgen's 1-based scheme, which is why
1021  // there is a +1 here
1022  Vector<std::set<unsigned>> node_on_edges(n_node + 1);
1023 
1024  // 0-based index scheme used to construct a global lookup for each
1025  // edge that will be used to uniquely construct interior edge nodes
1026  Nglobal_edge = 0;
1027 
1028  // Conversion from the edge numbers to the nodes at the end
1029  // of each each edge
1030  const unsigned first_local_edge_node[6] = {0, 0, 0, 1, 2, 1};
1031  const unsigned second_local_edge_node[6] = {1, 2, 3, 2, 3, 3};
1032 
1033  // Loop over the elements
1034  for (unsigned e = 0; e < n_element; e++)
1035  {
1036  // Each element has four faces
1037  Face_boundary[e].resize(4);
1038  Face_index[e].resize(4);
1039  // By default each face is NOT on a boundary
1040  for (unsigned i = 0; i < 4; i++)
1041  {
1042  Face_boundary[e][i] = 0;
1043  }
1044 
1045  Edge_index[e].resize(6);
1046 
1047  // Read out the global node numbers of the nodes from
1048  // tetgen's 1-based numbering.
1049  // The offset is to match the offset used above
1050  const unsigned element_offset = e * n_local_node;
1051  unsigned glob_num[4] = {0, 0, 0, 0};
1052  for (unsigned i = 0; i < 4; ++i)
1053  {
1054  glob_num[i] = Global_node[element_offset + ((i + 1) % 4)];
1055  }
1056 
1057  // Now we know the global node numbers of the elements' four nodes
1058  // in tetgen's 1-based numbering.
1059 
1060  // Determine whether any of the element faces have already been allocated
1061  // an index. This may be because they are located on boundaries, or have
1062  // already been visited The global index for the i-th face will be stored
1063  // in Face_index[i]
1064 
1065  // Loop over the local faces in the element
1066  for (unsigned i = 0; i < 4; ++i)
1067  {
1068  // On the i-th face, our numbering convention is such that
1069  // it is the (3-i)th node of the element that is omitted
1070  const unsigned omitted_node = 3 - i;
1071 
1072  // Start with the set of global faces associated with the i-th node
1073  // which is always on the i-th face
1074  std::set<unsigned> input = node_on_faces[glob_num[i]];
1075 
1076  // If there is no data yet, not point doing intersections
1077  // the face has not been set up
1078  if (input.size() > 0)
1079  {
1080  // Loop over the other nodes on the face
1081  unsigned local_node = i;
1082  for (unsigned i2 = 0; i2 < 3; ++i2)
1083  {
1084  // Create the next node index (mod 4)
1085  local_node = (local_node + 1) % 4;
1086  // Only take the intersection of nodes on the face
1087  // i.e. not 3-i
1088  if (local_node != omitted_node)
1089  {
1090  // Local storage for the intersection
1091  std::set<unsigned> intersection_result;
1092  // Calculate the intersection of the current input
1093  // and next node
1094  std::set_intersection(
1095  input.begin(),
1096  input.end(),
1097  node_on_faces[glob_num[local_node]].begin(),
1098  node_on_faces[glob_num[local_node]].end(),
1099  std::insert_iterator<std::set<unsigned>>(
1100  intersection_result, intersection_result.begin()));
1101  // Let the result be the next input
1102  input = intersection_result;
1103  }
1104  }
1105  }
1106 
1107  // If the nodes share more than one global face index, then
1108  // we have a problem
1109  if (input.size() > 1)
1110  {
1111  throw OomphLibError(
1112  "Nodes in scaffold mesh share more than one global face",
1113  OOMPH_CURRENT_FUNCTION,
1114  OOMPH_EXCEPTION_LOCATION);
1115  }
1116 
1117  // If the element's face is not already allocated, the intersection
1118  // will be empty
1119  if (input.size() == 0)
1120  {
1121  // Allocate the next global index
1122  Face_index[e][i] = Nglobal_face;
1123  // Associate the new face index with the nodes
1124  for (unsigned i2 = 0; i2 < 4; ++i2)
1125  {
1126  // Don't add the omitted node
1127  if (i2 != omitted_node)
1128  {
1129  node_on_faces[glob_num[i2]].insert(Nglobal_face);
1130  }
1131  }
1132  // Increment the global face index
1133  ++Nglobal_face;
1134  }
1135  // Otherwise we already have a face
1136  else if (input.size() == 1)
1137  {
1138  const unsigned global_face_index = *input.begin();
1139  // Set the face index
1140  Face_index[e][i] = global_face_index;
1141  // Allocate the boundary index, if it's a boundary
1142  if (global_face_index < n_face)
1143  {
1144  Face_boundary[e][i] = face_boundary[global_face_index];
1145  // Add the nodes to the boundary look-up scheme in
1146  // oomph-lib (0-based) index
1147  for (unsigned i2 = 0; i2 < 4; ++i2)
1148  {
1149  // Don't add the omitted node
1150  if (i2 != omitted_node)
1151  {
1152  add_boundary_node(face_boundary[global_face_index] - 1,
1153  Node_pt[glob_num[i2] - 1]);
1154  }
1155  }
1156  }
1157  }
1158  } // end of loop over the faces
1159 
1160 
1161  // Loop over the element edges and assign global edge numbers
1162  for (unsigned i = 0; i < 6; ++i)
1163  {
1164  // Storage for the result of the intersection
1165  std::vector<unsigned> local_edge_index;
1166 
1167  const unsigned first_global_num = glob_num[first_local_edge_node[i]];
1168  const unsigned second_global_num = glob_num[second_local_edge_node[i]];
1169 
1170  // Find the common global edge index for the nodes on
1171  // the i-th element edge
1172  std::set_intersection(node_on_edges[first_global_num].begin(),
1173  node_on_edges[first_global_num].end(),
1174  node_on_edges[second_global_num].begin(),
1175  node_on_edges[second_global_num].end(),
1176  std::insert_iterator<std::vector<unsigned>>(
1177  local_edge_index, local_edge_index.begin()));
1178 
1179  // If the nodes share more than one global edge index, then
1180  // we have a problem
1181  if (local_edge_index.size() > 1)
1182  {
1183  throw OomphLibError(
1184  "Nodes in scaffold mesh share more than one global edge",
1185  OOMPH_CURRENT_FUNCTION,
1186  OOMPH_EXCEPTION_LOCATION);
1187  }
1188 
1189  // If the element's edge is not already allocated, the intersection
1190  // will be empty
1191  if (local_edge_index.size() == 0)
1192  {
1193  // Allocate the next global index
1194  Edge_index[e][i] = Nglobal_edge;
1195  // Associate the new edge index with the nodes
1196  node_on_edges[first_global_num].insert(Nglobal_edge);
1197  node_on_edges[second_global_num].insert(Nglobal_edge);
1198  // Increment the global edge index
1199  ++Nglobal_edge;
1200  }
1201  // Otherwise we already have an edge
1202  else if (local_edge_index.size() == 1)
1203  {
1204  // Set the edge index from the result of the intersection
1205  Edge_index[e][i] = local_edge_index[0];
1206  }
1207  }
1208 
1209  } // end for e
1210 
1211 
1212  // Now determine whether any edges lie on boundaries by using the
1213  // face boundary scheme and
1214 
1215  // Resize the storage
1216  Edge_boundary.resize(Nglobal_edge, false);
1217 
1218  // Now loop over all the boundary faces and mark that all edges
1219  // must also lie on the bounadry
1220  for (unsigned i = 0; i < n_face; ++i)
1221  {
1222  {
1223  // Storage for the result of the intersection
1224  std::vector<unsigned> local_edge_index;
1225 
1226  // Find the common global edge index for first edge of the face
1227  std::set_intersection(node_on_edges[first_node[i]].begin(),
1228  node_on_edges[first_node[i]].end(),
1229  node_on_edges[second_node[i]].begin(),
1230  node_on_edges[second_node[i]].end(),
1231  std::insert_iterator<std::vector<unsigned>>(
1232  local_edge_index, local_edge_index.begin()));
1233 
1234  // If the nodes do not share exactly one global edge index, then
1235  // we have a problem
1236  if (local_edge_index.size() != 1)
1237  {
1238  throw OomphLibError(
1239  "Nodes in scaffold mesh face do not share exactly one global edge",
1240  OOMPH_CURRENT_FUNCTION,
1241  OOMPH_EXCEPTION_LOCATION);
1242  }
1243  else
1244  {
1245  Edge_boundary[local_edge_index[0]] = true;
1246  }
1247  }
1248 
1249  {
1250  // Storage for the result of the intersection
1251  std::vector<unsigned> local_edge_index;
1252 
1253  // Find the common global edge index for second edge of the face
1254  std::set_intersection(node_on_edges[second_node[i]].begin(),
1255  node_on_edges[second_node[i]].end(),
1256  node_on_edges[third_node[i]].begin(),
1257  node_on_edges[third_node[i]].end(),
1258  std::insert_iterator<std::vector<unsigned>>(
1259  local_edge_index, local_edge_index.begin()));
1260 
1261  // If the nodes do not share exactly one global edge index, then
1262  // we have a problem
1263  if (local_edge_index.size() != 1)
1264  {
1265  throw OomphLibError(
1266  "Nodes in scaffold mesh face do not share exactly one global edge",
1267  OOMPH_CURRENT_FUNCTION,
1268  OOMPH_EXCEPTION_LOCATION);
1269  }
1270  else
1271  {
1272  Edge_boundary[local_edge_index[0]] = true;
1273  }
1274  }
1275 
1276  {
1277  // Storage for the result of the intersection
1278  std::vector<unsigned> local_edge_index;
1279 
1280  // Find the common global edge index for third edge of the face
1281  std::set_intersection(node_on_edges[first_node[i]].begin(),
1282  node_on_edges[first_node[i]].end(),
1283  node_on_edges[third_node[i]].begin(),
1284  node_on_edges[third_node[i]].end(),
1285  std::insert_iterator<std::vector<unsigned>>(
1286  local_edge_index, local_edge_index.begin()));
1287 
1288  // If the nodes do not share exactly one global edge index, then
1289  // we have a problem
1290  if (local_edge_index.size() != 1)
1291  {
1292  throw OomphLibError(
1293  "Nodes in scaffold mesh face do not share exactly one global edge",
1294  OOMPH_CURRENT_FUNCTION,
1295  OOMPH_EXCEPTION_LOCATION);
1296  }
1297  else
1298  {
1299  Edge_boundary[local_edge_index[0]] = true;
1300  }
1301  }
1302  }
1303 
1304 
1305  } // end of constructor
tetgenio::trifacemarkerlist
int * trifacemarkerlist
Definition: tetgen.h:380
tetgenio::trifacelist
int * trifacelist
Definition: tetgen.h:378
tetgenio::numberofpoints
int numberofpoints
Definition: tetgen.h:306
tetgenio::pointlist
REAL * pointlist
Definition: tetgen.h:302
tetgenio::numberofcorners
int numberofcorners
Definition: tetgen.h:326
tetgenio::numberoftrifaces
int numberoftrifaces
Definition: tetgen.h:381
tetgenio::numberoftetrahedronattributes
int numberoftetrahedronattributes
Definition: tetgen.h:327
tetgenio::numberoftetrahedra
int numberoftetrahedra
Definition: tetgen.h:325
tetgenio::tetrahedronlist
int * tetrahedronlist
Definition: tetgen.h:321
tetgenio::tetrahedronattributelist
REAL * tetrahedronattributelist
Definition: tetgen.h:322
tetgenio::pointmarkerlist
int * pointmarkerlist
Definition: tetgen.h:305

References oomph::Mesh::add_boundary_node(), oomph::FiniteElement::construct_boundary_node(), oomph::FiniteElement::construct_node(), e(), Edge_boundary, Edge_index, Element_attribute, oomph::Mesh::Element_pt, Eigen::placeholders::end, face_boundary(), Face_boundary, Face_index, oomph::Mesh::finite_element_pt(), Global_node, global_node_number(), i, j, k, Nglobal_edge, Nglobal_face, oomph::FiniteElement::node_pt(), oomph::Mesh::Node_pt, tetgenio::numberofcorners, tetgenio::numberofpoints, tetgenio::numberoftetrahedra, tetgenio::numberoftetrahedronattributes, tetgenio::numberoftrifaces, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, tetgenio::pointlist, tetgenio::pointmarkerlist, oomph::Mesh::set_nboundary(), tetgenio::tetrahedronattributelist, tetgenio::tetrahedronlist, tetgenio::trifacelist, and tetgenio::trifacemarkerlist.

◆ ~TetgenScaffoldMesh()

oomph::TetgenScaffoldMesh::~TetgenScaffoldMesh ( )
inline

Empty destructor.

57 {}

Member Function Documentation

◆ edge_boundary()

bool oomph::TetgenScaffoldMesh::edge_boundary ( const unsigned &  i)
inline

Return a boolean indicating whether the i-th global edge is on a boundary

85  {
86  return Edge_boundary[i];
87  }

References Edge_boundary, and i.

◆ edge_index()

unsigned oomph::TetgenScaffoldMesh::edge_index ( const unsigned &  e,
const unsigned &  i 
) const
inline

Return the global index of the i-th edge in the e-th element: The global index starts from zero

92  {
93  return Edge_index[e][i];
94  }

References e(), Edge_index, and i.

◆ element_attribute()

double oomph::TetgenScaffoldMesh::element_attribute ( const unsigned &  e) const
inline

Return the attribute of the element e. NOTE: Attributes are doubles because tetgen forces us to! We only use them for region IDs

113  {
114  return Element_attribute[e];
115  }

References e(), and Element_attribute.

◆ face_boundary()

unsigned oomph::TetgenScaffoldMesh::face_boundary ( const unsigned &  e,
const unsigned &  i 
) const
inline

Return the boundary id of the i-th face in the e-th element: This is zero-based as in tetgen. Zero means the face is not on a boundary. Postive numbers identify the boundary. Will be reduced by one to identify the oomph-lib boundary.

72  {
73  return Face_boundary[e][i];
74  }

References e(), Face_boundary, and i.

Referenced by TetgenScaffoldMesh().

◆ face_index()

unsigned oomph::TetgenScaffoldMesh::face_index ( const unsigned &  e,
const unsigned &  i 
) const
inline

Return the global index of the i-th face in the e-th element: The global index starts from zero

105  {
106  return Face_index[e][i];
107  }

References e(), Face_index, and i.

◆ global_node_number()

unsigned oomph::TetgenScaffoldMesh::global_node_number ( const unsigned &  i)
inline

Return the global node of each local node listed element-by-element e*n_local_node + n_local Note that the node numbers are indexed from 1

63  {
64  return Global_node[i];
65  }

References Global_node, and i.

Referenced by TetgenScaffoldMesh().

◆ nglobal_edge()

unsigned oomph::TetgenScaffoldMesh::nglobal_edge ( )
inline

Return the number of internal edges.

78  {
79  return Nglobal_edge;
80  }

References Nglobal_edge.

◆ nglobal_face()

unsigned oomph::TetgenScaffoldMesh::nglobal_face ( )
inline

Return the number of internal face.

98  {
99  return Nglobal_face;
100  }

References Nglobal_face.

Member Data Documentation

◆ Edge_boundary

std::vector<bool> oomph::TetgenScaffoldMesh::Edge_boundary
protected

Vector of booleans to indicate whether a global edge lies on a boundary

Referenced by edge_boundary(), and TetgenScaffoldMesh().

◆ Edge_index

Vector<Vector<unsigned> > oomph::TetgenScaffoldMesh::Edge_index
protected

Vector of vectors containing the global edge index of.

Referenced by edge_index(), and TetgenScaffoldMesh().

◆ Element_attribute

Vector<double> oomph::TetgenScaffoldMesh::Element_attribute
protected

Vector of double attributes for each element. NOTE: This stores doubles because tetgen forces us to! We only use it for region IDs

Referenced by element_attribute(), and TetgenScaffoldMesh().

◆ Face_boundary

Vector<Vector<unsigned> > oomph::TetgenScaffoldMesh::Face_boundary
protected

Vector of vectors containing the boundary ids of the elements' faces

Referenced by face_boundary(), and TetgenScaffoldMesh().

◆ Face_index

Vector<Vector<unsigned> > oomph::TetgenScaffoldMesh::Face_index
protected

Vector of vectors containing the global edge index of.

Referenced by face_index(), and TetgenScaffoldMesh().

◆ Global_node

Vector<unsigned> oomph::TetgenScaffoldMesh::Global_node
protected

Storage for global node numbers listed element-by-element.

Referenced by global_node_number(), and TetgenScaffoldMesh().

◆ Nglobal_edge

unsigned oomph::TetgenScaffoldMesh::Nglobal_edge
protected

Storage for the number of global edges.

Referenced by nglobal_edge(), and TetgenScaffoldMesh().

◆ Nglobal_face

unsigned oomph::TetgenScaffoldMesh::Nglobal_face
protected

Storage for the number of global faces.

Referenced by nglobal_face(), and TetgenScaffoldMesh().

The documentation for this class was generated from the following files: