miscellaneous_ex6.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// <h1>Miscellaneous Example 6 - Meshing with LibMesh's TetGen and Triangle Interfaces</h1>
// \author John W. Peterson
// \date 2011
//
// 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/replicated_mesh.h"
#include "libmesh/utility.h" // libmesh_map_find()

// 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 libraries, like in example 2.
  LibMeshInit init (argc, argv);

  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");

  libMesh::out << "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_requires(3 <= LIBMESH_DIM, "3D support");

  libMesh::out << "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();

#ifdef LIBMESH_HAVE_EXODUS_API
  // Write the result to file
  mesh.write("delaunay_l_shaped_hole.e");
#endif

#else
  // Avoid compiler warnings
  libmesh_ignore(comm);
#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
  ReplicatedMesh 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"
  double 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
  double volume_constraint = 0.001;

  // Construct the Delaunay tetrahedralization
  TetGenMeshInterface t(mesh);

  // In debug mode, set the tetgen switches
  // -V (verbose) and
  // -CC (check consistency and constrained Delaunay)
  // In optimized mode, only switch Q (quiet) is set.
  // For more options, see tetgen website:
  // (http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html#cmd-m)
#ifdef DEBUG
  t.set_switches("VCC");
#endif
  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
#ifdef LIBMESH_HAVE_EXODUS_API
  mesh.write("hole_3D.e");
#endif

#else
  // Avoid compiler warnings
  libmesh_ignore(comm);
#endif // LIBMESH_HAVE_TETGEN
}













void add_cube_convex_hull_to_mesh(MeshBase & mesh,
                                  Point lower_limit,
                                  Point upper_limit)
{
#ifdef LIBMESH_HAVE_TETGEN
  ReplicatedMesh 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;

  for (auto & elem : cube_mesh.element_ptr_range())
    for (auto s : elem->side_index_range())
      if (elem->neighbor_ptr(s) == nullptr)
        {
          // Add the node IDs of this side to the set
          std::unique_ptr<Elem> side = elem->side_ptr(s);

          for (auto n : side->node_index_range())
            node_id_map.emplace(side->node_id(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_ref(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
  for (auto & old_elem : cube_mesh.element_ptr_range())
    if (old_elem->type() == TRI3)
      {
        Elem * new_elem = mesh.add_elem(Elem::build(TRI3));

        // Assign nodes in new elements.  Since this is an example,
        // we'll do it in several steps.
        for (auto i : old_elem->node_index_range())
          {
            // Mapping to node ID in input mesh
            unsigned int new_node_id = libmesh_map_find(node_id_map, old_elem->node_id(i));

            // Node pointer assigned from input mesh
            new_elem->set_node(i, mesh.node_ptr(new_node_id));
          }
      }
#else
  // Avoid compiler warnings
  libmesh_ignore(mesh, lower_limit, upper_limit);
#endif // LIBMESH_HAVE_TETGEN
}