mesh_triangle_interface.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
28 #include "libmesh/face_tri3.h"
29 #include "libmesh/face_tri6.h"
32 #include "libmesh/boundary_info.h"
35 
36 namespace libMesh
37 {
38  //
39  // Function definitions for the TriangleInterface class
40  //
41 
42  // Constructor
44  : _mesh(mesh),
45  _holes(NULL),
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 != NULL) && (!_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 (_mesh.n_nodes());
77 
79  const MeshBase::node_iterator node_end = _mesh.nodes_end();
80 
81  for (unsigned int ctr=0; node_it != node_end; ++node_it)
82  original_points[ctr++] = **node_it;
83 
84  // Clear out the mesh
85  _mesh.clear();
86 
87  // Make sure the new Mesh will be 2D
89 
90  // Insert a new point on each PSLG at some random location
91  // np=index into new points vector
92  // n =index into original points vector
93  for (unsigned int np=0, n=0; np<2*original_points.size(); ++np)
94  {
95  // the even entries are the original points
96  if (np%2==0)
97  _mesh.add_point(original_points[n++]);
98 
99  else // the odd entries are the midpoints of the original PSLG segments
100  _mesh.add_point ((original_points[n] + original_points[n-1])/2);
101  }
102  }
103 
104  // Regardless of whether we added additional points, the set of points to
105  // triangulate is now sitting in the mesh.
106 
107  // If the holes vector is non-NULL (and non-empty) we need to determine
108  // the number of additional points which the holes will add to the
109  // triangulation.
110  unsigned int n_hole_points = 0;
111 
112  if (have_holes)
113  {
114  for (unsigned int i=0; i<_holes->size(); ++i)
115  n_hole_points += (*_holes)[i]->n_points();
116  }
117 
118  // Triangle data structure for the mesh
119  TriangleWrapper::triangulateio initial;
120  TriangleWrapper::triangulateio final;
121 
122  // Pseudo-Constructor for the triangle io structs
123  TriangleWrapper::init(initial);
124  TriangleWrapper::init(final);
125 
126  initial.numberofpoints = _mesh.n_nodes() + n_hole_points;
127  initial.pointlist = static_cast<REAL*>(std::malloc(initial.numberofpoints * 2 * sizeof(REAL)));
128 
130  {
131  // Implicit segment ordering: One segment per point, including hole points
132  if (this->segments.empty())
133  initial.numberofsegments = initial.numberofpoints;
134 
135  // User-defined segment ordering: One segment per entry in the segments vector
136  else
137  initial.numberofsegments = this->segments.size();
138  }
139 
141  initial.numberofsegments = n_hole_points; // One segment for each hole point
142 
143  // Debugging
144  // libMesh::out << "Number of segments set to: " << initial.numberofsegments << std::endl;
145 
146  // Allocate space for the segments (2 int per segment)
147  if (initial.numberofsegments > 0)
148  {
149  initial.segmentlist = static_cast<int*> (std::malloc(initial.numberofsegments * 2 * sizeof(int)));
150  }
151 
152 
153  // Copy all the holes' points and segments into the triangle struct.
154 
155  // The hole_offset is a constant offset into the points vector which points
156  // past the end of the last hole point added.
157  unsigned int hole_offset=0;
158 
159  if (have_holes)
160  for (unsigned int i=0; i<_holes->size(); ++i)
161  {
162  for (unsigned int ctr=0, h=0; h<(*_holes)[i]->n_points(); ctr+=2, ++h)
163  {
164  Point p = (*_holes)[i]->point(h);
165 
166  const unsigned int index0 = 2*hole_offset+ctr;
167  const unsigned int index1 = 2*hole_offset+ctr+1;
168 
169  // Save the x,y locations in the triangle struct.
170  initial.pointlist[index0] = p(0);
171  initial.pointlist[index1] = p(1);
172 
173  // Set the points which define the segments
174  initial.segmentlist[index0] = hole_offset+h;
175  initial.segmentlist[index1] = (h==(*_holes)[i]->n_points()-1) ? hole_offset : hole_offset+h+1; // wrap around
176  }
177 
178  // Update the hole_offset for the next hole
179  hole_offset += (*_holes)[i]->n_points();
180  }
181 
182 
183  // Copy all the non-hole points and segments into the triangle struct.
184  {
187 
188  for (dof_id_type ctr=0; it != end; ctr+=2, ++it)
189  {
190  dof_id_type index = 2*hole_offset + ctr;
191 
192  // Get pointer to the current node
193  Node* node = *it;
194 
195  // Set x,y values in pointlist
196  initial.pointlist[index] = (*node)(0);
197  initial.pointlist[index+1] = (*node)(1);
198 
199  // If the user requested a PSLG, the non-hole points are also segments
201  {
202  // Use implicit ordering to define segments
203  if (this->segments.empty())
204  {
205  dof_id_type n = ctr/2; // ctr is always even
206  initial.segmentlist[index] = hole_offset+n;
207  initial.segmentlist[index+1] = (n==_mesh.n_nodes()-1) ? hole_offset : hole_offset+n+1; // wrap around
208  }
209  }
210  }
211  }
212 
213 
214  // If the user provided it, use his ordering to define the segments
215  for (unsigned int ctr=0, s=0; s<this->segments.size(); ctr+=2, ++s)
216  {
217  const unsigned int index0 = 2*hole_offset+ctr;
218  const unsigned int index1 = 2*hole_offset+ctr+1;
219 
220  initial.segmentlist[index0] = hole_offset + this->segments[s].first;
221  initial.segmentlist[index1] = hole_offset + this->segments[s].second;
222  }
223 
224 
225 
226  // Tell the input struct about the holes
227  if (have_holes)
228  {
229  initial.numberofholes = _holes->size();
230  initial.holelist = static_cast<REAL*>(std::malloc(initial.numberofholes * 2 * sizeof(REAL)));
231  for (unsigned int i=0, ctr=0; i<_holes->size(); ++i, ctr+=2)
232  {
233  Point inside_point = (*_holes)[i]->inside();
234  initial.holelist[ctr] = inside_point(0);
235  initial.holelist[ctr+1] = inside_point(1);
236  }
237  }
238 
239  // Set the triangulation flags.
240  // c ~ enclose convex hull with segments
241  // z ~ use zero indexing
242  // B ~ Suppresses boundary markers in the output
243  // Q ~ run in "quiet" mode
244  // p ~ Triangulates a Planar Straight Line Graph
245  // If the `p' switch is used, `segmentlist' must point to a list of
246  // segments, `numberofsegments' must be properly set, and
247  // `segmentmarkerlist' must either be set to NULL (in which case all
248  // markers default to zero), or must point to a list of markers.
249  // D ~ Conforming Delaunay: use this switch if you want all triangles
250  // in the mesh to be Delaunay, and not just constrained Delaunay
251  // q ~ Quality mesh generation with no angles smaller than 20 degrees.
252  // An alternate minimum angle may be specified after the q
253  // a ~ Imposes a maximum triangle area constraint.
254  // -P Suppresses the output .poly file. Saves disk space, but you lose the ability to maintain
255  // constraining segments on later refinements of the mesh.
256  // Create the flag strings, depends on element type
257  std::ostringstream flags;
258 
259  // Default flags always used
260  flags << "zBPQ";
261 
262  // Flags which are specific to the type of triangulation
263  switch (_triangulation_type)
264  {
266  {
267  flags << "c";
268  break;
269  }
270 
271  case PSLG:
272  {
273  flags << "p";
274  break;
275  }
276 
278  {
279  libmesh_error();
280  break;
281  }
282 
283  default:
284  {
285  libmesh_error();
286  }
287  }
288 
289 
290  // Flags specific to the type of element
291  switch (_elem_type)
292  {
293  case TRI3:
294  {
295  // do nothing.
296  break;
297  }
298 
299  case TRI6:
300  {
301  flags << "o2";
302  break;
303  }
304 
305  default:
306  {
307  libMesh::err << "ERROR: Unrecognized triangular element type." << std::endl;
308  libmesh_error();
309  }
310  }
311 
312 
313  // If we do have holes and the user asked to GENERATE_CONVEX_HULL,
314  // need to add the p flag so the triangulation respects those segments.
315  if ((_triangulation_type==GENERATE_CONVEX_HULL) && (have_holes))
316  flags << "p";
317 
318  // Finally, add the area constraint
319  if (_desired_area > TOLERANCE)
320  flags << "a" << std::fixed << _desired_area;
321 
322  // add minimum angle constraint
324  flags << "q" << std::fixed << _minimum_angle;
325 
326  // add user provided extra flags
327  if (_extra_flags.size() > 0)
328  flags << _extra_flags;
329 
330  // Refine the initial output to conform to the area constraint
331  TriangleWrapper::triangulate(const_cast<char*>(flags.str().c_str()),
332  &initial,
333  &final,
334  NULL); // voronoi ouput -- not used
335 
336 
337  // Send the information computed by Triangle to the Mesh.
339  _mesh,
340  _elem_type);
341 
342  // To the naked eye, a few smoothing iterations usually looks better,
343  // so we do this by default unless the user says not to.
344  if (this->_smooth_after_generating)
346 
347 
348  // Clean up.
351 
352  }
353 
354 } // namespace libMesh
355 
356 
357 
358 
359 
360 
361 
362 #endif // LIBMESH_HAVE_TRIANGLE

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo