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.
bool	pointOnLine (const Real px, const Real py, const Real param_1, const Real param_2, const Real param_3, const Real dis_tol=libMesh::TOLERANCE)
	Determines whether a point on XY-plane is on a given line, to within a tolerance.
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.
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.
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.
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.
std::vector< std::pair< Real, unsigned int > >	vertex_angles (const Elem &elem)
	Calculates the internal angles of a given 2D element.
std::vector< std::pair< Real, unsigned int > >	vertex_distances (const Elem &elem)
	Calculates the distances between the vertices of a given 2D element.
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.
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.
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.
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.
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.
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.
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.
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.
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.
void	lineRemoverMoveNode (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, const bool assign_ext_to_new, const bool side_to_remove)
bool	quasiTriElementsFixer (ReplicatedMesh &mesh, const std::set< subdomain_id_type > &subdomain_ids_set, const subdomain_id_type tri_elem_subdomain_shift, const SubdomainName tri_elem_subdomain_name_suffix)
void	triElemSplitter (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)
void	triElemSplitter (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)
void	quadElemSplitter (ReplicatedMesh &mesh, const dof_id_type elem_id, const subdomain_id_type tri_elem_subdomain_shift)
void	quadToTriOnLine (ReplicatedMesh &mesh, const std::vector< Real > &cut_line_params, const dof_id_type tri_subdomain_id_shift, const SubdomainName tri_elem_subdomain_name_suffix)
void	lineRemoverCutElemTri (ReplicatedMesh &mesh, const std::vector< Real > &cut_line_params, const subdomain_id_type block_id_to_remove, const boundary_id_type new_boundary_id)
void	lineRemoverCutElem (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)
void	boundaryTriElemImprover (ReplicatedMesh &mesh, const boundary_id_type boundary_to_improve)
void	makeImprovedTriElement (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, const std::vector< boundary_id_type > &boundary_ids_for_side_0, const std::vector< boundary_id_type > &boundary_ids_for_side_2)
bool	elemSideLocator (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)

Function Documentation

boundaryTriElemImprover() [1/2]

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

Referenced by lineRemoverCutElem().

boundaryTriElemImprover() [2/2]

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

Definition at line 1015 of file MooseMeshXYCuttingUtils.C.

{
  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 (const auto j : make_range(n_elem_extra_ids))
        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() [1/2]

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

Referenced by boundaryTriElemImprover().

elemSideLocator() [2/2]

bool MooseMeshXYCuttingUtils::elemSideLocator	(	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
	)

Definition at line 1249 of file MooseMeshXYCuttingUtils.C.

{
  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() [1/2]

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

Referenced by XYMeshLineCutter::generate().

lineRemoverCutElem() [2/2]

void MooseMeshXYCuttingUtils::lineRemoverCutElem	(	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
	)

Definition at line 996 of file MooseMeshXYCuttingUtils.C.

{
  // 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() [1/2]

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

Referenced by lineRemoverCutElem().

lineRemoverCutElemTri() [2/2]

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

Definition at line 774 of file MooseMeshXYCuttingUtils.C.

{
  // 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++)
  {
    const auto n_vertices = (*elem_it)->n_vertices();
    unsigned int n_points_on_line = 0;
    std::vector<unsigned short> node_side_rec(n_vertices, 0);
    for (const auto i : make_range(n_vertices))
    {
      // 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);
      if (pointOnLine(
              v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2]))
        ++n_points_on_line;
      else
        node_side_rec[i] = lineSideDeterminator(v_point(0),
                                                v_point(1),
                                                cut_line_params[0],
                                                cut_line_params[1],
                                                cut_line_params[2],
                                                true);
    }
    // This counts the booleans in node_side_rec, which does not include nodes
    // that are exactly on the line (these nodes are excluded from the
    // decision). In this case, num_nodes node lie on one side of the line and
    // node_side_rec.size() - n_nodes lie on the other side. In the case that
    // there are nodes on both sides of the line, we mark the element for
    // removal.
    const unsigned int num_nodes = std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0);
    if (num_nodes == node_side_rec.size() - n_points_on_line)
    {
      (*elem_it)->subdomain_id() = block_id_to_remove;
    }
    else if (num_nodes > 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 (const auto i : index_range(node_side_rec))
      {
        // 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 (const auto i : index_range(cross_elems))
  {
    // 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 (const auto j : make_range((*elem_it)->n_sides()))
      {
        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() [1/2]

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

Referenced by XYMeshLineCutter::generate().

lineRemoverMoveNode() [2/2]

void MooseMeshXYCuttingUtils::lineRemoverMoveNode	(	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,
		const bool	assign_ext_to_new,
		const bool	side_to_remove
	)

Definition at line 25 of file MooseMeshXYCuttingUtils.C.

{
  // 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 (const auto i : index_range(bdry_side_list))
    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 (const auto i : make_range((*elem_it)->n_vertices()))
    {
      // 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 (const auto i : make_range(mesh.elem_ptr(elem_id)->n_sides()))
    {
      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 (const auto i : make_range((*elem_it)->n_sides()))
    {
      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 (const auto i : index_range(slc_bdry_side_list))
  {
    // 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 (const auto i : index_range(node_list))
  {
    // 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 (const auto i : make_range((*elem_it)->n_sides()))
      {
        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 273 of file MooseMeshXYCuttingUtils.C.

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

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

makeImprovedTriElement() [1/2]

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

Referenced by boundaryTriElemImprover().

makeImprovedTriElement() [2/2]

void MooseMeshXYCuttingUtils::makeImprovedTriElement	(	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,
		const std::vector< boundary_id_type > &	boundary_ids_for_side_0,
		const std::vector< boundary_id_type > &	boundary_ids_for_side_2
	)

Definition at line 1219 of file MooseMeshXYCuttingUtils.C.

{
  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 (const auto j : index_range(extra_elem_ids))
  {
    elem_Tri3_new->set_extra_integer(j, extra_elem_ids[j]);
  }
}

pointOnLine()

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

Determines whether a point on XY-plane is on a given line, to within a tolerance.

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
dis_tol	tolerance used in determining whether the point is on the line

Returns

whether the point is on the line

Definition at line 262 of file MooseMeshXYCuttingUtils.C.

{
  return std::abs(px * param_1 + py * param_2 + param_3) <= dis_tol;
}

Referenced by lineRemoverCutElemTri(), and quadToTriOnLine().

quadElemSplitter() [1/2]

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

Referenced by quadToTriOnLine().

quadElemSplitter() [2/2]

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

Definition at line 611 of file MooseMeshXYCuttingUtils.C.

{
  // 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 (const auto i : index_range(bdry_side_list))
  {
    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 (const auto j : make_range(n_elem_extra_ids))
    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 (const auto j : make_range(n_elem_extra_ids))
    {
      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 (const auto j : make_range(n_elem_extra_ids))
    {
      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() [1/2]

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

Referenced by lineRemoverCutElem().

quadToTriOnLine() [2/2]

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

Definition at line 706 of file MooseMeshXYCuttingUtils.C.

{
  // 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 (const auto i : make_range(4))
      {
        const Point v_point = (*elem_it)->point(i);
        if (!pointOnLine(
                v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2]))
          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));
      }
      // This counts the booleans in node_side_rec, which does not include nodes
      // that are exactly on the line (these nodes are excluded from the
      // decision). In this case, num_nodes node lie on one side of the line and
      // node_side_rec.size() - n_nodes lie on the other side. In the case that
      // there are nodes on both sides of the line, we mark the element for
      // conversion.
      const auto num_nodes = std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0);
      if (num_nodes != (int)node_side_rec.size() && num_nodes > 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.set_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() [1/2]

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

Referenced by XYMeshLineCutter::generate().

quasiTriElementsFixer() [2/2]

bool MooseMeshXYCuttingUtils::quasiTriElementsFixer	(	ReplicatedMesh &	mesh,
		const std::set< subdomain_id_type > &	subdomain_ids_set,
		const subdomain_id_type	tri_elem_subdomain_shift,
		const SubdomainName	tri_elem_subdomain_name_suffix
	)

Definition at line 315 of file MooseMeshXYCuttingUtils.C.

{
  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 (const auto j : make_range(n_elem_extra_ids))
      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 (const auto j : make_range(n_elem_extra_ids))
      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.set_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/4]

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

Referenced by lineRemoverCutElemTri().

triElemSplitter() [2/4]

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

triElemSplitter() [3/4]

void MooseMeshXYCuttingUtils::triElemSplitter	(	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
	)

Definition at line 469 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();
 
  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 (const auto i : index_range(bdry_side_list))
  {
    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 (const auto j : make_range(n_elem_extra_ids))
    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 (const auto j : make_range(n_elem_extra_ids))
  {
    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() [4/4]

void MooseMeshXYCuttingUtils::triElemSplitter	(	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
	)

Definition at line 549 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 (const auto i : index_range(bdry_side_list))
  {
    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 (const auto j : make_range(n_elem_extra_ids))
    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 (const auto j : make_range(n_elem_extra_ids))
  {
    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 286 of file MooseMeshXYCuttingUtils.C.

{
  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);
}

Referenced by twoPointandLineIntersection().

twoPointandLineIntersection()

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

Definition at line 300 of file MooseMeshXYCuttingUtils.C.

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

Referenced by lineRemoverCutElemTri(), and lineRemoverMoveNode().

vertex_angles()

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

Definition at line 433 of file MooseMeshXYCuttingUtils.C.

{
  std::vector<std::pair<Real, unsigned int>> angles;
  const unsigned int n_vertices = elem.n_vertices();
 
  for (const auto i : make_range(n_vertices))
  {
    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;
}

Referenced by quasiTriElementsFixer().

vertex_distances()

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

Definition at line 454 of file MooseMeshXYCuttingUtils.C.

{
  std::vector<std::pair<Real, unsigned int>> distances;
  const unsigned int n_vertices = elem.n_vertices();
 
  for (const auto i : make_range(n_vertices))
  {
    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;
}

Referenced by quasiTriElementsFixer().