libMesh
miscellaneous_ex6.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Miscellaneous Example 6 - Meshing with LibMesh's TetGen and Triangle Interfaces</h1>
21 // \author John W. Peterson
22 // \date 2011
23 //
24 // LibMesh provides interfaces to both Triangle and TetGen for generating
25 // Delaunay triangulations and tetrahedralizations in two and three dimensions
26 // (respectively).
27 
28 // Local header files
29 #include "libmesh/elem.h"
30 #include "libmesh/face_tri3.h"
31 #include "libmesh/mesh.h"
32 #include "libmesh/mesh_generation.h"
33 #include "libmesh/mesh_tetgen_interface.h"
34 #include "libmesh/mesh_triangle_holes.h"
35 #include "libmesh/mesh_triangle_interface.h"
36 #include "libmesh/node.h"
37 #include "libmesh/replicated_mesh.h"
38 
39 // Bring in everything from the libMesh namespace
40 using namespace libMesh;
41 
42 // Major functions called by main
45 
46 // Helper routine for tetrahedralize_domain(). Adds the points and elements
47 // of a convex hull generated by TetGen to the input mesh
48 void add_cube_convex_hull_to_mesh(MeshBase & mesh, Point lower_limit, Point upper_limit);
49 
50 
51 
52 
53 // Begin the main program.
54 int main (int argc, char ** argv)
55 {
56  // Initialize libMesh and any dependent libraries, like in example 2.
57  LibMeshInit init (argc, argv);
58 
59  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
60 
61  libMesh::out << "Triangulating an L-shaped domain with holes" << std::endl;
62 
63  // 1.) 2D triangulation of L-shaped domain with three holes of different shape
64  triangulate_domain(init.comm());
65 
66  libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
67 
68  libMesh::out << "Tetrahedralizing a prismatic domain with a hole" << std::endl;
69 
70  // 2.) 3D tetrahedralization of rectangular domain with hole.
72 
73  return 0;
74 }
75 
76 
77 
78 
80 {
81 #ifdef LIBMESH_HAVE_TRIANGLE
82  // Use typedefs for slightly less typing.
83  typedef TriangleInterface::Hole Hole;
84  typedef TriangleInterface::PolygonHole PolygonHole;
85  typedef TriangleInterface::ArbitraryHole ArbitraryHole;
86 
87  // Libmesh mesh that will eventually be created.
88  Mesh mesh(comm, 2);
89 
90  // The points which make up the L-shape:
91  mesh.add_point(Point(0., 0.));
92  mesh.add_point(Point(0., -1.));
93  mesh.add_point(Point(-1., -1.));
94  mesh.add_point(Point(-1., 1.));
95  mesh.add_point(Point(1., 1.));
96  mesh.add_point(Point(1., 0.));
97 
98  // Declare the TriangleInterface object. This is where
99  // we can set parameters of the triangulation and where the
100  // actual triangulate function lives.
101  TriangleInterface t(mesh);
102 
103  // Customize the variables for the triangulation
104  t.desired_area() = .01;
105 
106  // A Planar Straight Line Graph (PSLG) is essentially a list
107  // of segments which have to exist in the final triangulation.
108  // For an L-shaped domain, Triangle will compute the convex
109  // hull of boundary points if we do not specify the PSLG.
110  // The PSLG algorithm is also required for triangulating domains
111  // containing holes
113 
114  // Turn on/off Laplacian mesh smoothing after generation.
115  // By default this is on.
116  t.smooth_after_generating() = true;
117 
118  // Define holes...
119 
120  // hole_1 is a circle (discretized by 50 points)
121  PolygonHole hole_1(Point(-0.5, 0.5), // center
122  0.25, // radius
123  50); // n. points
124 
125  // hole_2 is itself a triangle
126  PolygonHole hole_2(Point(0.5, 0.5), // center
127  0.1, // radius
128  3); // n. points
129 
130  // hole_3 is an ellipse of 100 points which we define here
131  Point ellipse_center(-0.5, -0.5);
132  const unsigned int n_ellipse_points=100;
133  std::vector<Point> ellipse_points(n_ellipse_points);
134  const Real
135  dtheta = 2*libMesh::pi / static_cast<Real>(n_ellipse_points),
136  a = .1,
137  b = .2;
138 
139  for (unsigned int i=0; i<n_ellipse_points; ++i)
140  ellipse_points[i]= Point(ellipse_center(0)+a*cos(i*dtheta),
141  ellipse_center(1)+b*sin(i*dtheta));
142 
143  ArbitraryHole hole_3(ellipse_center, ellipse_points);
144 
145  // Create the vector of Hole*'s ...
146  std::vector<Hole *> holes;
147  holes.push_back(&hole_1);
148  holes.push_back(&hole_2);
149  holes.push_back(&hole_3);
150 
151  // ... and attach it to the triangulator object
152  t.attach_hole_list(&holes);
153 
154  // Triangulate!
155  t.triangulate();
156 
157  // Write the result to file
158  mesh.write("delaunay_l_shaped_hole.e");
159 
160 #else
161  // Avoid compiler warnings
162  libmesh_ignore(comm);
163 #endif // LIBMESH_HAVE_TRIANGLE
164 }
165 
166 
167 
169 {
170 #ifdef LIBMESH_HAVE_TETGEN
171  // The algorithm is broken up into several steps:
172  // 1.) A convex hull is constructed for a rectangular hole.
173  // 2.) A convex hull is constructed for the domain exterior.
174  // 3.) Neighbor information is updated so TetGen knows there is a convex hull
175  // 4.) A vector of hole points is created.
176  // 5.) The domain is tetrahedralized, the mesh is written out, etc.
177 
178  // The mesh we will eventually generate
179  ReplicatedMesh mesh(comm, 3);
180 
181  // Lower and Upper bounding box limits for a rectangular hole within the unit cube.
182  Point hole_lower_limit(0.2, 0.2, 0.4);
183  Point hole_upper_limit(0.8, 0.8, 0.6);
184 
185  // 1.) Construct a convex hull for the hole
186  add_cube_convex_hull_to_mesh(mesh, hole_lower_limit, hole_upper_limit);
187 
188  // 2.) Generate elements comprising the outer boundary of the domain.
190  Point(0., 0., 0.),
191  Point(1., 1., 1.));
192 
193  // 3.) Update neighbor information so that TetGen can verify there is a convex hull.
194  mesh.find_neighbors();
195 
196  // 4.) Set up vector of hole points
197  std::vector<Point> hole(1);
198  hole[0] = Point(0.5*(hole_lower_limit + hole_upper_limit));
199 
200  // 5.) Set parameters and tetrahedralize the domain
201 
202  // 0 means "use TetGen default value"
203  Real quality_constraint = 2.0;
204 
205  // The volume constraint determines the max-allowed tetrahedral
206  // volume in the Mesh. TetGen will split cells which are larger than
207  // this size
208  Real volume_constraint = 0.001;
209 
210  // Construct the Delaunay tetrahedralization
211  TetGenMeshInterface t(mesh);
212 
213  // In debug mode, set the tetgen switches
214  // -V (verbose) and
215  // -CC (check consistency and constrained Delaunay)
216  // In optimized mode, only switch Q (quiet) is set.
217  // For more options, see tetgen website:
218  // (http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html#cmd-m)
219 #ifdef DEBUG
220  t.set_switches("VCC");
221 #endif
223  quality_constraint,
224  volume_constraint);
225 
226  // Find neighbors, etc in preparation for writing out the Mesh
227  mesh.prepare_for_use();
228 
229  // Finally, write out the result
230  mesh.write("hole_3D.e");
231 #else
232  // Avoid compiler warnings
233  libmesh_ignore(comm);
234 #endif // LIBMESH_HAVE_TETGEN
235 }
236 
237 
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 
248 
250  Point lower_limit,
251  Point upper_limit)
252 {
253 #ifdef LIBMESH_HAVE_TETGEN
254  ReplicatedMesh cube_mesh(mesh.comm(), 3);
255 
256  unsigned n_elem = 1;
257 
259  n_elem, n_elem, n_elem, // n. elements in each direction
260  lower_limit(0), upper_limit(0),
261  lower_limit(1), upper_limit(1),
262  lower_limit(2), upper_limit(2),
263  HEX8);
264 
265  // The pointset_convexhull() algorithm will ignore the Hex8s
266  // in the Mesh, and just construct the triangulation
267  // of the convex hull.
268  TetGenMeshInterface t(cube_mesh);
270 
271  // Now add all nodes from the boundary of the cube_mesh to the input mesh.
272 
273  // Map from "node id in cube_mesh" -> "node id in mesh". Initially inserted
274  // with a dummy value, later to be assigned a value by the input mesh.
275  std::map<unsigned, unsigned> node_id_map;
276  typedef std::map<unsigned, unsigned>::iterator iterator;
277 
278  for (auto & elem : cube_mesh.element_ptr_range())
279  for (auto s : elem->side_index_range())
280  if (elem->neighbor(s) == libmesh_nullptr)
281  {
282  // Add the node IDs of this side to the set
283  UniquePtr<Elem> side = elem->side(s);
284 
285  for (auto n : side->node_index_range())
286  node_id_map.insert(std::make_pair(side->node_id(n), /*dummy_value=*/0));
287  }
288 
289  // For each node in the map, insert it into the input mesh and keep
290  // track of the ID assigned.
291  for (iterator it=node_id_map.begin(); it != node_id_map.end(); ++it)
292  {
293  // Id of the node in the cube mesh
294  unsigned id = (*it).first;
295 
296  // Pointer to node in the cube mesh
297  Node & old_node = cube_mesh.node_ref(id);
298 
299  // Add geometric point to input mesh
300  Node * new_node = mesh.add_point (old_node);
301 
302  // Track ID value of new_node in map
303  (*it).second = new_node->id();
304  }
305 
306  // With the points added and the map data structure in place, we are
307  // ready to add each TRI3 element of the cube_mesh to the input Mesh
308  // with proper node assignments
309  for (auto & old_elem : cube_mesh.element_ptr_range())
310  if (old_elem->type() == TRI3)
311  {
312  Elem * new_elem = mesh.add_elem(new Tri3);
313 
314  // Assign nodes in new elements. Since this is an example,
315  // we'll do it in several steps.
316  for (auto i : old_elem->node_index_range())
317  {
318  // Locate old node ID in the map
319  iterator it = node_id_map.find(old_elem->node_id(i));
320 
321  // Check for not found
322  if (it == node_id_map.end())
323  libmesh_error_msg("Node id " << old_elem->node_id(i) << " not found in map!");
324 
325  // Mapping to node ID in input mesh
326  unsigned new_node_id = (*it).second;
327 
328  // Node pointer assigned from input mesh
329  new_elem->set_node(i) = mesh.node_ptr(new_node_id);
330  }
331  }
332 #else
333  // Avoid compiler warnings
334  libmesh_ignore(mesh);
335  libmesh_ignore(lower_limit);
336  libmesh_ignore(upper_limit);
337 #endif // LIBMESH_HAVE_TETGEN
338 }
Class TetGenMeshInterface provides an interface for tetrahedralization of meshes using the TetGen lib...
The Tri3 is an element in 2D composed of 3 nodes.
Definition: face_tri3.h:54
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
void set_switches(const std::string &)
Method to set switches to tetgen, allowing for different behaviours.
An abstract class for defining a 2-dimensional hole.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
bool & smooth_after_generating()
Sets/gets flag which tells whether to do Delaunay mesh smoothing after generating the grid...
void add_cube_convex_hull_to_mesh(MeshBase &mesh, Point lower_limit, Point upper_limit)
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
void attach_hole_list(const std::vector< Hole * > *holes)
Attaches a vector of Hole* pointers which will be meshed around.
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
virtual const Node * node_ptr(const dof_id_type i) const =0
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true) libmesh_override
Other functions from MeshBase requiring re-definition.
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
Use the triangle library to triangulate a Planar Straight Line Graph which is defined implicitly by t...
Real & desired_area()
Sets and/or gets the desired triangle area.
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void triangulate()
This is the main public interface for this function.
void triangulate_conformingDelaunayMesh_carvehole(const std::vector< Point > &holes, double quality_constraint=0., double volume_constraint=0.)
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...
int main(int argc, char **argv)
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual void write(const std::string &name) libmesh_override
Write the file specified by name.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
Another concrete instantiation of the hole, this one should be sufficiently general for most non-poly...
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
TriangulationType & triangulation_type()
Sets and/or gets the desired triangulation type.
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Parallel::Communicator & comm() const
OStreamProxy out
void tetrahedralize_domain(const Parallel::Communicator &comm)
void pointset_convexhull()
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...
dof_id_type id() const
Definition: dof_object.h:632
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
A C++ interface between LibMesh and the Triangle library written by J.R.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
void triangulate_domain(const Parallel::Communicator &comm)
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id) libmesh_override
functions for adding /deleting nodes elements.
const Real pi
.
Definition: libmesh.h:172
A concrete instantiation of the Hole class that describes polygonal (triangular, square, pentagonal, ...) holes.