The source file miscellaneous_ex6.C with comments:

Miscellaneous Example 6 - Meshing with LibMesh's TetGen and Triangle Interfaces



LibMesh provides interfaces to both Triangle and TetGen for generating Delaunay triangulations and tetrahedralizations in two and three dimensions (respectively).

Local header files
        #include "libmesh/elem.h"
        #include "libmesh/face_tri3.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/mesh_tetgen_interface.h"
        #include "libmesh/mesh_triangle_holes.h"
        #include "libmesh/mesh_triangle_interface.h"
        #include "libmesh/node.h"
        #include "libmesh/serial_mesh.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Major functions called by main
        void triangulate_domain(const Parallel::Communicator& comm);
        void tetrahedralize_domain(const Parallel::Communicator& comm);
        
Helper routine for tetrahedralize_domain(). Adds the points and elements of a convex hull generated by TetGen to the input mesh
        void add_cube_convex_hull_to_mesh(MeshBase& mesh, Point lower_limit, Point upper_limit);
        
        
        
        
Begin the main program.
        int main (int argc, char** argv)
        {
Initialize libMesh and any dependent libaries, like in example 2.
          LibMeshInit init (argc, argv);
        
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
          std::cout << "Triangulating an L-shaped domain with holes" << std::endl;
        
1.) 2D triangulation of L-shaped domain with three holes of different shape
          triangulate_domain(init.comm());
        
          libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
        
          std::cout << "Tetrahedralizing a prismatic domain with a hole" << std::endl;
        
2.) 3D tetrahedralization of rectangular domain with hole.
          tetrahedralize_domain(init.comm());
        
          return 0;
        }
        
        
        
        
        void triangulate_domain(const Parallel::Communicator& comm)
        {
        #ifdef LIBMESH_HAVE_TRIANGLE
Use typedefs for slightly less typing.
          typedef TriangleInterface::Hole Hole;
          typedef TriangleInterface::PolygonHole PolygonHole;
          typedef TriangleInterface::ArbitraryHole ArbitraryHole;
        
Libmesh mesh that will eventually be created.
          Mesh mesh(comm, 2);
        
The points which make up the L-shape:
          mesh.add_point(Point( 0. ,  0.));
          mesh.add_point(Point( 0. , -1.));
          mesh.add_point(Point(-1. , -1.));
          mesh.add_point(Point(-1. ,  1.));
          mesh.add_point(Point( 1. ,  1.));
          mesh.add_point(Point( 1. ,  0.));
        
Declare the TriangleInterface object. This is where we can set parameters of the triangulation and where the actual triangulate function lives.
          TriangleInterface t(mesh);
        
Customize the variables for the triangulation
          t.desired_area()       = .01;
        
A Planar Straight Line Graph (PSLG) is essentially a list of segments which have to exist in the final triangulation. For an L-shaped domain, Triangle will compute the convex hull of boundary points if we do not specify the PSLG. The PSLG algorithm is also required for triangulating domains containing holes
          t.triangulation_type() = TriangleInterface::PSLG;
        
Turn on/off Laplacian mesh smoothing after generation. By default this is on.
          t.smooth_after_generating() = true;
        
Define holes...

hole_1 is a circle (discretized by 50 points)
          PolygonHole hole_1(Point(-0.5,  0.5), // center
        		     0.25,              // radius
        		     50);               // n. points
        
hole_2 is itself a triangle
          PolygonHole hole_2(Point(0.5, 0.5),   // center
        		     0.1,               // radius
        		     3);                // n. points
        
hole_3 is an ellipse of 100 points which we define here
          Point ellipse_center(-0.5,  -0.5);
          const unsigned int n_ellipse_points=100;
          std::vector<Point> ellipse_points(n_ellipse_points);
          const Real
            dtheta = 2*libMesh::pi / static_cast<Real>(n_ellipse_points),
            a = .1,
            b = .2;
        
          for (unsigned int i=0; i<n_ellipse_points; ++i)
            ellipse_points[i]= Point(ellipse_center(0)+a*cos(i*dtheta),
        			     ellipse_center(1)+b*sin(i*dtheta));
        
          ArbitraryHole hole_3(ellipse_center, ellipse_points);
        
Create the vector of Hole*'s ...
          std::vector<Hole*> holes;
          holes.push_back(&hole_1);
          holes.push_back(&hole_2);
          holes.push_back(&hole_3);
        
... and attach it to the triangulator object
          t.attach_hole_list(&holes);
        
Triangulate!
          t.triangulate();
        
Write the result to file
          mesh.write("delaunay_l_shaped_hole.e");
        
        #endif // LIBMESH_HAVE_TRIANGLE
        }
        
        
        
        void tetrahedralize_domain(const Parallel::Communicator& comm)
        {
        #ifdef LIBMESH_HAVE_TETGEN
The algorithm is broken up into several steps: 1.) A convex hull is constructed for a rectangular hole. 2.) A convex hull is constructed for the domain exterior. 3.) Neighbor information is updated so TetGen knows there is a convex hull 4.) A vector of hole points is created. 5.) The domain is tetrahedralized, the mesh is written out, etc.

The mesh we will eventually generate
          SerialMesh mesh(comm,3);
        
Lower and Upper bounding box limits for a rectangular hole within the unit cube.
          Point hole_lower_limit(0.2, 0.2, 0.4);
          Point hole_upper_limit(0.8, 0.8, 0.6);
        
1.) Construct a convex hull for the hole
          add_cube_convex_hull_to_mesh(mesh, hole_lower_limit, hole_upper_limit);
        
2.) Generate elements comprising the outer boundary of the domain.
          add_cube_convex_hull_to_mesh(mesh, Point(0.,0.,0.), Point(1., 1., 1.));
        
3.) Update neighbor information so that TetGen can verify there is a convex hull.
          mesh.find_neighbors();
        
4.) Set up vector of hole points
          std::vector<Point> hole(1);
          hole[0] = Point( 0.5*(hole_lower_limit + hole_upper_limit) );
        
5.) Set parameters and tetrahedralize the domain

0 means "use TetGen default value"
          Real quality_constraint = 2.0;
        
The volume constraint determines the max-allowed tetrahedral volume in the Mesh. TetGen will split cells which are larger than this size
          Real volume_constraint = 0.001;
        
Construct the Delaunay tetrahedralization
          TetGenMeshInterface t(mesh);
          t.triangulate_conformingDelaunayMesh_carvehole(hole,
        						 quality_constraint,
        						 volume_constraint);
        
Find neighbors, etc in preparation for writing out the Mesh
          mesh.prepare_for_use();
        
Finally, write out the result
          mesh.write("hole_3D.e");
        #endif // LIBMESH_HAVE_TETGEN
        }
        
        
        
        
        
        
        
        
        
        
        
        
        
        void add_cube_convex_hull_to_mesh(MeshBase& mesh, Point lower_limit, Point upper_limit)
        {
        #ifdef LIBMESH_HAVE_TETGEN
          SerialMesh cube_mesh(mesh.comm(),3);
        
          unsigned n_elem = 1;
        
          MeshTools::Generation::build_cube(cube_mesh,
        				    n_elem,n_elem,n_elem, // n. elements in each direction
        				    lower_limit(0), upper_limit(0),
        				    lower_limit(1), upper_limit(1),
        				    lower_limit(2), upper_limit(2),
        				    HEX8);
        
The pointset_convexhull() algorithm will ignore the Hex8s in the Mesh, and just construct the triangulation of the convex hull.
          TetGenMeshInterface t(cube_mesh);
          t.pointset_convexhull();
        
Now add all nodes from the boundary of the cube_mesh to the input mesh.

Map from "node id in cube_mesh" -> "node id in mesh". Initially inserted with a dummy value, later to be assigned a value by the input mesh.
          std::map<unsigned,unsigned> node_id_map;
          typedef std::map<unsigned,unsigned>::iterator iterator;
        
          {
            MeshBase::element_iterator it = cube_mesh.elements_begin();
            const MeshBase::element_iterator end = cube_mesh.elements_end();
            for ( ; it != end; ++it)
              {
        	Elem* elem = *it;
        
        	for (unsigned s=0; s<elem->n_sides(); ++s)
        	  if (elem->neighbor(s) == NULL)
        	    {
Add the node IDs of this side to the set
                      AutoPtr<Elem> side = elem->side(s);
        
        	      for (unsigned n=0; n<side->n_nodes(); ++n)
        		node_id_map.insert( std::make_pair(side->node(n), /*dummy_value=*/0) );
        	    }
              }
          }
        
For each node in the map, insert it into the input mesh and keep track of the ID assigned.
          for (iterator it=node_id_map.begin(); it != node_id_map.end(); ++it)
            {
Id of the node in the cube mesh
              unsigned id = (*it).first;
        
Pointer to node in the cube mesh
              Node* old_node = cube_mesh.node_ptr(id);
        
Add geometric point to input mesh
              Node* new_node = mesh.add_point ( *old_node );
        
Track ID value of new_node in map
              (*it).second = new_node->id();
            }
        
With the points added and the map data structure in place, we are ready to add each TRI3 element of the cube_mesh to the input Mesh with proper node assignments
          {
            MeshBase::element_iterator       el     = cube_mesh.elements_begin();
            const MeshBase::element_iterator end_el = cube_mesh.elements_end();
        
            for (; el != end_el; ++el)
              {
        	Elem* old_elem = *el;
        
        	if (old_elem->type() == TRI3)
        	  {
        	    Elem* new_elem = mesh.add_elem(new Tri3);
        
Assign nodes in new elements. Since this is an example, we'll do it in several steps.
                    for (unsigned i=0; i<old_elem->n_nodes(); ++i)
        	      {
Locate old node ID in the map
                        iterator it = node_id_map.find(old_elem->node(i));
        
Check for not found
                        if (it == node_id_map.end())
        		  {
        		    libMesh::err << "Node id " << old_elem->node(i) << " not found in map!" << std::endl;
        		    libmesh_error();
        		  }
        
Mapping to node ID in input mesh
                        unsigned new_node_id = (*it).second;
        
Node pointer assigned from input mesh
                        new_elem->set_node(i) = mesh.node_ptr(new_node_id);
        	      }
        	  }
              }
          }
        #endif // LIBMESH_HAVE_TETGEN
        }



The source file miscellaneous_ex6.C without comments:

 
  
  #include "libmesh/elem.h"
  #include "libmesh/face_tri3.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/mesh_tetgen_interface.h"
  #include "libmesh/mesh_triangle_holes.h"
  #include "libmesh/mesh_triangle_interface.h"
  #include "libmesh/node.h"
  #include "libmesh/serial_mesh.h"
  
  using namespace libMesh;
  
  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)
  {
    LibMeshInit init (argc, argv);
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    std::cout << "Triangulating an L-shaped domain with holes" << std::endl;
  
    triangulate_domain(init.comm());
  
    libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
  
    std::cout << "Tetrahedralizing a prismatic domain with a hole" << std::endl;
  
    tetrahedralize_domain(init.comm());
  
    return 0;
  }
  
  
  
  
  void triangulate_domain(const Parallel::Communicator& comm)
  {
  #ifdef LIBMESH_HAVE_TRIANGLE
    typedef TriangleInterface::Hole Hole;
    typedef TriangleInterface::PolygonHole PolygonHole;
    typedef TriangleInterface::ArbitraryHole ArbitraryHole;
  
    Mesh mesh(comm, 2);
  
    mesh.add_point(Point( 0. ,  0.));
    mesh.add_point(Point( 0. , -1.));
    mesh.add_point(Point(-1. , -1.));
    mesh.add_point(Point(-1. ,  1.));
    mesh.add_point(Point( 1. ,  1.));
    mesh.add_point(Point( 1. ,  0.));
  
    TriangleInterface t(mesh);
  
    t.desired_area()       = .01;
  
    t.triangulation_type() = TriangleInterface::PSLG;
  
    t.smooth_after_generating() = true;
  
  
    PolygonHole hole_1(Point(-0.5,  0.5), // center
  		     0.25,              // radius
  		     50);               // n. points
  
    PolygonHole hole_2(Point(0.5, 0.5),   // center
  		     0.1,               // radius
  		     3);                // n. points
  
    Point ellipse_center(-0.5,  -0.5);
    const unsigned int n_ellipse_points=100;
    std::vector<Point> ellipse_points(n_ellipse_points);
    const Real
      dtheta = 2*libMesh::pi / static_cast<Real>(n_ellipse_points),
      a = .1,
      b = .2;
  
    for (unsigned int i=0; i<n_ellipse_points; ++i)
      ellipse_points[i]= Point(ellipse_center(0)+a*cos(i*dtheta),
  			     ellipse_center(1)+b*sin(i*dtheta));
  
    ArbitraryHole hole_3(ellipse_center, ellipse_points);
  
    std::vector<Hole*> holes;
    holes.push_back(&hole_1);
    holes.push_back(&hole_2);
    holes.push_back(&hole_3);
  
    t.attach_hole_list(&holes);
  
    t.triangulate();
  
    mesh.write("delaunay_l_shaped_hole.e");
  
  #endif // LIBMESH_HAVE_TRIANGLE
  }
  
  
  
  void tetrahedralize_domain(const Parallel::Communicator& comm)
  {
  #ifdef LIBMESH_HAVE_TETGEN
  
    SerialMesh mesh(comm,3);
  
    Point hole_lower_limit(0.2, 0.2, 0.4);
    Point hole_upper_limit(0.8, 0.8, 0.6);
  
    add_cube_convex_hull_to_mesh(mesh, hole_lower_limit, hole_upper_limit);
  
    add_cube_convex_hull_to_mesh(mesh, Point(0.,0.,0.), Point(1., 1., 1.));
  
    mesh.find_neighbors();
  
    std::vector<Point> hole(1);
    hole[0] = Point( 0.5*(hole_lower_limit + hole_upper_limit) );
  
  
    Real quality_constraint = 2.0;
  
    Real volume_constraint = 0.001;
  
    TetGenMeshInterface t(mesh);
    t.triangulate_conformingDelaunayMesh_carvehole(hole,
  						 quality_constraint,
  						 volume_constraint);
  
    mesh.prepare_for_use();
  
    mesh.write("hole_3D.e");
  #endif // LIBMESH_HAVE_TETGEN
  }
  
  
  
  
  
  
  
  
  
  
  
  
  
  void add_cube_convex_hull_to_mesh(MeshBase& mesh, Point lower_limit, Point upper_limit)
  {
  #ifdef LIBMESH_HAVE_TETGEN
    SerialMesh cube_mesh(mesh.comm(),3);
  
    unsigned n_elem = 1;
  
    MeshTools::Generation::build_cube(cube_mesh,
  				    n_elem,n_elem,n_elem, // n. elements in each direction
  				    lower_limit(0), upper_limit(0),
  				    lower_limit(1), upper_limit(1),
  				    lower_limit(2), upper_limit(2),
  				    HEX8);
  
    TetGenMeshInterface t(cube_mesh);
    t.pointset_convexhull();
  
  
    std::map<unsigned,unsigned> node_id_map;
    typedef std::map<unsigned,unsigned>::iterator iterator;
  
    {
      MeshBase::element_iterator it = cube_mesh.elements_begin();
      const MeshBase::element_iterator end = cube_mesh.elements_end();
      for ( ; it != end; ++it)
        {
  	Elem* elem = *it;
  
  	for (unsigned s=0; s<elem->n_sides(); ++s)
  	  if (elem->neighbor(s) == NULL)
  	    {
  	      AutoPtr<Elem> side = elem->side(s);
  
  	      for (unsigned n=0; n<side->n_nodes(); ++n)
  		node_id_map.insert( std::make_pair(side->node(n), /*dummy_value=*/0) );
  	    }
        }
    }
  
    for (iterator it=node_id_map.begin(); it != node_id_map.end(); ++it)
      {
        unsigned id = (*it).first;
  
        Node* old_node = cube_mesh.node_ptr(id);
  
        Node* new_node = mesh.add_point ( *old_node );
  
        (*it).second = new_node->id();
      }
  
    {
      MeshBase::element_iterator       el     = cube_mesh.elements_begin();
      const MeshBase::element_iterator end_el = cube_mesh.elements_end();
  
      for (; el != end_el; ++el)
        {
  	Elem* old_elem = *el;
  
  	if (old_elem->type() == TRI3)
  	  {
  	    Elem* new_elem = mesh.add_elem(new Tri3);
  
  	    for (unsigned i=0; i<old_elem->n_nodes(); ++i)
  	      {
  		iterator it = node_id_map.find(old_elem->node(i));
  
  		if (it == node_id_map.end())
  		  {
  		    libMesh::err << "Node id " << old_elem->node(i) << " not found in map!" << std::endl;
  		    libmesh_error();
  		  }
  
  		unsigned new_node_id = (*it).second;
  
  		new_elem->set_node(i) = mesh.node_ptr(new_node_id);
  	      }
  	  }
        }
    }
  #endif // LIBMESH_HAVE_TETGEN
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/miscellaneous/miscellaneous_ex6'
***************************************************************
* Running Example miscellaneous_ex6:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
Triangulating an L-shaped domain with holes
Tetrahedralizing a prismatic domain with a hole

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:51:25 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.269246, Active time=0.133998                                            |
 -----------------------------------------------------------------------------------------------------------
| Event                         nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                         w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|-----------------------------------------------------------------------------------------------------------|
|                                                                                                           |
|                                                                                                           |
| Mesh                                                                                                      |
|   find_neighbors()            5         0.0081      0.001623    0.0086      0.001712    6.06     6.39     |
|   renumber_nodes_and_elem()   6         0.0004      0.000072    0.0004      0.000072    0.32     0.32     |
|   write()                     2         0.0735      0.036749    0.0737      0.036837    54.85    54.98    |
|                                                                                                           |
| MeshCommunication                                                                                         |
|   compute_hilbert_indices()   2         0.0165      0.008250    0.0165      0.008250    12.31    12.31    |
|   find_global_indices()       2         0.0021      0.001046    0.0195      0.009739    1.56     14.54    |
|   parallel_sort()             2         0.0006      0.000277    0.0006      0.000316    0.41     0.47     |
|                                                                                                           |
| MeshTools::Generation                                                                                     |
|   build_cube()                2         0.0001      0.000048    0.0001      0.000048    0.07     0.07     |
|                                                                                                           |
| MetisPartitioner                                                                                          |
|   partition()                 1         0.0278      0.027838    0.0374      0.037440    20.77    27.94    |
|                                                                                                           |
| Parallel                                                                                                  |
|   allgather()                 6         0.0002      0.000026    0.0002      0.000037    0.11     0.17     |
|   broadcast()                 2         0.0000      0.000005    0.0000      0.000005    0.01     0.01     |
|   max(scalar)                 129       0.0006      0.000005    0.0006      0.000005    0.45     0.45     |
|   max(vector)                 31        0.0002      0.000007    0.0006      0.000018    0.16     0.42     |
|   min(bool)                   153       0.0005      0.000004    0.0005      0.000004    0.40     0.40     |
|   min(scalar)                 124       0.0010      0.000008    0.0010      0.000008    0.72     0.72     |
|   min(vector)                 31        0.0003      0.000008    0.0006      0.000020    0.19     0.46     |
|   probe()                     30        0.0001      0.000003    0.0001      0.000003    0.07     0.07     |
|   receive()                   30        0.0001      0.000003    0.0002      0.000006    0.07     0.14     |
|   send()                      30        0.0001      0.000002    0.0001      0.000002    0.05     0.05     |
|   send_receive()              34        0.0001      0.000004    0.0004      0.000013    0.10     0.33     |
|   sum()                       6         0.0001      0.000008    0.0001      0.000013    0.04     0.06     |
|                                                                                                           |
| Parallel::Request                                                                                         |
|   wait()                      30        0.0000      0.000002    0.0000      0.000002    0.03     0.03     |
|                                                                                                           |
| Partitioner                                                                                               |
|   set_node_processor_ids()    1         0.0011      0.001115    0.0013      0.001277    0.83     0.95     |
|   set_parent_processor_ids()  1         0.0005      0.000508    0.0005      0.000508    0.38     0.38     |
|   single_partition()          2         0.0000      0.000006    0.0000      0.000006    0.01     0.01     |
 -----------------------------------------------------------------------------------------------------------
| Totals:                       662       0.1340                                          100.00            |
 -----------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example miscellaneous_ex6:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/miscellaneous/miscellaneous_ex6'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo