MeshTriangulationUtils Namespace Reference

Classes
struct	XYDelaunayOptions
	Bundle of inputs for triangulateWithDelaunay. More...

Functions
std::set< std::size_t >	outerBoundaryIds (MeshGenerator &mg, MeshBase &boundary_mesh, const XYDelaunayOptions &opts)
	Resolves the outer-boundary selection of the options into the set of ids that define it: the ids of 'input_boundary_names' or of 'input_subdomain_names' in the boundary mesh.
std::unique_ptr< MeshBase >	triangulateWithDelaunay (MeshGenerator &mg, std::unique_ptr< MeshBase > boundary_mesh, std::vector< std::unique_ptr< MeshBase > > hole_meshes, const XYDelaunayOptions &xyd_opts)
	Performs a 2D Delaunay triangulation (via libMesh::Poly2TriTriangulator) inside a closed boundary mesh with optional holes, optionally stitches the resulting triangulation to each hole mesh that requests it, and applies subdomain / boundary renumbering.
void	finalizeTriangulation (MeshGenerator &mg, UnstructuredMesh &mesh, std::vector< std::unique_ptr< MeshBase > > &holes, const std::vector< bool > &holes_with_midpoints, const XYDelaunayOptions &opts)
	Performs the subdomain and boundary naming, the boundary id remapping (outer boundary to 0 and hole i to i + 1), and the hole stitching on a mesh that has already been triangulated.

Function Documentation

finalizeTriangulation()

void MeshTriangulationUtils::finalizeTriangulation	(	MeshGenerator &	mg,
		UnstructuredMesh &	mesh,
		std::vector< std::unique_ptr< MeshBase > > &	holes,
		const std::vector< bool > &	holes_with_midpoints,
		const XYDelaunayOptions &	opts
	)

Performs the subdomain and boundary naming, the boundary id remapping (outer boundary to 0 and hole i to i + 1), and the hole stitching on a mesh that has already been triangulated.

This is shared by the Poly2Tri-backed and the frontal triangulators.

Parameters

mg	The calling mesh generator (used for paramError reporting)
mesh	The triangulated mesh, modified in place
holes	The hole meshes used by the triangulation (empty vector if no holes)
holes_with_midpoints	Whether the boundary of each hole carries midpoints, as reported by the libMesh MeshedHole built from that hole before the triangulation
opts	Triangulation options (see XYDelaunayOptions)

Definition at line 186 of file MeshTriangulationUtils.C.

{
  SubdomainID output_subdomain_id = opts.has_output_subdomain_id ? opts.output_subdomain_id : 0;
 
  if (opts.has_output_subdomain_name)
  {
    auto id = MooseMeshUtils::getSubdomainID(opts.output_subdomain_name, mesh);
 
    if (id == Elem::invalid_subdomain_id)
    {
      if (!opts.has_output_subdomain_id)
      {
        // We'll probably need to make a new ID, then
        output_subdomain_id = MooseMeshUtils::getNextFreeSubdomainID(mesh);
 
        // But check the hole meshes for our output subdomain name too
        for (auto & hole_ptr : holes)
        {
          auto possible_sbdid =
              MooseMeshUtils::getSubdomainID(opts.output_subdomain_name, *hole_ptr);
          // Huh, it was in one of them
          if (possible_sbdid != Elem::invalid_subdomain_id)
          {
            output_subdomain_id = possible_sbdid;
            break;
          }
          output_subdomain_id =
              std::max(output_subdomain_id, MooseMeshUtils::getNextFreeSubdomainID(*hole_ptr));
        }
      }
    }
    else
    {
      if (opts.has_output_subdomain_id)
      {
        if (id != output_subdomain_id)
          mg.paramError("output_subdomain_name",
                        "name has been used by the input meshes and the corresponding id is not "
                        "equal to 'output_subdomain_id'");
      }
      else
        output_subdomain_id = id;
    }
    // We do not want to set an empty subdomain name
    if (opts.output_subdomain_name.size())
      mesh.set_subdomain_name(output_subdomain_id, opts.output_subdomain_name);
  }
 
  if (opts.smooth_tri || output_subdomain_id)
    for (auto elem : mesh.element_ptr_range())
    {
      mooseAssert(
          elem->type() ==
              (opts.tri_elem_type == "TRI6" ? TRI6 : (opts.tri_elem_type == "TRI7" ? TRI7 : TRI3)),
          "Unexpected element type " << libMesh::Utility::enum_to_string(elem->type())
                                     << " found in triangulation");
 
      elem->subdomain_id() = output_subdomain_id;
 
      // I do not trust Laplacian mesh smoothing not to invert
      // elements near reentrant corners.  Eventually we'll add better
      // smoothing options, but even those might have failure cases.
      // Better to always do extra tests here than to ever let users
      // try to run on a degenerate mesh.
      if (opts.smooth_tri)
      {
        auto cross_prod = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
 
        if (cross_prod(2) <= 0)
          mooseError("Inverted element found in triangulation.\n"
                     "Laplacian smoothing can create these at reentrant corners; disable it?");
      }
    }
 
  // The hole meshes are specified by the user, so they could have any
  // BCID or no BCID or any combination of BCIDs on their outer
  // boundary, so we'll have to set our own BCID to use for stitching
  // there.  We'll need to check all the holes for used BCIDs, if we
  // want to pick a new ID on hole N that doesn't conflict with any
  // IDs on hole M < N (or with the IDs on the new triangulation)
 
  // The new triangulation by default assigns BCID i+1 to hole i ...
  // but we can't even use this for mesh stitching, because we can't
  // be sure it isn't also already in use on the hole's mesh and so we
  // won't be able to safely clear it afterwards.
  const boundary_id_type end_bcid = holes.size() + 1;
 
  // If a hole has its boundary layer mesh, we need to move the hole bcid to the "real" hole
  // boundary in the boundary layer mesh. So we need to record them here.
  std::vector<BoundaryID> hole_boundary_rec(holes.size());
  std::iota(hole_boundary_rec.begin(), hole_boundary_rec.end(), 1);
 
  // For the hole meshes that need to be stitched, we would like to make sure the hole boundary ids
  // and output boundary id are not conflicting with the existing boundary ids of the hole meshes to
  // be stitched.
  BoundaryID free_boundary_id = 0;
  if (opts.stitch_holes.size())
  {
    for (auto hole_i : index_range(holes))
    {
      if (opts.stitch_holes[hole_i])
      {
        free_boundary_id =
            std::max(free_boundary_id, MooseMeshUtils::getNextFreeBoundaryID(*holes[hole_i]));
        holes[hole_i]->comm().max(free_boundary_id);
      }
    }
    for (auto h : index_range(holes))
    {
      libMesh::MeshTools::Modification::change_boundary_id(mesh, h + 1, h + 1 + free_boundary_id);
      hole_boundary_rec[h] = h + 1 + free_boundary_id;
    }
  }
  boundary_id_type new_hole_bcid = end_bcid + free_boundary_id;
 
  // We might be overriding the default bcid numbers.  We have to be
  // careful about how we renumber, though.  We pick unused temporary
  // numbers because e.g. "0->2, 2->0" is impossible to do
  // sequentially, but "0->N, 2->N+2, N->2, N+2->0" works.
  libMesh::MeshTools::Modification::change_boundary_id(
      mesh, 0, (opts.has_output_boundary ? end_bcid : 0) + free_boundary_id);
 
  if (!opts.hole_boundaries.empty())
  {
    auto hole_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, opts.hole_boundaries, true);
 
    for (auto h : index_range(holes))
      libMesh::MeshTools::Modification::change_boundary_id(
          mesh, h + 1 + free_boundary_id, h + 1 + free_boundary_id + end_bcid);
 
    for (auto h : index_range(holes))
    {
      libMesh::MeshTools::Modification::change_boundary_id(
          mesh, h + 1 + free_boundary_id + end_bcid, hole_boundary_ids[h]);
      hole_boundary_rec[h] = hole_boundary_ids[h];
      mesh.get_boundary_info().sideset_name(hole_boundary_ids[h]) = opts.hole_boundaries[h];
      new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(hole_boundary_ids[h] + 1));
    }
  }
 
  if (opts.has_output_boundary)
  {
    const std::vector<BoundaryID> output_boundary_id =
        MooseMeshUtils::getBoundaryIDs(mesh, {opts.output_boundary}, true);
 
    libMesh::MeshTools::Modification::change_boundary_id(
        mesh, end_bcid + free_boundary_id, output_boundary_id[0]);
    mesh.get_boundary_info().sideset_name(output_boundary_id[0]) = opts.output_boundary;
 
    new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(output_boundary_id[0] + 1));
  }
 
  bool doing_stitching = false;
 
  for (auto hole_i : index_range(holes))
  {
    const MeshBase & hole_mesh = *holes[hole_i];
    auto & hole_boundary_info = hole_mesh.get_boundary_info();
    const std::set<boundary_id_type> & local_hole_bcids = hole_boundary_info.get_boundary_ids();
 
    if (!local_hole_bcids.empty())
      new_hole_bcid = std::max(new_hole_bcid, boundary_id_type(*local_hole_bcids.rbegin() + 1));
    hole_mesh.comm().max(new_hole_bcid);
 
    if (hole_i < opts.stitch_holes.size() && opts.stitch_holes[hole_i])
      doing_stitching = true;
  }
 
  const boundary_id_type inner_bcid = new_hole_bcid + 1;
 
  // libMesh mesh stitching still requires a serialized mesh, and it's
  // cheaper to do that once than to do it once-per-hole
  libMesh::MeshSerializer serial(mesh, doing_stitching);
 
  // Define a reference map variable for subdomain map
  auto & main_subdomain_map = mesh.set_subdomain_name_map();
  for (auto hole_i : index_range(holes))
  {
    if (hole_i < opts.stitch_holes.size() && opts.stitch_holes[hole_i])
    {
      UnstructuredMesh & hole_mesh = dynamic_cast<UnstructuredMesh &>(*holes[hole_i]);
      // increase hole mesh order if the triangulation mesh has higher order
      if (!holes_with_midpoints[hole_i])
      {
        if (opts.tri_elem_type == "TRI6")
          hole_mesh.all_second_order();
        else if (opts.tri_elem_type == "TRI7")
          hole_mesh.all_complete_order();
      }
      auto & hole_boundary_info = hole_mesh.get_boundary_info();
 
      // Our algorithm here requires a serialized Mesh.  To avoid
      // redundant serialization and deserialization (libMesh
      // MeshedHole and stitch_meshes still also require
      // serialization) we'll do the serialization up front.
      libMesh::MeshSerializer serial_hole(hole_mesh);
 
      // It would have been nicer for MeshedHole to add the BCID
      // itself, but we want MeshedHole to work with a const mesh.
      // We'll still use MeshedHole, for its code distinguishing
      // outer boundaries from inner boundaries on a
      // hole-with-holes.
      const auto & hole_bdy_id_filter = (hole_i < opts.hole_boundary_id_filters.size())
                                            ? opts.hole_boundary_id_filters[hole_i]
                                            : std::set<std::size_t>();
      libMesh::TriangulatorInterface::MeshedHole mh{hole_mesh, hole_bdy_id_filter};
 
      // We have to translate from MeshedHole points to mesh
      // sides.
      std::unordered_map<Point, Point> next_hole_boundary_point;
      const int np = mh.n_points();
      for (auto pi : make_range(1, np))
        next_hole_boundary_point[mh.point(pi - 1)] = mh.point(pi);
      next_hole_boundary_point[mh.point(np - 1)] = mh.point(0);
 
#ifndef NDEBUG
      int found_hole_sides = 0;
#endif
      for (auto elem : hole_mesh.element_ptr_range())
      {
        if (elem->dim() != 2)
          mooseError("Non 2-D element found in hole; stitching is not supported.");
 
        auto ns = elem->n_sides();
        for (auto s : make_range(ns))
        {
          auto it_s = next_hole_boundary_point.find(elem->point(s));
          if (it_s != next_hole_boundary_point.end())
            if (it_s->second == elem->point((s + 1) % ns))
            {
              hole_boundary_info.add_side(elem, s, new_hole_bcid);
#ifndef NDEBUG
              ++found_hole_sides;
#endif
            }
        }
      }
      mooseAssert(found_hole_sides == np, "Failed to find full outer boundary of meshed hole");
 
      auto & mesh_boundary_info = mesh.get_boundary_info();
#ifndef NDEBUG
      int found_inner_sides = 0;
#endif
      for (auto elem : mesh.element_ptr_range())
      {
        auto ns = elem->n_sides();
        for (auto s : make_range(ns))
        {
          auto it_s = next_hole_boundary_point.find(elem->point((s + 1) % ns));
          if (it_s != next_hole_boundary_point.end())
            if (it_s->second == elem->point(s))
            {
              mesh_boundary_info.add_side(elem, s, inner_bcid);
#ifndef NDEBUG
              ++found_inner_sides;
#endif
            }
        }
      }
      mooseAssert(found_inner_sides == np, "Failed to find full boundary around meshed hole");
 
      // Retrieve subdomain name map from the mesh to be stitched and insert it into the main
      // subdomain map
      const auto & increment_subdomain_map = hole_mesh.get_subdomain_name_map();
      main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
 
      // We do not need the hole_bdy_id_filter anymore
      for (const auto & bcid : hole_bdy_id_filter)
        hole_boundary_info.remove_id(bcid);
      // If we are stitching a hole boundary layer mesh, we need to reassign the bcid
      if (hole_bdy_id_filter.size())
      {
        MooseMeshUtils::changeBoundaryId(
            hole_mesh,
            opts.hole_boundary_inner_id_defaults[hole_i].empty()
                ? 1
                : *opts.hole_boundary_inner_id_defaults[hole_i].begin(),
            hole_boundary_rec[hole_i],
            true);
        hole_mesh.get_boundary_info().sideset_name(hole_boundary_rec[hole_i]) =
            mesh.get_boundary_info().sideset_name(hole_boundary_rec[hole_i]);
        mesh.get_boundary_info().remove_id(hole_boundary_rec[hole_i]);
      }
 
      mesh.stitch_meshes(hole_mesh,
                         inner_bcid,
                         new_hole_bcid,
                         TOLERANCE,
                         /*clear_stitched_bcids*/ true,
                         opts.verbose_stitching,
                         opts.use_binary_search);
    }
  }
  // Check if one SubdomainName is shared by more than one subdomain ids
  std::set<SubdomainName> main_subdomain_map_name_list;
  for (auto const & id_name_pair : main_subdomain_map)
    main_subdomain_map_name_list.emplace(id_name_pair.second);
  if (main_subdomain_map.size() != main_subdomain_map_name_list.size())
    mg.paramError("holes", "The hole meshes contain subdomain name maps with conflicts.");
 
  mesh.unset_is_prepared();
}

Referenced by triangulateWithDelaunay().

outerBoundaryIds()

std::set< std::size_t > MeshTriangulationUtils::outerBoundaryIds	(	MeshGenerator &	mg,
		MeshBase &	boundary_mesh,
		const XYDelaunayOptions &	opts
	)

Resolves the outer-boundary selection of the options into the set of ids that define it: the ids of 'input_boundary_names' or of 'input_subdomain_names' in the boundary mesh.

This is shared by the Poly2Tri-backed and the frontal triangulators, so that both select the outer boundary the same way.

Parameters

mg	The calling mesh generator (used for paramError reporting)
boundary_mesh	The mesh the names are resolved against
opts	Triangulation options carrying the names (see XYDelaunayOptions)

Returns

The resolved ids, empty when neither parameter selects a boundary

Definition at line 29 of file MeshTriangulationUtils.C.

{
  std::set<std::size_t> bdy_ids;
 
  if (!opts.input_boundary_names.empty())
  {
    if (!opts.input_subdomain_names.empty())
      mg.paramError(
          "input_subdomain_names",
          "input_boundary_names and input_subdomain_names cannot both specify an outer boundary.");
 
    for (const auto & name : opts.input_boundary_names)
    {
      auto bcid = MooseMeshUtils::getBoundaryID(name, boundary_mesh);
      if (bcid == BoundaryInfo::invalid_id)
        mg.paramError("input_boundary_names", name, " is not a boundary name in the input mesh");
 
      bdy_ids.insert(bcid);
    }
  }
 
  if (!opts.input_subdomain_names.empty())
  {
    // Make sure subdomain info caches are up to date
    if (!boundary_mesh.preparation().has_cached_elem_data)
      boundary_mesh.cache_elem_data();
 
    const auto subdomain_ids =
        MooseMeshUtils::getSubdomainIDs(boundary_mesh, opts.input_subdomain_names);
 
    // Check that the requested subdomains exist in the mesh
    std::set<SubdomainID> subdomains;
    boundary_mesh.subdomain_ids(subdomains);
 
    for (auto i : index_range(subdomain_ids))
    {
      if (subdomain_ids[i] == Moose::INVALID_BLOCK_ID || !subdomains.count(subdomain_ids[i]))
        mg.paramError("input_subdomain_names",
                      opts.input_subdomain_names[i],
                      " was not found in the boundary mesh");
 
      bdy_ids.insert(subdomain_ids[i]);
    }
  }
 
  return bdy_ids;
}

Referenced by XYFrontalDelaunayGenerator::generate(), and triangulateWithDelaunay().

triangulateWithDelaunay()

std::unique_ptr< MeshBase > MeshTriangulationUtils::triangulateWithDelaunay	(	MeshGenerator &	mg,
		std::unique_ptr< MeshBase >	boundary_mesh,
		std::vector< std::unique_ptr< MeshBase > >	hole_meshes,
		const XYDelaunayOptions &	xyd_opts
	)

Performs a 2D Delaunay triangulation (via libMesh::Poly2TriTriangulator) inside a closed boundary mesh with optional holes, optionally stitches the resulting triangulation to each hole mesh that requests it, and applies subdomain / boundary renumbering.

This is the core algorithm shared by XYDelaunayGenerator and XYTriangleBoundaryLayerGenerator.

Parameters

mg	The calling mesh generator (used for paramError reporting and communicator)
boundary_mesh	The closed-loop boundary mesh
hole_meshes	Optional hole meshes (one mesh per hole; empty vector if no holes)
xyd_opts	Triangulation options (see XYDelaunayOptions)

Returns

The triangulated mesh with optionally stitched holes

Definition at line 78 of file MeshTriangulationUtils.C.

{
  // Put the boundary mesh in a local pointer
  std::unique_ptr<UnstructuredMesh> mesh =
      dynamic_pointer_cast<UnstructuredMesh>(std::move(boundary_mesh));
 
  // Get ready to triangulate the line segments we extract from it
  libMesh::Poly2TriTriangulator poly2tri(*mesh);
  poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
 
  // If we're using a user-requested subset of boundaries on that
  // mesh, get their ids.
  const std::set<std::size_t> bdy_ids = outerBoundaryIds(mg, *mesh, xyd_opts);
 
  if (!bdy_ids.empty())
    poly2tri.set_outer_boundary_ids(bdy_ids);
 
  poly2tri.set_interpolate_boundary_points(xyd_opts.add_nodes_per_boundary_segment);
  poly2tri.set_refine_boundary_allowed(xyd_opts.refine_bdy);
  poly2tri.set_verify_hole_boundaries(xyd_opts.verify_holes);
 
  poly2tri.desired_area() = xyd_opts.desired_area;
  poly2tri.minimum_angle() = 0; // Not yet supported
  poly2tri.smooth_after_generating() = xyd_opts.smooth_tri;
 
  std::vector<libMesh::TriangulatorInterface::MeshedHole> meshed_holes;
  std::vector<libMesh::TriangulatorInterface::Hole *> triangulator_hole_ptrs(hole_meshes.size());
  // This tells us the element orders of the hole meshes
  // For the boundary meshes, it can be access through poly2tri.segment_midpoints.
  std::vector<bool> holes_with_midpoints(hole_meshes.size());
  bool stitch_second_order_holes(false);
 
  // Make sure pointers here aren't invalidated by a resize
  meshed_holes.reserve(hole_meshes.size());
  for (auto hole_i : index_range(hole_meshes))
  {
    if (!hole_meshes[hole_i]->is_prepared())
      hole_meshes[hole_i]->prepare_for_use();
    if (hole_i < xyd_opts.hole_boundary_id_filters.size() &&
        !xyd_opts.hole_boundary_id_filters[hole_i].empty())
      meshed_holes.emplace_back(*hole_meshes[hole_i], xyd_opts.hole_boundary_id_filters[hole_i]);
    else
      meshed_holes.emplace_back(*hole_meshes[hole_i]);
    holes_with_midpoints[hole_i] = meshed_holes.back().n_midpoints();
    stitch_second_order_holes =
        xyd_opts.stitch_holes.empty()
            ? false
            : ((holes_with_midpoints[hole_i] && xyd_opts.stitch_holes[hole_i]) ||
               stitch_second_order_holes);
    if (hole_i < xyd_opts.refine_holes.size())
      meshed_holes.back().set_refine_boundary_allowed(xyd_opts.refine_holes[hole_i]);
 
    triangulator_hole_ptrs[hole_i] = &meshed_holes.back();
  }
  if (stitch_second_order_holes &&
      (xyd_opts.tri_elem_type == "TRI3" || xyd_opts.tri_elem_type == "DEFAULT"))
    mg.paramError(
        "tri_element_type",
        "Cannot use first order elements with stitched quadratic element holes. Please try "
        "to specify a higher-order tri_element_type or reduce the order of the hole inputs.");
 
  if (!triangulator_hole_ptrs.empty())
    poly2tri.attach_hole_list(&triangulator_hole_ptrs);
 
  if (xyd_opts.desired_area_func != "")
  {
    // poly2tri will clone this so it's fine going out of scope
    libMesh::ParsedFunction<Real> area_func{xyd_opts.desired_area_func};
    poly2tri.set_desired_area_function(&area_func);
  }
  else if (xyd_opts.use_auto_area_func)
  {
    poly2tri.set_auto_area_function(
        mg.comm(),
        xyd_opts.auto_area_function_num_points,
        xyd_opts.auto_area_function_power,
        xyd_opts.auto_area_func_default_size > 0.0 ? xyd_opts.auto_area_func_default_size : 0.0,
        xyd_opts.auto_area_func_default_size_dist > 0.0 ? xyd_opts.auto_area_func_default_size_dist
                                                        : -1.0);
  }
 
  if (xyd_opts.tri_elem_type == "TRI6")
    poly2tri.elem_type() = libMesh::ElemType::TRI6;
  else if (xyd_opts.tri_elem_type == "TRI7")
  {
    poly2tri.elem_type() = libMesh::ElemType::TRI7;
    // Snapping boundary mid-edge nodes onto a curved input boundary leaves the
    // TRI7 interior node at the stale straight-edge centroid. Ask libMesh to
    // move it to the curved-mapping centroid; this is also what makes
    // libMesh's Tri6-based tangling check valid for TRI7.
    poly2tri.set_fixup_tri7_center_nodes(true);
  }
  // Add interior points before triangulating. Only points inside the boundaries
  // will be meshed.
  for (const auto & point : xyd_opts.interior_points)
    mesh->add_point(point);
 
  poly2tri.triangulate();
 
  finalizeTriangulation(mg, *mesh, hole_meshes, holes_with_midpoints, xyd_opts);
 
  return mesh;
}

Referenced by XYFrontalDelaunayGenerator::buildBackgroundMesh(), BoundaryLayerUtils::buildBoundaryLayerRing(), and XYDelaunayGenerator::generate().