tetgen_mesh.template.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 
27 #ifndef OOMPH_TETGEN_MESH_HEADER
28 #define OOMPH_TETGEN_MESH_HEADER
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 
36 #ifdef OOMPH_HAS_MPI
37 // mpi headers
38 #include "mpi.h"
39 #endif
40 
41 #include "../generic/tetgen_scaffold_mesh.h"
42 #include "../generic/tet_mesh.h"
43 
44 namespace oomph
45 {
46  //=========start of TetgenMesh class======================================
47  /// Unstructured tet mesh based on output from Tetgen:
48  /// http://wias-berlin.de/software/tetgen/
49  //========================================================================
50  template<class ELEMENT>
51  class TetgenMesh : public virtual TetMeshBase
52  {
53  public:
54  /// Empty constructor
56  {
57  // Mesh can only be built with 3D Telements.
58  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
59  }
60 
61  /// Constructor with the input files
62  TetgenMesh(const std::string& node_file_name,
63  const std::string& element_file_name,
64  const std::string& face_file_name,
65  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
66  const bool& use_attributes = false)
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  }
97 
98 
99  /// Constructor with tetgenio data structure
100  TetgenMesh(tetgenio& tetgen_data,
101  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
102  const bool& use_attributes = false)
103 
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  }
136 
137 
138  /// Constructor with the input files. Setting the boolean
139  /// flag to true splits "corner" elements, i.e. elements that
140  /// that have at least three faces on a domain boundary. The
141  /// relevant elements are split without introducing hanging
142  /// nodes so the sons have a "worse" shape than their fathers.
143  /// However, this step avoids otherwise-hard-to-diagnose
144  /// problems in fluids problems where the application of
145  /// boundary conditions at such "corner" elements can
146  /// overconstrain the solution.
147  TetgenMesh(const std::string& node_file_name,
148  const std::string& element_file_name,
149  const std::string& face_file_name,
150  const bool& split_corner_elements,
151  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
152  const bool& use_attributes = false)
153 
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  }
193 
194  /// Constructor with tetgen data structure Setting the boolean
195  /// flag to true splits "corner" elements, i.e. elements that
196  /// that have at least three faces on a domain boundary. The
197  /// relevant elements are split without introducing hanging
198  /// nodes so the sons have a "worse" shape than their fathers.
199  /// However, this step avoids otherwise-hard-to-diagnose
200  /// problems in fluids problems where the application of
201  /// boundary conditions at such "corner" elements can
202  /// overconstrain the solution.
203  TetgenMesh(tetgenio& tetgen_data,
204  const bool& split_corner_elements,
205  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
206  const bool& use_attributes = false)
207 
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  }
246 
247 
248  /// Build mesh, based on a TetgenMeshFactedClosedSurface that
249  /// specifies the outer boundary of the domain and any number of internal
250  /// boundaries, specified by TetMeshFacetedSurfaces.
251  /// Also specify target size for uniform element size.
252  /// Optionally specify the target element volume in each region.
254  TetMeshFacetedClosedSurface* const& outer_boundary_pt,
255  Vector<TetMeshFacetedSurface*>& internal_surface_pt,
256  const double& element_volume,
257  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
258  const bool& use_attributes = false,
259  const bool& split_corner_elements = false,
260  Vector<double>* const& target_element_volume_in_region_pt = nullptr)
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
274  this->setup_reverse_lookup_schemes_for_faceted_surface(Outer_boundary_pt);
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  {
283  Internal_surface_pt[i]);
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!)
393  snap_nodes_onto_geometric_objects();
394  }
395 
396  /// Build tetgenio object from the TetMeshFacetedSurfaces
398  TetMeshFacetedSurface* const& outer_boundary_pt,
399  Vector<TetMeshFacetedSurface*>& internal_surface_pt,
400  Vector<double>* const& target_element_volume_in_region_pt,
401  tetgenio& tetgen_io)
402  {
403  // Pointer to Tetgen facet
404  tetgenio::facet* f;
405  // Pointer to Tetgen polygon
406  tetgenio::polygon* p;
407 
408  // Start all indices from 1 (it's a choice and we've made it
409  tetgen_io.firstnumber = 1;
410  /// ALH: This may not be needed
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  }
650 
651 
652  /// Empty destructor
654  {
655  if (Tetgenio_exists)
656  {
657  delete Tetgenio_pt;
658  }
659  }
660 
661 
662  /// Overload set_mesh_level_time_stepper so that the stored
663  /// time stepper now corresponds to the new timestepper
664  void set_mesh_level_time_stepper(TimeStepper* const& time_stepper_pt,
665  const bool& preserve_existing_data)
666  {
667  this->Time_stepper_pt = time_stepper_pt;
668  }
669 
670 
671  /// Boolen defining whether tetgenio object has been built or not
672  bool tetgenio_exists() const
673  {
674  return Tetgenio_exists;
675  }
676 
677 
678  /// Access to the triangulateio representation of the mesh
679  tetgenio*& tetgenio_pt()
680  {
681  return Tetgenio_pt;
682  }
683 
684  /// Set the tetgen pointer by a deep copy
685  void set_deep_copy_tetgenio_pt(tetgenio* const& tetgenio_pt)
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  }
697 
698  /// Transfer tetgenio data from the input to the output
699  /// The output is assumed to have been constructed and "empty"
700  void deep_copy_of_tetgenio(tetgenio* const& input_pt, tetgenio*& output_pt);
701 
702 
703  protected:
704  /// Build mesh from scaffold
705  void build_from_scaffold(TimeStepper* time_stepper_pt,
706  const bool& use_attributes);
707 
708  /// Function to setup the reverse look-up schemes
710  TetMeshFacetedSurface* const& faceted_surface_pt);
711 
712  /// Temporary scaffold mesh
713  TetgenScaffoldMesh* Tmp_mesh_pt;
714 
715  /// Boolean to indicate whether a tetgenio representation of the
716  /// mesh exists
718 
719  /// Tetgen representation of mesh
720  tetgenio* Tetgenio_pt;
721 
722  /// Boolean flag to indicate whether to use attributes or not
723  /// (required for multidomain meshes)
725 
726  }; // end class
727 
728 
729  /// ///////////////////////////////////////////////////////////////
730  /// ///////////////////////////////////////////////////////////////
731  /// ///////////////////////////////////////////////////////////////
732 
733 
734  //==============start_mesh=================================================
735  /// Tetgen-based mesh upgraded to become a solid mesh. Automatically
736  /// enumerates all boundaries.
737  //=========================================================================
738  template<class ELEMENT>
739  class SolidTetgenMesh : public virtual TetgenMesh<ELEMENT>,
740  public virtual SolidMesh
741  {
742  public:
743  /// Constructor. Boundary coordinates are setup
744  /// automatically.
745  SolidTetgenMesh(const std::string& node_file_name,
746  const std::string& element_file_name,
747  const std::string& face_file_name,
748  const bool& split_corner_elements,
749  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
750  const bool& use_attributes = false)
751  : TetgenMesh<ELEMENT>(node_file_name,
752  element_file_name,
753  face_file_name,
754  split_corner_elements,
755  time_stepper_pt,
756  use_attributes)
757  {
758  // Assign the Lagrangian coordinates
759  set_lagrangian_nodal_coordinates();
760  }
761 
762  /// Constructor. Boundary coordinates are re-setup
763  /// automatically, with the orientation of the outer unit
764  /// normal determined by switch_normal.
765  SolidTetgenMesh(const std::string& node_file_name,
766  const std::string& element_file_name,
767  const std::string& face_file_name,
768  const bool& split_corner_elements,
769  const bool& switch_normal,
770  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
771  const bool& use_attributes = false)
772  : TetgenMesh<ELEMENT>(node_file_name,
773  element_file_name,
774  face_file_name,
775  split_corner_elements,
776  time_stepper_pt,
777  use_attributes)
778  {
779  // Assign the Lagrangian coordinates
780  set_lagrangian_nodal_coordinates();
781 
782  // Re-setup boundary coordinates for all boundaries with specified
783  // orientation of nnormal
784  unsigned nb = this->nboundary();
785  for (unsigned b = 0; b < nb; b++)
786  {
787  this->template setup_boundary_coordinates<ELEMENT>(b, switch_normal);
788  }
789  }
790 
791  /// Empty Destructor
792  virtual ~SolidTetgenMesh() {}
793  };
794 
795 
796 } // namespace oomph
797 
798 #endif
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
SolidTetgenMesh(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)
Constructor. Boundary coordinates are setup automatically.
SolidTetgenMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &face_file_name, const bool &split_corner_elements, const bool &switch_normal, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
Constructor. Boundary coordinates are re-setup automatically, with the orientation of the outer unit ...
virtual ~SolidTetgenMesh()
Empty Destructor.
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen/.
tetgenio *& tetgenio_pt()
Access to the triangulateio representation of the mesh.
TetgenMesh(tetgenio &tetgen_data, const bool &split_corner_elements, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
Constructor with tetgen data structure Setting the boolean flag to true splits "corner" elements,...
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.
TetgenMesh()
Empty constructor.
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
tetgenio * Tetgenio_pt
Tetgen representation of mesh.
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)
Build mesh, based on a TetgenMeshFactedClosedSurface that specifies the outer boundary of the domain ...
void setup_reverse_lookup_schemes_for_faceted_surface(TetMeshFacetedSurface *const &faceted_surface_pt)
Function to setup the reverse look-up schemes.
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)
Constructor with the input files. Setting the boolean flag to true splits "corner" elements,...
void set_deep_copy_tetgenio_pt(tetgenio *const &tetgenio_pt)
Set the tetgen pointer by a deep copy.
void 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 an...
bool Use_attributes
Boolean flag to indicate whether to use attributes or not (required for multidomain meshes)
~TetgenMesh()
Empty destructor.
void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Overload set_mesh_level_time_stepper so that the stored time stepper now corresponds to the new times...
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.
bool Tetgenio_exists
Boolean to indicate whether a tetgenio representation of the mesh exists.
TetgenScaffoldMesh * Tmp_mesh_pt
Temporary scaffold mesh.
TetgenMesh(tetgenio &tetgen_data, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
Constructor with tetgenio data structure.
bool tetgenio_exists() const
Boolen defining whether tetgenio object has been built or not.
////////////////////////////////////////////////////////////////////// //////////////////////////////...