libMesh
Functions
miscellaneous_ex6.C File Reference

Go to the source code of this file.

Functions

void triangulate_domain (const Parallel::Communicator &comm)
 
void tetrahedralize_domain (const Parallel::Communicator &comm)
 
void add_cube_convex_hull_to_mesh (MeshBase &mesh, Point lower_limit, Point upper_limit)
 
int main (int argc, char **argv)
 

Function Documentation

void add_cube_convex_hull_to_mesh ( MeshBase mesh,
Point  lower_limit,
Point  upper_limit 
)

Definition at line 249 of file miscellaneous_ex6.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::MeshTools::Generation::build_cube(), libMesh::ParallelObject::comm(), libMesh::HEX8, libMesh::DofObject::id(), libMesh::libmesh_ignore(), libmesh_nullptr, libMesh::MeshTools::n_elem(), libMesh::MeshBase::node_ptr(), libMesh::TetGenMeshInterface::pointset_convexhull(), libMesh::Elem::set_node(), side, and libMesh::TRI3.

Referenced by tetrahedralize_domain().

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);
269  t.pointset_convexhull();
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...
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
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
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
const class libmesh_nullptr_t libmesh_nullptr
virtual const Node * node_ptr(const dof_id_type i) const =0
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.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
void libmesh_ignore(const T &)
const Parallel::Communicator & comm() const
dof_id_type id() const
Definition: dof_object.h:632
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.
int main ( int  argc,
char **  argv 
)

Definition at line 54 of file miscellaneous_ex6.C.

References libMesh::LibMeshInit::comm(), libMesh::TriangleWrapper::init(), libMesh::out, tetrahedralize_domain(), and triangulate_domain().

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 }
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
OStreamProxy out
void tetrahedralize_domain(const Parallel::Communicator &comm)
void triangulate_domain(const Parallel::Communicator &comm)
void tetrahedralize_domain ( const Parallel::Communicator comm)

Definition at line 168 of file miscellaneous_ex6.C.

References add_cube_convex_hull_to_mesh(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::libmesh_ignore(), mesh, libMesh::MeshBase::prepare_for_use(), libMesh::Real, libMesh::TetGenMeshInterface::set_switches(), libMesh::TetGenMeshInterface::triangulate_conformingDelaunayMesh_carvehole(), and libMesh::UnstructuredMesh::write().

Referenced by main().

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.
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
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
222  t.triangulate_conformingDelaunayMesh_carvehole(hole,
223  quality_constraint,
224  volume_constraint);
225 
226  // Find neighbors, etc in preparation for writing out the Mesh
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 }
Class TetGenMeshInterface provides an interface for tetrahedralization of meshes using the TetGen lib...
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
void add_cube_convex_hull_to_mesh(MeshBase &mesh, Point lower_limit, Point upper_limit)
MeshBase & mesh
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true)=0
Locate element face (edge in 2D) neighbors.
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
virtual void write(const std::string &name)=0
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void triangulate_domain ( const Parallel::Communicator comm)

Definition at line 79 of file miscellaneous_ex6.C.

References libMesh::DistributedMesh::add_point(), libMesh::TriangleInterface::attach_hole_list(), libMesh::TriangleInterface::desired_area(), libMesh::libmesh_ignore(), mesh, libMesh::pi, libMesh::TriangleInterface::PSLG, libMesh::Real, libMesh::TriangleInterface::smooth_after_generating(), libMesh::TriangleInterface::triangulate(), libMesh::TriangleInterface::triangulation_type(), and libMesh::UnstructuredMesh::write().

Referenced by main().

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.
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
112  t.triangulation_type() = TriangleInterface::PSLG;
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 }
An abstract class for defining a 2-dimensional hole.
MeshBase & mesh
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.
Another concrete instantiation of the hole, this one should be sufficiently general for most non-poly...
virtual void write(const std::string &name)=0
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
const Real pi
.
Definition: libmesh.h:172
A concrete instantiation of the Hole class that describes polygonal (triangular, square, pentagonal, ...) holes.