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// Version 1.0; svn revision $LastChangedRevision: 1174 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-11 10:03:56 +0100 (Wed, 11 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #ifndef OOMPH_TETGEN_MESH_HEADER
32 #define OOMPH_TETGEN_MESH_HEADER
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 
40 #ifdef OOMPH_HAS_MPI
41 //mpi headers
42 #include "mpi.h"
43 #endif
44 
45 #include "../generic/tetgen_scaffold_mesh.h"
46 #include "../generic/tet_mesh.h"
47 
48 namespace oomph
49 {
50 
51 //=========start of TetgenMesh class======================================
52 /// \short Unstructured tet mesh based on output from Tetgen:
53 /// http://wias-berlin.de/software/tetgen//
54 //========================================================================
55 template <class ELEMENT>
56 class TetgenMesh : public virtual TetMeshBase
57 {
58 
59 public:
60 
61  /// \short Empty constructor
63  {
64  // Mesh can only be built with 3D Telements.
65  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
66  }
67 
68  /// \short Constructor with the input files
69  TetgenMesh(const std::string& node_file_name,
70  const std::string& element_file_name,
71  const std::string& face_file_name,
72  TimeStepper* time_stepper_pt=
73  &Mesh::Default_TimeStepper,
74  const bool &use_attributes=false)
75  : Tetgenio_exists(false),
76  Tetgenio_pt(0)
77  {
78  // Mesh can only be built with 3D Telements.
79  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
80 
81  //Store the attributes
82  Use_attributes = use_attributes;
83 
84  // Store timestepper used to build elements
85  Time_stepper_pt = time_stepper_pt;
86 
87  // Build scaffold
88  Tmp_mesh_pt= new
89  TetgenScaffoldMesh(node_file_name,element_file_name,face_file_name);
90 
91  // Convert mesh from scaffold to actual mesh
92  build_from_scaffold(time_stepper_pt,use_attributes);
93 
94  // Kill the scaffold
95  delete Tmp_mesh_pt;
96  Tmp_mesh_pt=0;
97  }
98 
99 
100  /// \short Constructor with tetgenio data structure
101  TetgenMesh(tetgenio& tetgen_data,
102  TimeStepper* time_stepper_pt=
103  &Mesh::Default_TimeStepper,
104  const bool &use_attributes=false)
105 
106  {
107  // Mesh can only be built with 3D Telements.
108  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
109 
110  //Store the attributes
111  Use_attributes = use_attributes;
112 
113  // Store timestepper used to build elements
114  Time_stepper_pt = time_stepper_pt;
115 
116  //We do not have a tetgenio representation
117  Tetgenio_exists = false;
118  Tetgenio_pt = 0;
119 
120  // Build scaffold
121  Tmp_mesh_pt= new
122  TetgenScaffoldMesh(tetgen_data);
123 
124  // Convert mesh from scaffold to actual mesh
125  build_from_scaffold(time_stepper_pt,use_attributes);
126 
127  // Kill the scaffold
128  delete Tmp_mesh_pt;
129  Tmp_mesh_pt=0;
130  }
131 
132 
133  /// \short Constructor with the input files. Setting the boolean
134  /// flag to true splits "corner" elements, i.e. elements that
135  /// that have at least three faces on a domain boundary. The
136  /// relevant elements are split without introducing hanging
137  /// nodes so the sons have a "worse" shape than their fathers.
138  /// However, this step avoids otherwise-hard-to-diagnose
139  /// problems in fluids problems where the application of
140  /// boundary conditions at such "corner" elements can
141  /// overconstrain the solution.
142  TetgenMesh(const std::string& node_file_name,
143  const std::string& element_file_name,
144  const std::string& face_file_name,
145  const bool& split_corner_elements,
146  TimeStepper* time_stepper_pt=
147  &Mesh::Default_TimeStepper,
148  const bool &use_attributes=false) : Outer_boundary_pt(0)
149 
150  {
151  // Mesh can only be built with 3D Telements.
152  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
153 
154  //Throw an error if you try to split corner elements
155  //and use attributes
156  if(split_corner_elements && use_attributes)
157  {
158  std::ostringstream error_stream;
159  error_stream <<
160  "Using region attributes and split_corner_elements simultaneously\n"
161  <<
162  "is not guaranteed to work. The elements adjacent to boundaries\n"
163  <<
164  "accessed by region will not be set up correctly.\n"
165  <<
166  "\n Please fix this!\n\n";
167 
168  throw OomphLibError(error_stream.str(),
169  OOMPH_CURRENT_FUNCTION,
170  OOMPH_EXCEPTION_LOCATION);
171  }
172 
173  //Store the attributes
174  Use_attributes = use_attributes;
175 
176  // Store timestepper used to build elements
177  Time_stepper_pt = time_stepper_pt;
178 
179  // We do not have a tetgenio representation
180  this->Tetgenio_exists = false;;
181  this->Tetgenio_pt = 0;
182 
183  // Build scaffold
184  Tmp_mesh_pt= new
185  TetgenScaffoldMesh(node_file_name,element_file_name,face_file_name);
186 
187  // Convert mesh from scaffold to actual mesh
188  build_from_scaffold(time_stepper_pt,use_attributes);
189 
190  // Kill the scaffold
191  delete Tmp_mesh_pt;
192  Tmp_mesh_pt=0;
193 
194  // Split corner elements
195  if (split_corner_elements)
196  {
198  }
199 
200  // Setup boundary coordinates
201  unsigned nb=nboundary();
202  for (unsigned b=0;b<nb;b++)
203  {
204  bool switch_normal=false;
205  setup_boundary_coordinates(b,switch_normal);
206  }
207  }
208 
209  /// \short Constructor with tetgen data structure Setting the boolean
210  /// flag to true splits "corner" elements, i.e. elements that
211  /// that have at least three faces on a domain boundary. The
212  /// relevant elements are split without introducing hanging
213  /// nodes so the sons have a "worse" shape than their fathers.
214  /// However, this step avoids otherwise-hard-to-diagnose
215  /// problems in fluids problems where the application of
216  /// boundary conditions at such "corner" elements can
217  /// overconstrain the solution.
218  TetgenMesh(tetgenio &tetgen_data,
219  const bool& split_corner_elements,
220  TimeStepper* time_stepper_pt=
221  &Mesh::Default_TimeStepper,
222  const bool &use_attributes=false) : Outer_boundary_pt(0)
223 
224  {
225  // Mesh can only be built with 3D Telements.
226  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
227 
228  //Throw an error if you try to split corner elements
229  //and use attributes
230  if(split_corner_elements && use_attributes)
231  {
232  std::ostringstream error_stream;
233  error_stream <<
234  "Using region attributes and split_corner_elements simultaneously\n"
235  <<
236  "is not guaranteed to work. The elements adjacent to boundaries\n"
237  <<
238  "accessed by region will not be set up correctly.\n"
239  <<
240  "\n Please fix this!\n\n";
241 
242  throw OomphLibError(error_stream.str(),
243  OOMPH_CURRENT_FUNCTION,
244  OOMPH_EXCEPTION_LOCATION);
245  }
246 
247  //Store the attributes
248  Use_attributes = use_attributes;
249 
250  // Store timestepper used to build elements
251  Time_stepper_pt = time_stepper_pt;
252 
253  // We do not have a tetgenio representation
254  this->Tetgenio_exists = false;;
255  this->Tetgenio_pt = 0;
256 
257  // Build scaffold
258  Tmp_mesh_pt= new TetgenScaffoldMesh(tetgen_data);
259 
260  // Convert mesh from scaffold to actual mesh
261  build_from_scaffold(time_stepper_pt,use_attributes);
262 
263  // Kill the scaffold
264  delete Tmp_mesh_pt;
265  Tmp_mesh_pt=0;
266 
267  // Split corner elements
268  if (split_corner_elements)
269  {
271  }
272  }
273 
274 
275  /// \short Build mesh, based on a TetgenMeshClosedSurface that specifies
276  /// the outer boundary of the domain and any number of internal
277  /// closed curves, also specified by TriangleMeshClosedSurfaces.
278  /// Also specify target area for uniform element size.
279  TetgenMesh(TetgenMeshFacetedSurface* const &outer_boundary_pt,
280  Vector<TetgenMeshFacetedSurface*>&
281  internal_closed_surface_pt,
282  const double &element_volume,
283  TimeStepper* time_stepper_pt=
284  &Mesh::Default_TimeStepper,
285  const bool &use_attributes=false)
286  {
287  // Mesh can only be built with 3D Telements.
288  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(3);
289 
290  //Store the attributes
291  Use_attributes = use_attributes;
292 
293  // Store timestepper used to build elements
294  Time_stepper_pt = time_stepper_pt;
295 
296  //Store the boundary
297  Outer_boundary_pt = outer_boundary_pt;
298 
299  //Store the internal boundary
300  Internal_surface_pt = internal_closed_surface_pt;
301 
302  //Tetgen data structure for the input and output
303  tetgenio in;
304 
305  this->build_tetgenio(outer_boundary_pt,
306  internal_closed_surface_pt,
307  in);
308 
309  //Now tetrahedralise
310  std::stringstream input_string;
311  input_string << "pAa" << element_volume;
312 
313  //If any of the boundaries should not be split add the "Y" flag
314  bool can_boundaries_be_split = outer_boundary_pt->can_boundaries_be_split();
315  {
316  unsigned n_internal=internal_closed_surface_pt.size();
317  for(unsigned i=0;i<n_internal;i++)
318  {
319  can_boundaries_be_split &=
320  internal_closed_surface_pt[i]->can_boundaries_be_split();
321  }
322  }
323  //If we can't split the boundaries add the flag
324  if(can_boundaries_be_split==false) {input_string << "Y";}
325 
326  //Now convert to a C-style string
327  char tetswitches[100];
328  sprintf(tetswitches,"%s",input_string.str().c_str());
329 
330  //Make a new tetgen representation
331  this->Tetgenio_exists = true;
332  Tetgenio_pt = new tetgenio;
333 
334  tetrahedralize(tetswitches, &in, this->Tetgenio_pt);
335 
336  // Build scaffold
337  Tmp_mesh_pt= new TetgenScaffoldMesh(*this->Tetgenio_pt);
338 
339  //If any of the objects are different regions then we need to use
340  //the atributes
341  bool regions_exist = false;
342  {
343  unsigned n_internal=internal_closed_surface_pt.size();
344  for(unsigned i=0;i<n_internal;i++)
345  {
346  regions_exist |=
347  internal_closed_surface_pt[i]->is_region();
348  }
349  }
350  //If there are regions, use the attributes
351  if(regions_exist==true) {Use_attributes=true;}
352 
353  // Convert mesh from scaffold to actual mesh
354  build_from_scaffold(time_stepper_pt,Use_attributes);
355 
356  // Kill the scaffold
357  delete Tmp_mesh_pt;
358  Tmp_mesh_pt=0;
359  }
360 
361  ///\short Build tetgenio object from the TetgenMeshClosedSurfaces
362  void build_tetgenio(TetgenMeshFacetedSurface* const &outer_boundary_pt,
363  Vector<TetgenMeshFacetedSurface*>
364  &internal_closed_surface_pt,
365  tetgenio& tetgen_io)
366  {
367  //Pointer to Tetgen facet
368  tetgenio::facet *f;
369  //Pointer to Tetgen polygon
370  tetgenio::polygon *p;
371 
372  //Start all indices from 1 (it's a choice and we've made it
373  tetgen_io.firstnumber = 1;
374  ///ALH: This may not be needed
375  tetgen_io.useindex = true;
376 
377  //Find the number of internal surfaces
378  const unsigned n_internal = internal_closed_surface_pt.size();
379 
380  //Find the number of points on the outer boundary
381  const unsigned n_outer_vertex = outer_boundary_pt->nvertex();
382  tetgen_io.numberofpoints = n_outer_vertex;
383 
384  //Find the number of points on the inner boundaries and add to the totals
385  Vector<unsigned> n_internal_vertex(n_internal);
386  Vector<unsigned> internal_vertex_offset(n_internal);
387  for(unsigned h=0;h<n_internal;++h)
388  {
389  n_internal_vertex[h] = internal_closed_surface_pt[h]->nvertex();
390  internal_vertex_offset[h] = tetgen_io.numberofpoints;
391  tetgen_io.numberofpoints += n_internal_vertex[h];
392  }
393 
394  //Read the data into the point list
395  tetgen_io.pointlist = new double[tetgen_io.numberofpoints * 3];
396  tetgen_io.pointmarkerlist = new int[tetgen_io.numberofpoints];
397  unsigned counter=0;
398  for(unsigned n=0;n<n_outer_vertex;++n)
399  {
400  for(unsigned i=0;i<3;++i)
401  {
402  tetgen_io.pointlist[counter] =
403  outer_boundary_pt->vertex_coordinate(n,i);
404  ++counter;
405  }
406  }
407  for(unsigned h=0;h<n_internal;++h)
408  {
409  const unsigned n_inner = n_internal_vertex[h];
410  for(unsigned n=0;n<n_inner;++n)
411  {
412  for(unsigned i=0;i<3;++i)
413  {
414  tetgen_io.pointlist[counter] =
415  internal_closed_surface_pt[h]->vertex_coordinate(n,i);
416  ++counter;
417  }
418  }
419  }
420 
421 
422  //Set up the pointmarkers
423  counter=0;
424  for(unsigned n=0;n<n_outer_vertex;++n)
425  {
426  tetgen_io.pointmarkerlist[counter] =
427  outer_boundary_pt->vertex_boundary_id(n);
428  ++counter;
429  }
430  for(unsigned h=0;h<n_internal;++h)
431  {
432  const unsigned n_inner = n_internal_vertex[h];
433  for(unsigned n=0;n<n_inner;++n)
434  {
435  tetgen_io.pointmarkerlist[counter] =
436  internal_closed_surface_pt[h]->vertex_boundary_id(n);
437  ++counter;
438  }
439  }
440 
441 
442  //Now the facets
443  unsigned n_outer_facet = outer_boundary_pt->nfacet();
444  tetgen_io.numberoffacets = n_outer_facet;
445  Vector<unsigned> n_inner_facet(n_internal);
446  for(unsigned h=0;h<n_internal;++h)
447  {
448  n_inner_facet[h] = internal_closed_surface_pt[h]->nfacet();
449  tetgen_io.numberoffacets += n_inner_facet[h];
450  }
451 
452  tetgen_io.facetlist = new tetgenio::facet[tetgen_io.numberoffacets];
453  tetgen_io.facetmarkerlist = new int[tetgen_io.numberoffacets];
454 
455 
456  counter=0;
457  for(unsigned n=0;n<n_outer_facet;++n)
458  {
459  //Set pointer to facet
460  f = &tetgen_io.facetlist[counter];
461  f->numberofpolygons = 1;
462  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
463  f->numberofholes = 0;
464  f->holelist = NULL;
465  p = &f->polygonlist[0];
466 
467  Vector<unsigned> facet = outer_boundary_pt->facet(n);
468 
469  p->numberofvertices = facet.size();
470  p->vertexlist = new int[p->numberofvertices];
471  for(int i=0;i<p->numberofvertices;++i)
472  {
473  //The offset here is because we have insisted on 1-based indexing
474  p->vertexlist[i] = facet[i] + 1;
475  }
476 
477  //Set up the boundary markers
478  tetgen_io.facetmarkerlist[counter] =
479  outer_boundary_pt->facet_boundary_id(n);
480  //Increase the counter
481  ++counter;
482  }
483 
484  //Initialise the number of holes
485  tetgen_io.numberofholes=0;
486  //and the number of regions
487  tetgen_io.numberofregions=0;
488 
489  //Loop over the internal stuff
490  for(unsigned h=0;h<n_internal;++h)
491  {
492  for(unsigned n=0;n<n_inner_facet[h];++n)
493  {
494  //Set pointer to facet
495  f = &tetgen_io.facetlist[counter];
496  f->numberofpolygons = 1;
497  f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
498  f->numberofholes = 0;
499  f->holelist = NULL;
500  p = &f->polygonlist[0];
501 
502  Vector<unsigned> facet = internal_closed_surface_pt[h]->facet(n);
503 
504  p->numberofvertices = facet.size();
505  p->vertexlist = new int[p->numberofvertices];
506  for(int i=0;i<p->numberofvertices;++i)
507  {
508  //Add the 1-based and vertex offsets to get these number correct
509  p->vertexlist[i] = facet[i] + internal_vertex_offset[h] + 1;
510  }
511  //Set up the boundary markers
512  tetgen_io.facetmarkerlist[counter] =
513  internal_closed_surface_pt[h]->facet_boundary_id(n);
514  ++counter;
515  }
516 
517  //If it's a hole add it
518  if(internal_closed_surface_pt[h]->is_hole())
519  {
520  ++tetgen_io.numberofholes;
521  }
522  //Otherwise it may be region
523  else
524  {
525  if(internal_closed_surface_pt[h]->is_region())
526  {
527  ++tetgen_io.numberofregions;
528  }
529  }
530  }
531 
532  //Set storage for the holes
533  tetgen_io.holelist = new double[3*tetgen_io.numberofholes];
534  //Loop over all the internal boundaries again
535  counter=0;
536  for(unsigned h=0;h<n_internal;++h)
537  {
538  if(internal_closed_surface_pt[h]->is_hole())
539  {
540  for(unsigned i=0;i<3;++i)
541  {
542  tetgen_io.holelist[counter] =
543  internal_closed_surface_pt[h]->internal_point(i);
544  ++counter;
545  }
546  }
547  }
548 
549  //Set storage for the regions
550  tetgen_io.regionlist = new double[5*tetgen_io.numberofregions];
551  //Loop over all the internal boundaries again
552  counter=0;
553  for(unsigned h=0;h<n_internal;++h)
554  {
555  if(internal_closed_surface_pt[h]->is_region())
556  {
557  for(unsigned i=0;i<3;++i)
558  {
559  tetgen_io.regionlist[counter] =
560  internal_closed_surface_pt[h]->internal_point(i);
561  ++counter;
562  }
563  tetgen_io.regionlist[counter] =
564  static_cast<double>(internal_closed_surface_pt[h]->region());
565  ++counter;
566  //Area target
567  tetgen_io.regionlist[counter] = 0.0;
568  ++counter;
569  }
570  }
571  }
572 
573 
574  /// Empty destructor
576  {
577  if(Tetgenio_exists)
578  {
579  delete Tetgenio_pt;
580  }
581  }
582 
583 
584  /// \short Overload set_mesh_level_time_stepper so that the stored
585  /// time stepper now corresponds to the new timestepper
586  void set_mesh_level_time_stepper(TimeStepper* const &time_stepper_pt,
587  const bool &preserve_existing_data)
588  {this->Time_stepper_pt = time_stepper_pt;}
589 
590 
591 
592  /// Boolen defining whether tetgenio object has been built or not
593  bool tetgenio_exists() const {return Tetgenio_exists;}
594 
595 
596  /// Access to the triangulateio representation of the mesh
597  tetgenio* &tetgenio_pt() {return Tetgenio_pt;}
598 
599  /// Set the tetgen pointer by a deep copy
600  void set_deep_copy_tetgenio_pt(tetgenio* const &tetgenio_pt)
601  {
602  //Delete the existing data
603  if(Tetgenio_exists) {delete Tetgenio_pt;}
604  this->Tetgenio_pt= new tetgenio;
605  //Create a deep copy of tetgenio_pt and store the result in
606  //Tetgenio_pt
607  this->deep_copy_of_tetgenio(tetgenio_pt,this->Tetgenio_pt);
608  }
609 
610  /// Transfer tetgenio data from the input to the output
611  /// The output is assumed to have been constructed and "empty"
612  void deep_copy_of_tetgenio(tetgenio* const &input_pt,
613  tetgenio* &output_pt);
614 
615  /// \short Setup boundary coordinate on boundary b which is
616  /// assumed to be planar. Boundary coordinates are the
617  /// x-y coordinates in the plane of that boundary with the
618  /// x-axis along the line from the (lexicographically)
619  /// "lower left" to the "upper right" node. The y axis
620  /// is obtained by taking the cross-product of the positive
621  /// x direction with the outer unit normal computed by
622  /// the face elements.
623  void setup_boundary_coordinates(const unsigned& b)
624  {
625  // Dummy file
626  std::ofstream some_file;
627 
628  // Don't switch the normal
629  bool switch_normal=false;
630 
631  this->setup_boundary_coordinates_generic(b,switch_normal,some_file);
632  }
633 
634  /// \short Setup boundary coordinate on boundary b which is
635  /// assumed to be planar. Boundary coordinates are the
636  /// x-y coordinates in the plane of that boundary with the
637  /// x-axis along the line from the (lexicographically)
638  /// "lower left" to the "upper right" node. The y axis
639  /// is obtained by taking the cross-product of the positive
640  /// x direction with the outer unit normal computed by
641  /// the face elements.
642  void setup_boundary_coordinates(const unsigned& b, const bool& switch_normal)
643  {
644  // Dummy file
645  std::ofstream some_file;
646 
647  this->setup_boundary_coordinates_generic(b,switch_normal,some_file);
648  }
649 
650 
651  /// \short Setup boundary coordinate on boundary b which is
652  /// assumed to be planar. Boundary coordinates are the
653  /// x-y coordinates in the plane of that boundary with the
654  /// x-axis along the line from the (lexicographically)
655  /// "lower left" to the "upper right" node. The y axis
656  /// is obtained by taking the cross-product of the positive
657  /// x direction with the outer unit normal computed by
658  /// the face elements. Doc faces in output file.
659  void setup_boundary_coordinates(const unsigned& b,
660  std::ofstream& outfile)
661  {
662  // Don't switch the normal
663  bool switch_normal=false;
664 
665  this->setup_boundary_coordinates_generic(b,switch_normal,outfile);
666  }
667 
668  /// \short Setup boundary coordinate on boundary b which is
669  /// assumed to be planar. Boundary coordinates are the
670  /// x-y coordinates in the plane of that boundary with the
671  /// x-axis along the line from the (lexicographically)
672  /// "lower left" to the "upper right" node. The y axis
673  /// is obtained by taking the cross-product of the positive
674  /// x direction with the outer unit normal computed by
675  /// the face elements (or its negative if switch_normal is set
676  /// to true). Doc faces in output file.
677  void setup_boundary_coordinates(const unsigned& b,
678  const bool& switch_normal,
679  std::ofstream& outfile)
680  {
681  this->setup_boundary_coordinates_generic(b,switch_normal,outfile);
682  }
683 
684  /// \short Setup boundary coordinate on boundary b which is
685  /// assumed to be planar. Boundary coordinates are the
686  /// x-y coordinates in the plane of that boundary with the
687  /// x-axis along the line from the (lexicographically)
688  /// "lower left" to the "upper right" node. The y axis
689  /// is obtained by taking the cross-product of the positive
690  /// x direction with the outer unit normal computed by
691  /// the face elements (or its negative if switch_normal is set
692  /// to true). Doc faces in output file.
693  virtual void setup_boundary_coordinates_generic(const unsigned& b,
694  const bool& switch_normal,
695  std::ofstream& outfile);
696 
697 
698  /// Return the number of elements adjacent to boundary b in region r
699  inline unsigned nboundary_element_in_region(const unsigned &b,
700  const unsigned &r) const
701  {
702  //Need to use a constant iterator here to keep the function "const"
703  //Return an iterator to the appropriate entry, if we find it
704  std::map<unsigned,Vector<FiniteElement*> >::const_iterator it =
705  Boundary_region_element_pt[b].find(r);
706  if(it!=Boundary_region_element_pt[b].end())
707  {
708  return (it->second).size();
709  }
710  //Otherwise there are no elements adjacent to boundary b in the region r
711  else
712  {
713  return 0;
714  }
715  }
716 
717  /// Return pointer to the e-th element adjacent to boundary b in region r
718  FiniteElement* boundary_element_in_region_pt(const unsigned &b,
719  const unsigned &r,
720  const unsigned &e) const
721  {
722  //Use a constant iterator here to keep function "const" overall
723  std::map<unsigned,Vector<FiniteElement*> >::const_iterator it =
724  Boundary_region_element_pt[b].find(r);
725  if(it!=Boundary_region_element_pt[b].end())
726  {
727  return (it->second)[e];
728  }
729  else
730  {
731  return 0;
732  }
733  }
734 
735  /// Return face index of the e-th element adjacent to boundary b in region r
736  int face_index_at_boundary_in_region(const unsigned &b,
737  const unsigned &r,
738  const unsigned &e) const
739  {
740  //Use a constant iterator here to keep function "const" overall
741  std::map<unsigned,Vector<int> >::const_iterator it =
743  if(it!=Face_index_region_at_boundary[b].end())
744  {
745  return (it->second)[e];
746  }
747  else
748  {
749  std::ostringstream error_message;
750  error_message << "Face indices not set up for boundary "
751  << b << " in region " << r << "\n";
752  error_message
753  << "This probably means that the boundary is not adjacent to region\n";
754  throw OomphLibError(error_message.str(),
755  OOMPH_CURRENT_FUNCTION,
756  OOMPH_EXCEPTION_LOCATION);
757  }
758  }
759 
760 
761  /// Return the number of regions specified by attributes
762  unsigned nregion() {return Region_element_pt.size();}
763 
764  /// Return the number of elements in region i
765  unsigned nregion_element(const unsigned &i)
766  {return Region_element_pt[i].size();}
767 
768  /// Return the attribute associated with region i
769  double region_attribute(const unsigned &i)
770  {return Region_attribute[i];}
771 
772  /// Return the e-th element in the i-th region
773  FiniteElement* region_element_pt(const unsigned &i,
774  const unsigned &e)
775  {return Region_element_pt[i][e];}
776 
777  /// \short Snap boundaries specified by the IDs listed in boundary_id to
778  /// a quadratric surface, specified in the file
779  /// quadratic_surface_file_name. This is usually used with vmtk-based
780  /// meshes for which oomph-lib's xda to poly conversion code produces the files
781  /// "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat"
782  /// which specify the quadratic FSI boundary (for the fluid and the solid)
783  /// and the quadratic representation of the outer boundary of the solid.
784  /// When used with these files, the flag switch_normal should be
785  /// set to true when calling the function for the outer boundary of the
786  /// solid. The DocInfo object can be used to label optional output
787  /// files. (Uses directory and label).
788  void snap_to_quadratic_surface(const Vector<unsigned>& boundary_id,
789  const std::string& quadratic_surface_file_name,
790  const bool& switch_normal,
791  DocInfo& doc_info);
792 
793  /// \short Snap boundaries specified by the IDs listed in boundary_id to
794  /// a quadratric surface, specified in the file
795  /// quadratic_surface_file_name. This is usually used with vmtk-based
796  /// meshes for which oomph-lib's xda to poly conversion code produces the files
797  /// "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat"
798  /// which specify the quadratic FSI boundary (for the fluid and the solid)
799  /// and the quadratic representation of the outer boundary of the solid.
800  /// When used with these files, the flag switch_normal should be
801  /// set to true when calling the function for the outer boundary of the
802  /// solid.
803  void snap_to_quadratic_surface(const Vector<unsigned>& boundary_id,
804  const std::string& quadratic_surface_file_name,
805  const bool& switch_normal)
806  {
807  // Dummy doc info
808  DocInfo doc_info;
809  doc_info.disable_doc();
810  snap_to_quadratic_surface(boundary_id,quadratic_surface_file_name,
811  switch_normal,doc_info);
812 
813  }
814 
815 
816  /// \short Non-Delaunay split elements that have three faces on a boundary
817  /// into sons. Timestepper species timestepper for new nodes; defaults
818  /// to to steady timestepper.
819  void split_elements_in_corners(TimeStepper* time_stepper_pt=
820  &Mesh::Default_TimeStepper);
821 
822 
823  protected:
824 
825  /// Build mesh from scaffold
826  void build_from_scaffold(TimeStepper* time_stepper_pt,
827  const bool &use_attributes);
828 
829 
830  protected:
831 
832  /// \short Boolean to indicate whether a tetgenio representation of the
833  /// mesh exists
835 
836  /// \short Tetgen representation of mesh
837  tetgenio* Tetgenio_pt;
838 
839  /// Temporary scaffold mesh
840  TetgenScaffoldMesh* Tmp_mesh_pt;
841 
842  /// Vectors of elements in each region differentiated by attribute
843  Vector<Vector<FiniteElement*> > Region_element_pt;
844 
845  /// Vector of attributes associated with the elements in each region
846  Vector<double> Region_attribute;
847 
848  /// Storage for elements adjacent to a boundary in a particular region
849  Vector<std::map<unsigned,Vector<FiniteElement*> > >
851 
852  /// \short Storage for the face index adjacent to a boundary in a particular
853  /// region
854  Vector<std::map<unsigned,Vector<int> > > Face_index_region_at_boundary;
855 
856  /// Faceted surface that defines outer boundaries
857  TetgenMeshFacetedSurface* Outer_boundary_pt;
858 
859  /// Vector to faceted surfaces that define internal boundaries
860  Vector<TetgenMeshFacetedSurface*> Internal_surface_pt;
861 
862 
863  /// Timestepper used to build elements
864  TimeStepper* Time_stepper_pt;
865 
866  /// Boolean flag to indicate whether to use attributes or not
867  /// (required for multidomain meshes)
869 
870 }; //end class
871 
872 
873 
874 
875 //////////////////////////////////////////////////////////////////
876 //////////////////////////////////////////////////////////////////
877 //////////////////////////////////////////////////////////////////
878 
879 
880 //==============start_mesh=================================================
881 /// Tetgen-based mesh upgraded to become a solid mesh. Automatically
882 /// enumerates all boundaries
883 //=========================================================================
884 template<class ELEMENT>
885 class SolidTetMesh : public virtual TetgenMesh<ELEMENT>,
886  public virtual SolidMesh
887 {
888 
889 public:
890 
891  /// \short Constructor. Boundary coordinates are setup
892  /// automatically.
893  SolidTetMesh(const std::string& node_file_name,
894  const std::string& element_file_name,
895  const std::string& face_file_name,
896  const bool& split_corner_elements,
897  TimeStepper* time_stepper_pt=
898  &Mesh::Default_TimeStepper,
899  const bool &use_attributes=false) :
900  TetgenMesh<ELEMENT>(node_file_name, element_file_name,
901  face_file_name, split_corner_elements,
902  time_stepper_pt, use_attributes)
903  {
904  //Assign the Lagrangian coordinates
905  set_lagrangian_nodal_coordinates();
906 
907  // Find out elements next to boundary
908  setup_boundary_element_info();
909 
910  // Setup boundary coordinates for all boundaries without switching
911  // direction of normal
912  bool switch_normal=false;
913  unsigned nb=this->nboundary();
914  for (unsigned b=0;b<nb;b++)
915  {
916  this->setup_boundary_coordinates(b,switch_normal);
917  }
918  }
919 
920  /// \short Constructor. Boundary coordinates are setup
921  /// automatically, with the orientation of the outer unit
922  /// normal determined by switch_normal.
923  SolidTetMesh(const std::string& node_file_name,
924  const std::string& element_file_name,
925  const std::string& face_file_name,
926  const bool& split_corner_elements,
927  const bool& switch_normal,
928  TimeStepper* time_stepper_pt=
929  &Mesh::Default_TimeStepper,
930  const bool &use_attributes=false) :
931  TetgenMesh<ELEMENT>(node_file_name, element_file_name,
932  face_file_name, split_corner_elements,
933  time_stepper_pt, use_attributes)
934  {
935  //Assign the Lagrangian coordinates
936  set_lagrangian_nodal_coordinates();
937 
938  // Find out elements next to boundary
939  setup_boundary_element_info();
940 
941  // Setup boundary coordinates for all boundaries
942  unsigned nb=this->nboundary();
943  for (unsigned b=0;b<nb;b++)
944  {
945  this->setup_boundary_coordinates(b,switch_normal);
946  }
947  }
948 
949  /// Empty Destructor
950  virtual ~SolidTetMesh() { }
951 
952 };
953 
954 
955 
956 
957 
958 }
959 
960 #endif
void setup_boundary_coordinates(const unsigned &b, const bool &switch_normal, std::ofstream &outfile)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
void split_elements_in_corners(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Non-Delaunay split elements that have three faces on a boundary into sons. Timestepper species timest...
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region.
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
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...
TetgenMesh()
Empty constructor.
TetgenMesh(TetgenMeshFacetedSurface *const &outer_boundary_pt, Vector< TetgenMeshFacetedSurface * > &internal_closed_surface_pt, const double &element_volume, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
Build mesh, based on a TetgenMeshClosedSurface that specifies the outer boundary of the domain and an...
FiniteElement * boundary_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.
TetgenScaffoldMesh * Tmp_mesh_pt
Temporary scaffold mesh.
SolidTetMesh(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 setup automatically, with the orientation of the outer unit nor...
virtual ~SolidTetMesh()
Empty Destructor.
void snap_to_quadratic_surface(const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal, DocInfo &doc_info)
Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface, specified in the file quadratic_surface_file_name. This is usually used with vmtk-based meshes for which oomph-lib's xda to poly conversion code produces the files "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI boundary (for the fluid and the solid) and the quadratic representation of the outer boundary of the solid. When used with these files, the flag switch_normal should be set to true when calling the function for the outer boundary of the solid. The DocInfo object can be used to label optional output files. (Uses directory and label).
Vector< TetgenMeshFacetedSurface * > Internal_surface_pt
Vector to faceted surfaces that define internal boundaries.
TimeStepper * Time_stepper_pt
Timestepper used to build elements.
unsigned nregion()
Return the number of regions specified by attributes.
bool tetgenio_exists() const
Boolen defining whether tetgenio object has been built or not.
FiniteElement * region_element_pt(const unsigned &i, const unsigned &e)
Return the e-th element in the i-th region.
~TetgenMesh()
Empty destructor.
TetgenMesh(tetgenio &tetgen_data, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false)
Constructor with tetgenio data structure.
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
bool Tetgenio_exists
Boolean to indicate whether a tetgenio representation of the mesh exists.
void setup_boundary_coordinates(const unsigned &b, std::ofstream &outfile)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
Unstructured tet mesh based on output from Tetgen: http://wias-berlin.de/software/tetgen//.
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...
TetgenMeshFacetedSurface * Outer_boundary_pt
Faceted surface that defines outer boundaries.
void build_tetgenio(TetgenMeshFacetedSurface *const &outer_boundary_pt, Vector< TetgenMeshFacetedSurface * > &internal_closed_surface_pt, tetgenio &tetgen_io)
Build tetgenio object from the TetgenMeshClosedSurfaces.
void deep_copy_of_tetgenio(tetgenio *const &input_pt, tetgenio *&output_pt)
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...
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
tetgenio *& tetgenio_pt()
Access to the triangulateio representation of the mesh.
void setup_boundary_coordinates(const unsigned &b)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
unsigned nboundary_element_in_region(const unsigned &b, const unsigned &r) const
Return the number of elements adjacent to boundary b in region r.
double region_attribute(const unsigned &i)
Return the attribute associated with region i.
void setup_boundary_coordinates(const unsigned &b, const bool &switch_normal)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
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.
unsigned nregion_element(const unsigned &i)
Return the number of elements in region i.
void set_deep_copy_tetgenio_pt(tetgenio *const &tetgenio_pt)
Set the tetgen pointer by a deep copy.
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.
virtual void setup_boundary_coordinates_generic(const unsigned &b, const bool &switch_normal, std::ofstream &outfile)
Setup boundary coordinate on boundary b which is assumed to be planar. Boundary coordinates are the x...
void snap_to_quadratic_surface(const Vector< unsigned > &boundary_id, const std::string &quadratic_surface_file_name, const bool &switch_normal)
Snap boundaries specified by the IDs listed in boundary_id to a quadratric surface, specified in the file quadratic_surface_file_name. This is usually used with vmtk-based meshes for which oomph-lib's xda to poly conversion code produces the files "quadratic_fsi_boundary.dat" and "quadratic_outer_solid_boundary.dat" which specify the quadratic FSI boundary (for the fluid and the solid) and the quadratic representation of the outer boundary of the solid. When used with these files, the flag switch_normal should be set to true when calling the function for the outer boundary of the solid.
SolidTetMesh(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.
tetgenio * Tetgenio_pt
Tetgen representation of mesh.
Vector< Vector< FiniteElement * > > Region_element_pt
Vectors of elements in each region differentiated by attribute.