MooseMeshXYCuttingUtils Namespace Reference

Functions
void	lineRemoverMoveNode (libMesh::ReplicatedMesh &mesh, const std::vector< Real > &bdry_pars, const subdomain_id_type block_id_to_remove, const std::set< subdomain_id_type > &subdomain_ids_set, const boundary_id_type trimming_section_boundary_id, const boundary_id_type external_boundary_id, const std::vector< boundary_id_type > &other_boundaries_to_conform=std::vector< boundary_id_type >(), const bool assign_ext_to_new=false, const bool side_to_remove=true)
	Removes all the elements on one side of a given line and deforms the elements intercepted by the line to form a flat new boundary. More...
bool	lineSideDeterminator (const Real px, const Real py, const Real param_1, const Real param_2, const Real param_3, const bool direction_param, const Real dis_tol=libMesh::TOLERANCE)
	Determines whether a point on XY-plane is on the side of a given line that needs to be removed. More...
Point	twoLineIntersection (const Real param_11, const Real param_12, const Real param_13, const Real param_21, const Real param_22, const Real param_23)
	Calculates the intersection Point of two given straight lines. More...
Point	twoPointandLineIntersection (const Point &pt1, const Point &pt2, const Real param_1, const Real param_2, const Real param_3)
	Calculates the intersection Point of a straight line defined by two given points and another straight line. More...
bool	quasiTriElementsFixer (libMesh::ReplicatedMesh &mesh, const std::set< subdomain_id_type > &subdomain_ids_set, const subdomain_id_type tri_elem_subdomain_shift=Moose::INVALID_BLOCK_ID, const SubdomainName tri_elem_subdomain_name_suffix="tri")
	Fixes degenerate QUAD elements created by the hexagonal mesh trimming by converting them into TRI elements. More...
std::vector< std::pair< Real, unsigned int > >	vertex_angles (const Elem &elem)
	Calculates the internal angles of a given 2D element. More...
std::vector< std::pair< Real, unsigned int > >	vertex_distances (const Elem &elem)
	Calculates the distances between the vertices of a given 2D element. More...
void	triElemSplitter (libMesh::ReplicatedMesh &mesh, const dof_id_type elem_id, const unsigned short node_shift, const dof_id_type nid_3, const dof_id_type nid_4, const subdomain_id_type single_elem_side_id, const subdomain_id_type double_elem_side_id)
	Split a TRI3 element into three TRI3 elements based on two nodes on the two sides of the triangle. More...
void	triElemSplitter (libMesh::ReplicatedMesh &mesh, const dof_id_type elem_id, const unsigned short node_shift, const dof_id_type nid_m, const subdomain_id_type first_elem_side_id, const subdomain_id_type second_elem_side_id)
	Split a TRI3 element into two TRI3 elements based on one node on one side of the triangle. More...
void	quadElemSplitter (libMesh::ReplicatedMesh &mesh, const dof_id_type elem_id, const subdomain_id_type tri_elem_subdomain_shift)
	Split a QUAD4 element into two TRI3 elements. More...
void	quadToTriOnLine (libMesh::ReplicatedMesh &mesh, const std::vector< Real > &cut_line_params, const dof_id_type tri_subdomain_id_shift, const SubdomainName tri_elem_subdomain_name_suffix)
	Convert all the QUAD4 elements in the mesh that are crossed by the given line into TRI3 elements. More...
void	lineRemoverCutElemTri (libMesh::ReplicatedMesh &mesh, const std::vector< Real > &cut_line_params, const subdomain_id_type block_id_to_remove, const boundary_id_type new_boundary_id)
	Trim the 2D mesh by removing all the elements on one side of the given line. More...
void	lineRemoverCutElem (libMesh::ReplicatedMesh &mesh, const std::vector< Real > &cut_line_params, const dof_id_type tri_subdomain_id_shift, const SubdomainName tri_elem_subdomain_name_suffix, const subdomain_id_type block_id_to_remove, const boundary_id_type new_boundary_id, const bool improve_boundary_tri_elems=false)
	Trim the 2D mesh by removing all the elements on one side of the given line. More...
void	boundaryTriElemImprover (libMesh::ReplicatedMesh &mesh, const boundary_id_type boundary_to_improve)
	Improve the element quality of the boundary TRI3 elements of the given boundary. More...
void	makeImprovedTriElement (libMesh::ReplicatedMesh &mesh, const dof_id_type node_id_0, const dof_id_type node_id_1, const dof_id_type node_id_2, const subdomain_id_type subdomain_id, const std::vector< dof_id_type > &extra_elem_ids, const std::vector< boundary_id_type > &boundary_ids_for_side_1=std::vector< boundary_id_type >(), const std::vector< boundary_id_type > &boundary_ids_for_side_0=std::vector< boundary_id_type >(), const std::vector< boundary_id_type > &boundary_ids_for_side_2=std::vector< boundary_id_type >())
	Make a TRI3 element with the given node ids and subdomain id with boundary information. More...
bool	elemSideLocator (libMesh::ReplicatedMesh &mesh, const dof_id_type elem_id, const dof_id_type node_id_0, const dof_id_type node_id_1, unsigned short &side_id, bool &is_inverse)
	Check if there is a side in an element that contains the given pair of nodes; if yes, also find the side id and the direction of the two nodes in the side. More...
Point	twoPointandLineIntersection (const Point &pt1, const Point &pt2, const Real param_1, const Real param_2, const Real param_3)
std::vector< std::pair< Real, unsigned int > >	vertex_angles (const Elem &elem)
std::vector< std::pair< Real, unsigned int > >	vertex_distances (const Elem &elem)

Function Documentation

boundaryTriElemImprover()

void MooseMeshXYCuttingUtils::boundaryTriElemImprover	(	libMesh::ReplicatedMesh &	mesh,
		const boundary_id_type	boundary_to_improve
	)

Improve the element quality of the boundary TRI3 elements of the given boundary.

Parameters

mesh	input mesh with the boundary TRI3 elements that need to be improved
boundary_to_improve	boundary id of the boundary that needs to be improved

Definition at line 981 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElem().

{
  if (!MooseMeshUtils::hasBoundaryID(mesh, boundary_to_improve))
    mooseError(
        "MooseMeshXYCuttingUtils::boundaryTriElemImprover(): The boundary_to_improve provided "
        "does not exist in the given mesh.");
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto side_list = boundary_info.build_side_list();
  // Here we would like to collect the following information for all the TRI3 elements on the
  // boundary: Key: node id of the off-boundary node Value: a vector of tuples, each tuple contains
  // the following information:
  // 1. The element id of the element that is on the boundary to improve
  // 2. the one node id of that element that is on the boundary to improve
  // 3. the other node id of the element that is on the boundary to improve
  std::map<dof_id_type, std::vector<std::tuple<dof_id_type, dof_id_type, dof_id_type>>>
      tri3_elem_info;
  for (const auto & side : side_list)
  {
    if (std::get<2>(side) == boundary_to_improve)
    {
      Elem * elem = mesh.elem_ptr(std::get<0>(side));
      if (elem->type() == TRI3)
      {
        const auto key_node_id = elem->node_id((std::get<1>(side) + 2) % 3);
        const auto value_elem_id = elem->id();
        const auto value_node_id_1 = elem->node_id(std::get<1>(side));
        const auto value_node_id_2 = elem->node_id((std::get<1>(side) + 1) % 3);
        tri3_elem_info[key_node_id].push_back(
            std::make_tuple(value_elem_id, value_node_id_1, value_node_id_2));
      }
    }
  }
  // Elements that need to be removed
  std::vector<dof_id_type> elems_to_remove;
  // Now check if any group of TRI3 sharing an off-boundary node can be improved.
  for (const auto & tri_group : tri3_elem_info)
  {
    // It is possible to improve only when more than one TRI3 elements share the same off-boundary
    // node
    std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
    std::vector<dof_id_type> elem_id_list;
    for (const auto & tri : tri_group.second)
    {
      node_assm.push_back(std::make_pair(std::get<1>(tri), std::get<2>(tri)));
      elem_id_list.push_back(std::get<0>(tri));
    }
    std::vector<dof_id_type> ordered_node_list;
    std::vector<dof_id_type> ordered_elem_list;
    MooseMeshUtils::makeOrderedNodeList(
        node_assm, elem_id_list, ordered_node_list, ordered_elem_list);

    // For all the elements sharing the same off-boundary node, we need to know how many separated
    // subdomains are involved
    // If there are extra element ids defined on the mesh, they also want to retain their boundaries
    // Only triangle elements that share a side can be merged
    const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
    std::vector<std::tuple<subdomain_id_type, std::vector<dof_id_type>, unsigned int>> blocks_info;
    for (const auto & elem_id : ordered_elem_list)
    {
      std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
      // Record all the element extra integers of the original quad element
      for (unsigned int j = 0; j < n_elem_extra_ids; j++)
        exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
      if (!blocks_info.empty())
      {
        if (mesh.elem_ptr(elem_id)->subdomain_id() == std::get<0>(blocks_info.back()) &&
            exist_extra_ids == std::get<1>(blocks_info.back()))
        {
          std::get<2>(blocks_info.back())++;
          continue;
        }
      }
      blocks_info.push_back(
          std::make_tuple(mesh.elem_ptr(elem_id)->subdomain_id(), exist_extra_ids, 1));
    }
    // For each separated subdomain / set of extra ids, we try to improve the boundary elements
    unsigned int side_counter = 0;
    for (const auto & block_info : blocks_info)
    {
      const auto node_1 = mesh.node_ptr(ordered_node_list[side_counter]);
      // we do not need to subtract 1 for node_2
      const auto node_2 = mesh.node_ptr(ordered_node_list[side_counter + std::get<2>(block_info)]);
      const auto node_0 = mesh.node_ptr(tri_group.first);
      const Point v1 = *node_1 - *node_0;
      const Point v2 = *node_2 - *node_0;
      const Real angle = std::acos(v1 * v2 / v1.norm() / v2.norm()) / M_PI * 180.0;
      const std::vector<dof_id_type> block_elems(ordered_elem_list.begin() + side_counter,
                                                 ordered_elem_list.begin() + side_counter +
                                                     std::get<2>(block_info));
      // We assume that there are no sidesets defined inside a subdomain
      // For the first TRI3 element, we want to check if its side defined by node_0 and node_1 is
      // defined in any sidesets
      unsigned short side_id_0;
      unsigned short side_id_t;
      bool is_inverse_0;
      bool is_inverse_t;
      elemSideLocator(mesh,
                      block_elems.front(),
                      tri_group.first,
                      ordered_node_list[side_counter],
                      side_id_0,
                      is_inverse_0);
      elemSideLocator(mesh,
                      block_elems.back(),
                      ordered_node_list[side_counter + std::get<2>(block_info)],
                      tri_group.first,
                      side_id_t,
                      is_inverse_t);
      // Collect boundary information of the identified sides
      std::vector<boundary_id_type> side_0_boundary_ids;
      boundary_info.boundary_ids(
          mesh.elem_ptr(block_elems.front()), side_id_0, side_0_boundary_ids);
      std::vector<boundary_id_type> side_t_boundary_ids;
      boundary_info.boundary_ids(mesh.elem_ptr(block_elems.back()), side_id_t, side_t_boundary_ids);

      // Ideally we want this angle to be 60 degrees
      // In reality, we want one TRI3 element if the angle is less than 90 degrees;
      // we want two TRI3 elements if the angle is greater than 90 degrees and less than 135
      // degrees; we want three TRI3 elements if the angle is greater than 135 degrees and less than
      // 180 degrees.
      if (angle < 90.0)
      {
        if (std::get<2>(block_info) > 1)
        {
          makeImprovedTriElement(mesh,
                                 tri_group.first,
                                 ordered_node_list[side_counter],
                                 ordered_node_list[side_counter + std::get<2>(block_info)],
                                 std::get<0>(block_info),
                                 std::get<1>(block_info),
                                 {boundary_to_improve},
                                 side_0_boundary_ids,
                                 side_t_boundary_ids);
          elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
        }
      }
      else if (angle < 135.0)
      {
        // We can just add the middle node because there's nothing on the other side
        const auto node_m = mesh.add_point((*node_1 + *node_2) / 2.0);
        makeImprovedTriElement(mesh,
                               tri_group.first,
                               ordered_node_list[side_counter],
                               node_m->id(),
                               std::get<0>(block_info),
                               std::get<1>(block_info),
                               {boundary_to_improve},
                               side_0_boundary_ids,
                               std::vector<boundary_id_type>());
        makeImprovedTriElement(mesh,
                               tri_group.first,
                               node_m->id(),
                               ordered_node_list[side_counter + std::get<2>(block_info)],
                               std::get<0>(block_info),
                               std::get<1>(block_info),
                               {boundary_to_improve},
                               std::vector<boundary_id_type>(),
                               side_t_boundary_ids);
        elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
      }
      else
      {
        const auto node_m1 = mesh.add_point((*node_1 * 2.0 + *node_2) / 3.0);
        const auto node_m2 = mesh.add_point((*node_1 + *node_2 * 2.0) / 3.0);
        makeImprovedTriElement(mesh,
                               tri_group.first,
                               ordered_node_list[side_counter],
                               node_m1->id(),
                               std::get<0>(block_info),
                               std::get<1>(block_info),
                               {boundary_to_improve},
                               side_0_boundary_ids,
                               std::vector<boundary_id_type>());
        makeImprovedTriElement(mesh,
                               tri_group.first,
                               node_m1->id(),
                               node_m2->id(),
                               std::get<0>(block_info),
                               std::get<1>(block_info),
                               {boundary_to_improve},
                               std::vector<boundary_id_type>(),
                               std::vector<boundary_id_type>());
        makeImprovedTriElement(mesh,
                               tri_group.first,
                               node_m2->id(),
                               ordered_node_list[side_counter + std::get<2>(block_info)],
                               std::get<0>(block_info),
                               std::get<1>(block_info),
                               {boundary_to_improve},
                               std::vector<boundary_id_type>(),
                               side_t_boundary_ids);
        elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
      }
      side_counter += std::get<2>(block_info);
    }
    // TODO: Need to check if the new element is inverted?
  }
  // Delete the original elements
  for (const auto & elem_to_remove : elems_to_remove)
    mesh.delete_elem(mesh.elem_ptr(elem_to_remove));
  mesh.contract();
}

elemSideLocator()

bool MooseMeshXYCuttingUtils::elemSideLocator	(	libMesh::ReplicatedMesh &	mesh,
		const dof_id_type	elem_id,
		const dof_id_type	node_id_0,
		const dof_id_type	node_id_1,
		unsigned short &	side_id,
		bool &	is_inverse
	)

Check if there is a side in an element that contains the given pair of nodes; if yes, also find the side id and the direction of the two nodes in the side.

Parameters

mesh	input mesh with the element that needs to be checked
elem_id	id of the element that needs to be checked
node_id_0	id of the first node of the pair
node_id_1	id of the second node of the pair
side_id	id of the side that contains the pair of nodes
is_inverse	flag to indicate if the two nodes are in the same direction as the side

Returns

true if the element contains the side with the given pair of nodes

Definition at line 1215 of file MooseMeshXYCuttingUtils.C.

Referenced by boundaryTriElemImprover().

{
  Elem * elem = mesh.elem_ptr(elem_id);
  for (unsigned short i = 0; i < elem->n_sides(); i++)
  {
    if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_0 &&
        elem->side_ptr(i)->node_ptr(1)->id() == node_id_1)
    {
      side_id = i;
      is_inverse = false;
      return true;
    }
    else if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_1 &&
             elem->side_ptr(i)->node_ptr(1)->id() == node_id_0)
    {
      side_id = i;
      is_inverse = true;
      return true;
    }
  }
  return false;
}

lineRemoverCutElem()

void MooseMeshXYCuttingUtils::lineRemoverCutElem	(	libMesh::ReplicatedMesh &	mesh,
		const std::vector< Real > &	cut_line_params,
		const dof_id_type	tri_subdomain_id_shift,
		const SubdomainName	tri_elem_subdomain_name_suffix,
		const subdomain_id_type	block_id_to_remove,
		const boundary_id_type	new_boundary_id,
		const bool	improve_boundary_tri_elems = false
	)

Trim the 2D mesh by removing all the elements on one side of the given line.

Note that the mesh can only contain QUAD4 and TRI3 elements

Parameters

mesh	input mesh that need to be trimmed
cut_line_params	parameters of the line that cuts the input mesh
tri_subdomain_id_shift	subdomain id shift used to define the TRI element subdomains formed due to the trimming
tri_elem_subdomain_name_suffix	suffix used to name the TRI element subdomains formed due to the trimming
block_id_to_remove	a temporary subdomain id used to mark the elements that need to be removed
new_boundary_id	boundary id of the new boundary that forms due to the trimming
improve_boundary_tri_elems	flag to indicate whether the boundary TRI3 elements need to be improved

Definition at line 962 of file MooseMeshXYCuttingUtils.C.

Referenced by XYMeshLineCutter::generate().

{
  // Convert any quad elements crossed by the line into tri elements
  quadToTriOnLine(mesh, cut_line_params, tri_subdomain_id_shift, tri_elem_subdomain_name_suffix);
  // Then do the cutting for the preprocessed mesh that only contains tri elements crossed by the
  // cut line
  lineRemoverCutElemTri(mesh, cut_line_params, block_id_to_remove, new_boundary_id);

  if (improve_boundary_tri_elems)
    boundaryTriElemImprover(mesh, new_boundary_id);
}

lineRemoverCutElemTri()

void MooseMeshXYCuttingUtils::lineRemoverCutElemTri	(	libMesh::ReplicatedMesh &	mesh,
		const std::vector< Real > &	cut_line_params,
		const subdomain_id_type	block_id_to_remove,
		const boundary_id_type	new_boundary_id
	)

Trim the 2D mesh by removing all the elements on one side of the given line.

Note that the mesh needs to be pre-processed so that only TRI3 are crossed by the given line

Parameters

mesh	input mesh that need to be trimmed
cut_line_params	parameters of the line that cuts the input mesh
block_id_to_remove	a temporary subdomain id used to mark the elements that need to be removed
new_boundary_id	boundary id of the new boundary that forms due to the trimming

Definition at line 755 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElem().

{
  // Find all the elements that are across the cutting line
  std::vector<dof_id_type> cross_elems;
  // A vector for element specific information
  std::vector<std::vector<std::pair<dof_id_type, dof_id_type>>> node_pairs_vec;
  // A set for unique pairs
  std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs_unique_vec;
  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
       elem_it++)
  {
    std::vector<unsigned short> node_side_rec;
    const auto n_vertices = (*elem_it)->n_vertices();
    node_side_rec.resize(n_vertices);
    for (unsigned int i = 0; i < n_vertices; i++)
    {
      // First check if the vertex is in the XY Plane
      if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
        mooseError("MooseMeshXYCuttingUtils::lineRemoverCutElemTri() only works for 2D meshes in "
                   "XY Plane.");
      const Point v_point = (*elem_it)->point(i);
      node_side_rec[i] = lineSideDeterminator(
          v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2], true);
    }
    if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) == (int)node_side_rec.size())
    {
      (*elem_it)->subdomain_id() = block_id_to_remove;
    }
    else if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) > 0)
    {
      if ((*elem_it)->n_vertices() != 3 || (*elem_it)->n_nodes() != 3)
        mooseError("The element across the cutting line is not TRI3, which is not supported.");
      cross_elems.push_back((*elem_it)->id());
      // Then we need to check pairs of nodes that are on the different side
      std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs;
      for (unsigned int i = 0; i < node_side_rec.size(); i++)
      {
        // first node on removal side and second node on retaining side
        if (node_side_rec[i] > 0 && node_side_rec[(i + 1) % node_side_rec.size()] == 0)
        {
          // Removal side first
          node_pairs.push_back(
              std::make_pair((*elem_it)->node_ptr(i)->id(),
                             (*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id()));
          node_pairs_unique_vec.push_back(node_pairs.back());
        }
        // first node on retaining side and second node on removal side
        else if (node_side_rec[i] == 0 && node_side_rec[(i + 1) % node_side_rec.size()] > 0)
        {
          // Removal side first
          node_pairs.push_back(
              std::make_pair((*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id(),
                             (*elem_it)->node_ptr(i)->id()));
          node_pairs_unique_vec.push_back(node_pairs.back());
        }
      }
      node_pairs_vec.push_back(node_pairs);
    }
  }
  auto vec_ip = std::unique(node_pairs_unique_vec.begin(), node_pairs_unique_vec.end());
  node_pairs_unique_vec.resize(std::distance(node_pairs_unique_vec.begin(), vec_ip));

  // Loop over all the node pairs to define new nodes that sit on the cutting line
  std::vector<Node *> nodes_on_line;
  // whether the on-line node is overlapped with the node pairs or a brand new node
  std::vector<unsigned short> nodes_on_line_overlap;
  for (const auto & node_pair : node_pairs_unique_vec)
  {
    const Point pt1 = *mesh.node_ptr(node_pair.first);
    const Point pt2 = *mesh.node_ptr(node_pair.second);
    const Point pt_line = twoPointandLineIntersection(
        pt1, pt2, cut_line_params[0], cut_line_params[1], cut_line_params[2]);
    if ((pt_line - pt1).norm() < libMesh::TOLERANCE)
    {
      nodes_on_line.push_back(mesh.node_ptr(node_pair.first));
      nodes_on_line_overlap.push_back(1);
    }
    else if ((pt_line - pt2).norm() < libMesh::TOLERANCE)
    {
      nodes_on_line.push_back(mesh.node_ptr(node_pair.second));
      nodes_on_line_overlap.push_back(2);
    }
    else
    {
      nodes_on_line.push_back(mesh.add_point(pt_line));
      nodes_on_line_overlap.push_back(0);
    }
  }

  // make new elements
  for (unsigned int i = 0; i < cross_elems.size(); i++)
  {
    // Only TRI elements are involved after preprocessing
    auto cross_elem = mesh.elem_ptr(cross_elems[i]);
    auto node_0 = cross_elem->node_ptr(0);
    auto node_1 = cross_elem->node_ptr(1);
    auto node_2 = cross_elem->node_ptr(2);
    const std::vector<dof_id_type> tri_nodes = {node_0->id(), node_1->id(), node_2->id()};

    const auto online_node_index_1 = std::distance(node_pairs_unique_vec.begin(),
                                                   std::find(node_pairs_unique_vec.begin(),
                                                             node_pairs_unique_vec.end(),
                                                             node_pairs_vec[i][0]));
    const auto online_node_index_2 = std::distance(node_pairs_unique_vec.begin(),
                                                   std::find(node_pairs_unique_vec.begin(),
                                                             node_pairs_unique_vec.end(),
                                                             node_pairs_vec[i][1]));
    auto node_3 = nodes_on_line[online_node_index_1];
    auto node_4 = nodes_on_line[online_node_index_2];
    const auto node_3_overlap_flag = nodes_on_line_overlap[online_node_index_1];
    const auto node_4_overlap_flag = nodes_on_line_overlap[online_node_index_2];
    // Most common case, no overlapped nodes
    if (node_3_overlap_flag == 0 && node_4_overlap_flag == 0)
    {
      // True if the common node is on the removal side; false if on the retaining side
      const bool common_node_side = node_pairs_vec[i][0].first == node_pairs_vec[i][1].first;
      const subdomain_id_type block_id_to_assign_1 =
          common_node_side ? block_id_to_remove : cross_elem->subdomain_id();
      const subdomain_id_type block_id_to_assign_2 =
          common_node_side ? cross_elem->subdomain_id() : block_id_to_remove;
      // The reference node ids need to be adjusted according to the common node of the two cut
      // sides
      const dof_id_type common_node_id =
          common_node_side ? node_pairs_vec[i][0].first : node_pairs_vec[i][0].second;

      triElemSplitter(mesh,
                      cross_elem->id(),
                      std::distance(tri_nodes.begin(),
                                    std::find(tri_nodes.begin(), tri_nodes.end(), common_node_id)),
                      node_3->id(),
                      node_4->id(),
                      block_id_to_assign_1,
                      block_id_to_assign_2);
      mesh.delete_elem(cross_elem);
    }
    // both node_3 and node_4 are overlapped
    else if (node_3_overlap_flag > 0 && node_4_overlap_flag > 0)
    {
      // In this case, the entire element is on one side of the cutting line
      // No change needed just check which side the element is on
      cross_elem->subdomain_id() = lineSideDeterminator(cross_elem->vertex_average()(0),
                                                        cross_elem->vertex_average()(1),
                                                        cut_line_params[0],
                                                        cut_line_params[1],
                                                        cut_line_params[2],
                                                        true)
                                       ? block_id_to_remove
                                       : cross_elem->subdomain_id();
    }
    // node_3 or node_4 is overlapped
    else
    {
      const auto node_3_finder = std::distance(
          tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_3->id()));
      const auto node_4_finder = std::distance(
          tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_4->id()));
      // As only one of the two above values should be less than the three, the smaller one should
      // be used
      const dof_id_type node_id = node_3_finder < node_4_finder ? node_4->id() : node_3->id();
      const auto node_finder = std::min(node_3_finder, node_4_finder);

      triElemSplitter(
          mesh,
          cross_elem->id(),
          node_finder,
          node_id,
          tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
              ? block_id_to_remove
              : cross_elem->subdomain_id(),
          tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
              ? cross_elem->subdomain_id()
              : block_id_to_remove);
      mesh.delete_elem(cross_elem);
    }
  }
  mesh.contract();
  // Due to the complexity, we identify the new boundary here together instead of during cutting of
  // each element, because the preexisting element edges that are aligned with the cutting line also
  // need to be added to the new boundary.
  mesh.find_neighbors();
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
       elem_it++)
  {
    if ((*elem_it)->subdomain_id() != block_id_to_remove)
    {
      for (unsigned int j = 0; j < (*elem_it)->n_sides(); j++)
      {
        if ((*elem_it)->neighbor_ptr(j) != nullptr)
          if ((*elem_it)->neighbor_ptr(j)->subdomain_id() == block_id_to_remove)
            boundary_info.add_side(*elem_it, j, new_boundary_id);
      }
    }
  }

  // Delete the block to remove
  for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
       elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
       elem_it++)
    mesh.delete_elem(*elem_it);
  mesh.contract();
}

lineRemoverMoveNode()

void MooseMeshXYCuttingUtils::lineRemoverMoveNode	(	libMesh::ReplicatedMesh &	mesh,
		const std::vector< Real > &	bdry_pars,
		const subdomain_id_type	block_id_to_remove,
		const std::set< subdomain_id_type > &	subdomain_ids_set,
		const boundary_id_type	trimming_section_boundary_id,
		const boundary_id_type	external_boundary_id,
		const std::vector< boundary_id_type > &	other_boundaries_to_conform = std::vector<boundary_id_type>(),
		const bool	assign_ext_to_new = false,
		const bool	side_to_remove = true
	)

Removes all the elements on one side of a given line and deforms the elements intercepted by the line to form a flat new boundary.

Parameters

mesh	input mesh to perform line-based elements removing on
bdry_pars	line parameter sets {a, b, c} as in a*x+b*y+c=0
block_id_to_remove	subdomain id used to mark the elements that need to be removed
subdomain_ids_set	all the subdomain ids in the input mesh
trimming_section_boundary_id	ID of the new external boundary formed due to trimming
external_boundary_id	ID of the external boundary of the input mesh
other_boundaries_to_conform	IDs of the other boundaries that need to be conformed to during nodes moving
assign_ext_to_new	whether to assign external_boundary_id to the new boundary formed by removal
side_to_remove	which side of the mesh needs to be removed: true means ax+by+c>0 and false means ax+by+c<0

Definition at line 27 of file MooseMeshXYCuttingUtils.C.

Referenced by XYMeshLineCutter::generate().

{
  // Build boundary information of the mesh
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto bdry_side_list = boundary_info.build_side_list();
  // Only select the boundaries_to_conform
  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> slc_bdry_side_list;
  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
    if (std::get<2>(bdry_side_list[i]) == external_boundary_id ||
        std::find(other_boundaries_to_conform.begin(),
                  other_boundaries_to_conform.end(),
                  std::get<2>(bdry_side_list[i])) != other_boundaries_to_conform.end())
      slc_bdry_side_list.push_back(bdry_side_list[i]);

  // Assign block id for elements to be removed
  // Also record the elements crossed by the line and with its average vertices on the removal side
  std::vector<dof_id_type> crossed_elems_to_remove;
  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
       elem_it++)
  {
    // Check all the vertices of the element
    unsigned short removal_side_count = 0;
    for (unsigned int i = 0; i < (*elem_it)->n_vertices(); i++)
    {
      // First check if the vertex is on the XY-Plane
      if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
        mooseError(
            "MooseMeshXYCuttingUtils::lineRemoverMoveNode() only works for 2D meshes in XY plane.");
      if (lineSideDeterminator((*elem_it)->point(i)(0),
                               (*elem_it)->point(i)(1),
                               bdry_pars[0],
                               bdry_pars[1],
                               bdry_pars[2],
                               side_to_remove))
        removal_side_count++;
    }
    if (removal_side_count == (*elem_it)->n_vertices())
    {
      (*elem_it)->subdomain_id() = block_id_to_remove;
      continue;
    }
    // Check the average of the vertices of the element
    if (lineSideDeterminator((*elem_it)->vertex_average()(0),
                             (*elem_it)->vertex_average()(1),
                             bdry_pars[0],
                             bdry_pars[1],
                             bdry_pars[2],
                             side_to_remove))
      crossed_elems_to_remove.push_back((*elem_it)->id());
  }
  // Check each crossed element to see if removing it would lead to boundary moving
  for (const auto & elem_id : crossed_elems_to_remove)
  {
    bool remove_flag = true;
    for (unsigned int i = 0; i < mesh.elem_ptr(elem_id)->n_sides(); i++)
    {
      if (mesh.elem_ptr(elem_id)->neighbor_ptr(i) != nullptr)
        if (mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id() != block_id_to_remove &&
            std::find(crossed_elems_to_remove.begin(),
                      crossed_elems_to_remove.end(),
                      mesh.elem_ptr(elem_id)->neighbor_ptr(i)->id()) ==
                crossed_elems_to_remove.end())
        {
          if (mesh.elem_ptr(elem_id)->subdomain_id() !=
              mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id())
          {
            remove_flag = false;
            break;
          }
        }
    }
    if (remove_flag)
      mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
  }

  // Identify all the nodes that are on the interface between block_id_to_remove and other blocks
  // !!! We need a check here: if a node is on the retaining side, the removed element has a
  // different subdomain id
  std::vector<dof_id_type> node_list;
  for (auto elem_it = mesh.active_subdomain_set_elements_begin(subdomain_ids_set);
       elem_it != mesh.active_subdomain_set_elements_end(subdomain_ids_set);
       elem_it++)
  {
    for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
    {
      if ((*elem_it)->neighbor_ptr(i) != nullptr)
        if ((*elem_it)->neighbor_ptr(i)->subdomain_id() == block_id_to_remove)
        {
          node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(0)->id());
          node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(1)->id());
          boundary_info.add_side(*elem_it, i, trimming_section_boundary_id);
          if (assign_ext_to_new && trimming_section_boundary_id != external_boundary_id)
            boundary_info.add_side(*elem_it, i, external_boundary_id);
        }
    }
  }
  // Remove duplicate nodes
  const auto unique_it = std::unique(node_list.begin(), node_list.end());
  node_list.resize(std::distance(node_list.begin(), unique_it));
  // Mark those nodes that are on a boundary that requires conformality
  // If both nodes of a side are involved, we should only move one node
  std::vector<bool> node_list_flag(node_list.size(), false);
  std::vector<Point> node_list_point(node_list.size(), Point(0.0, 0.0, 0.0));
  // Loop over all the selected sides
  for (unsigned int i = 0; i < slc_bdry_side_list.size(); i++)
  {
    // Get the two node ids of the side
    dof_id_type side_id_0 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
                                ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
                                ->node_ptr(0)
                                ->id();
    dof_id_type side_id_1 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
                                ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
                                ->node_ptr(1)
                                ->id();
    // True means the selected bdry node is in the node list of the trimming interface
    bool side_id_0_in =
        !(std::find(node_list.begin(), node_list.end(), side_id_0) == node_list.end());
    bool side_id_1_in =
        !(std::find(node_list.begin(), node_list.end(), side_id_1) == node_list.end());

    // True means the selected bdry node is on the removal side of the trimming interface
    bool side_node_0_remove = lineSideDeterminator((*mesh.node_ptr(side_id_0))(0),
                                                   (*mesh.node_ptr(side_id_0))(1),
                                                   bdry_pars[0],
                                                   bdry_pars[1],
                                                   bdry_pars[2],
                                                   side_to_remove);
    bool side_node_1_remove = lineSideDeterminator((*mesh.node_ptr(side_id_1))(0),
                                                   (*mesh.node_ptr(side_id_1))(1),
                                                   bdry_pars[0],
                                                   bdry_pars[1],
                                                   bdry_pars[2],
                                                   side_to_remove);
    // If both nodes of that side are involved in the trimming interface
    if (side_id_0_in && side_id_1_in)
      // The side needs to be removed from the sideset because it is not longer an interface
      // The other node will be handled by other element's side
      boundary_info.remove_side(mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i])),
                                std::get<1>(slc_bdry_side_list[i]),
                                std::get<2>(slc_bdry_side_list[i]));
    // If node 0 is on the trimming interface, and the side is cut by the trimming line
    else if (side_id_0_in && (side_node_0_remove != side_node_1_remove))
    {
      // Use the intersection point as the destination of the node after moving
      node_list_flag[std::distance(
          node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_0))] = true;
      const Point p0 = *mesh.node_ptr(side_id_0);
      const Point p1 = *mesh.node_ptr(side_id_1);

      node_list_point[std::distance(node_list.begin(),
                                    std::find(node_list.begin(), node_list.end(), side_id_0))] =
          twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
    }
    // If node 1 is on the trimming interface, and the side is cut by the trimming line
    else if (side_id_1_in && (side_node_0_remove != side_node_1_remove))
    {
      // Use the intersection point as the destination of the node after moving
      node_list_flag[std::distance(
          node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_1))] = true;
      const Point p0 = *mesh.node_ptr(side_id_0);
      const Point p1 = *mesh.node_ptr(side_id_1);

      node_list_point[std::distance(node_list.begin(),
                                    std::find(node_list.begin(), node_list.end(), side_id_1))] =
          twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
    }
  }

  // move nodes
  for (unsigned int i = 0; i < node_list.size(); i++)
  {
    // This means the node is on both the trimming boundary and the original external
    // boundary/selected interface boundaries. In order to keep the shape of the original external
    // boundary, the node is moved along the original external boundary.
    if (node_list_flag[i])
      *(mesh.node_ptr(node_list[i])) = node_list_point[i];
    // This means the node does not need to conform to any boundaries.
    // Just move it along the normal direction of the trimming line.
    else
    {
      const Real x0 = (*(mesh.node_ptr(node_list[i])))(0);
      const Real y0 = (*(mesh.node_ptr(node_list[i])))(1);
      (*(mesh.node_ptr(node_list[i])))(0) =
          (bdry_pars[1] * (bdry_pars[1] * x0 - bdry_pars[0] * y0) - bdry_pars[0] * bdry_pars[2]) /
          (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
      (*(mesh.node_ptr(node_list[i])))(1) =
          (bdry_pars[0] * (-bdry_pars[1] * x0 + bdry_pars[0] * y0) - bdry_pars[1] * bdry_pars[2]) /
          (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
    }
  }

  // Delete the block
  for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
       elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
       elem_it++)
    mesh.delete_elem(*elem_it);
  mesh.contract();
  mesh.find_neighbors();
  // Delete zero volume elements
  std::vector<dof_id_type> zero_elems;
  for (auto elem_it = mesh.elements_begin(); elem_it != mesh.elements_end(); elem_it++)
  {
    if (MooseUtils::absoluteFuzzyEqual((*elem_it)->volume(), 0.0))
    {
      for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
      {
        if ((*elem_it)->neighbor_ptr(i) != nullptr)
        {
          boundary_info.add_side((*elem_it)->neighbor_ptr(i),
                                 ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
                                 external_boundary_id);
          boundary_info.add_side((*elem_it)->neighbor_ptr(i),
                                 ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
                                 trimming_section_boundary_id);
        }
      }
      zero_elems.push_back((*elem_it)->id());
    }
  }
  for (const auto & zero_elem : zero_elems)
    mesh.delete_elem(mesh.elem_ptr(zero_elem));
  mesh.contract();
  // As we modified the side_list, it is safer to clear the node_list
  boundary_info.clear_boundary_node_ids();
  mesh.prepare_for_use();
}

lineSideDeterminator()

bool MooseMeshXYCuttingUtils::lineSideDeterminator	(	const Real	px,
		const Real	py,
		const Real	param_1,
		const Real	param_2,
		const Real	param_3,
		const bool	direction_param,
		const Real	dis_tol = libMesh::TOLERANCE
	)

Determines whether a point on XY-plane is on the side of a given line that needs to be removed.

Parameters

px	x coordinate of the point
py	y coordinate of the point
param_1	parameter 1 (a) in line formula a*x+b*y+c=0
param_2	parameter 2 (b) in line formula a*x+b*y+c=0
param_3	parameter 3 (c) in line formula a*x+b*y+c=0
direction_param	which side is the side that needs to be removed
dis_tol	tolerance used in determining side

Returns

whether the point is on the side of the line that needed to be removed

Definition at line 264 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElemTri(), lineRemoverMoveNode(), and quadToTriOnLine().

{
  const Real tmp = px * param_1 + py * param_2 + param_3;
  return direction_param ? tmp >= dis_tol : tmp <= dis_tol;
}

makeImprovedTriElement()

void MooseMeshXYCuttingUtils::makeImprovedTriElement	(	libMesh::ReplicatedMesh &	mesh,
		const dof_id_type	node_id_0,
		const dof_id_type	node_id_1,
		const dof_id_type	node_id_2,
		const subdomain_id_type	subdomain_id,
		const std::vector< dof_id_type > &	extra_elem_ids,
		const std::vector< boundary_id_type > &	boundary_ids_for_side_1 = std::vector<boundary_id_type>(),
		const std::vector< boundary_id_type > &	boundary_ids_for_side_0 = std::vector<boundary_id_type>(),
		const std::vector< boundary_id_type > &	boundary_ids_for_side_2 = std::vector<boundary_id_type>()
	)

Make a TRI3 element with the given node ids and subdomain id with boundary information.

Parameters

mesh	input mesh where the TRI3 element needs to be added
node_id_0	id of the first node of the TRI3 element
node_id_1	id of the second node of the TRI3 element
node_id_2	id of the third node of the TRI3 element
subdomain_id	subdomain id of the TRI3 element
extra_elem_ids	extra element ids to be assigned to the TRI3 element
boundary_ids_for_side_1	boundary ids of the second side of the TRI3 element
boundary_ids_for_side_0	boundary ids of the first side of the TRI3 element
boundary_ids_for_side_2	boundary ids of the third side of the TRI3 element

Definition at line 1185 of file MooseMeshXYCuttingUtils.C.

Referenced by boundaryTriElemImprover().

{
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  Elem * elem_Tri3_new = mesh.add_elem(new Tri3);
  elem_Tri3_new->set_node(0, mesh.node_ptr(node_id_0));
  elem_Tri3_new->set_node(1, mesh.node_ptr(node_id_1));
  elem_Tri3_new->set_node(2, mesh.node_ptr(node_id_2));
  for (const auto & boundary_id_for_side_0 : boundary_ids_for_side_0)
    boundary_info.add_side(elem_Tri3_new, 0, boundary_id_for_side_0);
  for (const auto & boundary_id_for_side_1 : boundary_ids_for_side_1)
    boundary_info.add_side(elem_Tri3_new, 1, boundary_id_for_side_1);
  for (const auto & boundary_id_for_side_2 : boundary_ids_for_side_2)
    boundary_info.add_side(elem_Tri3_new, 2, boundary_id_for_side_2);
  elem_Tri3_new->subdomain_id() = subdomain_id;
  // Retain element extra integers
  for (unsigned int j = 0; j < extra_elem_ids.size(); j++)
  {
    elem_Tri3_new->set_extra_integer(j, extra_elem_ids[j]);
  }
}

quadElemSplitter()

void MooseMeshXYCuttingUtils::quadElemSplitter	(	libMesh::ReplicatedMesh &	mesh,
		const dof_id_type	elem_id,
		const subdomain_id_type	tri_elem_subdomain_shift
	)

Split a QUAD4 element into two TRI3 elements.

Parameters

mesh	input mesh with the QUAD4 element that needs to be split
elem_id	id of the QUAD4 element that needs to be split
tri_elem_subdomain_shift	subdomain id shift used to define the TRI element subdomains

Definition at line 601 of file MooseMeshXYCuttingUtils.C.

Referenced by quadToTriOnLine().

{
  // Build boundary information of the mesh
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto bdry_side_list = boundary_info.build_side_list();
  // Create a list of sidesets involving the element to be split
  std::vector<std::vector<boundary_id_type>> elem_side_list;
  elem_side_list.resize(4);
  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
  {
    if (std::get<0>(bdry_side_list[i]) == elem_id)
    {
      elem_side_list[std::get<1>(bdry_side_list[i])].push_back(std::get<2>(bdry_side_list[i]));
    }
  }

  auto node_0 = mesh.elem_ptr(elem_id)->node_ptr(0);
  auto node_1 = mesh.elem_ptr(elem_id)->node_ptr(1);
  auto node_2 = mesh.elem_ptr(elem_id)->node_ptr(2);
  auto node_3 = mesh.elem_ptr(elem_id)->node_ptr(3);

  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
  // Record all the element extra integers of the original quad element
  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
    exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);

  // There are two trivial ways to split a quad element
  // We prefer the way that leads to triangles with similar areas
  if (std::abs((*node_1 - *node_0).cross(*node_3 - *node_0).norm() -
               (*node_1 - *node_2).cross(*node_3 - *node_2).norm()) >
      std::abs((*node_0 - *node_1).cross(*node_2 - *node_1).norm() -
               (*node_0 - *node_3).cross(*node_2 - *node_3).norm()))
  {
    Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
    elem_Tri3_0->set_node(0, node_0);
    elem_Tri3_0->set_node(1, node_1);
    elem_Tri3_0->set_node(2, node_2);
    elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
    Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
    elem_Tri3_1->set_node(0, node_0);
    elem_Tri3_1->set_node(1, node_2);
    elem_Tri3_1->set_node(2, node_3);
    elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
    // Retain element extra integers
    for (unsigned int j = 0; j < n_elem_extra_ids; j++)
    {
      elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
      elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
    }

    // Add sideset information to the new elements
    for (const auto & side_info_0 : elem_side_list[0])
      boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
    for (const auto & side_info_1 : elem_side_list[1])
      boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
    for (const auto & side_info_2 : elem_side_list[2])
      boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
    for (const auto & side_info_3 : elem_side_list[3])
      boundary_info.add_side(elem_Tri3_1, 2, side_info_3);
  }
  else
  {
    Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
    elem_Tri3_0->set_node(0, node_0);
    elem_Tri3_0->set_node(1, node_1);
    elem_Tri3_0->set_node(2, node_3);
    elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
    Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
    elem_Tri3_1->set_node(0, node_1);
    elem_Tri3_1->set_node(1, node_2);
    elem_Tri3_1->set_node(2, node_3);
    elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
    // Retain element extra integers
    for (unsigned int j = 0; j < n_elem_extra_ids; j++)
    {
      elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
      elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
    }

    // Add sideset information to the new elements
    for (const auto & side_info_0 : elem_side_list[0])
      boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
    for (const auto & side_info_1 : elem_side_list[1])
      boundary_info.add_side(elem_Tri3_1, 0, side_info_1);
    for (const auto & side_info_2 : elem_side_list[2])
      boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
    for (const auto & side_info_3 : elem_side_list[3])
      boundary_info.add_side(elem_Tri3_0, 2, side_info_3);
  }
}

quadToTriOnLine()

void MooseMeshXYCuttingUtils::quadToTriOnLine	(	libMesh::ReplicatedMesh &	mesh,
		const std::vector< Real > &	cut_line_params,
		const dof_id_type	tri_subdomain_id_shift,
		const SubdomainName	tri_elem_subdomain_name_suffix
	)

Convert all the QUAD4 elements in the mesh that are crossed by the given line into TRI3 elements.

Parameters

mesh	input mesh with the QUAD4 elements that need to be converted
cut_line_params	parameters of the line that cuts the input mesh
tri_subdomain_id_shift	subdomain id shift used to define the TRI element subdomains generated due to the conversion
tri_elem_subdomain_name_suffix	suffix used to name the TRI element subdomains generated due to the conversion

Definition at line 696 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElem().

{
  // Preprocess: find all the quad elements that are across the cutting line
  std::vector<dof_id_type> cross_elems_quad;
  std::set<subdomain_id_type> new_subdomain_ids;
  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
       elem_it++)
  {
    if ((*elem_it)->n_vertices() == 4)
    {
      std::vector<unsigned short> node_side_rec;
      for (unsigned int i = 0; i < 4; i++)
      {
        const Point v_point = (*elem_it)->point(i);
        node_side_rec.push_back(lineSideDeterminator(v_point(0),
                                                     v_point(1),
                                                     cut_line_params[0],
                                                     cut_line_params[1],
                                                     cut_line_params[2],
                                                     true));
      }
      if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) != 4 &&
          std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) > 0)
      {
        cross_elems_quad.push_back((*elem_it)->id());
        new_subdomain_ids.emplace((*elem_it)->subdomain_id() + tri_subdomain_id_shift);
      }
    }
  }
  // Then convert these quad elements into tri elements
  for (const auto & cross_elem_quad : cross_elems_quad)
  {
    quadElemSplitter(mesh, cross_elem_quad, tri_subdomain_id_shift);
    mesh.delete_elem(mesh.elem_ptr(cross_elem_quad));
  }
  for (auto & nid : new_subdomain_ids)
  {
    const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
    if (MooseMeshUtils::getSubdomainID(
            (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
                              : old_name) +
                "_" + tri_elem_subdomain_name_suffix,
            mesh) != Moose::INVALID_BLOCK_ID)
      throw MooseException("The new subdomain name already exists in the mesh.");
    mesh.subdomain_name(nid) =
        (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
                          : old_name) +
        "_" + tri_elem_subdomain_name_suffix;
    mooseWarning("QUAD elements have been converted into TRI elements with a new "
                 "subdomain name: " +
                 mesh.subdomain_name(nid) + ".");
  }
  mesh.contract();
}

quasiTriElementsFixer()

bool MooseMeshXYCuttingUtils::quasiTriElementsFixer	(	libMesh::ReplicatedMesh &	mesh,
		const std::set< subdomain_id_type > &	subdomain_ids_set,
		const subdomain_id_type	tri_elem_subdomain_shift = Moose::INVALID_BLOCK_ID,
		const SubdomainName	tri_elem_subdomain_name_suffix = "tri"
	)

Fixes degenerate QUAD elements created by the hexagonal mesh trimming by converting them into TRI elements.

Parameters

mesh	input mesh with degenerate QUAD elements that need to be fixed
subdomain_ids_set	all the subdomain ids in the input mesh
tri_elem_subdomain_shift	subdomain id shift used to define the TRI element subdomains
tri_elem_subdomain_name_suffix	suffix used to name the TRI element subdomains

Returns

whether any elements have been fixed

Definition at line 306 of file MooseMeshXYCuttingUtils.C.

Referenced by XYMeshLineCutter::generate().

{
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  // Define the subdomain id shift for the new TRI3 element subdomain(s)
  const subdomain_id_type max_subdomain_id(*subdomain_ids_set.rbegin());
  const subdomain_id_type tri_subdomain_id_shift =
      tri_elem_subdomain_shift == Moose::INVALID_BLOCK_ID ? max_subdomain_id
                                                          : tri_elem_subdomain_shift;
  mooseAssert(std::numeric_limits<subdomain_id_type>::max() - max_subdomain_id >
                  tri_subdomain_id_shift,
              "The TRI elements subdomain id to be assigned may exceed the numeric limit.");
  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
  std::vector<std::tuple<Elem *, unsigned int, bool, bool>> bad_elems_rec;
  // Loop over all the active elements to find any degenerate QUAD elements
  for (auto & elem : as_range(mesh.active_elements_begin(), mesh.active_elements_end()))
  {
    // Two types of degenerate QUAD elements are identified:
    // (1) QUAD elements with three collinear vertices
    // (2) QUAD elements with two overlapped vertices
    const auto elem_angles = vertex_angles(*elem);
    const auto elem_distances = vertex_distances(*elem);
    // Type 1
    if (MooseUtils::absoluteFuzzyEqual(elem_angles.front().first, M_PI, 0.001))
    {
      bad_elems_rec.push_back(std::make_tuple(elem, elem_angles.front().second, false, true));
      continue;
    }
    // Type 2
    if (MooseUtils::absoluteFuzzyEqual(elem_distances.front().first, 0.0))
    {
      bad_elems_rec.push_back(std::make_tuple(elem, elem_distances.front().second, false, false));
    }
  }
  std::set<subdomain_id_type> new_subdomain_ids;
  // Loop over all the identified degenerate QUAD elements
  for (const auto & bad_elem : bad_elems_rec)
  {
    std::vector<boundary_id_type> elem_bdry_container_0;
    std::vector<boundary_id_type> elem_bdry_container_1;
    std::vector<boundary_id_type> elem_bdry_container_2;

    Elem * elem_0 = std::get<0>(bad_elem);
    if (std::get<3>(bad_elem))
    {
      // elems 1 and 2 are the neighboring elements of the degenerate element corresponding to the
      // two collinear sides.
      // For the degenerated element with three colinear vertices, if the elems 1 and 2 do not
      // exist, the two sides are on the external boundary formed by trimming.
      Elem * elem_1 = elem_0->neighbor_ptr(std::get<1>(bad_elem));
      Elem * elem_2 = elem_0->neighbor_ptr((std::get<1>(bad_elem) - 1) % elem_0->n_vertices());
      if ((elem_1 != nullptr || elem_2 != nullptr))
        throw MooseException("The input mesh has degenerate quad element before trimming.");
    }
    mesh.get_boundary_info().boundary_ids(elem_0, std::get<1>(bad_elem), elem_bdry_container_0);
    mesh.get_boundary_info().boundary_ids(
        elem_0, (std::get<1>(bad_elem) + 1) % elem_0->n_vertices(), elem_bdry_container_1);
    mesh.get_boundary_info().boundary_ids(
        elem_0, (std::get<1>(bad_elem) + 2) % elem_0->n_vertices(), elem_bdry_container_2);
    mesh.get_boundary_info().boundary_ids(
        elem_0, (std::get<1>(bad_elem) + 3) % elem_0->n_vertices(), elem_bdry_container_0);

    // Record subdomain id of the degenerate element
    auto elem_block_id = elem_0->subdomain_id();
    // Define the three of four nodes that will be used to generate the TRI element
    auto pt0 = elem_0->node_ptr((std::get<1>(bad_elem) + 1) % elem_0->n_vertices());
    auto pt1 = elem_0->node_ptr((std::get<1>(bad_elem) + 2) % elem_0->n_vertices());
    auto pt2 = elem_0->node_ptr((std::get<1>(bad_elem) + 3) % elem_0->n_vertices());
    // Record all the element extra integers of the degenerate element
    for (unsigned int j = 0; j < n_elem_extra_ids; j++)
      exist_extra_ids[j] = elem_0->get_extra_integer(j);
    // Delete the degenerate QUAD element
    mesh.delete_elem(elem_0);
    // Create the new TRI element
    Elem * elem_Tri3 = mesh.add_elem(new Tri3);
    elem_Tri3->set_node(0, pt0);
    elem_Tri3->set_node(1, pt1);
    elem_Tri3->set_node(2, pt2);
    // Retain the boundary information
    for (auto bdry_id : elem_bdry_container_0)
      boundary_info.add_side(elem_Tri3, 2, bdry_id);
    for (auto bdry_id : elem_bdry_container_1)
      boundary_info.add_side(elem_Tri3, 0, bdry_id);
    for (auto bdry_id : elem_bdry_container_2)
      boundary_info.add_side(elem_Tri3, 1, bdry_id);
    // Assign subdomain id for the TRI element by shifting its original subdomain id
    elem_Tri3->subdomain_id() = elem_block_id + tri_subdomain_id_shift;
    new_subdomain_ids.emplace(elem_block_id + tri_subdomain_id_shift);
    // Retain element extra integers
    for (unsigned int j = 0; j < n_elem_extra_ids; j++)
      elem_Tri3->set_extra_integer(j, exist_extra_ids[j]);
  }
  // Assign subdomain names for the new TRI elements
  for (auto & nid : new_subdomain_ids)
  {
    const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
    if (MooseMeshUtils::getSubdomainID(
            (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
                              : old_name) +
                "_" + tri_elem_subdomain_name_suffix,
            mesh) != Moose::INVALID_BLOCK_ID)
      throw MooseException("The new subdomain name already exists in the mesh.");
    mesh.subdomain_name(nid) =
        (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
                          : old_name) +
        "_" + tri_elem_subdomain_name_suffix;
    mooseWarning("Degenerate QUAD elements have been converted into TRI elements with a new "
                 "subdomain name: " +
                 mesh.subdomain_name(nid) + ".");
  }
  return bad_elems_rec.size();
}

triElemSplitter() [1/2]

void MooseMeshXYCuttingUtils::triElemSplitter	(	libMesh::ReplicatedMesh &	mesh,
		const dof_id_type	elem_id,
		const unsigned short	node_shift,
		const dof_id_type	nid_3,
		const dof_id_type	nid_4,
		const subdomain_id_type	single_elem_side_id,
		const subdomain_id_type	double_elem_side_id
	)

Split a TRI3 element into three TRI3 elements based on two nodes on the two sides of the triangle.

Parameters

mesh	input mesh with the TRI3 element that needs to be split
elem_id	id of the TRI3 element that needs to be split
node_shift	shift used to rotate the vertices to make sure the vertex that the two cut sides share is the first vertex
nid_3	id of the node on the first cut side of the triangle
nid_4	id of the node on the second cut side of the triangle
single_elem_side_id	subdomain id of the single element side
double_elem_side_id	subdomain id of the double element side

Definition at line 459 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElemTri().

{
  const auto elem_old = mesh.elem_ptr(elem_id);
  const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
  const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
  const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();

  const bool m1_side_flag =
      MooseUtils::absoluteFuzzyEqual((*(mesh.node_ptr(nid_3)) - *(mesh.node_ptr(nid_0)))
                                         .cross(*(mesh.node_ptr(nid_1)) - *(mesh.node_ptr(nid_0)))
                                         .norm(),
                                     0.0);
  const dof_id_type nid_m1 = m1_side_flag ? nid_3 : nid_4;
  const dof_id_type nid_m2 = m1_side_flag ? nid_4 : nid_3;
  // Build boundary information of the mesh
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto bdry_side_list = boundary_info.build_side_list();
  // Create a list of sidesets involving the element to be split
  std::vector<std::vector<boundary_id_type>> elem_side_list;
  elem_side_list.resize(3);
  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
  {
    if (std::get<0>(bdry_side_list[i]) == elem_id)
    {
      elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
          std::get<2>(bdry_side_list[i]));
    }
  }

  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
  // Record all the element extra integers of the original element
  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
    exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);

  Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
  elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
  elem_Tri3_0->set_node(1, mesh.node_ptr(nid_m1));
  elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m2));
  elem_Tri3_0->subdomain_id() = single_elem_side_id;
  Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
  elem_Tri3_1->set_node(0, mesh.node_ptr(nid_1));
  elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m2));
  elem_Tri3_1->set_node(2, mesh.node_ptr(nid_m1));
  elem_Tri3_1->subdomain_id() = double_elem_side_id;
  Elem * elem_Tri3_2 = mesh.add_elem(new Tri3);
  elem_Tri3_2->set_node(0, mesh.node_ptr(nid_2));
  elem_Tri3_2->set_node(1, mesh.node_ptr(nid_m2));
  elem_Tri3_2->set_node(2, mesh.node_ptr(nid_1));
  elem_Tri3_2->subdomain_id() = double_elem_side_id;
  // Retain element extra integers
  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
  {
    elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
    elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
    elem_Tri3_2->set_extra_integer(j, exist_extra_ids[j]);
  }

  // Add sideset information to the new elements
  for (const auto & side_info_0 : elem_side_list[0])
  {
    boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
    boundary_info.add_side(elem_Tri3_1, 2, side_info_0);
  }
  for (const auto & side_info_1 : elem_side_list[1])
    boundary_info.add_side(elem_Tri3_2, 2, side_info_1);
  for (const auto & side_info_2 : elem_side_list[2])
  {
    boundary_info.add_side(elem_Tri3_0, 2, side_info_2);
    boundary_info.add_side(elem_Tri3_2, 0, side_info_2);
  }
}

triElemSplitter() [2/2]

void MooseMeshXYCuttingUtils::triElemSplitter	(	libMesh::ReplicatedMesh &	mesh,
		const dof_id_type	elem_id,
		const unsigned short	node_shift,
		const dof_id_type	nid_m,
		const subdomain_id_type	first_elem_side_id,
		const subdomain_id_type	second_elem_side_id
	)

Split a TRI3 element into two TRI3 elements based on one node on one side of the triangle.

Parameters

mesh	input mesh with the TRI3 element that needs to be split
elem_id	id of the TRI3 element that needs to be split
node_shift	shift used to rotate the vertices to make sure the vertex corresponding to the cut side is the first vertex
nid_m	id of the node on the cut side of the triangle
first_elem_side_id	subdomain id of the first element side
second_elem_side_id	subdomain id of the second element side

Definition at line 539 of file MooseMeshXYCuttingUtils.C.

{
  const auto elem_old = mesh.elem_ptr(elem_id);
  const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
  const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
  const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
  // Build boundary information of the mesh
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto bdry_side_list = boundary_info.build_side_list();
  // Create a list of sidesets involving the element to be split
  std::vector<std::vector<boundary_id_type>> elem_side_list;
  elem_side_list.resize(3);
  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
  {
    if (std::get<0>(bdry_side_list[i]) == elem_id)
    {
      elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
          std::get<2>(bdry_side_list[i]));
    }
  }

  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
  // Record all the element extra integers of the original element
  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
    exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);

  Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
  elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
  elem_Tri3_0->set_node(1, mesh.node_ptr(nid_1));
  elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m));
  elem_Tri3_0->subdomain_id() = first_elem_side_id;
  Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
  elem_Tri3_1->set_node(0, mesh.node_ptr(nid_0));
  elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m));
  elem_Tri3_1->set_node(2, mesh.node_ptr(nid_2));
  elem_Tri3_1->subdomain_id() = second_elem_side_id;
  // Retain element extra integers
  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
  {
    elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
    elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
  }

  // Add sideset information to the new elements
  for (const auto & side_info_0 : elem_side_list[0])
    boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
  for (const auto & side_info_1 : elem_side_list[1])
  {
    boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
    boundary_info.add_side(elem_Tri3_1, 1, side_info_1);
  }
  for (const auto & side_info_2 : elem_side_list[2])
    boundary_info.add_side(elem_Tri3_1, 2, side_info_2);
}

twoLineIntersection()

Point MooseMeshXYCuttingUtils::twoLineIntersection	(	const Real	param_11,
		const Real	param_12,
		const Real	param_13,
		const Real	param_21,
		const Real	param_22,
		const Real	param_23
	)

Calculates the intersection Point of two given straight lines.

Parameters

param_11	parameter 1 (a) in line formula a*x+b*y+c=0 for the first line
param_12	parameter 2 (b) in line formula a*x+b*y+c=0 for the first line
param_13	parameter 3 (c) in line formula a*x+b*y+c=0 for the first line
param_21	parameter 1 (a) in line formula a*x+b*y+c=0 for the second line
param_22	parameter 2 (b) in line formula a*x+b*y+c=0 for the second line
param_23	parameter 3 (c) in line formula a*x+b*y+c=0 for the second line

Returns

intersect point of the two lines

Definition at line 277 of file MooseMeshXYCuttingUtils.C.

Referenced by twoPointandLineIntersection().

{
  return Point(
      (param_12 * param_23 - param_22 * param_13) / (param_11 * param_22 - param_21 * param_12),
      (param_13 * param_21 - param_23 * param_11) / (param_11 * param_22 - param_21 * param_12),
      0.0);
}

twoPointandLineIntersection() [1/2]

Point MooseMeshXYCuttingUtils::twoPointandLineIntersection	(	const Point &	pt1,
		const Point &	pt2,
		const Real	param_1,
		const Real	param_2,
		const Real	param_3
	)

Calculates the intersection Point of a straight line defined by two given points and another straight line.

Parameters

pt1	point 1 that defines the first straight line
pt2	point 2 that defines the first straight line
param_1	parameter 1 (a) in line formula a*x+b*y+c=0 for the second straight line
param_2	parameter 2 (b) in line formula a*x+b*y+c=0 for the second straight line
param_3	parameter 3 (c) in line formula a*x+b*y+c=0 for the second straight line

Returns

intersect point of the two lines

twoPointandLineIntersection() [2/2]

Point MooseMeshXYCuttingUtils::twoPointandLineIntersection	(	const Point &	pt1,
		const Point &	pt2,
		const Real	param_1,
		const Real	param_2,
		const Real	param_3
	)

Definition at line 291 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElemTri(), and lineRemoverMoveNode().

{
  return twoLineIntersection(param_1,
                             param_2,
                             param_3,
                             pt2(1) - pt1(1),
                             pt1(0) - pt2(0),
                             pt2(0) * pt1(1) - pt1(0) * pt2(1));
}

vertex_angles() [1/2]

std::vector<std::pair<Real, unsigned int> > MooseMeshXYCuttingUtils::vertex_angles

(

const Elem &

elem

)

Calculates the internal angles of a given 2D element.

Parameters

elem	the element that needs to be investigated

Returns

sizes of all the internal angles, sorted by their size

vertex_angles() [2/2]

std::vector<std::pair<Real, unsigned int> > MooseMeshXYCuttingUtils::vertex_angles

(

const Elem &

elem

)

Definition at line 423 of file MooseMeshXYCuttingUtils.C.

Referenced by quasiTriElementsFixer().

{
  std::vector<std::pair<Real, unsigned int>> angles;
  const unsigned int n_vertices = elem.n_vertices();

  for (unsigned int i = 0; i < n_vertices; i++)
  {
    Point v1 = (*elem.node_ptr((i - 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
    Point v2 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
    Real tmp = v1 * v2 / v1.norm() / v2.norm();
    if (tmp > 1.0)
      tmp = 1.0;
    else if (tmp < -1.0)
      tmp = -1.0;
    angles.push_back(std::make_pair(acos(tmp), i));
  }
  std::sort(angles.begin(), angles.end(), std::greater<>());
  return angles;
}

vertex_distances() [1/2]

std::vector<std::pair<Real, unsigned int> > MooseMeshXYCuttingUtils::vertex_distances

(

const Elem &

elem

)

Calculates the distances between the vertices of a given 2D element.

Parameters

elem	the element that needs to be investigated

Returns

values of all the distances, sorted by their value

vertex_distances() [2/2]

std::vector<std::pair<Real, unsigned int> > MooseMeshXYCuttingUtils::vertex_distances

(

const Elem &

elem

)

Definition at line 444 of file MooseMeshXYCuttingUtils.C.

Referenced by quasiTriElementsFixer().

{
  std::vector<std::pair<Real, unsigned int>> distances;
  const unsigned int n_vertices = elem.n_vertices();

  for (unsigned int i = 0; i < n_vertices; i++)
  {
    Point v1 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
    distances.push_back(std::make_pair(v1.norm(), i));
  }
  std::sort(distances.begin(), distances.end());
  return distances;
}