oomph::TetgenMesh< ELEMENT > Class Template Reference

#include <tetgen_mesh.template.h>

+ Inheritance diagram for oomph::TetgenMesh< ELEMENT >:

Public Member Functions

 TetgenMesh ()
 Empty constructor. More...
 
 TetgenMesh (const std::string &node_file_name, const std::string &element_file_name, const std::string &face_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
 Constructor with the input files. More...
 
 TetgenMesh (tetgenio &tetgen_data, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
 Constructor with tetgenio data structure. More...
 
 TetgenMesh (const std::string &node_file_name, const std::string &element_file_name, const std::string &face_file_name, const bool &split_corner_elements, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
 
 TetgenMesh (tetgenio &tetgen_data, const bool &split_corner_elements, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
 
 TetgenMesh (TetMeshFacetedClosedSurface *const &outer_boundary_pt, Vector< TetMeshFacetedSurface * > &internal_surface_pt, const double &element_volume, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false, const bool &split_corner_elements=false, Vector< double > *const &target_element_volume_in_region_pt=nullptr)
 
void build_tetgenio (TetMeshFacetedSurface *const &outer_boundary_pt, Vector< TetMeshFacetedSurface * > &internal_surface_pt, Vector< double > *const &target_element_volume_in_region_pt, tetgenio &tetgen_io)
 Build tetgenio object from the TetMeshFacetedSurfaces. More...
 
 ~TetgenMesh ()
 Empty destructor. More...
 
void set_mesh_level_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
bool tetgenio_exists () const
 Boolen defining whether tetgenio object has been built or not. More...
 
tetgenio *& tetgenio_pt ()
 Access to the triangulateio representation of the mesh. More...
 
void set_deep_copy_tetgenio_pt (tetgenio *const &tetgenio_pt)
 Set the tetgen pointer by a deep copy. More...
 
void deep_copy_of_tetgenio (tetgenio *const &input_pt, tetgenio *&output_pt)
 
- Public Member Functions inherited from oomph::TetMeshBase
 TetMeshBase ()
 Constructor. More...
 
 TetMeshBase (const TetMeshBase &node)=delete
 Broken copy constructor. More...
 
void operator= (const TetMeshBase &)=delete
 Broken assignment operator. More...
 
virtual ~TetMeshBase ()
 Destructor (empty) More...
 
void assess_mesh_quality (std::ofstream &some_file)
 
template<class ELEMENT >
void setup_boundary_coordinates (const unsigned &b)
 
template<class ELEMENT >
void setup_boundary_coordinates (const unsigned &b, const bool &switch_normal)
 
template<class ELEMENT >
void setup_boundary_coordinates (const unsigned &b, const bool &switch_normal, std::ofstream &outfile)
 
template<class ELEMENT >
void setup_boundary_coordinates (const unsigned &b, std::ofstream &outfile)
 
unsigned nboundary_element_in_region (const unsigned &b, const unsigned &r) const
 Return the number of elements adjacent to boundary b in region r. More...
 
FiniteElementboundary_element_in_region_pt (const unsigned &b, const unsigned &r, const unsigned &e) const
 Return pointer to the e-th element adjacent to boundary b in region r. More...
 
int face_index_at_boundary_in_region (const unsigned &b, const unsigned &r, const unsigned &e) const
 Return face index of the e-th element adjacent to boundary b in region r. More...
 
unsigned nregion ()
 Return the number of regions specified by attributes. More...
 
unsigned nregion_element (const unsigned &r)
 Return the number of elements in region r. More...
 
double region_attribute (const unsigned &i)
 
FiniteElementregion_element_pt (const unsigned &r, const unsigned &e)
 Return the e-th element in the r-th region. More...
 
template<class ELEMENT >
void snap_to_quadratic_surface (const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal, DocInfo &doc_info)
 
template<class ELEMENT >
void snap_to_quadratic_surface (const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal)
 
void snap_nodes_onto_geometric_objects ()
 
template<class ELEMENT >
void split_elements_in_corners (TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
void setup_boundary_element_info ()
 
void setup_boundary_element_info (std::ostream &outfile)
 
- 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 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...
 
Nodenode_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...
 
GeneralisedElementelement_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...
 
FiniteElementfinite_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...
 
Nodeboundary_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...
 
FiniteElementboundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Nodeget_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)
 
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...
 
OomphCommunicatorcommunicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Protected Member Functions

void build_from_scaffold (TimeStepper *time_stepper_pt, const bool &use_attributes)
 Build mesh from scaffold. More...
 
void setup_reverse_lookup_schemes_for_faceted_surface (TetMeshFacetedSurface *const &faceted_surface_pt)
 Function to setup the reverse look-up schemes. 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)
 

Protected Attributes

TetgenScaffoldMeshTmp_mesh_pt
 Temporary scaffold mesh. More...
 
bool Tetgenio_exists
 
tetgenioTetgenio_pt
 Tetgen representation of mesh. More...
 
bool Use_attributes
 
- Protected Attributes inherited from oomph::TetMeshBase
Vector< Vector< FiniteElement * > > Region_element_pt
 
Vector< doubleRegion_attribute
 
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
 Storage for elements adjacent to a boundary in a particular region. More...
 
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
 
TetMeshFacetedClosedSurfaceOuter_boundary_pt
 Faceted surface that defines outer boundaries. More...
 
Vector< TetMeshFacetedSurface * > Internal_surface_pt
 Vector to faceted surfaces that define internal boundaries. More...
 
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
 
std::map< unsigned, TetMeshFacet * > Tet_mesh_facet_pt
 
std::map< unsigned, Vector< Vector< double > > > Triangular_facet_vertex_boundary_coordinate
 
TimeStepperTime_stepper_pt
 Timestepper used to build nodes. More...
 
- 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< boolBoundary_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::TetMeshBase
static double Tolerance_for_boundary_finding = 1.0e-5
 
- 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...
 

Detailed Description

template<class ELEMENT>
class oomph::TetgenMesh< ELEMENT >

Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen/

Constructor & Destructor Documentation

◆ TetgenMesh() [1/6]

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::TetgenMesh ( )
inline

Empty constructor.

56  {
57  // Mesh can only be built with 3D Telements.
58  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
59  }

◆ TetgenMesh() [2/6]

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::TetgenMesh ( const std::string &  node_file_name,
const std::string &  element_file_name,
const std::string &  face_file_name,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper,
const bool use_attributes = false 
)
inline

Constructor with the input files.

67  : Tetgenio_exists(false), Tetgenio_pt(nullptr)
68  {
69  // Mesh can only be built with 3D Telements.
70  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
71 
72  // Store the attributes
73  Use_attributes = use_attributes;
74 
75  // Store timestepper used to build elements
76  Time_stepper_pt = time_stepper_pt;
77 
78  // Build scaffold
79  Tmp_mesh_pt = new TetgenScaffoldMesh(
80  node_file_name, element_file_name, face_file_name);
81 
82  // Convert mesh from scaffold to actual mesh
83  build_from_scaffold(time_stepper_pt, use_attributes);
84 
85  // Kill the scaffold
86  delete Tmp_mesh_pt;
87  Tmp_mesh_pt = 0;
88 
89  // Setup boundary coordinates
90  unsigned nb = nboundary();
91  for (unsigned b = 0; b < nb; b++)
92  {
93  bool switch_normal = false;
94  setup_boundary_coordinates<ELEMENT>(b, switch_normal);
95  }
96  }
Scalar * b
Definition: benchVecAdd.cpp:17
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
TimeStepper * Time_stepper_pt
Timestepper used to build nodes.
Definition: tet_mesh.h:1057
tetgenio * Tetgenio_pt
Tetgen representation of mesh.
Definition: tetgen_mesh.template.h:720
bool Use_attributes
Definition: tetgen_mesh.template.h:724
bool Tetgenio_exists
Definition: tetgen_mesh.template.h:717
TetgenScaffoldMesh * Tmp_mesh_pt
Temporary scaffold mesh.
Definition: tetgen_mesh.template.h:713
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
Definition: tetgen_mesh.template.cc:48
int nb
Definition: level2_impl.h:286

References b, oomph::TetgenMesh< ELEMENT >::build_from_scaffold(), nb, oomph::Mesh::nboundary(), oomph::TetMeshBase::Time_stepper_pt, oomph::TetgenMesh< ELEMENT >::Tmp_mesh_pt, and oomph::TetgenMesh< ELEMENT >::Use_attributes.

◆ TetgenMesh() [3/6]

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::TetgenMesh ( tetgenio tetgen_data,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper,
const bool use_attributes = false 
)
inline

Constructor with tetgenio data structure.

104  {
105  // Mesh can only be built with 3D Telements.
106  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
107 
108  // Store the attributes
109  Use_attributes = use_attributes;
110 
111  // Store timestepper used to build elements
112  Time_stepper_pt = time_stepper_pt;
113 
114  // We do not have a tetgenio representation
115  Tetgenio_exists = false;
116  Tetgenio_pt = 0;
117 
118  // Build scaffold
119  Tmp_mesh_pt = new TetgenScaffoldMesh(tetgen_data);
120 
121  // Convert mesh from scaffold to actual mesh
122  build_from_scaffold(time_stepper_pt, use_attributes);
123 
124  // Kill the scaffold
125  delete Tmp_mesh_pt;
126  Tmp_mesh_pt = 0;
127 
128  // Setup boundary coordinates
129  unsigned nb = nboundary();
130  for (unsigned b = 0; b < nb; b++)
131  {
132  bool switch_normal = false;
133  setup_boundary_coordinates<ELEMENT>(b, switch_normal);
134  }
135  }

References b, oomph::TetgenMesh< ELEMENT >::build_from_scaffold(), nb, oomph::Mesh::nboundary(), oomph::TetgenMesh< ELEMENT >::Tetgenio_exists, oomph::TetgenMesh< ELEMENT >::Tetgenio_pt, oomph::TetMeshBase::Time_stepper_pt, oomph::TetgenMesh< ELEMENT >::Tmp_mesh_pt, and oomph::TetgenMesh< ELEMENT >::Use_attributes.

◆ TetgenMesh() [4/6]

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::TetgenMesh ( const std::string &  node_file_name,
const std::string &  element_file_name,
const std::string &  face_file_name,
const bool split_corner_elements,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper,
const bool use_attributes = false 
)
inline

Constructor with the input files. Setting the boolean flag to true splits "corner" elements, i.e. elements that that have at least three faces on a domain boundary. The relevant elements are split without introducing hanging nodes so the sons have a "worse" shape than their fathers. However, this step avoids otherwise-hard-to-diagnose problems in fluids problems where the application of boundary conditions at such "corner" elements can overconstrain the solution.

154  {
155  // Mesh can only be built with 3D Telements.
156  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
157 
158  // Store the attributes
159  Use_attributes = use_attributes;
160 
161  // Store timestepper used to build elements
162  Time_stepper_pt = time_stepper_pt;
163 
164  // We do not have a tetgenio representation
165  this->Tetgenio_exists = false;
166  this->Tetgenio_pt = 0;
167 
168  // Build scaffold
169  Tmp_mesh_pt = new TetgenScaffoldMesh(
170  node_file_name, element_file_name, face_file_name);
171 
172  // Convert mesh from scaffold to actual mesh
173  build_from_scaffold(time_stepper_pt, use_attributes);
174 
175  // Kill the scaffold
176  delete Tmp_mesh_pt;
177  Tmp_mesh_pt = 0;
178 
179  // Split corner elements
180  if (split_corner_elements)
181  {
182  split_elements_in_corners<ELEMENT>();
183  }
184 
185  // Setup boundary coordinates
186  unsigned nb = nboundary();
187  for (unsigned b = 0; b < nb; b++)
188  {
189  bool switch_normal = false;
190  setup_boundary_coordinates<ELEMENT>(b, switch_normal);
191  }
192  }

References b, oomph::TetgenMesh< ELEMENT >::build_from_scaffold(), nb, oomph::Mesh::nboundary(), oomph::TetgenMesh< ELEMENT >::Tetgenio_exists, oomph::TetgenMesh< ELEMENT >::Tetgenio_pt, oomph::TetMeshBase::Time_stepper_pt, oomph::TetgenMesh< ELEMENT >::Tmp_mesh_pt, and oomph::TetgenMesh< ELEMENT >::Use_attributes.

◆ TetgenMesh() [5/6]

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::TetgenMesh ( tetgenio tetgen_data,
const bool split_corner_elements,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper,
const bool use_attributes = false 
)
inline

Constructor with tetgen data structure Setting the boolean flag to true splits "corner" elements, i.e. elements that that have at least three faces on a domain boundary. The relevant elements are split without introducing hanging nodes so the sons have a "worse" shape than their fathers. However, this step avoids otherwise-hard-to-diagnose problems in fluids problems where the application of boundary conditions at such "corner" elements can overconstrain the solution.

208  {
209  // Mesh can only be built with 3D Telements.
210  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
211 
212  // Store the attributes
213  Use_attributes = use_attributes;
214 
215  // Store timestepper used to build elements
216  Time_stepper_pt = time_stepper_pt;
217 
218  // We do not have a tetgenio representation
219  this->Tetgenio_exists = false;
220  this->Tetgenio_pt = nullptr;
221 
222  // Build scaffold
223  Tmp_mesh_pt = new TetgenScaffoldMesh(tetgen_data);
224 
225  // Convert mesh from scaffold to actual mesh
226  build_from_scaffold(time_stepper_pt, use_attributes);
227 
228  // Kill the scaffold
229  delete Tmp_mesh_pt;
230  Tmp_mesh_pt = nullptr;
231 
232  // Split corner elements
233  if (split_corner_elements)
234  {
235  split_elements_in_corners<ELEMENT>();
236  }
237 
238  // Setup boundary coordinates
239  unsigned nb = nboundary();
240  for (unsigned b = 0; b < nb; b++)
241  {
242  bool switch_normal = false;
243  setup_boundary_coordinates<ELEMENT>(b, switch_normal);
244  }
245  }

References b, oomph::TetgenMesh< ELEMENT >::build_from_scaffold(), nb, oomph::Mesh::nboundary(), oomph::TetgenMesh< ELEMENT >::Tetgenio_exists, oomph::TetgenMesh< ELEMENT >::Tetgenio_pt, oomph::TetMeshBase::Time_stepper_pt, oomph::TetgenMesh< ELEMENT >::Tmp_mesh_pt, and oomph::TetgenMesh< ELEMENT >::Use_attributes.

◆ TetgenMesh() [6/6]

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::TetgenMesh ( TetMeshFacetedClosedSurface *const &  outer_boundary_pt,
Vector< TetMeshFacetedSurface * > &  internal_surface_pt,
const double element_volume,
TimeStepper time_stepper_pt = &Mesh::Default_TimeStepper,
const bool use_attributes = false,
const bool split_corner_elements = false,
Vector< double > *const &  target_element_volume_in_region_pt = nullptr 
)
inline

Build mesh, based on a TetgenMeshFactedClosedSurface that specifies the outer boundary of the domain and any number of internal boundaries, specified by TetMeshFacetedSurfaces. Also specify target size for uniform element size. Optionally specify the target element volume in each region.

261  {
262  // Mesh can only be built with 3D Telements.
263  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
264 
265  // Store the attributes
266  Use_attributes = use_attributes;
267 
268  // Store timestepper used to build elements
269  Time_stepper_pt = time_stepper_pt;
270 
271  // Copy across
272  Outer_boundary_pt = outer_boundary_pt;
273  // Setup the reverse lookup scheme
275  // Store the internal boundary
276  Internal_surface_pt = internal_surface_pt;
277  // Setup the reverse lookup schemes
278  {
279  unsigned n = this->Internal_surface_pt.size();
280  for (unsigned i = 0; i < n; i++)
281  {
284  }
285  }
286 
287  // Tetgen data structure for the input and output
288  tetgenio in;
289  this->build_tetgenio(outer_boundary_pt,
290  internal_surface_pt,
291  target_element_volume_in_region_pt,
292  in);
293 
294 
295  // Now tetrahedralise
296 
297  // The 'p' switch reads a Piecewise Linear Complex, which generates a
298  // Delaunay tetrahedralization of the input. The 'q' switch prevents
299  // generation of high-aspect ratio tets (slivers), with the trailing float
300  // indicating the maximum allowed aspect ratio (default is 2.0). The 'a'
301  // switch without subsequent floating-point number switches on the
302  // specific element volume constraints for particular regions (the 5th
303  // index in the tetgenio.regionlist array). The 'A' enables region
304  // attributes, and the second 'a' switch with the subsequent float sets
305  // the global (non-region-specific) element volume constraint.
306  std::stringstream input_string;
307  input_string << "pq2.0aAa" << element_volume;
308 
309  // input_string << "Vpq1.414Aa" << element_volume;
310  // << "Q"; // Q for quiet!
311  // << "V"; // V for verbose incl. quality output!
312 
313  // If any of the boundaries should not be split add the "Y" flag
314  bool can_boundaries_be_split =
315  outer_boundary_pt->boundaries_can_be_split_in_tetgen();
316  {
317  unsigned n_internal = internal_surface_pt.size();
318  for (unsigned i = 0; i < n_internal; i++)
319  {
320  can_boundaries_be_split &=
321  internal_surface_pt[i]->boundaries_can_be_split_in_tetgen();
322  }
323  }
324 
325  // If we can't split the boundaries add the flag
326  if (can_boundaries_be_split == false)
327  {
328  input_string << "Y";
329  }
330 
331  // Now convert to a C-style string
332  char tetswitches[100];
333  sprintf(tetswitches, "%s", input_string.str().c_str());
334 
335  // Make a new tetgen representation
336  this->Tetgenio_exists = true;
337  Tetgenio_pt = new tetgenio;
338  tetrahedralize(tetswitches, &in, this->Tetgenio_pt);
339 
340  // Build scaffold
341  Tmp_mesh_pt = new TetgenScaffoldMesh(*this->Tetgenio_pt);
342 
343  // If any of the objects are different regions then we need to use
344  // the atributes
345  bool regions_exist = false;
346  {
347  unsigned n_internal = internal_surface_pt.size();
348  for (unsigned i = 0; i < n_internal; i++)
349  {
350  TetMeshFacetedClosedSurface* srf_pt =
351  dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[i]);
352  if (srf_pt != 0)
353  {
354  unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
355  for (unsigned j = 0; j < n_int_pts; j++)
356  {
357  regions_exist |=
358  srf_pt->internal_point_identifies_region_for_tetgen(j);
359  }
360  }
361  }
362  }
363 
364  // If there are regions, use the attributes
365  if (regions_exist)
366  {
367  Use_attributes = true;
368  }
369 
370  // Convert mesh from scaffold to actual mesh
371  build_from_scaffold(time_stepper_pt, Use_attributes);
372 
373  // Kill the scaffold
374  delete Tmp_mesh_pt;
375  Tmp_mesh_pt = 0;
376 
377  // Split corner elements
378  if (split_corner_elements)
379  {
380  split_elements_in_corners<ELEMENT>();
381  }
382 
383  // Setup boundary coordinates
384  unsigned nb = nboundary();
385  for (unsigned b = 0; b < nb; b++)
386  {
387  bool switch_normal = false;
388  setup_boundary_coordinates<ELEMENT>(b, switch_normal);
389  }
390 
391  // Now snap onto geometric objects associated with triangular facets
392  // (if any!)
394  }
int i
Definition: BiCGSTAB_step_by_step.cpp:9
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
void snap_nodes_onto_geometric_objects()
Definition: tet_mesh.cc:483
Vector< TetMeshFacetedSurface * > Internal_surface_pt
Vector to faceted surfaces that define internal boundaries.
Definition: tet_mesh.h:1041
TetMeshFacetedClosedSurface * Outer_boundary_pt
Faceted surface that defines outer boundaries.
Definition: tet_mesh.h:1038
void setup_reverse_lookup_schemes_for_faceted_surface(TetMeshFacetedSurface *const &faceted_surface_pt)
Function to setup the reverse look-up schemes.
Definition: tetgen_mesh.template.cc:828
void build_tetgenio(TetMeshFacetedSurface *const &outer_boundary_pt, Vector< TetMeshFacetedSurface * > &internal_surface_pt, Vector< double > *const &target_element_volume_in_region_pt, tetgenio &tetgen_io)
Build tetgenio object from the TetMeshFacetedSurfaces.
Definition: tetgen_mesh.template.h:397
Definition: tetgen.h:201
void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out, tetgenio *addin=NULL, tetgenio *bgmin=NULL)
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References b, oomph::TetMeshFacetedSurface::boundaries_can_be_split_in_tetgen(), oomph::TetgenMesh< ELEMENT >::build_from_scaffold(), oomph::TetgenMesh< ELEMENT >::build_tetgenio(), i, oomph::TetMeshFacetedClosedSurface::internal_point_identifies_region_for_tetgen(), oomph::TetMeshBase::Internal_surface_pt, j, n, nb, oomph::Mesh::nboundary(), oomph::TetMeshFacetedClosedSurface::ninternal_point_for_tetgen(), oomph::TetMeshBase::Outer_boundary_pt, oomph::TetgenMesh< ELEMENT >::setup_reverse_lookup_schemes_for_faceted_surface(), oomph::TetMeshBase::snap_nodes_onto_geometric_objects(), oomph::TetgenMesh< ELEMENT >::Tetgenio_exists, oomph::TetgenMesh< ELEMENT >::Tetgenio_pt, tetrahedralize(), oomph::TetMeshBase::Time_stepper_pt, oomph::TetgenMesh< ELEMENT >::Tmp_mesh_pt, and oomph::TetgenMesh< ELEMENT >::Use_attributes.

◆ ~TetgenMesh()

template<class ELEMENT >
oomph::TetgenMesh< ELEMENT >::~TetgenMesh ( )
inline

Empty destructor.

654  {
655  if (Tetgenio_exists)
656  {
657  delete Tetgenio_pt;
658  }
659  }

References oomph::TetgenMesh< ELEMENT >::Tetgenio_exists, and oomph::TetgenMesh< ELEMENT >::Tetgenio_pt.

Member Function Documentation

◆ build_from_scaffold()

template<class ELEMENT >
void TetgenMesh::build_from_scaffold ( TimeStepper time_stepper_pt,
const bool use_attributes 
)
protected

Build mesh from scaffold.

Build unstructured tet mesh based on output from scaffold.

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////

50  {
51  // Mesh can only be built with 3D Telements.
52  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
53 
54  // Create space for elements
55  unsigned nelem = Tmp_mesh_pt->nelement();
56  Element_pt.resize(nelem);
57 
58  // Create space for nodes
59  unsigned nnode_scaffold = Tmp_mesh_pt->nnode();
60  Node_pt.resize(nnode_scaffold);
61 
62  // Set number of boundaries
63  unsigned nbound = Tmp_mesh_pt->nboundary();
64  set_nboundary(nbound);
65  Boundary_element_pt.resize(nbound);
66  Face_index_at_boundary.resize(nbound);
67 
68  // If we have different regions, then resize the region
69  // information
70  if (use_attributes)
71  {
72  Boundary_region_element_pt.resize(nbound);
73  Face_index_region_at_boundary.resize(nbound);
74  }
75 
76  // Build elements
77  for (unsigned e = 0; e < nelem; e++)
78  {
79  Element_pt[e] = new ELEMENT;
80  }
81 
82  // Number of nodes per element
83  unsigned nnod_el = Tmp_mesh_pt->finite_element_pt(0)->nnode();
84 
85  // Setup map to check the (pseudo-)global node number
86  // Nodes whose number is zero haven't been copied across
87  // into the mesh yet.
88  std::map<Node*, unsigned> global_number;
89  unsigned global_count = 0;
90 
91  // Map of element attribute pairs
92  std::map<double, Vector<FiniteElement*>> element_attribute_map;
93 
94  // Loop over elements in scaffold mesh, visit their nodes
95  for (unsigned e = 0; e < nelem; e++)
96  {
97  // Loop over all nodes in element
98  for (unsigned j = 0; j < nnod_el; j++)
99  {
100  // Pointer to node in the scaffold mesh
101  Node* scaffold_node_pt = Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
102 
103  // Get the (pseudo-)global node number in scaffold mesh
104  // (It's zero [=default] if not visited this one yet)
105  unsigned j_global = global_number[scaffold_node_pt];
106 
107  // Haven't done this one yet
108  if (j_global == 0)
109  {
110  // Get pointer to set of mesh boundaries that this
111  // scaffold node occupies; NULL if the node is not on any boundary
112  std::set<unsigned>* boundaries_pt;
113  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
114 
115  // Is it on boundaries?
116  if (boundaries_pt != 0)
117  {
118  // Create new boundary node
119  Node* new_node_pt =
120  finite_element_pt(e)->construct_boundary_node(j, time_stepper_pt);
121 
122  // Give it a number (not necessarily the global node
123  // number in the scaffold mesh -- we just need something
124  // to keep track...)
125  global_count++;
126  global_number[scaffold_node_pt] = global_count;
127 
128  // Add to boundaries
129  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
130  it != boundaries_pt->end();
131  ++it)
132  {
133  add_boundary_node(*it, new_node_pt);
134  }
135  }
136  // Build normal node
137  else
138  {
139  // Create new normal node
140  finite_element_pt(e)->construct_node(j, time_stepper_pt);
141 
142  // Give it a number (not necessarily the global node
143  // number in the scaffold mesh -- we just need something
144  // to keep track...)
145  global_count++;
146  global_number[scaffold_node_pt] = global_count;
147  }
148 
149  // Copy new node, created using the NEW element's construct_node
150  // function into global storage, using the same global
151  // node number that we've just associated with the
152  // corresponding node in the scaffold mesh
153  Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
154 
155  // Assign coordinates
156  Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
157  Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
158  Node_pt[global_count - 1]->x(2) = scaffold_node_pt->x(2);
159  }
160  // This one has already been done: Copy across
161  else
162  {
163  finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
164  }
165  }
166 
167  // Store the attributes in the map
168  if (use_attributes)
169  {
170  element_attribute_map[Tmp_mesh_pt->element_attribute(e)].push_back(
172  }
173  }
174 
175  // Now let's construct lists
176  // Find the number of attributes
177  if (use_attributes)
178  {
179  unsigned n_attribute = element_attribute_map.size();
180  // There are n_attribute different regions
181  Region_element_pt.resize(n_attribute);
182  Region_attribute.resize(n_attribute);
183  // Copy the vectors in the map over to our internal storage
184  unsigned count = 0;
185  for (std::map<double, Vector<FiniteElement*>>::iterator it =
186  element_attribute_map.begin();
187  it != element_attribute_map.end();
188  ++it)
189  {
190  Region_attribute[count] = it->first;
191  Region_element_pt[count] = it->second;
192  ++count;
193  }
194  }
195 
196  // At this point we've created all the elements and
197  // created their vertex nodes. Now we need to create
198  // the additional (midside and internal) nodes!
199  unsigned boundary_id;
200 
201  // Get number of nodes along element edge and dimension of element (3)
202  // from first element
203  unsigned n_node_1d = finite_element_pt(0)->nnode_1d();
204 
205  // At the moment we're only able to deal with nnode_1d=2 or 3.
206  if ((n_node_1d != 2) && (n_node_1d != 3))
207  {
208  std::ostringstream error_message;
209  error_message << "Mesh generation from tetgen currently only works\n";
210  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
211  error_message << "for nnode_1d=" << n_node_1d << std::endl;
212 
213  throw OomphLibError(
214  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
215  }
216 
217  // Spatial dimension of element = number of local coordinates
218  unsigned dim = finite_element_pt(0)->dim();
219 
220  // Storage for the local coordinate of the new node
221  Vector<double> s(dim);
222 
223  // Get number of nodes in the element from first element
224  unsigned n_node = finite_element_pt(0)->nnode();
225 
226  // Storage for each global edge of the mesh
227  unsigned n_global_edge = Tmp_mesh_pt->nglobal_edge();
228  Vector<Vector<Node*>> nodes_on_global_edge(n_global_edge);
229 
230  // Storage for each global face of the mesh
231  unsigned n_global_face = Tmp_mesh_pt->nglobal_face();
232  Vector<Vector<Node*>> nodes_on_global_face(n_global_face);
233 
234  // Map storing the mid-side of an edge; edge identified by
235  // pointers to vertex nodes in scaffold mesh
236  // MapMatrix<Node*,Node*> central_edge_node_pt;
237  // Node* edge_node1_pt=0;
238  // Node* edge_node2_pt=0;
239 
240  // Map storing the mid point of a face; face identified by
241  // set of pointers to vertex nodes in scaffold mesh
242  // std::map<std::set<Node*>,Node*> central_face_node_pt;
243  // std::set<Node*> face_nodes_pt;
244 
245  // Mapping of Tetgen faces to face nodes in the enriched element
246  unsigned face_map[4] = {1, 0, 2, 3};
247 
248  // Storage for the faces shared by the edges
249  const unsigned faces_on_edge[6][2] = {
250  {0, 1}, {0, 2}, {1, 2}, {0, 3}, {2, 3}, {1, 3}};
251 
252  // Loop over all elements
253  for (unsigned e = 0; e < nelem; e++)
254  {
255  // Cache pointers to the elements
256  FiniteElement* const elem_pt = this->finite_element_pt(e);
257  FiniteElement* const tmp_elem_pt = Tmp_mesh_pt->finite_element_pt(e);
258 
259  // The number of edge nodes is 4 + 6*(n_node1d-2)
260  unsigned n_edge_node = 4 + 6 * (n_node_1d - 2);
261 
262  // Now loop over the edge nodes
263  // assuming that the numbering scheme is the same as the triangles
264  // which puts edge nodes in ascending order.
265  // We don't have any higher than quadratic at the moment, so it's
266  // a bit academic really.
267 
268  // Only bother if we actually have some edge nodes
269  if (n_edge_node > 4)
270  {
271  // Start from node number 4
272  unsigned n = 4;
273 
274  // Loop over the edges
275  for (unsigned j = 0; j < 6; ++j)
276  {
277  // Find the global edge index
278  unsigned edge_index = Tmp_mesh_pt->edge_index(e, j);
279 
280  // Use the intersection of the appropriate faces to determine
281  // whether the boundaries on which an edge lies
282  std::set<unsigned> edge_boundaries;
283  for (unsigned i = 0; i < 2; ++i)
284  {
285  unsigned face_boundary_id =
286  Tmp_mesh_pt->face_boundary(e, faces_on_edge[j][i]);
287  if (face_boundary_id > 0)
288  {
289  edge_boundaries.insert(face_boundary_id);
290  }
291  }
292 
293  // If the nodes on the edge have not been allocated, construct them
294  if (nodes_on_global_edge[edge_index].size() == 0)
295  {
296  // Now loop over the nodes on the edge
297  for (unsigned j2 = 0; j2 < n_node_1d - 2; ++j2)
298  {
299  // Storage for the new node
300  Node* new_node_pt = 0;
301 
302  // If the edge is on a boundary, determined from the
303  // scaffold mesh, construct a boundary node
304  // The use of the scaffold mesh's edge_boundary data structure
305  // ensures that a boundary node is created, even if the faces of
306  // the current element do not lie on boundaries
307  if (Tmp_mesh_pt->edge_boundary(edge_index) == true)
308  {
309  new_node_pt =
310  elem_pt->construct_boundary_node(n, time_stepper_pt);
311  // Add it to the boundaries in the set,
312  // remembering to subtract one to get to the oomph-lib numbering
313  // scheme
314  for (std::set<unsigned>::iterator it = edge_boundaries.begin();
315  it != edge_boundaries.end();
316  ++it)
317  {
318  this->add_boundary_node((*it) - 1, new_node_pt);
319  }
320  }
321  // Otherwise construct a normal node
322  else
323  {
324  new_node_pt = elem_pt->construct_node(n, time_stepper_pt);
325  }
326 
327  // Find the local coordinates of the node
328  elem_pt->local_coordinate_of_node(n, s);
329 
330  // Find the coordinates of the new node from the exiting and
331  // fully-functional element in the scaffold mesh
332  for (unsigned i = 0; i < dim; ++i)
333  {
334  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
335  }
336 
337  // Add the newly created node to the global node list
338  Node_pt.push_back(new_node_pt);
339  // Add to the edge index
340  nodes_on_global_edge[edge_index].push_back(new_node_pt);
341  // Increment the local node number
342  ++n;
343  } // end of loop over edge nodes
344  }
345  // Otherwise just set the pointers (orientation assumed the same)
346  else
347  {
348  for (unsigned j2 = 0; j2 < n_node_1d - 2; ++j2)
349  {
350  elem_pt->node_pt(n) = nodes_on_global_edge[edge_index][j2];
351  // It is possible that the edge may be on additional boundaries
352  // through another element
353  // So add again (note that this function will not add to
354  // boundaries twice)
355  for (std::set<unsigned>::iterator it = edge_boundaries.begin();
356  it != edge_boundaries.end();
357  ++it)
358  {
359  this->add_boundary_node((*it) - 1, elem_pt->node_pt(n));
360  }
361  ++n;
362  }
363  }
364  } // End of loop over edges
365 
366  // Deal with enriched elements
367  if (n_node == 15)
368  {
369  // Now loop over the faces
370  for (unsigned j = 0; j < 4; ++j)
371  {
372  // Find the boundary id of the face (need to map between node
373  // numbers and the face)
374  boundary_id = Tmp_mesh_pt->face_boundary(e, face_map[j]);
375 
376  // Find the global face index (need to map between node numbers and
377  // the face)
378  unsigned face_index = Tmp_mesh_pt->face_index(e, face_map[j]);
379 
380  // If the nodes on the face have not been allocated
381  if (nodes_on_global_face[face_index].size() == 0)
382  {
383  // Storage for the new node
384  Node* new_node_pt = 0;
385 
386  // If the face is on a boundary, construct a boundary node
387  if (boundary_id > 0)
388  {
389  new_node_pt =
390  elem_pt->construct_boundary_node(n, time_stepper_pt);
391  // Add it to the boundary
392  this->add_boundary_node(boundary_id - 1, new_node_pt);
393  }
394  // Otherwise construct a normal node
395  else
396  {
397  new_node_pt = elem_pt->construct_node(n, time_stepper_pt);
398  }
399 
400  // Find the local coordinates of the node
401  elem_pt->local_coordinate_of_node(n, s);
402 
403  // Find the coordinates of the new node from the exiting and
404  // fully-functional element in the scaffold mesh
405  for (unsigned i = 0; i < dim; ++i)
406  {
407  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
408  }
409 
410  // Add the newly created node to the global node list
411  Node_pt.push_back(new_node_pt);
412  // Add to the face index
413  nodes_on_global_face[face_index].push_back(new_node_pt);
414  // Increment the local node number
415  ++n;
416  }
417  // Otherwise just set the single node from the face element
418  else
419  {
420  elem_pt->node_pt(n) = nodes_on_global_face[face_index][0];
421  ++n;
422  }
423  } // end of loop over faces
424 
425  // Construct the element's central node, which is not on a boundary
426  {
427  Node* new_node_pt =
428  finite_element_pt(e)->construct_node(n, time_stepper_pt);
429  Node_pt.push_back(new_node_pt);
430 
431  // Find the local coordinates of the node
432  elem_pt->local_coordinate_of_node(n, s);
433 
434  // Find the coordinates of the new node from the existing
435  // and fully-functional element in the scaffold mesh
436  for (unsigned i = 0; i < dim; i++)
437  {
438  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
439  }
440  }
441  } // End of enriched case
442 
443  } // End of case for edge nodes
444 
445 
446  // Now loop over the faces to setup the information about elements
447  // adjacent to the boundary
448  for (unsigned j = 0; j < 4; ++j)
449  {
450  // Find the boundary id of the face
451  boundary_id = Tmp_mesh_pt->face_boundary(e, j);
452 
453  if (boundary_id > 0)
454  {
455  Boundary_element_pt[boundary_id - 1].push_back(elem_pt);
456  // Need to put a shift in here because of an inconsistent naming
457  // convention between tetgen and our faces
458  // Tetgen Face 0 is our Face 3
459  // Tetgen Face 1 is our Face 2
460  // Tetgen Face 2 is our Face 1
461  // Tetgen Face 3 is our Face 0
462  Face_index_at_boundary[boundary_id - 1].push_back(3 - j);
463 
464  // If using regions set up the boundary information
465  if (use_attributes)
466  {
467  // Element adjacent to boundary
468  Boundary_region_element_pt[boundary_id - 1]
469  [static_cast<unsigned>(
471  .push_back(elem_pt);
472  // Need to put a shift in here because of an inconsistent naming
473  // convention between triangle and face elements
474  Face_index_region_at_boundary[boundary_id - 1]
475  [static_cast<unsigned>(
477  .push_back(3 - j);
478  }
479  }
480  } // End of loop over faces
481 
482 
483  // Lookup scheme has now been setup
485 
486 
487  // /*
488 
489  // //For standard quadratic elements all nodes are edge nodes
490  // unsigned n_vertex_and_edge_node = n_node;
491  // //If we have enriched elements, there are only 10 vertex and edge
492  // nodes if(n_node==15)
493  // {
494  // //There are only 10 vertex and edge nodes
495  // n_vertex_and_edge_node = 10;
496  // }
497 
498  // // Loop over the new (edge) nodes in the element and create them.
499  // for(unsigned j=4;j<n_vertex_and_edge_node;j++)
500  // {
501 
502  // // Figure out edge nodes
503  // switch(j)
504  // {
505 
506  // // Node 4 is between nodes 0 and 1
507  // case 4:
508 
509  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
510  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
511  // break;
512 
513 
514  // // Node 5 is between nodes 0 and 2
515  // case 5:
516 
517  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
518  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
519  // break;
520 
521  // // Node 6 is between nodes 0 and 3
522  // case 6:
523 
524  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
525  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
526  // break;
527 
528  // // Node 7 is between nodes 1 and 2
529  // case 7:
530 
531  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
532  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
533  // break;
534 
535  // // Node 8 is between nodes 2 and 3
536  // case 8:
537 
538  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
539  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
540  // break;
541 
542  // // Node 9 is between nodes 1 and 3
543  // case 9:
544 
545  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
546  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
547  // break;
548 
549  // break;
550 
551  // //Error
552  // default:
553 
554  // throw OomphLibError("More than ten edge nodes in Tet Element",
555  // OOMPH_CURRENT_FUNCTION,
556  // OOMPH_EXCEPTION_LOCATION);
557  // }
558 
559 
560  // // Do we need a boundary node?
561  // bool need_boundary_node=false;
562 
563  // // hierher At some point fine tune via intersection of boundary
564  // sets if
565  // (edge_node1_pt->is_on_boundary()&&edge_node2_pt->is_on_boundary())
566  // {
567  // need_boundary_node=true;
568  // }
569 
570  // // Do we need a new node?
571  // if (central_edge_node_pt(edge_node1_pt,edge_node2_pt)==0)
572  // {
573  // Node* new_node_pt=0;
574 
575  // // Create new boundary node
576  // if (need_boundary_node)
577  // {
578  // new_node_pt=finite_element_pt(e)->
579  // construct_boundary_node(j,time_stepper_pt);
580  // }
581  // // Create new normal node
582  // else
583  // {
584  // new_node_pt=finite_element_pt(e)->
585  // construct_node(j,time_stepper_pt);
586  // }
587  // Node_pt.push_back(new_node_pt);
588 
589  // // Now indicate existence of newly created mideside node in map
590  // central_edge_node_pt(edge_node1_pt,edge_node2_pt)=new_node_pt;
591  // central_edge_node_pt(edge_node2_pt,edge_node1_pt)=new_node_pt;
592 
593  // // What are the node's local coordinates?
594  // finite_element_pt(e)->local_coordinate_of_node(j,s);
595 
596  // // Find the coordinates of the new node from the existing
597  // // and fully-functional element in the scaffold mesh
598  // for(unsigned i=0;i<dim;i++)
599  // {
600  // new_node_pt->x(i)=
601  // Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
602  // }
603  // }
604  // else
605  // {
606  // // Set pointer to the existing node
607  // finite_element_pt(e)->node_pt(j)=
608  // central_edge_node_pt(edge_node1_pt,edge_node2_pt);
609  // }
610 
611  // } // end of loop over new nodes
612 
613  // //Need to sort out the face nodes
614  // if(n_node==15)
615  // {
616  // // Loop over the new (face) nodes in the element and create them.
617  // for(unsigned j=10;j<14;j++)
618  // {
619  // //Empty the set of face nodes
620  // face_nodes_pt.clear();
621  // // Figure out face nodes
622  // switch(j)
623  // {
624 
625  // // Node 10 is between nodes 0 and 1 and 3
626  // case 10:
627 
628  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
629  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
630  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
631  // break;
632 
633  // // Node 11 is between nodes 0 and 1 and 2
634  // case 11:
635 
636  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
637  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
638  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
639  // break;
640 
641  // // Node 12 is between nodes 0 and 2 and 3
642  // case 12:
643 
644  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
645  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
646  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
647  // break;
648 
649  // // Node 13 is between nodes 1 and 2 and 3
650  // case 13:
651 
652  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
653  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
654  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
655  // break;
656 
657 
658  // //Error
659  // default:
660 
661  // throw OomphLibError("More than four face nodes in Tet Element",
662  // OOMPH_CURRENT_FUNCTION,
663  // OOMPH_EXCEPTION_LOCATION);
664  // }
665 
666  // // Do we need a boundary node?
667  // bool need_boundary_node=false;
668 
669  // //Work it out from the face boundary
670  // boundary_id = Tmp_mesh_pt->face_boundary(e,face_map[j-10]);
671  // //If it's non-zero then we do need to create a boundary node
672  // if(boundary_id!=0) {need_boundary_node=true;}
673 
674  // // Do we need a new node?
675  // if (central_face_node_pt[face_nodes_pt]==0)
676  // {
677  // Node* new_node_pt=0;
678 
679  // // Create a new boundary node
680  // if (need_boundary_node)
681  // {
682  // new_node_pt=finite_element_pt(e)->
683  // construct_boundary_node(j,time_stepper_pt);
684  // //Add it to the boundary
685  // add_boundary_node(boundary_id-1,new_node_pt);
686  // }
687  // // Create new normal node
688  // else
689  // {
690  // new_node_pt=finite_element_pt(e)->
691  // construct_node(j,time_stepper_pt);
692  // }
693  // Node_pt.push_back(new_node_pt);
694 
695  // // Now indicate existence of newly created mideside node in map
696  // central_face_node_pt[face_nodes_pt]=new_node_pt;
697 
698  // // What are the node's local coordinates?
699  // finite_element_pt(e)->local_coordinate_of_node(j,s);
700 
701  // // Find the coordinates of the new node from the existing
702  // // and fully-functional element in the scaffold mesh
703  // for(unsigned i=0;i<dim;i++)
704  // {
705  // new_node_pt->x(i)=
706  // Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
707  // }
708  // }
709  // else
710  // {
711  // // Set pointer to the existing node
712  // finite_element_pt(e)->node_pt(j)=
713  // central_face_node_pt[face_nodes_pt];
714  // }
715  // } //End of loop over face nodes
716 
717  // //Construct the element's central node, which is not on a boundary
718  // {
719  // Node* new_node_pt=
720  // finite_element_pt(e)->construct_node(14,time_stepper_pt);
721  // Node_pt.push_back(new_node_pt);
722 
723  // // What are the node's local coordinates?
724  // finite_element_pt(e)->local_coordinate_of_node(14,s);
725  // // Find the coordinates of the new node from the existing
726  // // and fully-functional element in the scaffold mesh
727  // for(unsigned i=0;i<dim;i++)
728  // {
729  // new_node_pt->x(i)=
730  // Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
731  // }
732  // }
733  // } //End of enriched case
734 
735  // } //end of loop over elements
736 
737 
738  // //Boundary conditions
739 
740  // // Matrix map to check if a node has already been added to
741  // // the boundary number b
742  // MapMatrixMixed<Node*,unsigned,bool> bound_node_done;
743 
744  // // Loop over elements
745  // for (unsigned e=0;e<nelem;e++)
746  // {
747  // // Loop over new local nodes
748  // for(unsigned j=4;j<n_node;j++)
749  // {
750  // // Pointer to the element's local node
751  // Node* loc_node_pt=finite_element_pt(e)->node_pt(j);
752 
753  // // This will have to be changed for higher-order elements
754  // //=======================================================
755 
756  // // These are the face nodes on the element's face 0:
757  // if ( (j==4) || (j==5) || (j==7) )
758  // {
759  // boundary_id=Tmp_mesh_pt->face_boundary(e,0);
760  // if (boundary_id!=0)
761  // {
762  // if (!bound_node_done(loc_node_pt,boundary_id-1))
763  // {
764  // add_boundary_node(boundary_id-1,loc_node_pt);
765  // bound_node_done(loc_node_pt,boundary_id-1)=true;
766  // }
767  // }
768  // }
769 
770 
771  // // These are the face nodes on the element's face 1:
772  // if ( (j==4) || (j==6) || (j==9) )
773  // {
774  // boundary_id=Tmp_mesh_pt->face_boundary(e,1);
775  // if (boundary_id!=0)
776  // {
777  // if (!bound_node_done(loc_node_pt,boundary_id-1))
778  // {
779  // add_boundary_node(boundary_id-1,loc_node_pt);
780  // bound_node_done(loc_node_pt,boundary_id-1)=true;
781  // }
782  // }
783  // }
784 
785  // // These are the face nodes on the element's face 2:
786  // if ( (j==5) || (j==6) || (j==8) )
787  // {
788  // boundary_id=Tmp_mesh_pt->face_boundary(e,2);
789  // if (boundary_id!=0)
790  // {
791  // if (!bound_node_done(loc_node_pt,boundary_id-1))
792  // {
793  // add_boundary_node(boundary_id-1,loc_node_pt);
794  // bound_node_done(loc_node_pt,boundary_id-1)=true;
795  // }
796  // }
797  // }
798 
799 
800  // // These are the face nodes on the element's face 3:
801  // if ( (j==7) || (j==8) || (j==9) )
802  // {
803  // boundary_id=Tmp_mesh_pt->face_boundary(e,3);
804  // if (boundary_id!=0)
805  // {
806  // if (!bound_node_done(loc_node_pt,boundary_id-1))
807  // {
808  // add_boundary_node(boundary_id-1,loc_node_pt);
809  // bound_node_done(loc_node_pt,boundary_id-1)=true;
810  // }
811  // }
812  // }
813 
814  // } //end for j
815 
816  // */
817 
818  } // end for e
819 
820 
821  } // end function
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
virtual Node * construct_node(const unsigned &n)
Definition: elements.h:2509
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned dim() const
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual unsigned nnode_1d() const
Definition: elements.h:2218
virtual Node * construct_boundary_node(const unsigned &n)
Definition: elements.h:2538
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
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
Vector< Vector< FiniteElement * > > Boundary_element_pt
Definition: mesh.h:91
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:87
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
Vector< Vector< int > > Face_index_at_boundary
Definition: mesh.h:95
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Definition: nodes.h:1365
Vector< double > Region_attribute
Definition: tet_mesh.h:1027
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
Definition: tet_mesh.h:1031
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Definition: tet_mesh.h:1035
Vector< Vector< FiniteElement * > > Region_element_pt
Definition: tet_mesh.h:1022
unsigned face_boundary(const unsigned &e, const unsigned &i) const
Definition: tetgen_scaffold_mesh.h:71
unsigned nglobal_edge()
Return the number of internal edges.
Definition: tetgen_scaffold_mesh.h:77
bool edge_boundary(const unsigned &i)
Definition: tetgen_scaffold_mesh.h:84
double element_attribute(const unsigned &e) const
Definition: tetgen_scaffold_mesh.h:112
unsigned edge_index(const unsigned &e, const unsigned &i) const
Definition: tetgen_scaffold_mesh.h:91
unsigned nglobal_face()
Return the number of internal face.
Definition: tetgen_scaffold_mesh.h:97
unsigned face_index(const unsigned &e, const unsigned &i) const
Definition: tetgen_scaffold_mesh.h:104
RealScalar s
Definition: level1_cplx_impl.h:130
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References oomph::FiniteElement::construct_boundary_node(), oomph::FiniteElement::construct_node(), e(), oomph::Node::get_boundaries_pt(), i, oomph::FiniteElement::interpolated_x(), j, oomph::FiniteElement::local_coordinate_of_node(), n, oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, size, and oomph::Node::x().

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), and oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ build_tetgenio()

template<class ELEMENT >
void oomph::TetgenMesh< ELEMENT >::build_tetgenio ( TetMeshFacetedSurface *const &  outer_boundary_pt,
Vector< TetMeshFacetedSurface * > &  internal_surface_pt,
Vector< double > *const &  target_element_volume_in_region_pt,
tetgenio tetgen_io 
)
inline

Build tetgenio object from the TetMeshFacetedSurfaces.

ALH: This may not be needed

402  {
403  // Pointer to Tetgen facet
405  // Pointer to Tetgen polygon
407 
408  // Start all indices from 1 (it's a choice and we've made it
409  tetgen_io.firstnumber = 1;
411  tetgen_io.useindex = true;
412 
413  // Find the number of internal surfaces
414  const unsigned n_internal = internal_surface_pt.size();
415 
416  // Find the number of points on the outer boundary
417  const unsigned n_outer_vertex = outer_boundary_pt->nvertex();
418  tetgen_io.numberofpoints = n_outer_vertex;
419 
420  // Find the number of points on the inner boundaries and add to the totals
421  Vector<unsigned> n_internal_vertex(n_internal);
422  Vector<unsigned> internal_vertex_offset(n_internal);
423  for (unsigned h = 0; h < n_internal; ++h)
424  {
425  n_internal_vertex[h] = internal_surface_pt[h]->nvertex();
426  internal_vertex_offset[h] = tetgen_io.numberofpoints;
427  tetgen_io.numberofpoints += n_internal_vertex[h];
428  }
429 
430  // Read the data into the point list
431  tetgen_io.pointlist = new double[tetgen_io.numberofpoints * 3];
432  tetgen_io.pointmarkerlist = new int[tetgen_io.numberofpoints];
433  unsigned counter = 0;
434  for (unsigned n = 0; n < n_outer_vertex; ++n)
435  {
436  for (unsigned i = 0; i < 3; ++i)
437  {
438  tetgen_io.pointlist[counter] =
439  outer_boundary_pt->vertex_coordinate(n, i);
440  ++counter;
441  }
442  }
443  for (unsigned h = 0; h < n_internal; ++h)
444  {
445  const unsigned n_inner = n_internal_vertex[h];
446  for (unsigned n = 0; n < n_inner; ++n)
447  {
448  for (unsigned i = 0; i < 3; ++i)
449  {
450  tetgen_io.pointlist[counter] =
451  internal_surface_pt[h]->vertex_coordinate(n, i);
452  ++counter;
453  }
454  }
455  }
456 
457 
458  // Set up the pointmarkers
459  counter = 0;
460  for (unsigned n = 0; n < n_outer_vertex; ++n)
461  {
462  tetgen_io.pointmarkerlist[counter] =
463  outer_boundary_pt->one_based_vertex_boundary_id(n);
464  ++counter;
465  }
466  for (unsigned h = 0; h < n_internal; ++h)
467  {
468  const unsigned n_inner = n_internal_vertex[h];
469  for (unsigned n = 0; n < n_inner; ++n)
470  {
471  tetgen_io.pointmarkerlist[counter] =
472  internal_surface_pt[h]->one_based_vertex_boundary_id(n);
473  ++counter;
474  }
475  }
476 
477 
478  // Now the facets
479  unsigned n_outer_facet = outer_boundary_pt->nfacet();
480  tetgen_io.numberoffacets = n_outer_facet;
481  Vector<unsigned> n_inner_facet(n_internal);
482  for (unsigned h = 0; h < n_internal; ++h)
483  {
484  n_inner_facet[h] = internal_surface_pt[h]->nfacet();
485  tetgen_io.numberoffacets += n_inner_facet[h];
486  }
487 
488  tetgen_io.facetlist = new tetgenio::facet[tetgen_io.numberoffacets];
489  tetgen_io.facetmarkerlist = new int[tetgen_io.numberoffacets];
490 
491 
492  counter = 0;
493  for (unsigned n = 0; n < n_outer_facet; ++n)
494  {
495  // Set pointer to facet
496  f = &tetgen_io.facetlist[counter];
497  f->numberofpolygons = 1;
498  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
499  f->numberofholes = 0;
500  f->holelist = NULL;
501  p = &f->polygonlist[0];
502 
503  Vector<unsigned> facet = outer_boundary_pt->vertex_index_in_tetgen(n);
504 
505  p->numberofvertices = facet.size();
506  p->vertexlist = new int[p->numberofvertices];
507  for (int i = 0; i < p->numberofvertices; ++i)
508  {
509  // The offset here is because we have insisted on 1-based indexing
510  p->vertexlist[i] = facet[i] + 1;
511  }
512 
513  // Set up the boundary markers
514  tetgen_io.facetmarkerlist[counter] =
515  outer_boundary_pt->one_based_facet_boundary_id(n);
516  // Increase the counter
517  ++counter;
518  }
519 
520  // Initialise the number of holes
521  tetgen_io.numberofholes = 0;
522  // and the number of regions
523  tetgen_io.numberofregions = 0;
524 
525  // Loop over the internal stuff
526  for (unsigned h = 0; h < n_internal; ++h)
527  {
528  for (unsigned n = 0; n < n_inner_facet[h]; ++n)
529  {
530  // Set pointer to facet
531  f = &tetgen_io.facetlist[counter];
532  f->numberofpolygons = 1;
533  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
534  f->numberofholes = 0;
535  f->holelist = NULL;
536  p = &f->polygonlist[0];
537 
538  Vector<unsigned> facet =
539  internal_surface_pt[h]->vertex_index_in_tetgen(n);
540 
541  p->numberofvertices = facet.size();
542  p->vertexlist = new int[p->numberofvertices];
543  for (int i = 0; i < p->numberofvertices; ++i)
544  {
545  // Add the 1-based and vertex offsets to get these number correct
546  p->vertexlist[i] = facet[i] + internal_vertex_offset[h] + 1;
547  }
548  // Set up the boundary markers
549  tetgen_io.facetmarkerlist[counter] =
550  internal_surface_pt[h]->one_based_facet_boundary_id(n);
551  ++counter;
552  }
553 
554  // If it's a hole/region add it
555  TetMeshFacetedClosedSurface* srf_pt =
556  dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[h]);
557  if (srf_pt != 0)
558  {
559  unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
560  for (unsigned j = 0; j < n_int_pts; j++)
561  {
562  if (srf_pt->internal_point_identifies_hole_for_tetgen(j))
563  {
564  ++tetgen_io.numberofholes;
565  }
566  // Otherwise it may be region
567  else
568  {
569  if (srf_pt->internal_point_identifies_region_for_tetgen(j))
570  {
571  ++tetgen_io.numberofregions;
572  }
573  }
574  }
575  }
576  }
577 
578  // Set storage for the holes
579  tetgen_io.holelist = new double[3 * tetgen_io.numberofholes];
580 
581  // Loop over all the internal boundaries again
582  counter = 0;
583  for (unsigned h = 0; h < n_internal; ++h)
584  {
585  TetMeshFacetedClosedSurface* srf_pt =
586  dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[h]);
587  if (srf_pt != 0)
588  {
589  unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
590  for (unsigned j = 0; j < n_int_pts; j++)
591  {
592  if (srf_pt->internal_point_identifies_hole_for_tetgen(j))
593  {
594  for (unsigned i = 0; i < 3; ++i)
595  {
596  tetgen_io.holelist[counter] =
597  srf_pt->internal_point_for_tetgen(j, i);
598  ++counter;
599  }
600  }
601  }
602  }
603  }
604 
605  // Set storage for the regions
606  tetgen_io.regionlist = new double[5 * tetgen_io.numberofregions];
607 
608  // Loop over all the internal boundaries again
609  counter = 0;
610  for (unsigned h = 0; h < n_internal; ++h)
611  {
612  TetMeshFacetedClosedSurface* srf_pt =
613  dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[h]);
614  if (srf_pt != 0)
615  {
616  unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
617  for (unsigned j = 0; j < n_int_pts; j++)
618  {
619  if (srf_pt->internal_point_identifies_region_for_tetgen(j))
620  {
621  for (unsigned i = 0; i < 3; ++i)
622  {
623  tetgen_io.regionlist[counter] =
624  srf_pt->internal_point_for_tetgen(j, i);
625  ++counter;
626  }
627  tetgen_io.regionlist[counter] =
628  static_cast<double>(srf_pt->region_id_for_tetgen(j));
629  ++counter;
630 
631  // if there's no target volumes specified, default to zero
632  if (target_element_volume_in_region_pt == nullptr)
633  {
634  tetgen_io.regionlist[counter] = 0;
635  }
636  else
637  {
638  // deliberate integer division here to round down to the region
639  // number (five doubles per region)
640  tetgen_io.regionlist[counter] =
641  (*target_element_volume_in_region_pt)[unsigned(counter / 5)];
642  }
643 
644  ++counter;
645  }
646  }
647  }
648  }
649  }
float * p
Definition: Tutorial_Map_using.cpp:9
int numberofregions
Definition: tetgen.h:349
bool useindex
Definition: tetgen.h:291
int * facetmarkerlist
Definition: tetgen.h:332
REAL * holelist
Definition: tetgen.h:338
int numberofpoints
Definition: tetgen.h:306
REAL * pointlist
Definition: tetgen.h:302
REAL * regionlist
Definition: tetgen.h:348
int firstnumber
Definition: tetgen.h:285
int numberofholes
Definition: tetgen.h:339
facet * facetlist
Definition: tetgen.h:331
int numberoffacets
Definition: tetgen.h:333
int * pointmarkerlist
Definition: tetgen.h:305
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
Definition: tetgen.h:229
Definition: tetgen.h:217

References f(), tetgenio::facetlist, tetgenio::facetmarkerlist, tetgenio::firstnumber, tetgenio::holelist, i, oomph::TetMeshFacetedClosedSurface::internal_point_for_tetgen(), oomph::TetMeshFacetedClosedSurface::internal_point_identifies_hole_for_tetgen(), oomph::TetMeshFacetedClosedSurface::internal_point_identifies_region_for_tetgen(), j, n, oomph::TetMeshFacetedSurface::nfacet(), oomph::TetMeshFacetedClosedSurface::ninternal_point_for_tetgen(), tetgenio::numberoffacets, tetgenio::numberofholes, tetgenio::numberofpoints, tetgenio::numberofregions, oomph::TetMeshFacetedSurface::nvertex(), oomph::TetMeshFacetedSurface::one_based_facet_boundary_id(), oomph::TetMeshFacetedSurface::one_based_vertex_boundary_id(), p, tetgenio::pointlist, tetgenio::pointmarkerlist, oomph::TetMeshFacetedClosedSurface::region_id_for_tetgen(), tetgenio::regionlist, tetgenio::useindex, oomph::TetMeshFacetedSurface::vertex_coordinate(), and oomph::TetMeshFacetedSurface::vertex_index_in_tetgen().

Referenced by oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ deep_copy_of_tetgenio()

template<class ELEMENT >
void TetgenMesh::deep_copy_of_tetgenio ( tetgenio *const &  input_pt,
tetgenio *&  output_pt 
)

Transfer tetgenio data from the input to the output The output is assumed to have been constructed and "empty"

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// Transfer tetgenio data from the input to the output The output is assumed to have been constructed and "empty"

867  {
868  // Copy data values
869  output_pt->firstnumber = input_pt->firstnumber;
870  output_pt->mesh_dim = input_pt->mesh_dim;
871  output_pt->useindex = input_pt->useindex;
872 
873  // Copy the number of points
874  output_pt->numberofpoints = input_pt->numberofpoints;
875  output_pt->numberofpointattributes = input_pt->numberofpointattributes;
876  output_pt->numberofpointmtrs = input_pt->numberofpointmtrs;
877 
878  const int n_point = output_pt->numberofpoints;
879  if (n_point > 0)
880  {
881  output_pt->pointlist = new double[n_point * 3];
882  // Copy the values
883  for (int n = 0; n < (n_point * 3); ++n)
884  {
885  output_pt->pointlist[n] = input_pt->pointlist[n];
886  }
887 
888  // If there are point markers copy as well
889  if (input_pt->pointmarkerlist != (int*)NULL)
890  {
891  output_pt->pointmarkerlist = new int[n_point];
892  for (int n = 0; n < n_point; ++n)
893  {
894  output_pt->pointmarkerlist[n] = input_pt->pointmarkerlist[n];
895  }
896  }
897 
898  const int n_attr = output_pt->numberofpointattributes;
899  if (n_attr > 0)
900  {
901  output_pt->pointattributelist = new double[n_point * n_attr];
902  for (int n = 0; n < (n_point * n_attr); ++n)
903  {
904  output_pt->pointattributelist[n] = input_pt->pointattributelist[n];
905  }
906  }
907  } // End of point case
908 
909  // Now add in metric tensors (if there are any)
910  const int n_point_mtr = output_pt->numberofpointmtrs;
911  if (n_point_mtr > 0)
912  {
913  output_pt->pointmtrlist = new double[n_point * n_point_mtr];
914  for (int n = 0; n < (n_point * n_point_mtr); ++n)
915  {
916  output_pt->pointmtrlist[n] = input_pt->pointmtrlist[n];
917  }
918  }
919 
920  // Next tetrahedrons
921  output_pt->numberoftetrahedra = input_pt->numberoftetrahedra;
922  output_pt->numberofcorners = input_pt->numberofcorners;
923  output_pt->numberoftetrahedronattributes =
925 
926  const int n_tetra = output_pt->numberoftetrahedra;
927  const int n_corner = output_pt->numberofcorners;
928  if (n_tetra > 0)
929  {
930  output_pt->tetrahedronlist = new int[n_tetra * n_corner];
931  for (int n = 0; n < (n_tetra * n_corner); ++n)
932  {
933  output_pt->tetrahedronlist[n] = input_pt->tetrahedronlist[n];
934  }
935 
936  // Add in the volume constriants
937  if (input_pt->tetrahedronvolumelist != (double*)NULL)
938  {
939  output_pt->tetrahedronvolumelist = new double[n_tetra];
940  for (int n = 0; n < n_tetra; ++n)
941  {
942  output_pt->tetrahedronvolumelist[n] =
943  input_pt->tetrahedronvolumelist[n];
944  }
945  }
946 
947  // Add in the attributes
948  const int n_tetra_attr = output_pt->numberoftetrahedronattributes;
949  if (n_tetra_attr > 0)
950  {
951  output_pt->tetrahedronattributelist =
952  new double[n_tetra * n_tetra_attr];
953  for (int n = 0; n < (n_tetra * n_tetra_attr); ++n)
954  {
955  output_pt->tetrahedronattributelist[n] =
956  input_pt->tetrahedronattributelist[n];
957  }
958  }
959 
960  // Add in the neighbour list
961  if (input_pt->neighborlist != (int*)NULL)
962  {
963  output_pt->neighborlist = new int[n_tetra * 4];
964  for (int n = 0; n < (n_tetra * 4); ++n)
965  {
966  output_pt->neighborlist = input_pt->neighborlist;
967  }
968  }
969  } // End of tetrahedron section
970 
971  // Now arrange the facet list
972  output_pt->numberoffacets = input_pt->numberoffacets;
973  const int n_facet = output_pt->numberoffacets;
974  if (n_facet > 0)
975  {
976  output_pt->facetlist = new tetgenio::facet[n_facet];
977  for (int n = 0; n < n_facet; ++n)
978  {
979  tetgenio::facet* input_f_pt = &input_pt->facetlist[n];
980  tetgenio::facet* output_f_pt = &output_pt->facetlist[n];
981 
982  // Copy polygons and holes from the facets
983  output_f_pt->numberofpolygons = input_f_pt->numberofpolygons;
984 
985  // Loop over polygons
986  const int n_poly = output_f_pt->numberofpolygons;
987  if (n_poly > 0)
988  {
989  output_f_pt->polygonlist = new tetgenio::polygon[n_poly];
990  for (int p = 0; p < n_poly; ++p)
991  {
992  tetgenio::polygon* output_p_pt = &output_f_pt->polygonlist[p];
993  tetgenio::polygon* input_p_pt = &input_f_pt->polygonlist[p];
994  // Find out how many vertices each polygon has
995  output_p_pt->numberofvertices = input_p_pt->numberofvertices;
996  // Now copy of the vertices
997  const int n_vertex = output_p_pt->numberofvertices;
998  if (n_vertex > 0)
999  {
1000  output_p_pt->vertexlist = new int[n_vertex];
1001  for (int v = 0; v < n_vertex; ++v)
1002  {
1003  output_p_pt->vertexlist[v] = input_p_pt->vertexlist[v];
1004  }
1005  }
1006  }
1007  }
1008 
1009  // Hole information
1010  output_f_pt->numberofholes = input_f_pt->numberofholes;
1011  const int n_hole = output_f_pt->numberofholes;
1012  if (n_hole > 0)
1013  {
1014  throw OomphLibError("Don't know how to handle holes yet\n",
1017  }
1018  } // end of loop over facets
1019 
1020  // Add the facetmarkers if there are any
1021  if (input_pt->facetmarkerlist != (int*)NULL)
1022  {
1023  output_pt->facetmarkerlist = new int[n_facet];
1024  for (int n = 0; n < n_facet; ++n)
1025  {
1026  output_pt->facetmarkerlist[n] = input_pt->facetmarkerlist[n];
1027  }
1028  }
1029  }
1030 
1031  // Now the holes
1032  output_pt->numberofholes = input_pt->numberofholes;
1033  const int n_hole = output_pt->numberofholes;
1034  if (n_hole > 0)
1035  {
1036  output_pt->holelist = new double[n_hole * 3];
1037  for (int n = 0; n < (n_hole * 3); ++n)
1038  {
1039  output_pt->holelist[n] = input_pt->holelist[n];
1040  }
1041  }
1042 
1043  // Now the regions
1044  output_pt->numberofregions = input_pt->numberofregions;
1045  const int n_region = output_pt->numberofregions;
1046  if (n_region > 0)
1047  {
1048  output_pt->regionlist = new double[n_region * 5];
1049  for (int n = 0; n < (n_region * 5); ++n)
1050  {
1051  output_pt->regionlist[n] = input_pt->regionlist[n];
1052  }
1053  }
1054 
1055  // Now the facet constraints
1056  output_pt->numberoffacetconstraints = input_pt->numberoffacetconstraints;
1057  const int n_facet_const = output_pt->numberoffacetconstraints;
1058  if (n_facet_const > 0)
1059  {
1060  output_pt->facetconstraintlist = new double[n_facet_const * 2];
1061  for (int n = 0; n < (n_facet_const * 2); ++n)
1062  {
1063  output_pt->facetconstraintlist[n] = input_pt->facetconstraintlist[n];
1064  }
1065  }
1066 
1067  // Now the segment constraints
1068  output_pt->numberofsegmentconstraints =
1069  input_pt->numberofsegmentconstraints;
1070  const int n_seg_const = output_pt->numberofsegmentconstraints;
1071  if (n_seg_const > 0)
1072  {
1073  output_pt->segmentconstraintlist = new double[n_seg_const * 2];
1074  for (int n = 0; n < (n_seg_const * 2); ++n)
1075  {
1076  output_pt->segmentconstraintlist[n] =
1077  input_pt->segmentconstraintlist[n];
1078  }
1079  }
1080 
1081  // Now the face list
1082  output_pt->numberoftrifaces = input_pt->numberoftrifaces;
1083  const int n_tri_face = output_pt->numberoftrifaces;
1084  if (n_tri_face > 0)
1085  {
1086  output_pt->trifacelist = new int[n_tri_face * 3];
1087  for (int n = 0; n < (n_tri_face * 3); ++n)
1088  {
1089  output_pt->trifacelist[n] = input_pt->trifacelist[n];
1090  }
1091 
1092  output_pt->trifacemarkerlist = new int[n_tri_face];
1093  for (int n = 0; n < n_tri_face; ++n)
1094  {
1095  output_pt->trifacemarkerlist[n] = input_pt->trifacemarkerlist[n];
1096  }
1097  }
1098 
1099  // Now the edge list
1100  output_pt->numberofedges = input_pt->numberofedges;
1101  const int n_edge = output_pt->numberofedges;
1102  if (n_edge > 0)
1103  {
1104  output_pt->edgelist = new int[n_edge * 2];
1105  for (int n = 0; n < (n_edge * 2); ++n)
1106  {
1107  output_pt->edgelist[n] = input_pt->edgelist[n];
1108  }
1109 
1110  output_pt->edgemarkerlist = new int[n_edge];
1111  for (int n = 0; n < n_edge; ++n)
1112  {
1113  output_pt->edgemarkerlist[n] = input_pt->edgemarkerlist[n];
1114  }
1115  }
1116  }
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
int numberofpointmtrs
Definition: tetgen.h:308
REAL * tetrahedronvolumelist
Definition: tetgen.h:323
int * trifacemarkerlist
Definition: tetgen.h:380
int * trifacelist
Definition: tetgen.h:378
int * edgemarkerlist
Definition: tetgen.h:388
int numberofpointattributes
Definition: tetgen.h:307
REAL * pointattributelist
Definition: tetgen.h:303
int numberofsegmentconstraints
Definition: tetgen.h:363
int numberofedges
Definition: tetgen.h:389
int mesh_dim
Definition: tetgen.h:288
int numberofcorners
Definition: tetgen.h:326
REAL * facetconstraintlist
Definition: tetgen.h:355
REAL * segmentconstraintlist
Definition: tetgen.h:362
int * neighborlist
Definition: tetgen.h:324
int numberoftrifaces
Definition: tetgen.h:381
int numberoftetrahedronattributes
Definition: tetgen.h:327
int numberoffacetconstraints
Definition: tetgen.h:356
int * edgelist
Definition: tetgen.h:387
int numberoftetrahedra
Definition: tetgen.h:325
REAL * pointmtrlist
Definition: tetgen.h:304
int * tetrahedronlist
Definition: tetgen.h:321
REAL * tetrahedronattributelist
Definition: tetgen.h:322
int numberofpolygons
Definition: tetgen.h:231
int numberofholes
Definition: tetgen.h:233
polygon * polygonlist
Definition: tetgen.h:230
int numberofvertices
Definition: tetgen.h:219
int * vertexlist
Definition: tetgen.h:218

References tetgenio::edgelist, tetgenio::edgemarkerlist, tetgenio::facetconstraintlist, tetgenio::facetlist, tetgenio::facetmarkerlist, tetgenio::firstnumber, tetgenio::holelist, tetgenio::mesh_dim, n, tetgenio::neighborlist, tetgenio::numberofcorners, tetgenio::numberofedges, tetgenio::numberoffacetconstraints, tetgenio::numberoffacets, tetgenio::facet::numberofholes, tetgenio::numberofholes, tetgenio::numberofpointattributes, tetgenio::numberofpointmtrs, tetgenio::numberofpoints, tetgenio::facet::numberofpolygons, tetgenio::numberofregions, tetgenio::numberofsegmentconstraints, tetgenio::numberoftetrahedra, tetgenio::numberoftetrahedronattributes, tetgenio::numberoftrifaces, tetgenio::polygon::numberofvertices, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, tetgenio::pointattributelist, tetgenio::pointlist, tetgenio::pointmarkerlist, tetgenio::pointmtrlist, tetgenio::facet::polygonlist, tetgenio::regionlist, tetgenio::segmentconstraintlist, tetgenio::tetrahedronattributelist, tetgenio::tetrahedronlist, tetgenio::tetrahedronvolumelist, tetgenio::trifacelist, tetgenio::trifacemarkerlist, tetgenio::useindex, v, and tetgenio::polygon::vertexlist.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), and oomph::TetgenMesh< ELEMENT >::set_deep_copy_tetgenio_pt().

◆ set_deep_copy_tetgenio_pt()

template<class ELEMENT >
void oomph::TetgenMesh< ELEMENT >::set_deep_copy_tetgenio_pt ( tetgenio *const &  tetgenio_pt)
inline

Set the tetgen pointer by a deep copy.

686  {
687  // Delete the existing data
688  if (Tetgenio_exists)
689  {
690  delete Tetgenio_pt;
691  }
692  this->Tetgenio_pt = new tetgenio;
693  // Create a deep copy of tetgenio_pt and store the result in
694  // Tetgenio_pt
695  this->deep_copy_of_tetgenio(tetgenio_pt, this->Tetgenio_pt);
696  }
void deep_copy_of_tetgenio(tetgenio *const &input_pt, tetgenio *&output_pt)
Definition: tetgen_mesh.template.cc:865

References oomph::TetgenMesh< ELEMENT >::deep_copy_of_tetgenio(), oomph::TetgenMesh< ELEMENT >::Tetgenio_exists, and oomph::TetgenMesh< ELEMENT >::Tetgenio_pt.

◆ set_mesh_level_time_stepper()

template<class ELEMENT >
void oomph::TetgenMesh< ELEMENT >::set_mesh_level_time_stepper ( TimeStepper *const &  time_stepper_pt,
const bool preserve_existing_data 
)
inlinevirtual

Overload set_mesh_level_time_stepper so that the stored time stepper now corresponds to the new timestepper

Reimplemented from oomph::Mesh.

666  {
667  this->Time_stepper_pt = time_stepper_pt;
668  }

References oomph::TetMeshBase::Time_stepper_pt.

◆ setup_reverse_lookup_schemes_for_faceted_surface()

template<class ELEMENT >
void TetgenMesh::setup_reverse_lookup_schemes_for_faceted_surface ( TetMeshFacetedSurface *const &  faceted_surface_pt)
protected

Function to setup the reverse look-up schemes.

Helper function to set up the reverse look up scheme for facets. This is used to set up boundary coordinates.

830  {
831  // Set up reverse lookup scheme for a given faceted surface
832  // Assumption is that there there is one boundary ID per facet.
833  unsigned n_facet = faceted_surface_pt->nfacet();
834  for (unsigned f = 0; f < n_facet; f++)
835  {
836  unsigned b = faceted_surface_pt->one_based_facet_boundary_id(f);
837  if (b != 0)
838  {
839  this->Tet_mesh_faceted_surface_pt[b - 1] = faceted_surface_pt;
840  this->Tet_mesh_facet_pt[b - 1] = faceted_surface_pt->facet_pt(f);
841  }
842  else
843  {
844  std::ostringstream error_message;
845  error_message << "Boundary IDs have to be one-based. Yours is " << b
846  << "\n";
847  throw OomphLibError(error_message.str(),
850  }
851  }
852  }
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
Definition: tet_mesh.h:1045
std::map< unsigned, TetMeshFacet * > Tet_mesh_facet_pt
Definition: tet_mesh.h:1049

References b, f(), oomph::TetMeshFacetedSurface::facet_pt(), oomph::TetMeshFacetedSurface::nfacet(), oomph::TetMeshFacetedSurface::one_based_facet_boundary_id(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::RefineableTetgenMesh(), and oomph::TetgenMesh< ELEMENT >::TetgenMesh().

◆ tetgenio_exists()

template<class ELEMENT >
bool oomph::TetgenMesh< ELEMENT >::tetgenio_exists ( ) const
inline

Boolen defining whether tetgenio object has been built or not.

673  {
674  return Tetgenio_exists;
675  }

References oomph::TetgenMesh< ELEMENT >::Tetgenio_exists.

◆ tetgenio_pt()

template<class ELEMENT >
tetgenio*& oomph::TetgenMesh< ELEMENT >::tetgenio_pt ( )
inline

Access to the triangulateio representation of the mesh.

680  {
681  return Tetgenio_pt;
682  }

References oomph::TetgenMesh< ELEMENT >::Tetgenio_pt.

Referenced by oomph::RefineableTetgenMesh< ELEMENT >::adapt().

Member Data Documentation

◆ Tetgenio_exists

◆ Tetgenio_pt

◆ Tmp_mesh_pt

template<class ELEMENT >
TetgenScaffoldMesh* oomph::TetgenMesh< ELEMENT >::Tmp_mesh_pt
protected

◆ Use_attributes

template<class ELEMENT >
bool oomph::TetgenMesh< ELEMENT >::Use_attributes
protected

Boolean flag to indicate whether to use attributes or not (required for multidomain meshes)

Referenced by oomph::TetgenMesh< ELEMENT >::TetgenMesh().


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