BoundaryLayerUtils Namespace Reference

Functions
std::unique_ptr< MeshBase >	buildBoundaryLayerRing (MeshGenerator &mg, MeshBase &input_mesh, const std::vector< BoundaryName > &boundary_names, unsigned int num_layers, Real thickness, Real layer_bias, bool outward, const MooseEnum &tri_elem_type, SubdomainID output_subdomain_id, const SubdomainName &output_subdomain_name)
	Builds a conformal boundary-layer ring of triangulated annuli along a boundary of an input 2D mesh (or a 1D loop).
std::vector< Point >	generateOffsetPolyline (MeshGenerator *mg, std::unique_ptr< libMesh::UnstructuredMesh > &ply_mesh_u, std::vector< Point > &points, std::vector< Point > &mid_points, const bool outward, const Real thickness)
	Generates a list of points offset from the input boundary polyline by a specified thickness in either outward/inward direction.
void	collectExteriorVertexPointsFromMesh (libMesh::TriangulatorInterface::MeshedHole &bdry_mh, std::vector< Point > &points, std::vector< Point > &mid_points, const bool skip_node_reduction=false)
	Collects key vertex points (and optional midpoints) from a meshed hole, optionally discarding colinear vertices.
Point	getKeyNormal (const Elem *elem, const unsigned int s, const unsigned int node_index)
	Extracts the normal vector of the EDGE3 side of a quadratic element at a given node index.

Function Documentation

buildBoundaryLayerRing()

std::unique_ptr< MeshBase > BoundaryLayerUtils::buildBoundaryLayerRing	(	MeshGenerator &	mg,
		MeshBase &	input_mesh,
		const std::vector< BoundaryName > &	boundary_names,
		unsigned int	num_layers,
		Real	thickness,
		Real	layer_bias,
		bool	outward,
		const MooseEnum &	tri_elem_type,
		SubdomainID	output_subdomain_id,
		const SubdomainName &	output_subdomain_name
	)

Builds a conformal boundary-layer ring of triangulated annuli along a boundary of an input 2D mesh (or a 1D loop).

Generates num_layers + 1 parallel polylines (geometric progression of thicknesses with the given bias), triangulates each annulus with triangulateWithDelaunay, and sequentially stitches them. The output mesh carries 2 * num_layers boundary ids, with the innermost being id 1 and the outermost being id (num_layers - 1) * 2.

Parameters

mg	The calling mesh generator (for paramError reporting + buildMeshBaseObject)
input_mesh	The 2D-XY input mesh or 1D closed loop providing the seed boundary
boundary_names	Subset of boundary names on input_mesh defining the seed boundary; if empty, the external boundary of input_mesh is auto-detected via MeshedHole
num_layers	Number of element layers to generate
thickness	Total boundary-layer thickness
layer_bias	Geometric growth factor between successive layer thicknesses (1.0 = uniform)
outward	If true, the layer grows outward from the seed boundary; else inward
tri_elem_type	Triangle element type ("TRI3", "TRI6", "TRI7", or "DEFAULT")
output_subdomain_id	Subdomain id assigned to all generated triangles (0 = default)
output_subdomain_name	Subdomain name assigned to output_subdomain_id (empty = unnamed)

Returns

The stitched boundary-layer ring mesh

Definition at line 30 of file BoundaryLayerUtils.C.

{
  // Extract seed boundary polyline points (vertices + optional midpoints) from input_mesh.
  std::set<std::size_t> bdry_id_set;
  if (!boundary_names.empty())
  {
    auto ids =
        MooseMeshUtils::getBoundaryIDs(input_mesh, boundary_names, /*generate unknown*/ false);
    bdry_id_set.insert(ids.begin(), ids.end());
  }
  TriangulatorInterface::MeshedHole bdry_mh(input_mesh, bdry_id_set);
  std::vector<Point> cur_pts;
  std::vector<Point> cur_mids;
  // Preserve every input boundary node (including colinear interior nodes on straight edges) so
  // the ring's innermost polyline can be stitched back to the input mesh's boundary exactly when
  // keep_input is requested by a downstream caller.
  collectExteriorVertexPointsFromMesh(bdry_mh,
                                      cur_pts,
                                      cur_mids,
                                      /*skip_node_reduction=*/true);
 
  // Geometric progression of incremental thicknesses. layer_thicknesses[i] is the offset
  // distance from polyline i-1 to polyline i (for i >= 1); layer_thicknesses[0] is unused.
  std::vector<Real> layer_thicknesses(num_layers + 1);
  const Real unit_thickness =
      (layer_bias == 1.0)
          ? (thickness / num_layers)
          : (thickness / (std::pow(layer_bias, num_layers) - 1.0) * (layer_bias - 1.0));
  layer_thicknesses[0] = 0.0;
  for (auto i : make_range(std::vector<Real>::size_type(1), layer_thicknesses.size()))
    layer_thicknesses[i] = unit_thickness * std::pow(layer_bias, i - 1);
 
  // Generate the N+1 polyline meshes. Polylines are indexed so that polyline 0 is geometrically
  // innermost (smallest) and polyline N is outermost (largest). Iteration order depends on
  // direction: for outward, build polyline 0 first (= input boundary); for inward, build polyline
  // N first.
  std::vector<std::unique_ptr<MeshBase>> polylines(num_layers + 1);
  for (auto layer_i : make_range(num_layers + 1))
  {
    const unsigned int layer_index = outward ? layer_i : (num_layers - layer_i);
 
    // Build the polyline mesh for storage as polylines[layer_index].
    auto ply = std::make_unique<ReplicatedMesh>(mg.comm());
    MooseMeshUtils::buildPolyLineMesh(*ply,
                                      cur_pts,
                                      cur_mids,
                                      /*loop=*/true,
                                      BoundaryName(),
                                      BoundaryName(),
                                      std::vector<unsigned int>({1}));
    polylines[layer_index] = std::move(ply);
 
    if (layer_i + 1 < num_layers + 1)
    {
      // Build a sacrificial polyline mesh to feed generateOffsetPolyline (it triangulates the
      // input mesh in place to compute side normals; we discard it afterwards).
      auto ply_for_offset = std::make_unique<ReplicatedMesh>(mg.comm());
      MooseMeshUtils::buildPolyLineMesh(*ply_for_offset,
                                        cur_pts,
                                        cur_mids,
                                        /*loop=*/true,
                                        BoundaryName(),
                                        BoundaryName(),
                                        std::vector<unsigned int>({1}));
      std::unique_ptr<UnstructuredMesh> ply_for_offset_u =
          dynamic_pointer_cast<UnstructuredMesh>(std::move(ply_for_offset));
 
      std::vector<Point> next_combined = generateOffsetPolyline(
          &mg, ply_for_offset_u, cur_pts, cur_mids, outward, layer_thicknesses[layer_i + 1]);
      if (cur_mids.empty())
        cur_pts = std::move(next_combined);
      else
      {
        const auto n_vert = cur_pts.size();
        cur_pts.assign(next_combined.begin(), next_combined.begin() + n_vert);
        cur_mids.assign(next_combined.begin() + n_vert, next_combined.end());
      }
    }
  }
 
  // Triangulate each annulus (between polylines[i] and polylines[i+1]). Stitch each annulus to
  // the accumulating ring (except the first one).
  std::unique_ptr<MeshBase> ring;
  for (auto i : make_range(num_layers))
  {
    MeshTriangulationUtils::XYDelaunayOptions xyd_opts;
    xyd_opts.refine_bdy = false;
    xyd_opts.verify_holes = false;
    xyd_opts.stitch_holes = {i > 0};
    xyd_opts.refine_holes = {false};
    xyd_opts.tri_elem_type = std::string(tri_elem_type);
    if (output_subdomain_id != 0)
    {
      xyd_opts.has_output_subdomain_id = true;
      xyd_opts.output_subdomain_id = output_subdomain_id;
    }
    if (output_subdomain_name.size())
    {
      xyd_opts.has_output_subdomain_name = true;
      xyd_opts.output_subdomain_name = output_subdomain_name;
    }
 
    std::vector<std::unique_ptr<MeshBase>> holes;
    holes.reserve(1);
    if (i == 0)
      holes.push_back(std::move(polylines[0]));
    else
      holes.push_back(std::move(ring));
 
    ring = MeshTriangulationUtils::triangulateWithDelaunay(
        mg, std::move(polylines[i + 1]), std::move(holes), xyd_opts);
  }
  // We now have 2 * num_layers boundaries
  // Let's only keep the innermost (1) and outermost (2 * num_layers) boundaries, and remove all
  // intermediate ring bcids
  std::vector<BoundaryID> bids_to_delete;
  for (const auto b : make_range(2 * num_layers))
  {
    if (b == 1 || b == (num_layers - 1) * 2)
      continue;
    bids_to_delete.push_back(b);
  }
  auto & bi = ring->get_boundary_info();
  for (auto b : bids_to_delete)
    bi.remove_id(b);
 
  return ring;
}

Referenced by XYDelaunayGenerator::generate(), and XYTriangleBoundaryLayerGenerator::generate().

collectExteriorVertexPointsFromMesh()

void BoundaryLayerUtils::collectExteriorVertexPointsFromMesh	(	libMesh::TriangulatorInterface::MeshedHole &	bdry_mh,
		std::vector< Point > &	points,
		std::vector< Point > &	mid_points,
		const bool	skip_node_reduction = false
	)

Collects key vertex points (and optional midpoints) from a meshed hole, optionally discarding colinear vertices.

Parameters

bdry_mh	The 2D MeshedHole object from which to collect key points
points	The vector to which collected vertex points are appended (in order)
mid_points	The vector to which collected midpoints are appended; only populated when the MeshedHole reports a midpoint per side AND skip_node_reduction is true
skip_node_reduction	If true, retain every vertex and midpoint; if false, drop vertices that are colinear with their two neighbors and do not collect midpoints

Definition at line 321 of file BoundaryLayerUtils.C.

{
  for (const auto i : make_range(bdry_mh.n_points()))
  {
    if (skip_node_reduction || !geom_utils::arePointsColinear(
                                   bdry_mh.point((i - 1 + bdry_mh.n_points()) % bdry_mh.n_points()),
                                   bdry_mh.point(i),
                                   bdry_mh.point((i + 1) % bdry_mh.n_points())))
    {
      points.push_back(bdry_mh.point(i));
      if (bdry_mh.n_midpoints() == 1 && skip_node_reduction)
        mid_points.push_back(bdry_mh.midpoint(0, i));
    }
  }
}

Referenced by buildBoundaryLayerRing(), and generateOffsetPolyline().

generateOffsetPolyline()

std::vector< Point > BoundaryLayerUtils::generateOffsetPolyline	(	MeshGenerator *	mg,
		std::unique_ptr< libMesh::UnstructuredMesh > &	ply_mesh_u,
		std::vector< Point > &	points,
		std::vector< Point > &	mid_points,
		const bool	outward,
		const Real	thickness
	)

Generates a list of points offset from the input boundary polyline by a specified thickness in either outward/inward direction.

Parameters

mg	The mesh generator calling this function, used for paramError reporting
ply_mesh_u	The 1D loop polyline mesh of the original boundary; the volume it encloses will be triangulated in-place to compute the normal and define the inward/outward directions
points	The vertex points of the original polyline. If empty, populated from the input mesh via collectExteriorVertexPointsFromMesh.
mid_points	The midpoints of the original polyline (optional to enable quadratic elements). If both this and points are empty, populated from the input mesh.
outward	Whether to offset in the outward (true) or inward (false) direction
thickness	The offset distance (>= 0)

Returns

Offset points: vertices first, then midpoints (length = points.size() + mid_points.size())

Definition at line 169 of file BoundaryLayerUtils.C.

{
  // If the input points are empty, we will extract them from the input 1D mesh
  if (points.empty())
  {
    mooseAssert(mid_points.empty(),
                "If the input points are empty, the input mid_points must be also empty.");
 
    TriangulatorInterface::MeshedHole bdry_mh(*ply_mesh_u);
    collectExteriorVertexPointsFromMesh(bdry_mh, points, mid_points);
  }
  // Generate a very simple triangulation mesh so that we can get the outward normal vectors
  libMesh::Poly2TriTriangulator poly2tri(*ply_mesh_u);
  poly2tri.triangulation_type() = libMesh::TriangulatorInterface::PSLG;
 
  poly2tri.set_interpolate_boundary_points(0);
  poly2tri.set_refine_boundary_allowed(false);
  poly2tri.set_verify_hole_boundaries(false);
  poly2tri.desired_area() = 0;
  poly2tri.minimum_angle() = 0; // Not yet supported
  poly2tri.smooth_after_generating() = false;
  if (mid_points.size())
    poly2tri.elem_type() = libMesh::ElemType::TRI6;
  poly2tri.triangulate();
 
  // We need to serialize the mesh for next steps
  libMesh::MeshSerializer serial(*ply_mesh_u);
  // The mesh now only contains one side set that corresponds to the outer boundary with an ID of 0
  auto bdry_list(ply_mesh_u->get_boundary_info().build_side_list());
 
  // For each vertex, the shifting direction to form the offset is defined by the normal vectors of
  // the two sides that contain the vertex.
  // We gather the normals for each node on the boundary here, which are all vertices because of
  // the pre-selection of nodes.
  std::map<dof_id_type, std::vector<Point>> node_normal_map;
  std::map<dof_id_type, Point> mid_node_normal_map;
  for (const auto & bside : bdry_list)
  {
    const auto & side = ply_mesh_u->elem_ptr(std::get<0>(bside))->side_ptr(std::get<1>(bside));
    // For linear elements, the side normal is constant and can be obtained straightforwardly;
    // For quadratic elements, we need the normal at the shared vertices
    const Point side_normal_0 =
        mid_points.size()
            ? getKeyNormal(ply_mesh_u->elem_ptr(std::get<0>(bside)), std::get<1>(bside), 0)
            : ply_mesh_u->elem_ptr(std::get<0>(bside))
                  ->side_vertex_average_normal(std::get<1>(bside));
    const Point side_normal_1 =
        mid_points.size()
            ? getKeyNormal(ply_mesh_u->elem_ptr(std::get<0>(bside)), std::get<1>(bside), 1)
            : side_normal_0;
 
    if (node_normal_map.count(side->node_ptr(0)->id()))
      node_normal_map[side->node_ptr(0)->id()].push_back(side_normal_0);
    else
      node_normal_map[side->node_ptr(0)->id()] = {side_normal_0};
    if (node_normal_map.count(side->node_ptr(1)->id()))
      node_normal_map[side->node_ptr(1)->id()].push_back(side_normal_1);
    else
      node_normal_map[side->node_ptr(1)->id()] = {side_normal_1};
 
    if (mid_points.size())
    {
      const Point mid_node_normal =
          getKeyNormal(ply_mesh_u->elem_ptr(std::get<0>(bside)), std::get<1>(bside), 2);
      mid_node_normal_map
          [ply_mesh_u->elem_ptr(std::get<0>(bside))->node_ptr(std::get<1>(bside) + 3)->id()] =
              mid_node_normal;
    }
  }
 
  std::vector<Point> mod_reduced_pts_list(points);
  for (const auto & [node_id, normal_vecs] : node_normal_map)
  {
    mooseAssert(normal_vecs.size() == 2,
                "Each vertex should be connected to exactly two sides in a polygon.");
 
    const Point original_pt = *(ply_mesh_u->node_ptr(node_id));
    // Form an average normal at the vertex from the two connected sides' normals
    const Point move_dir =
        (normal_vecs.front() + normal_vecs.back()).unit() * (outward ? 1.0 : -1.0);
    // Consider four points of interest to determine the moving distance
    // 1. the vertex point
    // 2. point along normal_vecs.front() from the vertex with a distance thickness
    // 3. point along normal_vecs.back() from the vertex with a distance thickness
    // 4. point along move_dir from the vertex with a distance mov_dist
    // The four points form a kite shape with its symmetry axis along 1-4 direction
    // Angle 1-2-4 and angle 1-3-4 are right angles
    // Angle 2-1-3 (i.e., theta) can be calculated using the interior product of
    // v1 = normal_vecs.front() and v2 = normal_vecs.back():
    // cos(theta) = (v1 . v2) / (|v1| * |v2|)
    // Because of the right angles, the distance between 1 and 4 can be calculated as:
    // mov_dist = thickness / cos(theta/2)
    // and cos(theta/2) = sqrt((1 + cos(theta)) / 2)
    const Real mov_dist =
        thickness / std::sqrt((1.0 + (normal_vecs.front() * normal_vecs.back()) /
                                         (normal_vecs.front().norm() * normal_vecs.back().norm())) /
                              2.0);
    mooseAssert(std::count(points.begin(), points.end(), original_pt) == 1,
                "The original point should be found exactly once in the reduced points list.");
    mod_reduced_pts_list[std::distance(points.begin(),
                                       std::find(points.begin(), points.end(), original_pt))] =
        original_pt + move_dir * mov_dist;
  }
 
  // To ensure no overlapping, we need to check set of four points
  // p1 and p2 should be a pair of points before and after shifting
  // p3 and p4 should be a pair of adjacent shifted points that are neither not p2
  for (const auto & i_node_1 : make_range(mod_reduced_pts_list.size()))
  {
    const Point & p1 = points[i_node_1];
    const Point & p2 = mod_reduced_pts_list[i_node_1];
    for (const auto & i_node_2 : make_range(mod_reduced_pts_list.size()))
    {
      if (i_node_2 == i_node_1 || (i_node_2 + 1) % mod_reduced_pts_list.size() == i_node_1)
        continue;
      const Point & p3 = mod_reduced_pts_list[i_node_2];
      const Point & p4 = mod_reduced_pts_list[(i_node_2 + 1) % mod_reduced_pts_list.size()];
      if (thickness > 0)
        if (geom_utils::segmentsIntersect(p1, p2, p3, p4))
          mg->paramError(
              "thickness",
              "The thickness is so large that the mesh is tangled because the offset nodes "
              "are no longer in the same order when following the original boundary. Please "
              "reduce the thickness value.");
    }
  }
 
  std::vector<Point> mid_mod_reduced_pts_list(mid_points.size());
  for (const auto & [node_id, normal_vec] : mid_node_normal_map)
  {
    const Point original_pt = *(ply_mesh_u->node_ptr(node_id));
    const Point move_dir = normal_vec.unit() * (outward ? 1.0 : -1.0);
    mid_mod_reduced_pts_list[std::distance(
        mid_points.begin(), std::find(mid_points.begin(), mid_points.end(), original_pt))] =
        original_pt + move_dir * thickness;
  }
 
  // combine mod_reduced_pts_list and mid_mod_reduced_pts_list to get the final list of points for
  // the layer mesh
  std::vector<Point> layer_pts_list(mod_reduced_pts_list);
  layer_pts_list.insert(
      layer_pts_list.end(), mid_mod_reduced_pts_list.begin(), mid_mod_reduced_pts_list.end());
 
  return layer_pts_list;
}

Referenced by buildBoundaryLayerRing().

getKeyNormal()

Point BoundaryLayerUtils::getKeyNormal	(	const Elem *	elem,
		const unsigned int	s,
		const unsigned int	node_index
	)

Extracts the normal vector of the EDGE3 side of a quadratic element at a given node index.

Parameters

elem	The element containing the side of interest
s	The side index of interest
node_index	The index of the node on the side at which to extract the normal vector

Returns

The normal vector at the node index on the side of interest

Definition at line 341 of file BoundaryLayerUtils.C.

{
  const std::unique_ptr<const Elem> face = elem->build_side_ptr(s);
  mooseAssert(face->type() == ElemType::EDGE3,
              "Only elements with EDGE3 sides are supported in this function.");
  mooseAssert(node_index < 3,
              "The node index for an EDGE3 side should be 0, 1, or 2 (for the two vertices and the "
              "midpoint).");
  std::unique_ptr<libMesh::FEBase> fe(
      libMesh::FEBase::build(2, libMesh::FEType(elem->default_order())));
  const std::vector<Point> & normals = fe->get_normals();
  std::vector<Point> ref_pts = {face->reference_elem()->point(node_index)};
  fe->reinit(elem, s, TOLERANCE, &ref_pts);
  return normals[0];
}

Referenced by generateOffsetPolyline().