libMesh
mesh_triangle_interface.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 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_HAVE_TRIANGLE
22 
23 // C/C++ includes
24 #include <sstream>
25 
26 #include "libmesh/mesh_triangle_interface.h"
27 #include "libmesh/unstructured_mesh.h"
28 #include "libmesh/face_tri3.h"
29 #include "libmesh/face_tri6.h"
30 #include "libmesh/mesh_generation.h"
31 #include "libmesh/mesh_smoother_laplace.h"
32 #include "libmesh/boundary_info.h"
33 #include "libmesh/mesh_triangle_holes.h"
34 #include "libmesh/mesh_triangle_wrapper.h"
35 
36 namespace libMesh
37 {
38 //
39 // Function definitions for the TriangleInterface class
40 //
41 
42 // Constructor
44  : _mesh(mesh),
45  _holes(libmesh_nullptr),
46  _elem_type(TRI3),
47  _desired_area(0.1),
48  _minimum_angle(20.0),
49  _extra_flags(""),
50  _triangulation_type(GENERATE_CONVEX_HULL),
51  _insert_extra_points(false),
52  _smooth_after_generating(true),
53  _serializer(_mesh)
54 {}
55 
56 
57 
58 // Primary function responsible for performing the triangulation
60 {
61  // Will the triangulation have holes?
62  const bool have_holes = ((_holes != libmesh_nullptr) && (!_holes->empty()));
63 
64  // If the initial PSLG is really simple, e.g. an L-shaped domain or
65  // a square/rectangle, the resulting triangulation may be very
66  // "structured" looking. Sometimes this is a problem if your
67  // intention is to work with an "unstructured" looking grid. We can
68  // attempt to work around this limitation by inserting midpoints
69  // into the original PSLG. Inserting additional points into a
70  // set of points meant to be a convex hull usually makes less sense.
71 
72  // May or may not need to insert new points ...
74  {
75  // Make a copy of the original points from the Mesh
76  std::vector<Point> original_points;
77  original_points.reserve (_mesh.n_nodes());
78  for (auto & node : _mesh.node_ptr_range())
79  original_points.push_back(*node);
80 
81  // Clear out the mesh
82  _mesh.clear();
83 
84  // Make sure the new Mesh will be 2D
86 
87  // Insert a new point on each PSLG at some random location
88  // np=index into new points vector
89  // n =index into original points vector
90  for (std::size_t np=0, n=0; np<2*original_points.size(); ++np)
91  {
92  // the even entries are the original points
93  if (np%2==0)
94  _mesh.add_point(original_points[n++]);
95 
96  else // the odd entries are the midpoints of the original PSLG segments
97  _mesh.add_point ((original_points[n] + original_points[n-1])/2);
98  }
99  }
100 
101  // Regardless of whether we added additional points, the set of points to
102  // triangulate is now sitting in the mesh.
103 
104  // If the holes vector is non-NULL (and non-empty) we need to determine
105  // the number of additional points which the holes will add to the
106  // triangulation.
107  unsigned int n_hole_points = 0;
108 
109  if (have_holes)
110  {
111  for (std::size_t i=0; i<_holes->size(); ++i)
112  n_hole_points += (*_holes)[i]->n_points();
113  }
114 
115  // Triangle data structure for the mesh
116  TriangleWrapper::triangulateio initial;
117  TriangleWrapper::triangulateio final;
118 
119  // Pseudo-Constructor for the triangle io structs
120  TriangleWrapper::init(initial);
121  TriangleWrapper::init(final);
122 
123  initial.numberofpoints = _mesh.n_nodes() + n_hole_points;
124  initial.pointlist = static_cast<REAL*>(std::malloc(initial.numberofpoints * 2 * sizeof(REAL)));
125 
127  {
128  // Implicit segment ordering: One segment per point, including hole points
129  if (this->segments.empty())
130  initial.numberofsegments = initial.numberofpoints;
131 
132  // User-defined segment ordering: One segment per entry in the segments vector
133  else
134  initial.numberofsegments = this->segments.size();
135  }
136 
138  initial.numberofsegments = n_hole_points; // One segment for each hole point
139 
140  // Debugging
141  // libMesh::out << "Number of segments set to: " << initial.numberofsegments << std::endl;
142 
143  // Allocate space for the segments (2 int per segment)
144  if (initial.numberofsegments > 0)
145  {
146  initial.segmentlist = static_cast<int *> (std::malloc(initial.numberofsegments * 2 * sizeof(int)));
147  }
148 
149 
150  // Copy all the holes' points and segments into the triangle struct.
151 
152  // The hole_offset is a constant offset into the points vector which points
153  // past the end of the last hole point added.
154  unsigned int hole_offset=0;
155 
156  if (have_holes)
157  for (std::size_t i=0; i<_holes->size(); ++i)
158  {
159  for (unsigned int ctr=0, h=0; h<(*_holes)[i]->n_points(); ctr+=2, ++h)
160  {
161  Point p = (*_holes)[i]->point(h);
162 
163  const unsigned int index0 = 2*hole_offset+ctr;
164  const unsigned int index1 = 2*hole_offset+ctr+1;
165 
166  // Save the x,y locations in the triangle struct.
167  initial.pointlist[index0] = p(0);
168  initial.pointlist[index1] = p(1);
169 
170  // Set the points which define the segments
171  initial.segmentlist[index0] = hole_offset+h;
172  initial.segmentlist[index1] = (h==(*_holes)[i]->n_points()-1) ? hole_offset : hole_offset+h+1; // wrap around
173  }
174 
175  // Update the hole_offset for the next hole
176  hole_offset += (*_holes)[i]->n_points();
177  }
178 
179 
180  // Copy all the non-hole points and segments into the triangle struct.
181  {
182  dof_id_type ctr=0;
183  for (auto & node : _mesh.node_ptr_range())
184  {
185  dof_id_type index = 2*hole_offset + ctr;
186 
187  // Set x,y values in pointlist
188  initial.pointlist[index] = (*node)(0);
189  initial.pointlist[index+1] = (*node)(1);
190 
191  // If the user requested a PSLG, the non-hole points are also segments
193  {
194  // Use implicit ordering to define segments
195  if (this->segments.empty())
196  {
197  dof_id_type n = ctr/2; // ctr is always even
198  initial.segmentlist[index] = hole_offset+n;
199  initial.segmentlist[index+1] = (n==_mesh.n_nodes()-1) ? hole_offset : hole_offset+n+1; // wrap around
200  }
201  }
202 
203  ctr +=2;
204  }
205  }
206 
207 
208  // If the user provided it, use his ordering to define the segments
209  for (std::size_t ctr=0, s=0; s<this->segments.size(); ctr+=2, ++s)
210  {
211  const unsigned int index0 = 2*hole_offset+ctr;
212  const unsigned int index1 = 2*hole_offset+ctr+1;
213 
214  initial.segmentlist[index0] = hole_offset + this->segments[s].first;
215  initial.segmentlist[index1] = hole_offset + this->segments[s].second;
216  }
217 
218 
219 
220  // Tell the input struct about the holes
221  if (have_holes)
222  {
223  initial.numberofholes = _holes->size();
224  initial.holelist = static_cast<REAL*>(std::malloc(initial.numberofholes * 2 * sizeof(REAL)));
225  for (std::size_t i=0, ctr=0; i<_holes->size(); ++i, ctr+=2)
226  {
227  Point inside_point = (*_holes)[i]->inside();
228  initial.holelist[ctr] = inside_point(0);
229  initial.holelist[ctr+1] = inside_point(1);
230  }
231  }
232 
233  // Set the triangulation flags.
234  // c ~ enclose convex hull with segments
235  // z ~ use zero indexing
236  // B ~ Suppresses boundary markers in the output
237  // Q ~ run in "quiet" mode
238  // p ~ Triangulates a Planar Straight Line Graph
239  // If the `p' switch is used, `segmentlist' must point to a list of
240  // segments, `numberofsegments' must be properly set, and
241  // `segmentmarkerlist' must either be set to NULL (in which case all
242  // markers default to zero), or must point to a list of markers.
243  // D ~ Conforming Delaunay: use this switch if you want all triangles
244  // in the mesh to be Delaunay, and not just constrained Delaunay
245  // q ~ Quality mesh generation with no angles smaller than 20 degrees.
246  // An alternate minimum angle may be specified after the q
247  // a ~ Imposes a maximum triangle area constraint.
248  // -P Suppresses the output .poly file. Saves disk space, but you lose the ability to maintain
249  // constraining segments on later refinements of the mesh.
250  // Create the flag strings, depends on element type
251  std::ostringstream flags;
252 
253  // Default flags always used
254  flags << "zBPQ";
255 
256  // Flags which are specific to the type of triangulation
257  switch (_triangulation_type)
258  {
260  {
261  flags << "c";
262  break;
263  }
264 
265  case PSLG:
266  {
267  flags << "p";
268  break;
269  }
270 
272  libmesh_error_msg("ERROR: INVALID_TRIANGULATION_TYPE selected!");
273 
274  default:
275  libmesh_error_msg("Unrecognized _triangulation_type");
276  }
277 
278 
279  // Flags specific to the type of element
280  switch (_elem_type)
281  {
282  case TRI3:
283  {
284  // do nothing.
285  break;
286  }
287 
288  case TRI6:
289  {
290  flags << "o2";
291  break;
292  }
293 
294  default:
295  libmesh_error_msg("ERROR: Unrecognized triangular element type.");
296  }
297 
298 
299  // If we do have holes and the user asked to GENERATE_CONVEX_HULL,
300  // need to add the p flag so the triangulation respects those segments.
301  if ((_triangulation_type==GENERATE_CONVEX_HULL) && (have_holes))
302  flags << "p";
303 
304  // Finally, add the area constraint
305  if (_desired_area > TOLERANCE)
306  flags << "a" << std::fixed << _desired_area;
307 
308  // add minimum angle constraint
310  flags << "q" << std::fixed << _minimum_angle;
311 
312  // add user provided extra flags
313  if (_extra_flags.size() > 0)
314  flags << _extra_flags;
315 
316  // Refine the initial output to conform to the area constraint
317  TriangleWrapper::triangulate(const_cast<char *>(flags.str().c_str()),
318  &initial,
319  &final,
320  libmesh_nullptr); // voronoi ouput -- not used
321 
322 
323  // Send the information computed by Triangle to the Mesh.
325  _mesh,
326  _elem_type);
327 
328  // To the naked eye, a few smoothing iterations usually looks better,
329  // so we do this by default unless the user says not to.
330  if (this->_smooth_after_generating)
332 
333 
334  // Clean up.
337 
338  // Prepare the mesh for use before returning. This ensures (among
339  // other things) that it is partitioned and therefore users can
340  // iterate over local elements, etc.
342 }
343 
344 } // namespace libMesh
345 
346 
347 
348 
349 
350 
351 
352 #endif // LIBMESH_HAVE_TRIANGLE
Real _desired_area
The desired area for the elements in the resulting mesh.
Uses the triangle library to first generate a convex hull from the set of points passed in...
UnstructuredMesh & _mesh
Reference to the mesh which is to be created by triangle.
std::string _extra_flags
Additional flags to be passed to triangle.
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
TriangulationType _triangulation_type
The type of triangulation to perform: choices are: convex hull PSLG.
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...
This class defines the data structures necessary for Laplace smoothing.
const std::vector< Hole * > * _holes
A pointer to a vector of Hole*s.
void triangulate()
This is the main public interface for this function.
The UnstructuredMesh class is derived from the MeshBase class.
void copy_tri_to_mesh(const triangulateio &triangle_data_input, UnstructuredMesh &mesh_output, const ElemType type)
Copies triangulation data computed by triangle from a triangulateio object to a LibMesh mesh...
bool _insert_extra_points
Flag which tells whether or not to insert additional nodes before triangulation.
virtual SimpleRange< node_iterator > node_ptr_range()=0
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
Real _minimum_angle
Minimum angle in triangles.
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
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:199
std::vector< std::pair< unsigned int, unsigned int > > segments
When constructing a PSLG, it is often not possible to do so implicitly through the ordering of the po...
virtual void clear()
Deletes all the data that are currently stored.
Definition: mesh_base.C:285
bool _smooth_after_generating
Flag which tells whether we should smooth the mesh after it is generated.
void destroy(triangulateio &t, IO_Type)
Frees any memory which has been dynamically allocated by Triangle.
virtual dof_id_type n_nodes() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void smooth() libmesh_override
Redefinition of the smooth function from the base class.
TriangleInterface(UnstructuredMesh &mesh)
The constructor.
uint8_t dof_id_type
Definition: id_types.h:64
ElemType _elem_type
The type of elements to generate.