https://mooseframework.inl.gov
Functions
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 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. 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
meshinput mesh with the boundary TRI3 elements that need to be improved
boundary_to_improveboundary id of the boundary that needs to be improved

Definition at line 1015 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElem().

1016 {
1017  if (!MooseMeshUtils::hasBoundaryID(mesh, boundary_to_improve))
1018  mooseError(
1019  "MooseMeshXYCuttingUtils::boundaryTriElemImprover(): The boundary_to_improve provided "
1020  "does not exist in the given mesh.");
1021  BoundaryInfo & boundary_info = mesh.get_boundary_info();
1022  auto side_list = boundary_info.build_side_list();
1023  // Here we would like to collect the following information for all the TRI3 elements on the
1024  // boundary: Key: node id of the off-boundary node Value: a vector of tuples, each tuple contains
1025  // the following information:
1026  // 1. The element id of the element that is on the boundary to improve
1027  // 2. the one node id of that element that is on the boundary to improve
1028  // 3. the other node id of the element that is on the boundary to improve
1029  std::map<dof_id_type, std::vector<std::tuple<dof_id_type, dof_id_type, dof_id_type>>>
1030  tri3_elem_info;
1031  for (const auto & side : side_list)
1032  {
1033  if (std::get<2>(side) == boundary_to_improve)
1034  {
1035  Elem * elem = mesh.elem_ptr(std::get<0>(side));
1036  if (elem->type() == TRI3)
1037  {
1038  const auto key_node_id = elem->node_id((std::get<1>(side) + 2) % 3);
1039  const auto value_elem_id = elem->id();
1040  const auto value_node_id_1 = elem->node_id(std::get<1>(side));
1041  const auto value_node_id_2 = elem->node_id((std::get<1>(side) + 1) % 3);
1042  tri3_elem_info[key_node_id].push_back(
1043  std::make_tuple(value_elem_id, value_node_id_1, value_node_id_2));
1044  }
1045  }
1046  }
1047  // Elements that need to be removed
1048  std::vector<dof_id_type> elems_to_remove;
1049  // Now check if any group of TRI3 sharing an off-boundary node can be improved.
1050  for (const auto & tri_group : tri3_elem_info)
1051  {
1052  // It is possible to improve only when more than one TRI3 elements share the same off-boundary
1053  // node
1054  std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
1055  std::vector<dof_id_type> elem_id_list;
1056  for (const auto & tri : tri_group.second)
1057  {
1058  node_assm.push_back(std::make_pair(std::get<1>(tri), std::get<2>(tri)));
1059  elem_id_list.push_back(std::get<0>(tri));
1060  }
1061  std::vector<dof_id_type> ordered_node_list;
1062  std::vector<dof_id_type> ordered_elem_list;
1064  node_assm, elem_id_list, ordered_node_list, ordered_elem_list);
1065 
1066  // For all the elements sharing the same off-boundary node, we need to know how many separated
1067  // subdomains are involved
1068  // If there are extra element ids defined on the mesh, they also want to retain their boundaries
1069  // Only triangle elements that share a side can be merged
1070  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
1071  std::vector<std::tuple<subdomain_id_type, std::vector<dof_id_type>, unsigned int>> blocks_info;
1072  for (const auto & elem_id : ordered_elem_list)
1073  {
1074  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
1075  // Record all the element extra integers of the original quad element
1076  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
1077  exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
1078  if (!blocks_info.empty())
1079  {
1080  if (mesh.elem_ptr(elem_id)->subdomain_id() == std::get<0>(blocks_info.back()) &&
1081  exist_extra_ids == std::get<1>(blocks_info.back()))
1082  {
1083  std::get<2>(blocks_info.back())++;
1084  continue;
1085  }
1086  }
1087  blocks_info.push_back(
1088  std::make_tuple(mesh.elem_ptr(elem_id)->subdomain_id(), exist_extra_ids, 1));
1089  }
1090  // For each separated subdomain / set of extra ids, we try to improve the boundary elements
1091  unsigned int side_counter = 0;
1092  for (const auto & block_info : blocks_info)
1093  {
1094  const auto node_1 = mesh.node_ptr(ordered_node_list[side_counter]);
1095  // we do not need to subtract 1 for node_2
1096  const auto node_2 = mesh.node_ptr(ordered_node_list[side_counter + std::get<2>(block_info)]);
1097  const auto node_0 = mesh.node_ptr(tri_group.first);
1098  const Point v1 = *node_1 - *node_0;
1099  const Point v2 = *node_2 - *node_0;
1100  const Real angle = std::acos(v1 * v2 / v1.norm() / v2.norm()) / M_PI * 180.0;
1101  const std::vector<dof_id_type> block_elems(ordered_elem_list.begin() + side_counter,
1102  ordered_elem_list.begin() + side_counter +
1103  std::get<2>(block_info));
1104  // We assume that there are no sidesets defined inside a subdomain
1105  // For the first TRI3 element, we want to check if its side defined by node_0 and node_1 is
1106  // defined in any sidesets
1107  unsigned short side_id_0;
1108  unsigned short side_id_t;
1109  bool is_inverse_0;
1110  bool is_inverse_t;
1111  elemSideLocator(mesh,
1112  block_elems.front(),
1113  tri_group.first,
1114  ordered_node_list[side_counter],
1115  side_id_0,
1116  is_inverse_0);
1117  elemSideLocator(mesh,
1118  block_elems.back(),
1119  ordered_node_list[side_counter + std::get<2>(block_info)],
1120  tri_group.first,
1121  side_id_t,
1122  is_inverse_t);
1123  // Collect boundary information of the identified sides
1124  std::vector<boundary_id_type> side_0_boundary_ids;
1125  boundary_info.boundary_ids(
1126  mesh.elem_ptr(block_elems.front()), side_id_0, side_0_boundary_ids);
1127  std::vector<boundary_id_type> side_t_boundary_ids;
1128  boundary_info.boundary_ids(mesh.elem_ptr(block_elems.back()), side_id_t, side_t_boundary_ids);
1129 
1130  // Ideally we want this angle to be 60 degrees
1131  // In reality, we want one TRI3 element if the angle is less than 90 degrees;
1132  // we want two TRI3 elements if the angle is greater than 90 degrees and less than 135
1133  // degrees; we want three TRI3 elements if the angle is greater than 135 degrees and less than
1134  // 180 degrees.
1135  if (angle < 90.0)
1136  {
1137  if (std::get<2>(block_info) > 1)
1138  {
1140  tri_group.first,
1141  ordered_node_list[side_counter],
1142  ordered_node_list[side_counter + std::get<2>(block_info)],
1143  std::get<0>(block_info),
1144  std::get<1>(block_info),
1145  {boundary_to_improve},
1146  side_0_boundary_ids,
1147  side_t_boundary_ids);
1148  elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1149  }
1150  }
1151  else if (angle < 135.0)
1152  {
1153  // We can just add the middle node because there's nothing on the other side
1154  const auto node_m = mesh.add_point((*node_1 + *node_2) / 2.0);
1156  tri_group.first,
1157  ordered_node_list[side_counter],
1158  node_m->id(),
1159  std::get<0>(block_info),
1160  std::get<1>(block_info),
1161  {boundary_to_improve},
1162  side_0_boundary_ids,
1163  std::vector<boundary_id_type>());
1165  tri_group.first,
1166  node_m->id(),
1167  ordered_node_list[side_counter + std::get<2>(block_info)],
1168  std::get<0>(block_info),
1169  std::get<1>(block_info),
1170  {boundary_to_improve},
1171  std::vector<boundary_id_type>(),
1172  side_t_boundary_ids);
1173  elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1174  }
1175  else
1176  {
1177  const auto node_m1 = mesh.add_point((*node_1 * 2.0 + *node_2) / 3.0);
1178  const auto node_m2 = mesh.add_point((*node_1 + *node_2 * 2.0) / 3.0);
1180  tri_group.first,
1181  ordered_node_list[side_counter],
1182  node_m1->id(),
1183  std::get<0>(block_info),
1184  std::get<1>(block_info),
1185  {boundary_to_improve},
1186  side_0_boundary_ids,
1187  std::vector<boundary_id_type>());
1189  tri_group.first,
1190  node_m1->id(),
1191  node_m2->id(),
1192  std::get<0>(block_info),
1193  std::get<1>(block_info),
1194  {boundary_to_improve},
1195  std::vector<boundary_id_type>(),
1196  std::vector<boundary_id_type>());
1198  tri_group.first,
1199  node_m2->id(),
1200  ordered_node_list[side_counter + std::get<2>(block_info)],
1201  std::get<0>(block_info),
1202  std::get<1>(block_info),
1203  {boundary_to_improve},
1204  std::vector<boundary_id_type>(),
1205  side_t_boundary_ids);
1206  elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1207  }
1208  side_counter += std::get<2>(block_info);
1209  }
1210  // TODO: Need to check if the new element is inverted?
1211  }
1212  // Delete the original elements
1213  for (const auto & elem_to_remove : elems_to_remove)
1214  mesh.delete_elem(mesh.elem_ptr(elem_to_remove));
1215  mesh.contract();
1216 }
void makeOrderedNodeList(std::vector< std::pair< dof_id_type, dof_id_type >> &node_assm, std::vector< dof_id_type > &elem_id_list, std::vector< dof_id_type > &midpoint_node_list, std::vector< dof_id_type > &ordered_node_list, std::vector< dof_id_type > &ordered_elem_id_list)
Convert a list of sides in the form of a vector of pairs of node ids into a list of ordered nodes bas...
auto norm() const -> decltype(std::norm(Real()))
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:323
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.
unsigned int n_elem_integers() const
MeshBase & mesh
bool hasBoundaryID(const MeshBase &input_mesh, const BoundaryID id)
Whether a particular boundary ID exists in the mesh.
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
const BoundaryInfo & get_boundary_info() const
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
TRI3
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
virtual void delete_elem(Elem *e)=0
dof_id_type id() const
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.
virtual const Elem * elem_ptr(const dof_id_type i) const=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool contract()=0
subdomain_id_type subdomain_id() const
virtual const Node * node_ptr(const dof_id_type i) const=0
virtual ElemType type() const=0
dof_id_type node_id(const unsigned int i) const
dof_id_type get_extra_integer(const unsigned int index) const

◆ 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
meshinput mesh with the element that needs to be checked
elem_idid of the element that needs to be checked
node_id_0id of the first node of the pair
node_id_1id of the second node of the pair
side_idid of the side that contains the pair of nodes
is_inverseflag 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 1249 of file MooseMeshXYCuttingUtils.C.

Referenced by boundaryTriElemImprover().

1255 {
1256  Elem * elem = mesh.elem_ptr(elem_id);
1257  for (unsigned short i = 0; i < elem->n_sides(); i++)
1258  {
1259  if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_0 &&
1260  elem->side_ptr(i)->node_ptr(1)->id() == node_id_1)
1261  {
1262  side_id = i;
1263  is_inverse = false;
1264  return true;
1265  }
1266  else if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_1 &&
1267  elem->side_ptr(i)->node_ptr(1)->id() == node_id_0)
1268  {
1269  side_id = i;
1270  is_inverse = true;
1271  return true;
1272  }
1273  }
1274  return false;
1275 }
const boundary_id_type side_id
MeshBase & mesh
virtual const Elem * elem_ptr(const dof_id_type i) const=0
virtual unsigned int n_sides() const=0
virtual std::unique_ptr< Elem > side_ptr(unsigned int i)=0

◆ 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
meshinput mesh that need to be trimmed
cut_line_paramsparameters of the line that cuts the input mesh
tri_subdomain_id_shiftsubdomain id shift used to define the TRI element subdomains formed due to the trimming
tri_elem_subdomain_name_suffixsuffix used to name the TRI element subdomains formed due to the trimming
block_id_to_removea temporary subdomain id used to mark the elements that need to be removed
new_boundary_idboundary id of the new boundary that forms due to the trimming
improve_boundary_tri_elemsflag to indicate whether the boundary TRI3 elements need to be improved

Definition at line 996 of file MooseMeshXYCuttingUtils.C.

Referenced by XYMeshLineCutter::generate().

1003 {
1004  // Convert any quad elements crossed by the line into tri elements
1005  quadToTriOnLine(mesh, cut_line_params, tri_subdomain_id_shift, tri_elem_subdomain_name_suffix);
1006  // Then do the cutting for the preprocessed mesh that only contains tri elements crossed by the
1007  // cut line
1008  lineRemoverCutElemTri(mesh, cut_line_params, block_id_to_remove, new_boundary_id);
1009 
1010  if (improve_boundary_tri_elems)
1011  boundaryTriElemImprover(mesh, new_boundary_id);
1012 }
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 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.

◆ 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
meshinput mesh that need to be trimmed
cut_line_paramsparameters of the line that cuts the input mesh
block_id_to_removea temporary subdomain id used to mark the elements that need to be removed
new_boundary_idboundary id of the new boundary that forms due to the trimming

Definition at line 774 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElem().

778 {
779  // Find all the elements that are across the cutting line
780  std::vector<dof_id_type> cross_elems;
781  // A vector for element specific information
782  std::vector<std::vector<std::pair<dof_id_type, dof_id_type>>> node_pairs_vec;
783  // A set for unique pairs
784  std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs_unique_vec;
785  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
786  elem_it++)
787  {
788  const auto n_vertices = (*elem_it)->n_vertices();
789  unsigned int n_points_on_line = 0;
790  std::vector<unsigned short> node_side_rec(n_vertices, 0);
791  for (unsigned int i = 0; i < n_vertices; i++)
792  {
793  // First check if the vertex is in the XY Plane
794  if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
795  mooseError("MooseMeshXYCuttingUtils::lineRemoverCutElemTri() only works for 2D meshes in "
796  "XY Plane.");
797  const Point v_point = (*elem_it)->point(i);
798  if (pointOnLine(
799  v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2]))
800  ++n_points_on_line;
801  else
802  node_side_rec[i] = lineSideDeterminator(v_point(0),
803  v_point(1),
804  cut_line_params[0],
805  cut_line_params[1],
806  cut_line_params[2],
807  true);
808  }
809  // This counts the booleans in node_side_rec, which does not include nodes
810  // that are exactly on the line (these nodes are excluded from the
811  // decision). In this case, num_nodes node lie on one side of the line and
812  // node_side_rec.size() - n_nodes lie on the other side. In the case that
813  // there are nodes on both sides of the line, we mark the element for
814  // removal.
815  const unsigned int num_nodes = std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0);
816  if (num_nodes == node_side_rec.size() - n_points_on_line)
817  {
818  (*elem_it)->subdomain_id() = block_id_to_remove;
819  }
820  else if (num_nodes > 0)
821  {
822  if ((*elem_it)->n_vertices() != 3 || (*elem_it)->n_nodes() != 3)
823  mooseError("The element across the cutting line is not TRI3, which is not supported.");
824  cross_elems.push_back((*elem_it)->id());
825  // Then we need to check pairs of nodes that are on the different side
826  std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs;
827  for (unsigned int i = 0; i < node_side_rec.size(); i++)
828  {
829  // first node on removal side and second node on retaining side
830  if (node_side_rec[i] > 0 && node_side_rec[(i + 1) % node_side_rec.size()] == 0)
831  {
832  // Removal side first
833  node_pairs.push_back(
834  std::make_pair((*elem_it)->node_ptr(i)->id(),
835  (*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id()));
836  node_pairs_unique_vec.push_back(node_pairs.back());
837  }
838  // first node on retaining side and second node on removal side
839  else if (node_side_rec[i] == 0 && node_side_rec[(i + 1) % node_side_rec.size()] > 0)
840  {
841  // Removal side first
842  node_pairs.push_back(
843  std::make_pair((*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id(),
844  (*elem_it)->node_ptr(i)->id()));
845  node_pairs_unique_vec.push_back(node_pairs.back());
846  }
847  }
848  node_pairs_vec.push_back(node_pairs);
849  }
850  }
851  auto vec_ip = std::unique(node_pairs_unique_vec.begin(), node_pairs_unique_vec.end());
852  node_pairs_unique_vec.resize(std::distance(node_pairs_unique_vec.begin(), vec_ip));
853 
854  // Loop over all the node pairs to define new nodes that sit on the cutting line
855  std::vector<Node *> nodes_on_line;
856  // whether the on-line node is overlapped with the node pairs or a brand new node
857  std::vector<unsigned short> nodes_on_line_overlap;
858  for (const auto & node_pair : node_pairs_unique_vec)
859  {
860  const Point pt1 = *mesh.node_ptr(node_pair.first);
861  const Point pt2 = *mesh.node_ptr(node_pair.second);
862  const Point pt_line = twoPointandLineIntersection(
863  pt1, pt2, cut_line_params[0], cut_line_params[1], cut_line_params[2]);
864  if ((pt_line - pt1).norm() < libMesh::TOLERANCE)
865  {
866  nodes_on_line.push_back(mesh.node_ptr(node_pair.first));
867  nodes_on_line_overlap.push_back(1);
868  }
869  else if ((pt_line - pt2).norm() < libMesh::TOLERANCE)
870  {
871  nodes_on_line.push_back(mesh.node_ptr(node_pair.second));
872  nodes_on_line_overlap.push_back(2);
873  }
874  else
875  {
876  nodes_on_line.push_back(mesh.add_point(pt_line));
877  nodes_on_line_overlap.push_back(0);
878  }
879  }
880 
881  // make new elements
882  for (unsigned int i = 0; i < cross_elems.size(); i++)
883  {
884  // Only TRI elements are involved after preprocessing
885  auto cross_elem = mesh.elem_ptr(cross_elems[i]);
886  auto node_0 = cross_elem->node_ptr(0);
887  auto node_1 = cross_elem->node_ptr(1);
888  auto node_2 = cross_elem->node_ptr(2);
889  const std::vector<dof_id_type> tri_nodes = {node_0->id(), node_1->id(), node_2->id()};
890 
891  const auto online_node_index_1 = std::distance(node_pairs_unique_vec.begin(),
892  std::find(node_pairs_unique_vec.begin(),
893  node_pairs_unique_vec.end(),
894  node_pairs_vec[i][0]));
895  const auto online_node_index_2 = std::distance(node_pairs_unique_vec.begin(),
896  std::find(node_pairs_unique_vec.begin(),
897  node_pairs_unique_vec.end(),
898  node_pairs_vec[i][1]));
899  auto node_3 = nodes_on_line[online_node_index_1];
900  auto node_4 = nodes_on_line[online_node_index_2];
901  const auto node_3_overlap_flag = nodes_on_line_overlap[online_node_index_1];
902  const auto node_4_overlap_flag = nodes_on_line_overlap[online_node_index_2];
903  // Most common case, no overlapped nodes
904  if (node_3_overlap_flag == 0 && node_4_overlap_flag == 0)
905  {
906  // True if the common node is on the removal side; false if on the retaining side
907  const bool common_node_side = node_pairs_vec[i][0].first == node_pairs_vec[i][1].first;
908  const subdomain_id_type block_id_to_assign_1 =
909  common_node_side ? block_id_to_remove : cross_elem->subdomain_id();
910  const subdomain_id_type block_id_to_assign_2 =
911  common_node_side ? cross_elem->subdomain_id() : block_id_to_remove;
912  // The reference node ids need to be adjusted according to the common node of the two cut
913  // sides
914  const dof_id_type common_node_id =
915  common_node_side ? node_pairs_vec[i][0].first : node_pairs_vec[i][0].second;
916 
917  triElemSplitter(mesh,
918  cross_elem->id(),
919  std::distance(tri_nodes.begin(),
920  std::find(tri_nodes.begin(), tri_nodes.end(), common_node_id)),
921  node_3->id(),
922  node_4->id(),
923  block_id_to_assign_1,
924  block_id_to_assign_2);
925  mesh.delete_elem(cross_elem);
926  }
927  // both node_3 and node_4 are overlapped
928  else if (node_3_overlap_flag > 0 && node_4_overlap_flag > 0)
929  {
930  // In this case, the entire element is on one side of the cutting line
931  // No change needed just check which side the element is on
932  cross_elem->subdomain_id() = lineSideDeterminator(cross_elem->vertex_average()(0),
933  cross_elem->vertex_average()(1),
934  cut_line_params[0],
935  cut_line_params[1],
936  cut_line_params[2],
937  true)
938  ? block_id_to_remove
939  : cross_elem->subdomain_id();
940  }
941  // node_3 or node_4 is overlapped
942  else
943  {
944  const auto node_3_finder = std::distance(
945  tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_3->id()));
946  const auto node_4_finder = std::distance(
947  tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_4->id()));
948  // As only one of the two above values should be less than the three, the smaller one should
949  // be used
950  const dof_id_type node_id = node_3_finder < node_4_finder ? node_4->id() : node_3->id();
951  const auto node_finder = std::min(node_3_finder, node_4_finder);
952 
954  mesh,
955  cross_elem->id(),
956  node_finder,
957  node_id,
958  tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
959  ? block_id_to_remove
960  : cross_elem->subdomain_id(),
961  tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
962  ? cross_elem->subdomain_id()
963  : block_id_to_remove);
964  mesh.delete_elem(cross_elem);
965  }
966  }
967  mesh.contract();
968  // Due to the complexity, we identify the new boundary here together instead of during cutting of
969  // each element, because the preexisting element edges that are aligned with the cutting line also
970  // need to be added to the new boundary.
972  BoundaryInfo & boundary_info = mesh.get_boundary_info();
973  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
974  elem_it++)
975  {
976  if ((*elem_it)->subdomain_id() != block_id_to_remove)
977  {
978  for (unsigned int j = 0; j < (*elem_it)->n_sides(); j++)
979  {
980  if ((*elem_it)->neighbor_ptr(j) != nullptr)
981  if ((*elem_it)->neighbor_ptr(j)->subdomain_id() == block_id_to_remove)
982  boundary_info.add_side(*elem_it, j, new_boundary_id);
983  }
984  }
985  }
986 
987  // Delete the block to remove
988  for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
989  elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
990  elem_it++)
991  mesh.delete_elem(*elem_it);
992  mesh.contract();
993 }
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...
KOKKOS_INLINE_FUNCTION const T * find(const T &target, const T *const begin, const T *const end)
Find a value in an array.
Definition: KokkosUtils.h:42
Point twoPointandLineIntersection(const Point &pt1, const Point &pt2, const Real param_1, const Real param_2, const Real param_3)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:323
static constexpr Real TOLERANCE
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true)=0
virtual void delete_elem(Elem *e)=0
dof_id_type id() const
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. ...
auto norm(const T &a) -> decltype(std::abs(a))
virtual const Elem * elem_ptr(const dof_id_type i) const=0
virtual bool contract()=0
const Node * node_ptr(const unsigned int i) const
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...
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual const Node * node_ptr(const dof_id_type i) const=0
auto min(const L &left, const R &right)
uint8_t dof_id_type

◆ 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
meshinput mesh to perform line-based elements removing on
bdry_parsline parameter sets {a, b, c} as in a*x+b*y+c=0
block_id_to_removesubdomain id used to mark the elements that need to be removed
subdomain_ids_setall the subdomain ids in the input mesh
trimming_section_boundary_idID of the new external boundary formed due to trimming
external_boundary_idID of the external boundary of the input mesh
other_boundaries_to_conformIDs of the other boundaries that need to be conformed to during nodes moving
assign_ext_to_newwhether to assign external_boundary_id to the new boundary formed by removal
side_to_removewhich 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().

36 {
37  // Build boundary information of the mesh
38  BoundaryInfo & boundary_info = mesh.get_boundary_info();
39  auto bdry_side_list = boundary_info.build_side_list();
40  // Only select the boundaries_to_conform
41  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> slc_bdry_side_list;
42  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
43  if (std::get<2>(bdry_side_list[i]) == external_boundary_id ||
44  std::find(other_boundaries_to_conform.begin(),
45  other_boundaries_to_conform.end(),
46  std::get<2>(bdry_side_list[i])) != other_boundaries_to_conform.end())
47  slc_bdry_side_list.push_back(bdry_side_list[i]);
48 
49  // Assign block id for elements to be removed
50  // Also record the elements crossed by the line and with its average vertices on the removal side
51  std::vector<dof_id_type> crossed_elems_to_remove;
52  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
53  elem_it++)
54  {
55  // Check all the vertices of the element
56  unsigned short removal_side_count = 0;
57  for (unsigned int i = 0; i < (*elem_it)->n_vertices(); i++)
58  {
59  // First check if the vertex is on the XY-Plane
60  if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
61  mooseError(
62  "MooseMeshXYCuttingUtils::lineRemoverMoveNode() only works for 2D meshes in XY plane.");
63  if (lineSideDeterminator((*elem_it)->point(i)(0),
64  (*elem_it)->point(i)(1),
65  bdry_pars[0],
66  bdry_pars[1],
67  bdry_pars[2],
68  side_to_remove))
69  removal_side_count++;
70  }
71  if (removal_side_count == (*elem_it)->n_vertices())
72  {
73  (*elem_it)->subdomain_id() = block_id_to_remove;
74  continue;
75  }
76  // Check the average of the vertices of the element
77  if (lineSideDeterminator((*elem_it)->vertex_average()(0),
78  (*elem_it)->vertex_average()(1),
79  bdry_pars[0],
80  bdry_pars[1],
81  bdry_pars[2],
82  side_to_remove))
83  crossed_elems_to_remove.push_back((*elem_it)->id());
84  }
85  // Check each crossed element to see if removing it would lead to boundary moving
86  for (const auto & elem_id : crossed_elems_to_remove)
87  {
88  bool remove_flag = true;
89  for (unsigned int i = 0; i < mesh.elem_ptr(elem_id)->n_sides(); i++)
90  {
91  if (mesh.elem_ptr(elem_id)->neighbor_ptr(i) != nullptr)
92  if (mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id() != block_id_to_remove &&
93  std::find(crossed_elems_to_remove.begin(),
94  crossed_elems_to_remove.end(),
95  mesh.elem_ptr(elem_id)->neighbor_ptr(i)->id()) ==
96  crossed_elems_to_remove.end())
97  {
98  if (mesh.elem_ptr(elem_id)->subdomain_id() !=
99  mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id())
100  {
101  remove_flag = false;
102  break;
103  }
104  }
105  }
106  if (remove_flag)
107  mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
108  }
109 
110  // Identify all the nodes that are on the interface between block_id_to_remove and other blocks
111  // !!! We need a check here: if a node is on the retaining side, the removed element has a
112  // different subdomain id
113  std::vector<dof_id_type> node_list;
114  for (auto elem_it = mesh.active_subdomain_set_elements_begin(subdomain_ids_set);
115  elem_it != mesh.active_subdomain_set_elements_end(subdomain_ids_set);
116  elem_it++)
117  {
118  for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
119  {
120  if ((*elem_it)->neighbor_ptr(i) != nullptr)
121  if ((*elem_it)->neighbor_ptr(i)->subdomain_id() == block_id_to_remove)
122  {
123  node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(0)->id());
124  node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(1)->id());
125  boundary_info.add_side(*elem_it, i, trimming_section_boundary_id);
126  if (assign_ext_to_new && trimming_section_boundary_id != external_boundary_id)
127  boundary_info.add_side(*elem_it, i, external_boundary_id);
128  }
129  }
130  }
131  // Remove duplicate nodes
132  const auto unique_it = std::unique(node_list.begin(), node_list.end());
133  node_list.resize(std::distance(node_list.begin(), unique_it));
134  // Mark those nodes that are on a boundary that requires conformality
135  // If both nodes of a side are involved, we should only move one node
136  std::vector<bool> node_list_flag(node_list.size(), false);
137  std::vector<Point> node_list_point(node_list.size(), Point(0.0, 0.0, 0.0));
138  // Loop over all the selected sides
139  for (unsigned int i = 0; i < slc_bdry_side_list.size(); i++)
140  {
141  // Get the two node ids of the side
142  dof_id_type side_id_0 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
143  ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
144  ->node_ptr(0)
145  ->id();
146  dof_id_type side_id_1 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
147  ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
148  ->node_ptr(1)
149  ->id();
150  // True means the selected bdry node is in the node list of the trimming interface
151  bool side_id_0_in =
152  !(std::find(node_list.begin(), node_list.end(), side_id_0) == node_list.end());
153  bool side_id_1_in =
154  !(std::find(node_list.begin(), node_list.end(), side_id_1) == node_list.end());
155 
156  // True means the selected bdry node is on the removal side of the trimming interface
157  bool side_node_0_remove = lineSideDeterminator((*mesh.node_ptr(side_id_0))(0),
158  (*mesh.node_ptr(side_id_0))(1),
159  bdry_pars[0],
160  bdry_pars[1],
161  bdry_pars[2],
162  side_to_remove);
163  bool side_node_1_remove = lineSideDeterminator((*mesh.node_ptr(side_id_1))(0),
164  (*mesh.node_ptr(side_id_1))(1),
165  bdry_pars[0],
166  bdry_pars[1],
167  bdry_pars[2],
168  side_to_remove);
169  // If both nodes of that side are involved in the trimming interface
170  if (side_id_0_in && side_id_1_in)
171  // The side needs to be removed from the sideset because it is not longer an interface
172  // The other node will be handled by other element's side
173  boundary_info.remove_side(mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i])),
174  std::get<1>(slc_bdry_side_list[i]),
175  std::get<2>(slc_bdry_side_list[i]));
176  // If node 0 is on the trimming interface, and the side is cut by the trimming line
177  else if (side_id_0_in && (side_node_0_remove != side_node_1_remove))
178  {
179  // Use the intersection point as the destination of the node after moving
180  node_list_flag[std::distance(
181  node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_0))] = true;
182  const Point p0 = *mesh.node_ptr(side_id_0);
183  const Point p1 = *mesh.node_ptr(side_id_1);
184 
185  node_list_point[std::distance(node_list.begin(),
186  std::find(node_list.begin(), node_list.end(), side_id_0))] =
187  twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
188  }
189  // If node 1 is on the trimming interface, and the side is cut by the trimming line
190  else if (side_id_1_in && (side_node_0_remove != side_node_1_remove))
191  {
192  // Use the intersection point as the destination of the node after moving
193  node_list_flag[std::distance(
194  node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_1))] = true;
195  const Point p0 = *mesh.node_ptr(side_id_0);
196  const Point p1 = *mesh.node_ptr(side_id_1);
197 
198  node_list_point[std::distance(node_list.begin(),
199  std::find(node_list.begin(), node_list.end(), side_id_1))] =
200  twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
201  }
202  }
203 
204  // move nodes
205  for (unsigned int i = 0; i < node_list.size(); i++)
206  {
207  // This means the node is on both the trimming boundary and the original external
208  // boundary/selected interface boundaries. In order to keep the shape of the original external
209  // boundary, the node is moved along the original external boundary.
210  if (node_list_flag[i])
211  *(mesh.node_ptr(node_list[i])) = node_list_point[i];
212  // This means the node does not need to conform to any boundaries.
213  // Just move it along the normal direction of the trimming line.
214  else
215  {
216  const Real x0 = (*(mesh.node_ptr(node_list[i])))(0);
217  const Real y0 = (*(mesh.node_ptr(node_list[i])))(1);
218  (*(mesh.node_ptr(node_list[i])))(0) =
219  (bdry_pars[1] * (bdry_pars[1] * x0 - bdry_pars[0] * y0) - bdry_pars[0] * bdry_pars[2]) /
220  (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
221  (*(mesh.node_ptr(node_list[i])))(1) =
222  (bdry_pars[0] * (-bdry_pars[1] * x0 + bdry_pars[0] * y0) - bdry_pars[1] * bdry_pars[2]) /
223  (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
224  }
225  }
226 
227  // Delete the block
228  for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
229  elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
230  elem_it++)
231  mesh.delete_elem(*elem_it);
232  mesh.contract();
234  // Delete zero volume elements
235  std::vector<dof_id_type> zero_elems;
236  for (auto elem_it = mesh.elements_begin(); elem_it != mesh.elements_end(); elem_it++)
237  {
238  if (MooseUtils::absoluteFuzzyEqual((*elem_it)->volume(), 0.0))
239  {
240  for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
241  {
242  if ((*elem_it)->neighbor_ptr(i) != nullptr)
243  {
244  boundary_info.add_side((*elem_it)->neighbor_ptr(i),
245  ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
246  external_boundary_id);
247  boundary_info.add_side((*elem_it)->neighbor_ptr(i),
248  ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
249  trimming_section_boundary_id);
250  }
251  }
252  zero_elems.push_back((*elem_it)->id());
253  }
254  }
255  for (const auto & zero_elem : zero_elems)
256  mesh.delete_elem(mesh.elem_ptr(zero_elem));
257  mesh.contract();
258  // As we modified the side_list, it is safer to clear the node_list
259  boundary_info.clear_boundary_node_ids();
261 }
KOKKOS_INLINE_FUNCTION const T * find(const T &target, const T *const begin, const T *const end)
Find a value in an array.
Definition: KokkosUtils.h:42
Point twoPointandLineIntersection(const Point &pt1, const Point &pt2, const Real param_1, const Real param_2, const Real param_3)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:323
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true)=0
virtual void delete_elem(Elem *e)=0
dof_id_type id() const
virtual const Elem * elem_ptr(const dof_id_type i) const=0
virtual unsigned int n_sides() const=0
const Elem * neighbor_ptr(unsigned int i) const
void remove_side(const Elem *elem, const unsigned short int side)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool contract()=0
virtual std::unique_ptr< Elem > side_ptr(unsigned int i)=0
subdomain_id_type subdomain_id() const
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...
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual const Node * node_ptr(const dof_id_type i) const=0
uint8_t dof_id_type

◆ 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
pxx coordinate of the point
pyy coordinate of the point
param_1parameter 1 (a) in line formula a*x+b*y+c=0
param_2parameter 2 (b) in line formula a*x+b*y+c=0
param_3parameter 3 (c) in line formula a*x+b*y+c=0
direction_paramwhich side is the side that needs to be removed
dis_toltolerance used in determining side
Returns
whether the point is on the side of the line that needed to be removed

Definition at line 275 of file MooseMeshXYCuttingUtils.C.

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

282 {
283  const Real tmp = px * param_1 + py * param_2 + param_3;
284  return direction_param ? tmp >= dis_tol : tmp <= dis_tol;
285 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ 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
meshinput mesh where the TRI3 element needs to be added
node_id_0id of the first node of the TRI3 element
node_id_1id of the second node of the TRI3 element
node_id_2id of the third node of the TRI3 element
subdomain_idsubdomain id of the TRI3 element
extra_elem_idsextra element ids to be assigned to the TRI3 element
boundary_ids_for_side_1boundary ids of the second side of the TRI3 element
boundary_ids_for_side_0boundary ids of the first side of the TRI3 element
boundary_ids_for_side_2boundary ids of the third side of the TRI3 element

Definition at line 1219 of file MooseMeshXYCuttingUtils.C.

Referenced by boundaryTriElemImprover().

1228 {
1229  BoundaryInfo & boundary_info = mesh.get_boundary_info();
1230  Elem * elem_Tri3_new = mesh.add_elem(new Tri3);
1231  elem_Tri3_new->set_node(0, mesh.node_ptr(node_id_0));
1232  elem_Tri3_new->set_node(1, mesh.node_ptr(node_id_1));
1233  elem_Tri3_new->set_node(2, mesh.node_ptr(node_id_2));
1234  for (const auto & boundary_id_for_side_0 : boundary_ids_for_side_0)
1235  boundary_info.add_side(elem_Tri3_new, 0, boundary_id_for_side_0);
1236  for (const auto & boundary_id_for_side_1 : boundary_ids_for_side_1)
1237  boundary_info.add_side(elem_Tri3_new, 1, boundary_id_for_side_1);
1238  for (const auto & boundary_id_for_side_2 : boundary_ids_for_side_2)
1239  boundary_info.add_side(elem_Tri3_new, 2, boundary_id_for_side_2);
1240  elem_Tri3_new->subdomain_id() = subdomain_id;
1241  // Retain element extra integers
1242  for (unsigned int j = 0; j < extra_elem_ids.size(); j++)
1243  {
1244  elem_Tri3_new->set_extra_integer(j, extra_elem_ids[j]);
1245  }
1246 }
virtual Node *& set_node(const unsigned int i)
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
virtual Elem * add_elem(Elem *e)=0
subdomain_id_type subdomain_id() const
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual const Node * node_ptr(const dof_id_type i) const=0
void set_extra_integer(const unsigned int index, const dof_id_type value)

◆ 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
pxx coordinate of the point
pyy coordinate of the point
param_1parameter 1 (a) in line formula a*x+b*y+c=0
param_2parameter 2 (b) in line formula a*x+b*y+c=0
param_3parameter 3 (c) in line formula a*x+b*y+c=0
dis_toltolerance used in determining whether the point is on the line
Returns
whether the point is on the line

Definition at line 264 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElemTri(), and quadToTriOnLine().

270 {
271  return std::abs(px * param_1 + py * param_2 + param_3) <= dis_tol;
272 }
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:50

◆ 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
meshinput mesh with the QUAD4 element that needs to be split
elem_idid of the QUAD4 element that needs to be split
tri_elem_subdomain_shiftsubdomain id shift used to define the TRI element subdomains

Definition at line 612 of file MooseMeshXYCuttingUtils.C.

Referenced by quadToTriOnLine().

615 {
616  // Build boundary information of the mesh
617  BoundaryInfo & boundary_info = mesh.get_boundary_info();
618  auto bdry_side_list = boundary_info.build_side_list();
619  // Create a list of sidesets involving the element to be split
620  std::vector<std::vector<boundary_id_type>> elem_side_list;
621  elem_side_list.resize(4);
622  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
623  {
624  if (std::get<0>(bdry_side_list[i]) == elem_id)
625  {
626  elem_side_list[std::get<1>(bdry_side_list[i])].push_back(std::get<2>(bdry_side_list[i]));
627  }
628  }
629 
630  auto node_0 = mesh.elem_ptr(elem_id)->node_ptr(0);
631  auto node_1 = mesh.elem_ptr(elem_id)->node_ptr(1);
632  auto node_2 = mesh.elem_ptr(elem_id)->node_ptr(2);
633  auto node_3 = mesh.elem_ptr(elem_id)->node_ptr(3);
634 
635  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
636  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
637  // Record all the element extra integers of the original quad element
638  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
639  exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
640 
641  // There are two trivial ways to split a quad element
642  // We prefer the way that leads to triangles with similar areas
643  if (std::abs((*node_1 - *node_0).cross(*node_3 - *node_0).norm() -
644  (*node_1 - *node_2).cross(*node_3 - *node_2).norm()) >
645  std::abs((*node_0 - *node_1).cross(*node_2 - *node_1).norm() -
646  (*node_0 - *node_3).cross(*node_2 - *node_3).norm()))
647  {
648  Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
649  elem_Tri3_0->set_node(0, node_0);
650  elem_Tri3_0->set_node(1, node_1);
651  elem_Tri3_0->set_node(2, node_2);
652  elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
653  Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
654  elem_Tri3_1->set_node(0, node_0);
655  elem_Tri3_1->set_node(1, node_2);
656  elem_Tri3_1->set_node(2, node_3);
657  elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
658  // Retain element extra integers
659  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
660  {
661  elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
662  elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
663  }
664 
665  // Add sideset information to the new elements
666  for (const auto & side_info_0 : elem_side_list[0])
667  boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
668  for (const auto & side_info_1 : elem_side_list[1])
669  boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
670  for (const auto & side_info_2 : elem_side_list[2])
671  boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
672  for (const auto & side_info_3 : elem_side_list[3])
673  boundary_info.add_side(elem_Tri3_1, 2, side_info_3);
674  }
675  else
676  {
677  Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
678  elem_Tri3_0->set_node(0, node_0);
679  elem_Tri3_0->set_node(1, node_1);
680  elem_Tri3_0->set_node(2, node_3);
681  elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
682  Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
683  elem_Tri3_1->set_node(0, node_1);
684  elem_Tri3_1->set_node(1, node_2);
685  elem_Tri3_1->set_node(2, node_3);
686  elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
687  // Retain element extra integers
688  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
689  {
690  elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
691  elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
692  }
693 
694  // Add sideset information to the new elements
695  for (const auto & side_info_0 : elem_side_list[0])
696  boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
697  for (const auto & side_info_1 : elem_side_list[1])
698  boundary_info.add_side(elem_Tri3_1, 0, side_info_1);
699  for (const auto & side_info_2 : elem_side_list[2])
700  boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
701  for (const auto & side_info_3 : elem_side_list[3])
702  boundary_info.add_side(elem_Tri3_0, 2, side_info_3);
703  }
704 }
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:50
virtual Node *& set_node(const unsigned int i)
unsigned int n_elem_integers() const
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
virtual Elem * add_elem(Elem *e)=0
auto norm(const T &a) -> decltype(std::abs(a))
virtual const Elem * elem_ptr(const dof_id_type i) const=0
subdomain_id_type subdomain_id() const
const Node * node_ptr(const unsigned int i) const
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
void set_extra_integer(const unsigned int index, const dof_id_type value)
dof_id_type get_extra_integer(const unsigned int index) const

◆ 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
meshinput mesh with the QUAD4 elements that need to be converted
cut_line_paramsparameters of the line that cuts the input mesh
tri_subdomain_id_shiftsubdomain id shift used to define the TRI element subdomains generated due to the conversion
tri_elem_subdomain_name_suffixsuffix used to name the TRI element subdomains generated due to the conversion

Definition at line 707 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElem().

711 {
712  // Preprocess: find all the quad elements that are across the cutting line
713  std::vector<dof_id_type> cross_elems_quad;
714  std::set<subdomain_id_type> new_subdomain_ids;
715  for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
716  elem_it++)
717  {
718  if ((*elem_it)->n_vertices() == 4)
719  {
720  std::vector<unsigned short> node_side_rec;
721  for (unsigned int i = 0; i < 4; i++)
722  {
723  const Point v_point = (*elem_it)->point(i);
724  if (!pointOnLine(
725  v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2]))
726  node_side_rec.push_back(lineSideDeterminator(v_point(0),
727  v_point(1),
728  cut_line_params[0],
729  cut_line_params[1],
730  cut_line_params[2],
731  true));
732  }
733  // This counts the booleans in node_side_rec, which does not include nodes
734  // that are exactly on the line (these nodes are excluded from the
735  // decision). In this case, num_nodes node lie on one side of the line and
736  // node_side_rec.size() - n_nodes lie on the other side. In the case that
737  // there are nodes on both sides of the line, we mark the element for
738  // conversion.
739  const auto num_nodes = std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0);
740  if (num_nodes != (int)node_side_rec.size() && num_nodes > 0)
741  {
742  cross_elems_quad.push_back((*elem_it)->id());
743  new_subdomain_ids.emplace((*elem_it)->subdomain_id() + tri_subdomain_id_shift);
744  }
745  }
746  }
747  // Then convert these quad elements into tri elements
748  for (const auto & cross_elem_quad : cross_elems_quad)
749  {
750  quadElemSplitter(mesh, cross_elem_quad, tri_subdomain_id_shift);
751  mesh.delete_elem(mesh.elem_ptr(cross_elem_quad));
752  }
753  for (auto & nid : new_subdomain_ids)
754  {
755  const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
757  (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
758  : old_name) +
759  "_" + tri_elem_subdomain_name_suffix,
761  throw MooseException("The new subdomain name already exists in the mesh.");
762  mesh.subdomain_name(nid) =
763  (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
764  : old_name) +
765  "_" + tri_elem_subdomain_name_suffix;
766  mooseWarning("QUAD elements have been converted into TRI elements with a new "
767  "subdomain name: " +
768  mesh.subdomain_name(nid) + ".");
769  }
770  mesh.contract();
771 }
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:357
MeshBase & mesh
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.
SubdomainID getSubdomainID(const SubdomainName &subdomain_name, const MeshBase &mesh)
Gets the subdomain ID associated with the given SubdomainName.
const SubdomainID INVALID_BLOCK_ID
Definition: MooseTypes.C:20
virtual void delete_elem(Elem *e)=0
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. ...
std::string & subdomain_name(subdomain_id_type id)
virtual const Elem * elem_ptr(const dof_id_type i) const=0
Provides a way for users to bail out of the current solve.
virtual bool contract()=0
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...

◆ 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
meshinput mesh with degenerate QUAD elements that need to be fixed
subdomain_ids_setall the subdomain ids in the input mesh
tri_elem_subdomain_shiftsubdomain id shift used to define the TRI element subdomains
tri_elem_subdomain_name_suffixsuffix used to name the TRI element subdomains
Returns
whether any elements have been fixed

Definition at line 317 of file MooseMeshXYCuttingUtils.C.

Referenced by XYMeshLineCutter::generate().

321 {
322  BoundaryInfo & boundary_info = mesh.get_boundary_info();
323  // Define the subdomain id shift for the new TRI3 element subdomain(s)
324  const subdomain_id_type max_subdomain_id(*subdomain_ids_set.rbegin());
325  const subdomain_id_type tri_subdomain_id_shift =
326  tri_elem_subdomain_shift == Moose::INVALID_BLOCK_ID ? max_subdomain_id
327  : tri_elem_subdomain_shift;
328  mooseAssert(std::numeric_limits<subdomain_id_type>::max() - max_subdomain_id >
329  tri_subdomain_id_shift,
330  "The TRI elements subdomain id to be assigned may exceed the numeric limit.");
331  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
332  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
333  std::vector<std::tuple<Elem *, unsigned int, bool, bool>> bad_elems_rec;
334  // Loop over all the active elements to find any degenerate QUAD elements
335  for (auto & elem : as_range(mesh.active_elements_begin(), mesh.active_elements_end()))
336  {
337  // Two types of degenerate QUAD elements are identified:
338  // (1) QUAD elements with three collinear vertices
339  // (2) QUAD elements with two overlapped vertices
340  const auto elem_angles = vertex_angles(*elem);
341  const auto elem_distances = vertex_distances(*elem);
342  // Type 1
343  if (MooseUtils::absoluteFuzzyEqual(elem_angles.front().first, M_PI, 0.001))
344  {
345  bad_elems_rec.push_back(std::make_tuple(elem, elem_angles.front().second, false, true));
346  continue;
347  }
348  // Type 2
349  if (MooseUtils::absoluteFuzzyEqual(elem_distances.front().first, 0.0))
350  {
351  bad_elems_rec.push_back(std::make_tuple(elem, elem_distances.front().second, false, false));
352  }
353  }
354  std::set<subdomain_id_type> new_subdomain_ids;
355  // Loop over all the identified degenerate QUAD elements
356  for (const auto & bad_elem : bad_elems_rec)
357  {
358  std::vector<boundary_id_type> elem_bdry_container_0;
359  std::vector<boundary_id_type> elem_bdry_container_1;
360  std::vector<boundary_id_type> elem_bdry_container_2;
361 
362  Elem * elem_0 = std::get<0>(bad_elem);
363  if (std::get<3>(bad_elem))
364  {
365  // elems 1 and 2 are the neighboring elements of the degenerate element corresponding to the
366  // two collinear sides.
367  // For the degenerated element with three colinear vertices, if the elems 1 and 2 do not
368  // exist, the two sides are on the external boundary formed by trimming.
369  Elem * elem_1 = elem_0->neighbor_ptr(std::get<1>(bad_elem));
370  Elem * elem_2 = elem_0->neighbor_ptr((std::get<1>(bad_elem) - 1) % elem_0->n_vertices());
371  if ((elem_1 != nullptr || elem_2 != nullptr))
372  throw MooseException("The input mesh has degenerate quad element before trimming.");
373  }
374  mesh.get_boundary_info().boundary_ids(elem_0, std::get<1>(bad_elem), elem_bdry_container_0);
376  elem_0, (std::get<1>(bad_elem) + 1) % elem_0->n_vertices(), elem_bdry_container_1);
378  elem_0, (std::get<1>(bad_elem) + 2) % elem_0->n_vertices(), elem_bdry_container_2);
380  elem_0, (std::get<1>(bad_elem) + 3) % elem_0->n_vertices(), elem_bdry_container_0);
381 
382  // Record subdomain id of the degenerate element
383  auto elem_block_id = elem_0->subdomain_id();
384  // Define the three of four nodes that will be used to generate the TRI element
385  auto pt0 = elem_0->node_ptr((std::get<1>(bad_elem) + 1) % elem_0->n_vertices());
386  auto pt1 = elem_0->node_ptr((std::get<1>(bad_elem) + 2) % elem_0->n_vertices());
387  auto pt2 = elem_0->node_ptr((std::get<1>(bad_elem) + 3) % elem_0->n_vertices());
388  // Record all the element extra integers of the degenerate element
389  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
390  exist_extra_ids[j] = elem_0->get_extra_integer(j);
391  // Delete the degenerate QUAD element
392  mesh.delete_elem(elem_0);
393  // Create the new TRI element
394  Elem * elem_Tri3 = mesh.add_elem(new Tri3);
395  elem_Tri3->set_node(0, pt0);
396  elem_Tri3->set_node(1, pt1);
397  elem_Tri3->set_node(2, pt2);
398  // Retain the boundary information
399  for (auto bdry_id : elem_bdry_container_0)
400  boundary_info.add_side(elem_Tri3, 2, bdry_id);
401  for (auto bdry_id : elem_bdry_container_1)
402  boundary_info.add_side(elem_Tri3, 0, bdry_id);
403  for (auto bdry_id : elem_bdry_container_2)
404  boundary_info.add_side(elem_Tri3, 1, bdry_id);
405  // Assign subdomain id for the TRI element by shifting its original subdomain id
406  elem_Tri3->subdomain_id() = elem_block_id + tri_subdomain_id_shift;
407  new_subdomain_ids.emplace(elem_block_id + tri_subdomain_id_shift);
408  // Retain element extra integers
409  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
410  elem_Tri3->set_extra_integer(j, exist_extra_ids[j]);
411  }
412  // Assign subdomain names for the new TRI elements
413  for (auto & nid : new_subdomain_ids)
414  {
415  const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
417  (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
418  : old_name) +
419  "_" + tri_elem_subdomain_name_suffix,
421  throw MooseException("The new subdomain name already exists in the mesh.");
422  mesh.subdomain_name(nid) =
423  (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
424  : old_name) +
425  "_" + tri_elem_subdomain_name_suffix;
426  mooseWarning("Degenerate QUAD elements have been converted into TRI elements with a new "
427  "subdomain name: " +
428  mesh.subdomain_name(nid) + ".");
429  }
430  return bad_elems_rec.size();
431 }
std::vector< std::pair< Real, unsigned int > > vertex_angles(const Elem &elem)
virtual Node *& set_node(const unsigned int i)
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:357
unsigned int n_elem_integers() const
MeshBase & mesh
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
const BoundaryInfo & get_boundary_info() const
SubdomainID getSubdomainID(const SubdomainName &subdomain_name, const MeshBase &mesh)
Gets the subdomain ID associated with the given SubdomainName.
auto max(const L &left, const R &right)
const SubdomainID INVALID_BLOCK_ID
Definition: MooseTypes.C:20
virtual void delete_elem(Elem *e)=0
std::vector< std::pair< Real, unsigned int > > vertex_distances(const Elem &elem)
virtual Elem * add_elem(Elem *e)=0
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
std::string & subdomain_name(subdomain_id_type id)
const Elem * neighbor_ptr(unsigned int i) const
Provides a way for users to bail out of the current solve.
virtual unsigned int n_vertices() const=0
subdomain_id_type subdomain_id() const
const Node * node_ptr(const unsigned int i) const
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
void set_extra_integer(const unsigned int index, const dof_id_type value)
dof_id_type get_extra_integer(const unsigned int index) const

◆ 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
meshinput mesh with the TRI3 element that needs to be split
elem_idid of the TRI3 element that needs to be split
node_shiftshift used to rotate the vertices to make sure the vertex that the two cut sides share is the first vertex
nid_3id of the node on the first cut side of the triangle
nid_4id of the node on the second cut side of the triangle
single_elem_side_idsubdomain id of the single element side
double_elem_side_idsubdomain id of the double element side

Definition at line 470 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElemTri().

477 {
478  const auto elem_old = mesh.elem_ptr(elem_id);
479  const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
480  const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
481  const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
482 
483  const bool m1_side_flag =
484  MooseUtils::absoluteFuzzyEqual((*(mesh.node_ptr(nid_3)) - *(mesh.node_ptr(nid_0)))
485  .cross(*(mesh.node_ptr(nid_1)) - *(mesh.node_ptr(nid_0)))
486  .norm(),
487  0.0);
488  const dof_id_type nid_m1 = m1_side_flag ? nid_3 : nid_4;
489  const dof_id_type nid_m2 = m1_side_flag ? nid_4 : nid_3;
490  // Build boundary information of the mesh
491  BoundaryInfo & boundary_info = mesh.get_boundary_info();
492  auto bdry_side_list = boundary_info.build_side_list();
493  // Create a list of sidesets involving the element to be split
494  std::vector<std::vector<boundary_id_type>> elem_side_list;
495  elem_side_list.resize(3);
496  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
497  {
498  if (std::get<0>(bdry_side_list[i]) == elem_id)
499  {
500  elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
501  std::get<2>(bdry_side_list[i]));
502  }
503  }
504 
505  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
506  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
507  // Record all the element extra integers of the original element
508  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
509  exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
510 
511  Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
512  elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
513  elem_Tri3_0->set_node(1, mesh.node_ptr(nid_m1));
514  elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m2));
515  elem_Tri3_0->subdomain_id() = single_elem_side_id;
516  Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
517  elem_Tri3_1->set_node(0, mesh.node_ptr(nid_1));
518  elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m2));
519  elem_Tri3_1->set_node(2, mesh.node_ptr(nid_m1));
520  elem_Tri3_1->subdomain_id() = double_elem_side_id;
521  Elem * elem_Tri3_2 = mesh.add_elem(new Tri3);
522  elem_Tri3_2->set_node(0, mesh.node_ptr(nid_2));
523  elem_Tri3_2->set_node(1, mesh.node_ptr(nid_m2));
524  elem_Tri3_2->set_node(2, mesh.node_ptr(nid_1));
525  elem_Tri3_2->subdomain_id() = double_elem_side_id;
526  // Retain element extra integers
527  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
528  {
529  elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
530  elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
531  elem_Tri3_2->set_extra_integer(j, exist_extra_ids[j]);
532  }
533 
534  // Add sideset information to the new elements
535  for (const auto & side_info_0 : elem_side_list[0])
536  {
537  boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
538  boundary_info.add_side(elem_Tri3_1, 2, side_info_0);
539  }
540  for (const auto & side_info_1 : elem_side_list[1])
541  boundary_info.add_side(elem_Tri3_2, 2, side_info_1);
542  for (const auto & side_info_2 : elem_side_list[2])
543  {
544  boundary_info.add_side(elem_Tri3_0, 2, side_info_2);
545  boundary_info.add_side(elem_Tri3_2, 0, side_info_2);
546  }
547 }
virtual Node *& set_node(const unsigned int i)
unsigned int n_elem_integers() const
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
dof_id_type id() const
virtual Elem * add_elem(Elem *e)=0
virtual const Elem * elem_ptr(const dof_id_type i) const=0
subdomain_id_type subdomain_id() const
const Node * node_ptr(const unsigned int i) const
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual const Node * node_ptr(const dof_id_type i) const=0
void set_extra_integer(const unsigned int index, const dof_id_type value)
dof_id_type get_extra_integer(const unsigned int index) const
uint8_t dof_id_type

◆ 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
meshinput mesh with the TRI3 element that needs to be split
elem_idid of the TRI3 element that needs to be split
node_shiftshift used to rotate the vertices to make sure the vertex corresponding to the cut side is the first vertex
nid_mid of the node on the cut side of the triangle
first_elem_side_idsubdomain id of the first element side
second_elem_side_idsubdomain id of the second element side

Definition at line 550 of file MooseMeshXYCuttingUtils.C.

556 {
557  const auto elem_old = mesh.elem_ptr(elem_id);
558  const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
559  const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
560  const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
561  // Build boundary information of the mesh
562  BoundaryInfo & boundary_info = mesh.get_boundary_info();
563  auto bdry_side_list = boundary_info.build_side_list();
564  // Create a list of sidesets involving the element to be split
565  std::vector<std::vector<boundary_id_type>> elem_side_list;
566  elem_side_list.resize(3);
567  for (unsigned int i = 0; i < bdry_side_list.size(); i++)
568  {
569  if (std::get<0>(bdry_side_list[i]) == elem_id)
570  {
571  elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
572  std::get<2>(bdry_side_list[i]));
573  }
574  }
575 
576  const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
577  std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
578  // Record all the element extra integers of the original element
579  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
580  exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
581 
582  Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
583  elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
584  elem_Tri3_0->set_node(1, mesh.node_ptr(nid_1));
585  elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m));
586  elem_Tri3_0->subdomain_id() = first_elem_side_id;
587  Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
588  elem_Tri3_1->set_node(0, mesh.node_ptr(nid_0));
589  elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m));
590  elem_Tri3_1->set_node(2, mesh.node_ptr(nid_2));
591  elem_Tri3_1->subdomain_id() = second_elem_side_id;
592  // Retain element extra integers
593  for (unsigned int j = 0; j < n_elem_extra_ids; j++)
594  {
595  elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
596  elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
597  }
598 
599  // Add sideset information to the new elements
600  for (const auto & side_info_0 : elem_side_list[0])
601  boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
602  for (const auto & side_info_1 : elem_side_list[1])
603  {
604  boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
605  boundary_info.add_side(elem_Tri3_1, 1, side_info_1);
606  }
607  for (const auto & side_info_2 : elem_side_list[2])
608  boundary_info.add_side(elem_Tri3_1, 2, side_info_2);
609 }
virtual Node *& set_node(const unsigned int i)
unsigned int n_elem_integers() const
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
dof_id_type id() const
virtual Elem * add_elem(Elem *e)=0
virtual const Elem * elem_ptr(const dof_id_type i) const=0
subdomain_id_type subdomain_id() const
const Node * node_ptr(const unsigned int i) const
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual const Node * node_ptr(const dof_id_type i) const=0
void set_extra_integer(const unsigned int index, const dof_id_type value)
dof_id_type get_extra_integer(const unsigned int index) const
uint8_t dof_id_type

◆ 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_11parameter 1 (a) in line formula a*x+b*y+c=0 for the first line
param_12parameter 2 (b) in line formula a*x+b*y+c=0 for the first line
param_13parameter 3 (c) in line formula a*x+b*y+c=0 for the first line
param_21parameter 1 (a) in line formula a*x+b*y+c=0 for the second line
param_22parameter 2 (b) in line formula a*x+b*y+c=0 for the second line
param_23parameter 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 288 of file MooseMeshXYCuttingUtils.C.

Referenced by twoPointandLineIntersection().

294 {
295  return Point(
296  (param_12 * param_23 - param_22 * param_13) / (param_11 * param_22 - param_21 * param_12),
297  (param_13 * param_21 - param_23 * param_11) / (param_11 * param_22 - param_21 * param_12),
298  0.0);
299 }

◆ 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
pt1point 1 that defines the first straight line
pt2point 2 that defines the first straight line
param_1parameter 1 (a) in line formula a*x+b*y+c=0 for the second straight line
param_2parameter 2 (b) in line formula a*x+b*y+c=0 for the second straight line
param_3parameter 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 302 of file MooseMeshXYCuttingUtils.C.

Referenced by lineRemoverCutElemTri(), and lineRemoverMoveNode().

307 {
308  return twoLineIntersection(param_1,
309  param_2,
310  param_3,
311  pt2(1) - pt1(1),
312  pt1(0) - pt2(0),
313  pt2(0) * pt1(1) - pt1(0) * pt2(1));
314 }
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.

◆ 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
elemthe 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 434 of file MooseMeshXYCuttingUtils.C.

Referenced by quasiTriElementsFixer().

435 {
436  std::vector<std::pair<Real, unsigned int>> angles;
437  const unsigned int n_vertices = elem.n_vertices();
438 
439  for (unsigned int i = 0; i < n_vertices; i++)
440  {
441  Point v1 = (*elem.node_ptr((i - 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
442  Point v2 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
443  Real tmp = v1 * v2 / v1.norm() / v2.norm();
444  if (tmp > 1.0)
445  tmp = 1.0;
446  else if (tmp < -1.0)
447  tmp = -1.0;
448  angles.push_back(std::make_pair(acos(tmp), i));
449  }
450  std::sort(angles.begin(), angles.end(), std::greater<>());
451  return angles;
452 }
auto norm() const -> decltype(std::norm(Real()))
virtual unsigned int n_vertices() const=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Node * node_ptr(const unsigned int i) const

◆ 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
elemthe 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 455 of file MooseMeshXYCuttingUtils.C.

Referenced by quasiTriElementsFixer().

456 {
457  std::vector<std::pair<Real, unsigned int>> distances;
458  const unsigned int n_vertices = elem.n_vertices();
459 
460  for (unsigned int i = 0; i < n_vertices; i++)
461  {
462  Point v1 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
463  distances.push_back(std::make_pair(v1.norm(), i));
464  }
465  std::sort(distances.begin(), distances.end());
466  return distances;
467 }
auto norm() const -> decltype(std::norm(Real()))
virtual unsigned int n_vertices() const=0
const Node * node_ptr(const unsigned int i) const