mesh_tetgen_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 #include "libmesh/libmesh_config.h"
19 #ifdef LIBMESH_HAVE_TETGEN
20 
21 
22 // C++ includes
23 #include <sstream>
24 
25 // Local includes
26 #include "libmesh/cell_tet4.h"
27 #include "libmesh/face_tri3.h"
30 #include "libmesh/utility.h" // binary_find
32 
33 namespace libMesh
34 {
35 
36  //----------------------------------------------------------------------
37  // TetGenMeshInterface class members
39  _mesh(mesh),
40  _serializer(_mesh)
41  {
42  }
43 
44 
45 
47  {
48  // class tetgen_wrapper allows library access on a basic level:
49  TetGenWrapper tetgen_wrapper;
50 
51  // fill input structure with point set data:
52  this->fill_pointlist(tetgen_wrapper);
53 
54  // Run TetGen triangulation method:
55  // Q = quiet, no terminal output
56  // V = verbose, more terminal output
57  // Note: if no switch is used, the input must be a list of 3D points
58  // (.node file) and the Delaunay tetrahedralization of this point set
59  // will be generated.
60 
61  // Can we apply quality and volume constraints in
62  // triangulate_pointset()?. On at least one test problem,
63  // specifying any quality or volume constraints here causes tetgen
64  // to segfault down in the insphere method: a NULL pointer is passed
65  // to the routine.
66  std::ostringstream oss;
67  oss << "Q"; // quiet operation
68  // oss << "V"; // verbose operation
69  //oss << "q" << std::fixed << 2.0; // quality constraint
70  //oss << "a" << std::fixed << 100.; // volume constraint
71  tetgen_wrapper.set_switches(oss.str());
72 
73  // Run tetgen
74  tetgen_wrapper.run_tetgen();
75 
76  // save elements to mesh structure, nodes will not be changed:
77  const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();
78 
79  // Vector that temporarily holds the node labels defining element.
80  unsigned int node_labels[4];
81 
82  for (unsigned int i=0; i<num_elements; ++i)
83  {
84  Elem* elem = new Tet4;
85 
86  // Get the nodes associated with this element
87  for (unsigned int j=0; j<elem->n_nodes(); ++j)
88  node_labels[j] = tetgen_wrapper.get_element_node(i,j);
89 
90  // Associate the nodes with this element
91  this->assign_nodes_to_elem(node_labels, elem);
92 
93  // Finally, add this element to the mesh.
94  this->_mesh.add_elem(elem);
95  }
96  }
97 
98 
99 
101  {
102  // class tetgen_wrapper allows library access on a basic level
103  TetGenWrapper tetgen_wrapper;
104 
105  // Copy Mesh's node points into TetGen data structure
106  this->fill_pointlist(tetgen_wrapper);
107 
108  // Run TetGen triangulation method:
109  // Q = quiet, no terminal output
110  // Note: if no switch is used, the input must be a list of 3D points
111  // (.node file) and the Delaunay tetrahedralization of this point set
112  // will be generated. In this particular function, we are throwing
113  // away the tetrahedra generated by TetGen, and keeping only the
114  // convex hull...
115  tetgen_wrapper.set_switches("Q");
116  tetgen_wrapper.run_tetgen();
117  unsigned int num_elements = tetgen_wrapper.get_numberoftrifaces();
118 
119  // Delete *all* old elements. Yes, we legally delete elements while
120  // iterating over them because no entries from the underlying container
121  // are actually erased.
122  {
125  for ( ; it != end; ++it)
126  this->_mesh.delete_elem (*it);
127  }
128 
129 
130  // Add the 2D elements which comprise the convex hull back to the mesh.
131  // Vector that temporarily holds the node labels defining element.
132  unsigned int node_labels[3];
133 
134  for (unsigned int i=0; i<num_elements; ++i)
135  {
136  Elem* elem = new Tri3;
137 
138  // Get node labels associated with this element
139  for (unsigned int j=0; j<elem->n_nodes(); ++j)
140  node_labels[j] = tetgen_wrapper.get_triface_node(i,j);
141 
142  this->assign_nodes_to_elem(node_labels, elem);
143 
144  // Finally, add this element to the mesh.
145  this->_mesh.add_elem(elem);
146  }
147  }
148 
149 
150 
151 
152 
154  double volume_constraint)
155  {
156  // start triangulation method with empty holes list:
157  std::vector<Point> noholes;
158  triangulate_conformingDelaunayMesh_carvehole(noholes, quality_constraint, volume_constraint);
159  }
160 
161 
162 
164  double quality_constraint,
165  double volume_constraint)
166  {
167  // Before calling this function, the Mesh must contain a convex hull
168  // of TRI3 elements which define the boundary.
169  unsigned hull_integrity_check = check_hull_integrity();
170 
171  // Possibly die if hull integrity check failed
172  this->process_hull_integrity_result(hull_integrity_check);
173 
174  // class tetgen_wrapper allows library access on a basic level
175  TetGenWrapper tetgen_wrapper;
176 
177  // Copy Mesh's node points into TetGen data structure
178  this->fill_pointlist(tetgen_wrapper);
179 
180  // >>> fill input structure "tetgenio" with facet data:
181  int facet_num = this->_mesh.n_elem();
182 
183  // allocate memory in "tetgenio" structure:
184  tetgen_wrapper.allocate_facetlist(facet_num, holes.size());
185 
186 
187  // Set up tetgen data structures with existing facet information
188  // from the convex hull.
189  {
190  int insertnum = 0;
193  for (; it != end ; ++it)
194  {
195  tetgen_wrapper.allocate_facet_polygonlist(insertnum, 1);
196  tetgen_wrapper.allocate_polygon_vertexlist(insertnum, 0, 3);
197 
198  Elem* elem = *it;
199 
200  for (unsigned int j=0; j<elem->n_nodes(); ++j)
201  {
202  // We need to get the sequential index of elem->get_node(j), but
203  // it should already be stored in _sequential_to_libmesh_node_map...
204  unsigned libmesh_node_id = elem->node(j);
205 
206  // The libmesh node IDs may not be sequential, but can we assume
207  // they are at least in order??? We will do so here.
208  std::vector<unsigned>::iterator node_iter =
211  libmesh_node_id);
212 
213  // Check to see if not found: this could also indicate the sequential
214  // node map is not sorted...
215  if (node_iter == _sequential_to_libmesh_node_map.end())
216  {
217  libMesh::err << "Global node " << libmesh_node_id << " not found in sequential node map!" << std::endl;
218  libmesh_error();
219  }
220 
221  std::vector<unsigned>::difference_type
222  sequential_index = std::distance(_sequential_to_libmesh_node_map.begin(), node_iter);
223 
224  // Debugging:
225  // libMesh::out << "libmesh_node_id=" << libmesh_node_id
226  // << ", sequential_index=" << sequential_index
227  // << std::endl;
228 
229  tetgen_wrapper.set_vertex(insertnum, // facet number
230  0, // polygon (always 0)
231  j, // local vertex index in tetgen input
232  sequential_index);
233  }
234 
235  // Go to next facet in polygonlist
236  insertnum++;
237  }
238  }
239 
240 
241 
242  // fill hole list (if there are holes):
243  if (holes.size() > 0)
244  {
245  std::vector<Point>::const_iterator ihole;
246  unsigned hole_index = 0;
247  for (ihole=holes.begin(); ihole!=holes.end(); ++ihole)
248  tetgen_wrapper.set_hole(hole_index++, (*ihole)(0), (*ihole)(1), (*ihole)(2));
249  }
250 
251 
252  // Run TetGen triangulation method:
253  // p = tetrahedralizes a piecewise linear complex (see definition in user manual)
254  // Q = quiet, no terminal output
255  // q = specify a minimum radius/edge ratio
256  // a = tetrahedron volume constraint
257 
258  // assemble switches:
259  std::ostringstream oss; // string holding switches
260  oss << "pQ";
261 
262  if (quality_constraint != 0)
263  oss << "q" << std::fixed << quality_constraint;
264 
265  if (volume_constraint != 0)
266  oss << "a" << std::fixed << volume_constraint;
267 
268  std::string params = oss.str();
269 
270  tetgen_wrapper.set_switches(params); // TetGen switches: Piecewise linear complex, Quiet mode
271  tetgen_wrapper.run_tetgen();
272 
273  // => nodes:
274  unsigned int old_nodesnum = this->_mesh.n_nodes();
275  REAL x=0., y=0., z=0.;
276  const unsigned int num_nodes = tetgen_wrapper.get_numberofpoints();
277 
278  // Debugging:
279  // libMesh::out << "Original mesh had " << old_nodesnum << " nodes." << std::endl;
280  // libMesh::out << "Reserving space for " << num_nodes << " total nodes." << std::endl;
281 
282  // Reserve space for additional nodes in the node map
283  _sequential_to_libmesh_node_map.reserve(num_nodes);
284 
285  // Add additional nodes to the Mesh.
286  // Original code had i<=num_nodes here (Note: the indexing is:
287  // foo[3*i], [3*i+1], [3*i+2]) But according to the TetGen docs, "In
288  // all cases, the first item in any array is stored starting at
289  // index [0]."
290  for (unsigned int i=old_nodesnum; i<num_nodes; i++)
291  {
292  // Fill in x, y, z values
293  tetgen_wrapper.get_output_node(i, x,y,z);
294 
295  // Catch the node returned by add_point()... this will tell us the ID
296  // assigned by the Mesh.
297  Node* new_node = this->_mesh.add_point ( Point(x,y,z) );
298 
299  // Store this new ID in our sequential-to-libmesh node mapping array
300  _sequential_to_libmesh_node_map.push_back( new_node->id() );
301  }
302 
303  // Debugging:
304  // std::copy(_sequential_to_libmesh_node_map.begin(),
305  // _sequential_to_libmesh_node_map.end(),
306  // std::ostream_iterator<unsigned>(std::cout, " "));
307  // std::cout << std::endl;
308 
309 
310  // => tetrahedra:
311  const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();
312 
313  // Vector that temporarily holds the node labels defining element connectivity.
314  unsigned int node_labels[4];
315 
316  for (unsigned int i=0; i<num_elements; i++)
317  {
318  // TetGen only supports Tet4 elements.
319  Elem* elem = new Tet4;
320 
321  // Fill up the the node_labels vector
322  for (unsigned int j=0; j<elem->n_nodes(); j++)
323  node_labels[j] = tetgen_wrapper.get_element_node(i,j);
324 
325  // Associate nodes with this element
326  this->assign_nodes_to_elem(node_labels, elem);
327 
328  // Finally, add this element to the mesh
329  this->_mesh.add_elem(elem);
330  }
331 
332  // Delete original convex hull elements. Is there ever a case where
333  // we should not do this?
334  this->delete_2D_hull_elements();
335  }
336 
337 
338 
339 
340 
342  {
343  // fill input structure with point set data:
344  wrapper.allocate_pointlist( this->_mesh.n_nodes() );
345 
346  // Make enough space to store a mapping between the implied sequential
347  // node numbering used in tetgen and libmesh's (possibly) non-sequential
348  // numbering scheme.
351 
352  {
353  unsigned index = 0;
355  const MeshBase::node_iterator end = this->_mesh.nodes_end();
356  for ( ; it != end; ++it)
357  {
358  _sequential_to_libmesh_node_map[index] = (*it)->id();
359  wrapper.set_node(index++, (**it)(0), (**it)(1), (**it)(2));
360  }
361  }
362  }
363 
364 
365 
366 
367 
368  void TetGenMeshInterface::assign_nodes_to_elem(unsigned* node_labels, Elem* elem)
369  {
370  for (unsigned int j=0; j<elem->n_nodes(); ++j)
371  {
372  // Get the mapped node index to ask the Mesh for
373  unsigned mapped_node_id = _sequential_to_libmesh_node_map[ node_labels[j] ];
374 
375  // Parallel mesh can return NULL pointers, this is bad...
376  Node* current_node = this->_mesh.node_ptr( mapped_node_id );
377 
378  if (current_node == NULL)
379  {
380  libMesh::err << "Error! Mesh returned NULL node pointer!" << std::endl;
381  libmesh_error();
382  }
383 
384  elem->set_node(j) = current_node;
385  }
386  }
387 
388 
389 
390 
391 
393  {
394  // Check for easy return: if the Mesh is empty (i.e. if
395  // somebody called triangulate_conformingDelaunayMesh on
396  // a Mesh with no elements, then hull integrity check must
397  // fail...
398  if (_mesh.n_elem() == 0)
399  return 3;
400 
403 
404  for (; it != end ; ++it)
405  {
406  Elem* elem = *it;
407 
408  // Check for proper element type
409  if (elem->type() != TRI3)
410  {
411  //libMesh::err << "ERROR: Some of the elements in the original mesh were not TRI3!" << std::endl;
412  //libmesh_error();
413  return 1;
414  }
415 
416  for (unsigned int i=0; i<elem->n_neighbors(); ++i)
417  {
418  if (elem->neighbor(i) == NULL)
419  {
420  // libMesh::err << "ERROR: Non-convex hull, cannot be tetrahedralized." << std::endl;
421  // libmesh_error();
422  return 2;
423  }
424  }
425  }
426 
427  // If we made it here, return success!
428  return 0;
429  }
430 
431 
432 
433 
434 
436  {
437  if (result != 0)
438  {
439  libMesh::err << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
440 
441  if (result==1)
442  {
443  libMesh::err << "Non-TRI3 elements were found in the input Mesh. ";
444  libMesh::err << "A constrained Delaunay triangulation requires a convex hull of TRI3 elements." << std::endl;
445  }
446 
447  libMesh::err << "Consider calling TetGenMeshInterface::pointset_convexhull() followed " << std::endl;
448  libMesh::err << "by Mesh::find_neighbors() first." << std::endl;
449 
450  libmesh_error();
451  }
452  }
453 
454 
455 
456 
458  {
461 
462  for (; it != end ; ++it)
463  {
464  Elem* elem = *it;
465 
466  // Check for proper element type. Yes, we legally delete elements while
467  // iterating over them because no entries from the underlying container
468  // are actually erased.
469  if (elem->type() == TRI3)
470  _mesh.delete_elem(elem);
471  }
472  }
473 
474 
475 
476 } // namespace libMesh
477 
478 
479 #endif // #ifdef LIBMESH_HAVE_TETGEN

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

Hosted By:
SourceForge.net Logo