Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // MOOSE includes
11 : #include "MooseMeshXYCuttingUtils.h"
12 : #include "MooseMeshUtils.h"
13 :
14 : #include "libmesh/elem.h"
15 : #include "libmesh/boundary_info.h"
16 : #include "libmesh/mesh_base.h"
17 : #include "libmesh/parallel.h"
18 : #include "libmesh/parallel_algebra.h"
19 : #include "libmesh/face_tri3.h"
20 :
21 : using namespace libMesh;
22 :
23 : namespace MooseMeshXYCuttingUtils
24 : {
25 :
26 : void
27 28 : lineRemoverMoveNode(ReplicatedMesh & mesh,
28 : const std::vector<Real> & bdry_pars,
29 : const subdomain_id_type block_id_to_remove,
30 : const std::set<subdomain_id_type> & subdomain_ids_set,
31 : const boundary_id_type trimming_section_boundary_id,
32 : const boundary_id_type external_boundary_id,
33 : const std::vector<boundary_id_type> & other_boundaries_to_conform,
34 : const bool assign_ext_to_new,
35 : const bool side_to_remove)
36 : {
37 : // Build boundary information of the mesh
38 28 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
39 28 : auto bdry_side_list = boundary_info.build_side_list();
40 : // Only select the boundaries_to_conform
41 28 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> slc_bdry_side_list;
42 2780 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
43 4560 : if (std::get<2>(bdry_side_list[i]) == external_boundary_id ||
44 1808 : std::find(other_boundaries_to_conform.begin(),
45 : other_boundaries_to_conform.end(),
46 6368 : std::get<2>(bdry_side_list[i])) != other_boundaries_to_conform.end())
47 1712 : 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 28 : std::vector<dof_id_type> crossed_elems_to_remove;
52 6876 : 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 6848 : unsigned short removal_side_count = 0;
57 34240 : for (unsigned int i = 0; i < (*elem_it)->n_vertices(); i++)
58 : {
59 : // First check if the vertex is on the XY-Plane
60 27392 : if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
61 0 : mooseError(
62 : "MooseMeshXYCuttingUtils::lineRemoverMoveNode() only works for 2D meshes in XY plane.");
63 27392 : if (lineSideDeterminator((*elem_it)->point(i)(0),
64 27392 : (*elem_it)->point(i)(1),
65 27392 : bdry_pars[0],
66 27392 : bdry_pars[1],
67 27392 : bdry_pars[2],
68 : side_to_remove))
69 16880 : removal_side_count++;
70 : }
71 6848 : if (removal_side_count == (*elem_it)->n_vertices())
72 : {
73 3848 : (*elem_it)->subdomain_id() = block_id_to_remove;
74 3848 : continue;
75 : }
76 : // Check the average of the vertices of the element
77 3000 : if (lineSideDeterminator((*elem_it)->vertex_average()(0),
78 6000 : (*elem_it)->vertex_average()(1),
79 3000 : bdry_pars[0],
80 3000 : bdry_pars[1],
81 3000 : bdry_pars[2],
82 : side_to_remove))
83 440 : crossed_elems_to_remove.push_back((*elem_it)->id());
84 28 : }
85 : // Check each crossed element to see if removing it would lead to boundary moving
86 468 : for (const auto & elem_id : crossed_elems_to_remove)
87 : {
88 440 : bool remove_flag = true;
89 1848 : for (unsigned int i = 0; i < mesh.elem_ptr(elem_id)->n_sides(); i++)
90 : {
91 1520 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i) != nullptr)
92 2260 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id() != block_id_to_remove &&
93 804 : std::find(crossed_elems_to_remove.begin(),
94 : crossed_elems_to_remove.end(),
95 1608 : mesh.elem_ptr(elem_id)->neighbor_ptr(i)->id()) ==
96 2260 : crossed_elems_to_remove.end())
97 : {
98 576 : if (mesh.elem_ptr(elem_id)->subdomain_id() !=
99 576 : mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id())
100 : {
101 112 : remove_flag = false;
102 112 : break;
103 : }
104 : }
105 : }
106 440 : if (remove_flag)
107 328 : 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 28 : std::vector<dof_id_type> node_list;
114 2700 : for (auto elem_it = mesh.active_subdomain_set_elements_begin(subdomain_ids_set);
115 2700 : elem_it != mesh.active_subdomain_set_elements_end(subdomain_ids_set);
116 : elem_it++)
117 : {
118 13360 : for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
119 : {
120 10688 : if ((*elem_it)->neighbor_ptr(i) != nullptr)
121 10176 : if ((*elem_it)->neighbor_ptr(i)->subdomain_id() == block_id_to_remove)
122 : {
123 760 : node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(0)->id());
124 760 : node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(1)->id());
125 760 : boundary_info.add_side(*elem_it, i, trimming_section_boundary_id);
126 760 : if (assign_ext_to_new && trimming_section_boundary_id != external_boundary_id)
127 0 : boundary_info.add_side(*elem_it, i, external_boundary_id);
128 : }
129 : }
130 28 : }
131 : // Remove duplicate nodes
132 28 : const auto unique_it = std::unique(node_list.begin(), node_list.end());
133 28 : 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 28 : std::vector<bool> node_list_flag(node_list.size(), false);
137 28 : std::vector<Point> node_list_point(node_list.size(), Point(0.0, 0.0, 0.0));
138 : // Loop over all the selected sides
139 1740 : for (unsigned int i = 0; i < slc_bdry_side_list.size(); i++)
140 : {
141 : // Get the two node ids of the side
142 1712 : dof_id_type side_id_0 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
143 1712 : ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
144 1712 : ->node_ptr(0)
145 1712 : ->id();
146 1712 : dof_id_type side_id_1 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
147 1712 : ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
148 1712 : ->node_ptr(1)
149 1712 : ->id();
150 : // True means the selected bdry node is in the node list of the trimming interface
151 : bool side_id_0_in =
152 1712 : !(std::find(node_list.begin(), node_list.end(), side_id_0) == node_list.end());
153 : bool side_id_1_in =
154 1712 : !(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 1712 : bool side_node_0_remove = lineSideDeterminator((*mesh.node_ptr(side_id_0))(0),
158 1712 : (*mesh.node_ptr(side_id_0))(1),
159 1712 : bdry_pars[0],
160 1712 : bdry_pars[1],
161 1712 : bdry_pars[2],
162 : side_to_remove);
163 1712 : bool side_node_1_remove = lineSideDeterminator((*mesh.node_ptr(side_id_1))(0),
164 1712 : (*mesh.node_ptr(side_id_1))(1),
165 1712 : bdry_pars[0],
166 1712 : bdry_pars[1],
167 1712 : bdry_pars[2],
168 : side_to_remove);
169 : // If both nodes of that side are involved in the trimming interface
170 1712 : 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 8 : boundary_info.remove_side(mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i])),
174 8 : std::get<1>(slc_bdry_side_list[i]),
175 8 : 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 1704 : 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 68 : node_list_flag[std::distance(
181 68 : node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_0))] = true;
182 68 : const Point p0 = *mesh.node_ptr(side_id_0);
183 68 : const Point p1 = *mesh.node_ptr(side_id_1);
184 :
185 68 : node_list_point[std::distance(node_list.begin(),
186 68 : std::find(node_list.begin(), node_list.end(), side_id_0))] =
187 68 : twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
188 68 : }
189 : // If node 1 is on the trimming interface, and the side is cut by the trimming line
190 1636 : 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 48 : node_list_flag[std::distance(
194 48 : node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_1))] = true;
195 48 : const Point p0 = *mesh.node_ptr(side_id_0);
196 48 : const Point p1 = *mesh.node_ptr(side_id_1);
197 :
198 48 : node_list_point[std::distance(node_list.begin(),
199 48 : std::find(node_list.begin(), node_list.end(), side_id_1))] =
200 48 : twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
201 : }
202 : }
203 :
204 : // move nodes
205 1268 : 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 1240 : if (node_list_flag[i])
211 116 : *(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 1124 : const Real x0 = (*(mesh.node_ptr(node_list[i])))(0);
217 1124 : const Real y0 = (*(mesh.node_ptr(node_list[i])))(1);
218 1124 : (*(mesh.node_ptr(node_list[i])))(0) =
219 1124 : (bdry_pars[1] * (bdry_pars[1] * x0 - bdry_pars[0] * y0) - bdry_pars[0] * bdry_pars[2]) /
220 1124 : (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
221 1124 : (*(mesh.node_ptr(node_list[i])))(1) =
222 1124 : (bdry_pars[0] * (-bdry_pars[1] * x0 + bdry_pars[0] * y0) - bdry_pars[1] * bdry_pars[2]) /
223 1124 : (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
224 : }
225 : }
226 :
227 : // Delete the block
228 4204 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
229 4204 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
230 : elem_it++)
231 4204 : mesh.delete_elem(*elem_it);
232 28 : mesh.contract();
233 28 : mesh.find_neighbors();
234 : // Delete zero volume elements
235 28 : std::vector<dof_id_type> zero_elems;
236 2700 : for (auto elem_it = mesh.elements_begin(); elem_it != mesh.elements_end(); elem_it++)
237 : {
238 2672 : if (MooseUtils::absoluteFuzzyEqual((*elem_it)->volume(), 0.0))
239 : {
240 80 : for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
241 : {
242 64 : if ((*elem_it)->neighbor_ptr(i) != nullptr)
243 : {
244 32 : boundary_info.add_side((*elem_it)->neighbor_ptr(i),
245 32 : ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
246 : external_boundary_id);
247 32 : boundary_info.add_side((*elem_it)->neighbor_ptr(i),
248 32 : ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
249 : trimming_section_boundary_id);
250 : }
251 : }
252 16 : zero_elems.push_back((*elem_it)->id());
253 : }
254 28 : }
255 44 : for (const auto & zero_elem : zero_elems)
256 16 : mesh.delete_elem(mesh.elem_ptr(zero_elem));
257 28 : mesh.contract();
258 : // As we modified the side_list, it is safer to clear the node_list
259 28 : boundary_info.clear_boundary_node_ids();
260 28 : mesh.prepare_for_use();
261 28 : }
262 :
263 : bool
264 86864 : lineSideDeterminator(const Real px,
265 : const Real py,
266 : const Real param_1,
267 : const Real param_2,
268 : const Real param_3,
269 : const bool direction_param,
270 : const Real dis_tol)
271 : {
272 86864 : const Real tmp = px * param_1 + py * param_2 + param_3;
273 86864 : return direction_param ? tmp >= dis_tol : tmp <= dis_tol;
274 : }
275 :
276 : Point
277 2212 : twoLineIntersection(const Real param_11,
278 : const Real param_12,
279 : const Real param_13,
280 : const Real param_21,
281 : const Real param_22,
282 : const Real param_23)
283 : {
284 : return Point(
285 2212 : (param_12 * param_23 - param_22 * param_13) / (param_11 * param_22 - param_21 * param_12),
286 2212 : (param_13 * param_21 - param_23 * param_11) / (param_11 * param_22 - param_21 * param_12),
287 2212 : 0.0);
288 : }
289 :
290 : Point
291 2212 : twoPointandLineIntersection(const Point & pt1,
292 : const Point & pt2,
293 : const Real param_1,
294 : const Real param_2,
295 : const Real param_3)
296 : {
297 6636 : return twoLineIntersection(param_1,
298 : param_2,
299 : param_3,
300 2212 : pt2(1) - pt1(1),
301 2212 : pt1(0) - pt2(0),
302 4424 : pt2(0) * pt1(1) - pt1(0) * pt2(1));
303 : }
304 :
305 : bool
306 28 : quasiTriElementsFixer(ReplicatedMesh & mesh,
307 : const std::set<subdomain_id_type> & subdomain_ids_set,
308 : const subdomain_id_type tri_elem_subdomain_shift,
309 : const SubdomainName tri_elem_subdomain_name_suffix)
310 : {
311 28 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
312 : // Define the subdomain id shift for the new TRI3 element subdomain(s)
313 28 : const subdomain_id_type max_subdomain_id(*subdomain_ids_set.rbegin());
314 28 : const subdomain_id_type tri_subdomain_id_shift =
315 28 : tri_elem_subdomain_shift == Moose::INVALID_BLOCK_ID ? max_subdomain_id
316 : : tri_elem_subdomain_shift;
317 : mooseAssert(std::numeric_limits<subdomain_id_type>::max() - max_subdomain_id >
318 : tri_subdomain_id_shift,
319 : "The TRI elements subdomain id to be assigned may exceed the numeric limit.");
320 28 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
321 28 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
322 28 : std::vector<std::tuple<Elem *, unsigned int, bool, bool>> bad_elems_rec;
323 : // Loop over all the active elements to find any degenerate QUAD elements
324 5340 : for (auto & elem : as_range(mesh.active_elements_begin(), mesh.active_elements_end()))
325 : {
326 : // Two types of degenerate QUAD elements are identified:
327 : // (1) QUAD elements with three collinear vertices
328 : // (2) QUAD elements with two overlapped vertices
329 2656 : const auto elem_angles = vertex_angles(*elem);
330 2656 : const auto elem_distances = vertex_distances(*elem);
331 : // Type 1
332 2656 : if (MooseUtils::absoluteFuzzyEqual(elem_angles.front().first, M_PI, 0.001))
333 : {
334 208 : bad_elems_rec.push_back(std::make_tuple(elem, elem_angles.front().second, false, true));
335 208 : continue;
336 : }
337 : // Type 2
338 2448 : if (MooseUtils::absoluteFuzzyEqual(elem_distances.front().first, 0.0))
339 : {
340 8 : bad_elems_rec.push_back(std::make_tuple(elem, elem_distances.front().second, false, false));
341 : }
342 2892 : }
343 28 : std::set<subdomain_id_type> new_subdomain_ids;
344 : // Loop over all the identified degenerate QUAD elements
345 244 : for (const auto & bad_elem : bad_elems_rec)
346 : {
347 216 : std::vector<boundary_id_type> elem_bdry_container_0;
348 216 : std::vector<boundary_id_type> elem_bdry_container_1;
349 216 : std::vector<boundary_id_type> elem_bdry_container_2;
350 :
351 216 : Elem * elem_0 = std::get<0>(bad_elem);
352 216 : if (std::get<3>(bad_elem))
353 : {
354 : // elems 1 and 2 are the neighboring elements of the degenerate element corresponding to the
355 : // two collinear sides.
356 : // For the degenerated element with three colinear vertices, if the elems 1 and 2 do not
357 : // exist, the two sides are on the external boundary formed by trimming.
358 208 : Elem * elem_1 = elem_0->neighbor_ptr(std::get<1>(bad_elem));
359 208 : Elem * elem_2 = elem_0->neighbor_ptr((std::get<1>(bad_elem) - 1) % elem_0->n_vertices());
360 208 : if ((elem_1 != nullptr || elem_2 != nullptr))
361 0 : throw MooseException("The input mesh has degenerate quad element before trimming.");
362 : }
363 216 : mesh.get_boundary_info().boundary_ids(elem_0, std::get<1>(bad_elem), elem_bdry_container_0);
364 216 : mesh.get_boundary_info().boundary_ids(
365 216 : elem_0, (std::get<1>(bad_elem) + 1) % elem_0->n_vertices(), elem_bdry_container_1);
366 216 : mesh.get_boundary_info().boundary_ids(
367 216 : elem_0, (std::get<1>(bad_elem) + 2) % elem_0->n_vertices(), elem_bdry_container_2);
368 216 : mesh.get_boundary_info().boundary_ids(
369 216 : elem_0, (std::get<1>(bad_elem) + 3) % elem_0->n_vertices(), elem_bdry_container_0);
370 :
371 : // Record subdomain id of the degenerate element
372 216 : auto elem_block_id = elem_0->subdomain_id();
373 : // Define the three of four nodes that will be used to generate the TRI element
374 216 : auto pt0 = elem_0->node_ptr((std::get<1>(bad_elem) + 1) % elem_0->n_vertices());
375 216 : auto pt1 = elem_0->node_ptr((std::get<1>(bad_elem) + 2) % elem_0->n_vertices());
376 216 : auto pt2 = elem_0->node_ptr((std::get<1>(bad_elem) + 3) % elem_0->n_vertices());
377 : // Record all the element extra integers of the degenerate element
378 216 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
379 0 : exist_extra_ids[j] = elem_0->get_extra_integer(j);
380 : // Delete the degenerate QUAD element
381 216 : mesh.delete_elem(elem_0);
382 : // Create the new TRI element
383 216 : Elem * elem_Tri3 = mesh.add_elem(new Tri3);
384 216 : elem_Tri3->set_node(0, pt0);
385 216 : elem_Tri3->set_node(1, pt1);
386 216 : elem_Tri3->set_node(2, pt2);
387 : // Retain the boundary information
388 432 : for (auto bdry_id : elem_bdry_container_0)
389 216 : boundary_info.add_side(elem_Tri3, 2, bdry_id);
390 256 : for (auto bdry_id : elem_bdry_container_1)
391 40 : boundary_info.add_side(elem_Tri3, 0, bdry_id);
392 284 : for (auto bdry_id : elem_bdry_container_2)
393 68 : boundary_info.add_side(elem_Tri3, 1, bdry_id);
394 : // Assign subdomain id for the TRI element by shifting its original subdomain id
395 216 : elem_Tri3->subdomain_id() = elem_block_id + tri_subdomain_id_shift;
396 216 : new_subdomain_ids.emplace(elem_block_id + tri_subdomain_id_shift);
397 : // Retain element extra integers
398 216 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
399 0 : elem_Tri3->set_extra_integer(j, exist_extra_ids[j]);
400 216 : }
401 : // Assign subdomain names for the new TRI elements
402 116 : for (auto & nid : new_subdomain_ids)
403 : {
404 92 : const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
405 92 : if (MooseMeshUtils::getSubdomainID(
406 184 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
407 184 : : old_name) +
408 368 : "_" + tri_elem_subdomain_name_suffix,
409 92 : mesh) != Moose::INVALID_BLOCK_ID)
410 4 : throw MooseException("The new subdomain name already exists in the mesh.");
411 88 : mesh.subdomain_name(nid) =
412 176 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
413 176 : : old_name) +
414 264 : "_" + tri_elem_subdomain_name_suffix;
415 88 : mooseWarning("Degenerate QUAD elements have been converted into TRI elements with a new "
416 88 : "subdomain name: " +
417 176 : mesh.subdomain_name(nid) + ".");
418 92 : }
419 48 : return bad_elems_rec.size();
420 36 : }
421 :
422 : std::vector<std::pair<Real, unsigned int>>
423 2656 : vertex_angles(const Elem & elem)
424 : {
425 2656 : std::vector<std::pair<Real, unsigned int>> angles;
426 2656 : const unsigned int n_vertices = elem.n_vertices();
427 :
428 13280 : for (unsigned int i = 0; i < n_vertices; i++)
429 : {
430 10624 : Point v1 = (*elem.node_ptr((i - 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
431 10624 : Point v2 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
432 10624 : Real tmp = v1 * v2 / v1.norm() / v2.norm();
433 10624 : if (tmp > 1.0)
434 0 : tmp = 1.0;
435 10624 : else if (tmp < -1.0)
436 44 : tmp = -1.0;
437 10624 : angles.push_back(std::make_pair(acos(tmp), i));
438 : }
439 2656 : std::sort(angles.begin(), angles.end(), std::greater<>());
440 2656 : return angles;
441 0 : }
442 :
443 : std::vector<std::pair<Real, unsigned int>>
444 2656 : vertex_distances(const Elem & elem)
445 : {
446 2656 : std::vector<std::pair<Real, unsigned int>> distances;
447 2656 : const unsigned int n_vertices = elem.n_vertices();
448 :
449 13280 : for (unsigned int i = 0; i < n_vertices; i++)
450 : {
451 10624 : Point v1 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
452 10624 : distances.push_back(std::make_pair(v1.norm(), i));
453 : }
454 2656 : std::sort(distances.begin(), distances.end());
455 2656 : return distances;
456 0 : }
457 :
458 : void
459 800 : triElemSplitter(ReplicatedMesh & mesh,
460 : const dof_id_type elem_id,
461 : const unsigned short node_shift,
462 : const dof_id_type nid_3,
463 : const dof_id_type nid_4,
464 : const subdomain_id_type single_elem_side_id,
465 : const subdomain_id_type double_elem_side_id)
466 : {
467 800 : const auto elem_old = mesh.elem_ptr(elem_id);
468 800 : const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
469 800 : const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
470 800 : const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
471 :
472 : const bool m1_side_flag =
473 800 : MooseUtils::absoluteFuzzyEqual((*(mesh.node_ptr(nid_3)) - *(mesh.node_ptr(nid_0)))
474 800 : .cross(*(mesh.node_ptr(nid_1)) - *(mesh.node_ptr(nid_0)))
475 800 : .norm(),
476 1600 : 0.0);
477 800 : const dof_id_type nid_m1 = m1_side_flag ? nid_3 : nid_4;
478 800 : const dof_id_type nid_m2 = m1_side_flag ? nid_4 : nid_3;
479 : // Build boundary information of the mesh
480 800 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
481 800 : auto bdry_side_list = boundary_info.build_side_list();
482 : // Create a list of sidesets involving the element to be split
483 800 : std::vector<std::vector<boundary_id_type>> elem_side_list;
484 800 : elem_side_list.resize(3);
485 115680 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
486 : {
487 114880 : if (std::get<0>(bdry_side_list[i]) == elem_id)
488 : {
489 352 : elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
490 176 : std::get<2>(bdry_side_list[i]));
491 : }
492 : }
493 :
494 800 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
495 800 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
496 : // Record all the element extra integers of the original element
497 800 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
498 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
499 :
500 800 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
501 800 : elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
502 800 : elem_Tri3_0->set_node(1, mesh.node_ptr(nid_m1));
503 800 : elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m2));
504 800 : elem_Tri3_0->subdomain_id() = single_elem_side_id;
505 800 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
506 800 : elem_Tri3_1->set_node(0, mesh.node_ptr(nid_1));
507 800 : elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m2));
508 800 : elem_Tri3_1->set_node(2, mesh.node_ptr(nid_m1));
509 800 : elem_Tri3_1->subdomain_id() = double_elem_side_id;
510 800 : Elem * elem_Tri3_2 = mesh.add_elem(new Tri3);
511 800 : elem_Tri3_2->set_node(0, mesh.node_ptr(nid_2));
512 800 : elem_Tri3_2->set_node(1, mesh.node_ptr(nid_m2));
513 800 : elem_Tri3_2->set_node(2, mesh.node_ptr(nid_1));
514 800 : elem_Tri3_2->subdomain_id() = double_elem_side_id;
515 : // Retain element extra integers
516 800 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
517 : {
518 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
519 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
520 0 : elem_Tri3_2->set_extra_integer(j, exist_extra_ids[j]);
521 : }
522 :
523 : // Add sideset information to the new elements
524 848 : for (const auto & side_info_0 : elem_side_list[0])
525 : {
526 48 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
527 48 : boundary_info.add_side(elem_Tri3_1, 2, side_info_0);
528 : }
529 848 : for (const auto & side_info_1 : elem_side_list[1])
530 48 : boundary_info.add_side(elem_Tri3_2, 2, side_info_1);
531 880 : for (const auto & side_info_2 : elem_side_list[2])
532 : {
533 80 : boundary_info.add_side(elem_Tri3_0, 2, side_info_2);
534 80 : boundary_info.add_side(elem_Tri3_2, 0, side_info_2);
535 : }
536 800 : }
537 :
538 : void
539 96 : triElemSplitter(ReplicatedMesh & mesh,
540 : const dof_id_type elem_id,
541 : const unsigned short node_shift,
542 : const dof_id_type nid_m,
543 : const subdomain_id_type first_elem_side_id,
544 : const subdomain_id_type second_elem_side_id)
545 : {
546 96 : const auto elem_old = mesh.elem_ptr(elem_id);
547 96 : const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
548 96 : const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
549 96 : const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
550 : // Build boundary information of the mesh
551 96 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
552 96 : auto bdry_side_list = boundary_info.build_side_list();
553 : // Create a list of sidesets involving the element to be split
554 96 : std::vector<std::vector<boundary_id_type>> elem_side_list;
555 96 : elem_side_list.resize(3);
556 10560 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
557 : {
558 10464 : if (std::get<0>(bdry_side_list[i]) == elem_id)
559 : {
560 96 : elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
561 48 : std::get<2>(bdry_side_list[i]));
562 : }
563 : }
564 :
565 96 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
566 96 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
567 : // Record all the element extra integers of the original element
568 96 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
569 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
570 :
571 96 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
572 96 : elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
573 96 : elem_Tri3_0->set_node(1, mesh.node_ptr(nid_1));
574 96 : elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m));
575 96 : elem_Tri3_0->subdomain_id() = first_elem_side_id;
576 96 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
577 96 : elem_Tri3_1->set_node(0, mesh.node_ptr(nid_0));
578 96 : elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m));
579 96 : elem_Tri3_1->set_node(2, mesh.node_ptr(nid_2));
580 96 : elem_Tri3_1->subdomain_id() = second_elem_side_id;
581 : // Retain element extra integers
582 96 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
583 : {
584 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
585 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
586 : }
587 :
588 : // Add sideset information to the new elements
589 112 : for (const auto & side_info_0 : elem_side_list[0])
590 16 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
591 96 : for (const auto & side_info_1 : elem_side_list[1])
592 : {
593 0 : boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
594 0 : boundary_info.add_side(elem_Tri3_1, 1, side_info_1);
595 : }
596 128 : for (const auto & side_info_2 : elem_side_list[2])
597 32 : boundary_info.add_side(elem_Tri3_1, 2, side_info_2);
598 96 : }
599 :
600 : void
601 816 : quadElemSplitter(ReplicatedMesh & mesh,
602 : const dof_id_type elem_id,
603 : const subdomain_id_type tri_elem_subdomain_shift)
604 : {
605 : // Build boundary information of the mesh
606 816 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
607 816 : auto bdry_side_list = boundary_info.build_side_list();
608 : // Create a list of sidesets involving the element to be split
609 816 : std::vector<std::vector<boundary_id_type>> elem_side_list;
610 816 : elem_side_list.resize(4);
611 90032 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
612 : {
613 89216 : if (std::get<0>(bdry_side_list[i]) == elem_id)
614 : {
615 460 : elem_side_list[std::get<1>(bdry_side_list[i])].push_back(std::get<2>(bdry_side_list[i]));
616 : }
617 : }
618 :
619 816 : auto node_0 = mesh.elem_ptr(elem_id)->node_ptr(0);
620 816 : auto node_1 = mesh.elem_ptr(elem_id)->node_ptr(1);
621 816 : auto node_2 = mesh.elem_ptr(elem_id)->node_ptr(2);
622 816 : auto node_3 = mesh.elem_ptr(elem_id)->node_ptr(3);
623 :
624 816 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
625 816 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
626 : // Record all the element extra integers of the original quad element
627 816 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
628 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
629 :
630 : // There are two trivial ways to split a quad element
631 : // We prefer the way that leads to triangles with similar areas
632 816 : if (std::abs((*node_1 - *node_0).cross(*node_3 - *node_0).norm() -
633 816 : (*node_1 - *node_2).cross(*node_3 - *node_2).norm()) >
634 816 : std::abs((*node_0 - *node_1).cross(*node_2 - *node_1).norm() -
635 816 : (*node_0 - *node_3).cross(*node_2 - *node_3).norm()))
636 : {
637 232 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
638 232 : elem_Tri3_0->set_node(0, node_0);
639 232 : elem_Tri3_0->set_node(1, node_1);
640 232 : elem_Tri3_0->set_node(2, node_2);
641 232 : elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
642 232 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
643 232 : elem_Tri3_1->set_node(0, node_0);
644 232 : elem_Tri3_1->set_node(1, node_2);
645 232 : elem_Tri3_1->set_node(2, node_3);
646 232 : elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
647 : // Retain element extra integers
648 232 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
649 : {
650 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
651 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
652 : }
653 :
654 : // Add sideset information to the new elements
655 232 : for (const auto & side_info_0 : elem_side_list[0])
656 0 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
657 348 : for (const auto & side_info_1 : elem_side_list[1])
658 116 : boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
659 232 : for (const auto & side_info_2 : elem_side_list[2])
660 0 : boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
661 232 : for (const auto & side_info_3 : elem_side_list[3])
662 0 : boundary_info.add_side(elem_Tri3_1, 2, side_info_3);
663 : }
664 : else
665 : {
666 584 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
667 584 : elem_Tri3_0->set_node(0, node_0);
668 584 : elem_Tri3_0->set_node(1, node_1);
669 584 : elem_Tri3_0->set_node(2, node_3);
670 584 : elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
671 584 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
672 584 : elem_Tri3_1->set_node(0, node_1);
673 584 : elem_Tri3_1->set_node(1, node_2);
674 584 : elem_Tri3_1->set_node(2, node_3);
675 584 : elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
676 : // Retain element extra integers
677 584 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
678 : {
679 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
680 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
681 : }
682 :
683 : // Add sideset information to the new elements
684 616 : for (const auto & side_info_0 : elem_side_list[0])
685 32 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
686 760 : for (const auto & side_info_1 : elem_side_list[1])
687 176 : boundary_info.add_side(elem_Tri3_1, 0, side_info_1);
688 640 : for (const auto & side_info_2 : elem_side_list[2])
689 56 : boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
690 664 : for (const auto & side_info_3 : elem_side_list[3])
691 80 : boundary_info.add_side(elem_Tri3_0, 2, side_info_3);
692 : }
693 816 : }
694 :
695 : void
696 44 : quadToTriOnLine(ReplicatedMesh & mesh,
697 : const std::vector<Real> & cut_line_params,
698 : const dof_id_type tri_subdomain_id_shift,
699 : const SubdomainName tri_elem_subdomain_name_suffix)
700 : {
701 : // Preprocess: find all the quad elements that are across the cutting line
702 44 : std::vector<dof_id_type> cross_elems_quad;
703 44 : std::set<subdomain_id_type> new_subdomain_ids;
704 7148 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
705 : elem_it++)
706 : {
707 7104 : if ((*elem_it)->n_vertices() == 4)
708 : {
709 7104 : std::vector<unsigned short> node_side_rec;
710 35520 : for (unsigned int i = 0; i < 4; i++)
711 : {
712 28416 : const Point v_point = (*elem_it)->point(i);
713 28416 : node_side_rec.push_back(lineSideDeterminator(v_point(0),
714 28416 : v_point(1),
715 28416 : cut_line_params[0],
716 28416 : cut_line_params[1],
717 28416 : cut_line_params[2],
718 : true));
719 : }
720 10304 : if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) != 4 &&
721 3200 : std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) > 0)
722 : {
723 816 : cross_elems_quad.push_back((*elem_it)->id());
724 816 : new_subdomain_ids.emplace((*elem_it)->subdomain_id() + tri_subdomain_id_shift);
725 : }
726 7104 : }
727 44 : }
728 : // Then convert these quad elements into tri elements
729 860 : for (const auto & cross_elem_quad : cross_elems_quad)
730 : {
731 816 : quadElemSplitter(mesh, cross_elem_quad, tri_subdomain_id_shift);
732 816 : mesh.delete_elem(mesh.elem_ptr(cross_elem_quad));
733 : }
734 148 : for (auto & nid : new_subdomain_ids)
735 : {
736 108 : const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
737 108 : if (MooseMeshUtils::getSubdomainID(
738 216 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
739 216 : : old_name) +
740 432 : "_" + tri_elem_subdomain_name_suffix,
741 108 : mesh) != Moose::INVALID_BLOCK_ID)
742 4 : throw MooseException("The new subdomain name already exists in the mesh.");
743 104 : mesh.subdomain_name(nid) =
744 208 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
745 208 : : old_name) +
746 312 : "_" + tri_elem_subdomain_name_suffix;
747 104 : mooseWarning("QUAD elements have been converted into TRI elements with a new "
748 104 : "subdomain name: " +
749 208 : mesh.subdomain_name(nid) + ".");
750 108 : }
751 40 : mesh.contract();
752 48 : }
753 :
754 : void
755 40 : lineRemoverCutElemTri(ReplicatedMesh & mesh,
756 : const std::vector<Real> & cut_line_params,
757 : const subdomain_id_type block_id_to_remove,
758 : const boundary_id_type new_boundary_id)
759 : {
760 : // Find all the elements that are across the cutting line
761 40 : std::vector<dof_id_type> cross_elems;
762 : // A vector for element specific information
763 40 : std::vector<std::vector<std::pair<dof_id_type, dof_id_type>>> node_pairs_vec;
764 : // A set for unique pairs
765 40 : std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs_unique_vec;
766 6480 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
767 : elem_it++)
768 : {
769 6440 : std::vector<unsigned short> node_side_rec;
770 6440 : const auto n_vertices = (*elem_it)->n_vertices();
771 6440 : node_side_rec.resize(n_vertices);
772 30840 : for (unsigned int i = 0; i < n_vertices; i++)
773 : {
774 : // First check if the vertex is in the XY Plane
775 24400 : if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
776 0 : mooseError("MooseMeshXYCuttingUtils::lineRemoverCutElemTri() only works for 2D meshes in "
777 : "XY Plane.");
778 24400 : const Point v_point = (*elem_it)->point(i);
779 24400 : node_side_rec[i] = lineSideDeterminator(
780 24400 : v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2], true);
781 : }
782 6440 : if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) == (int)node_side_rec.size())
783 : {
784 3424 : (*elem_it)->subdomain_id() = block_id_to_remove;
785 : }
786 3016 : else if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) > 0)
787 : {
788 1128 : if ((*elem_it)->n_vertices() != 3 || (*elem_it)->n_nodes() != 3)
789 0 : mooseError("The element across the cutting line is not TRI3, which is not supported.");
790 1128 : cross_elems.push_back((*elem_it)->id());
791 : // Then we need to check pairs of nodes that are on the different side
792 1128 : std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs;
793 4512 : for (unsigned int i = 0; i < node_side_rec.size(); i++)
794 : {
795 : // first node on removal side and second node on retaining side
796 3384 : if (node_side_rec[i] > 0 && node_side_rec[(i + 1) % node_side_rec.size()] == 0)
797 : {
798 : // Removal side first
799 1128 : node_pairs.push_back(
800 1128 : std::make_pair((*elem_it)->node_ptr(i)->id(),
801 1128 : (*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id()));
802 1128 : node_pairs_unique_vec.push_back(node_pairs.back());
803 : }
804 : // first node on retaining side and second node on removal side
805 2256 : else if (node_side_rec[i] == 0 && node_side_rec[(i + 1) % node_side_rec.size()] > 0)
806 : {
807 : // Removal side first
808 1128 : node_pairs.push_back(
809 1128 : std::make_pair((*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id(),
810 1128 : (*elem_it)->node_ptr(i)->id()));
811 1128 : node_pairs_unique_vec.push_back(node_pairs.back());
812 : }
813 : }
814 1128 : node_pairs_vec.push_back(node_pairs);
815 1128 : }
816 6480 : }
817 40 : auto vec_ip = std::unique(node_pairs_unique_vec.begin(), node_pairs_unique_vec.end());
818 40 : node_pairs_unique_vec.resize(std::distance(node_pairs_unique_vec.begin(), vec_ip));
819 :
820 : // Loop over all the node pairs to define new nodes that sit on the cutting line
821 40 : std::vector<Node *> nodes_on_line;
822 : // whether the on-line node is overlapped with the node pairs or a brand new node
823 40 : std::vector<unsigned short> nodes_on_line_overlap;
824 2136 : for (const auto & node_pair : node_pairs_unique_vec)
825 : {
826 2096 : const Point pt1 = *mesh.node_ptr(node_pair.first);
827 2096 : const Point pt2 = *mesh.node_ptr(node_pair.second);
828 2096 : const Point pt_line = twoPointandLineIntersection(
829 2096 : pt1, pt2, cut_line_params[0], cut_line_params[1], cut_line_params[2]);
830 2096 : if ((pt_line - pt1).norm() < libMesh::TOLERANCE)
831 : {
832 0 : nodes_on_line.push_back(mesh.node_ptr(node_pair.first));
833 0 : nodes_on_line_overlap.push_back(1);
834 : }
835 2096 : else if ((pt_line - pt2).norm() < libMesh::TOLERANCE)
836 : {
837 544 : nodes_on_line.push_back(mesh.node_ptr(node_pair.second));
838 544 : nodes_on_line_overlap.push_back(2);
839 : }
840 : else
841 : {
842 1552 : nodes_on_line.push_back(mesh.add_point(pt_line));
843 1552 : nodes_on_line_overlap.push_back(0);
844 : }
845 : }
846 :
847 : // make new elements
848 1168 : for (unsigned int i = 0; i < cross_elems.size(); i++)
849 : {
850 : // Only TRI elements are involved after preprocessing
851 1128 : auto cross_elem = mesh.elem_ptr(cross_elems[i]);
852 1128 : auto node_0 = cross_elem->node_ptr(0);
853 1128 : auto node_1 = cross_elem->node_ptr(1);
854 1128 : auto node_2 = cross_elem->node_ptr(2);
855 1128 : const std::vector<dof_id_type> tri_nodes = {node_0->id(), node_1->id(), node_2->id()};
856 :
857 1128 : const auto online_node_index_1 = std::distance(node_pairs_unique_vec.begin(),
858 : std::find(node_pairs_unique_vec.begin(),
859 : node_pairs_unique_vec.end(),
860 1128 : node_pairs_vec[i][0]));
861 1128 : const auto online_node_index_2 = std::distance(node_pairs_unique_vec.begin(),
862 : std::find(node_pairs_unique_vec.begin(),
863 : node_pairs_unique_vec.end(),
864 1128 : node_pairs_vec[i][1]));
865 1128 : auto node_3 = nodes_on_line[online_node_index_1];
866 1128 : auto node_4 = nodes_on_line[online_node_index_2];
867 1128 : const auto node_3_overlap_flag = nodes_on_line_overlap[online_node_index_1];
868 1128 : const auto node_4_overlap_flag = nodes_on_line_overlap[online_node_index_2];
869 : // Most common case, no overlapped nodes
870 1128 : if (node_3_overlap_flag == 0 && node_4_overlap_flag == 0)
871 : {
872 : // True if the common node is on the removal side; false if on the retaining side
873 800 : const bool common_node_side = node_pairs_vec[i][0].first == node_pairs_vec[i][1].first;
874 : const subdomain_id_type block_id_to_assign_1 =
875 800 : common_node_side ? block_id_to_remove : cross_elem->subdomain_id();
876 : const subdomain_id_type block_id_to_assign_2 =
877 800 : common_node_side ? cross_elem->subdomain_id() : block_id_to_remove;
878 : // The reference node ids need to be adjusted according to the common node of the two cut
879 : // sides
880 : const dof_id_type common_node_id =
881 800 : common_node_side ? node_pairs_vec[i][0].first : node_pairs_vec[i][0].second;
882 :
883 1600 : triElemSplitter(mesh,
884 : cross_elem->id(),
885 800 : std::distance(tri_nodes.begin(),
886 : std::find(tri_nodes.begin(), tri_nodes.end(), common_node_id)),
887 : node_3->id(),
888 : node_4->id(),
889 : block_id_to_assign_1,
890 : block_id_to_assign_2);
891 800 : mesh.delete_elem(cross_elem);
892 800 : }
893 : // both node_3 and node_4 are overlapped
894 328 : else if (node_3_overlap_flag > 0 && node_4_overlap_flag > 0)
895 : {
896 : // In this case, the entire element is on one side of the cutting line
897 : // No change needed just check which side the element is on
898 464 : cross_elem->subdomain_id() = lineSideDeterminator(cross_elem->vertex_average()(0),
899 464 : cross_elem->vertex_average()(1),
900 232 : cut_line_params[0],
901 232 : cut_line_params[1],
902 232 : cut_line_params[2],
903 : true)
904 232 : ? block_id_to_remove
905 0 : : cross_elem->subdomain_id();
906 : }
907 : // node_3 or node_4 is overlapped
908 : else
909 : {
910 96 : const auto node_3_finder = std::distance(
911 96 : tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_3->id()));
912 96 : const auto node_4_finder = std::distance(
913 96 : tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_4->id()));
914 : // As only one of the two above values should be less than the three, the smaller one should
915 : // be used
916 96 : const dof_id_type node_id = node_3_finder < node_4_finder ? node_4->id() : node_3->id();
917 96 : const auto node_finder = std::min(node_3_finder, node_4_finder);
918 :
919 144 : triElemSplitter(
920 : mesh,
921 : cross_elem->id(),
922 : node_finder,
923 : node_id,
924 96 : tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
925 48 : ? block_id_to_remove
926 48 : : cross_elem->subdomain_id(),
927 96 : tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
928 48 : ? cross_elem->subdomain_id()
929 : : block_id_to_remove);
930 96 : mesh.delete_elem(cross_elem);
931 : }
932 1128 : }
933 40 : mesh.contract();
934 : // Due to the complexity, we identify the new boundary here together instead of during cutting of
935 : // each element, because the preexisting element edges that are aligned with the cutting line also
936 : // need to be added to the new boundary.
937 40 : mesh.find_neighbors();
938 40 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
939 8176 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
940 : elem_it++)
941 : {
942 8136 : if ((*elem_it)->subdomain_id() != block_id_to_remove)
943 : {
944 14448 : for (unsigned int j = 0; j < (*elem_it)->n_sides(); j++)
945 : {
946 11280 : if ((*elem_it)->neighbor_ptr(j) != nullptr)
947 10704 : if ((*elem_it)->neighbor_ptr(j)->subdomain_id() == block_id_to_remove)
948 960 : boundary_info.add_side(*elem_it, j, new_boundary_id);
949 : }
950 : }
951 40 : }
952 :
953 : // Delete the block to remove
954 5008 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
955 5008 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
956 : elem_it++)
957 5008 : mesh.delete_elem(*elem_it);
958 40 : mesh.contract();
959 40 : }
960 :
961 : void
962 44 : lineRemoverCutElem(ReplicatedMesh & mesh,
963 : const std::vector<Real> & cut_line_params,
964 : const dof_id_type tri_subdomain_id_shift,
965 : const SubdomainName tri_elem_subdomain_name_suffix,
966 : const subdomain_id_type block_id_to_remove,
967 : const boundary_id_type new_boundary_id,
968 : const bool improve_boundary_tri_elems)
969 : {
970 : // Convert any quad elements crossed by the line into tri elements
971 44 : quadToTriOnLine(mesh, cut_line_params, tri_subdomain_id_shift, tri_elem_subdomain_name_suffix);
972 : // Then do the cutting for the preprocessed mesh that only contains tri elements crossed by the
973 : // cut line
974 40 : lineRemoverCutElemTri(mesh, cut_line_params, block_id_to_remove, new_boundary_id);
975 :
976 40 : if (improve_boundary_tri_elems)
977 8 : boundaryTriElemImprover(mesh, new_boundary_id);
978 40 : }
979 :
980 : void
981 8 : boundaryTriElemImprover(ReplicatedMesh & mesh, const boundary_id_type boundary_to_improve)
982 : {
983 8 : if (!MooseMeshUtils::hasBoundaryID(mesh, boundary_to_improve))
984 0 : mooseError(
985 : "MooseMeshXYCuttingUtils::boundaryTriElemImprover(): The boundary_to_improve provided "
986 : "does not exist in the given mesh.");
987 8 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
988 8 : auto side_list = boundary_info.build_side_list();
989 : // Here we would like to collect the following information for all the TRI3 elements on the
990 : // boundary: Key: node id of the off-boundary node Value: a vector of tuples, each tuple contains
991 : // the following information:
992 : // 1. The element id of the element that is on the boundary to improve
993 : // 2. the one node id of that element that is on the boundary to improve
994 : // 3. the other node id of the element that is on the boundary to improve
995 : std::map<dof_id_type, std::vector<std::tuple<dof_id_type, dof_id_type, dof_id_type>>>
996 8 : tri3_elem_info;
997 904 : for (const auto & side : side_list)
998 : {
999 896 : if (std::get<2>(side) == boundary_to_improve)
1000 : {
1001 416 : Elem * elem = mesh.elem_ptr(std::get<0>(side));
1002 416 : if (elem->type() == TRI3)
1003 : {
1004 416 : const auto key_node_id = elem->node_id((std::get<1>(side) + 2) % 3);
1005 416 : const auto value_elem_id = elem->id();
1006 416 : const auto value_node_id_1 = elem->node_id(std::get<1>(side));
1007 416 : const auto value_node_id_2 = elem->node_id((std::get<1>(side) + 1) % 3);
1008 832 : tri3_elem_info[key_node_id].push_back(
1009 832 : std::make_tuple(value_elem_id, value_node_id_1, value_node_id_2));
1010 : }
1011 : }
1012 : }
1013 : // Elements that need to be removed
1014 8 : std::vector<dof_id_type> elems_to_remove;
1015 : // Now check if any group of TRI3 sharing an off-boundary node can be improved.
1016 216 : for (const auto & tri_group : tri3_elem_info)
1017 : {
1018 : // It is possible to improve only when more than one TRI3 elements share the same off-boundary
1019 : // node
1020 208 : std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
1021 208 : std::vector<dof_id_type> elem_id_list;
1022 624 : for (const auto & tri : tri_group.second)
1023 : {
1024 416 : node_assm.push_back(std::make_pair(std::get<1>(tri), std::get<2>(tri)));
1025 416 : elem_id_list.push_back(std::get<0>(tri));
1026 : }
1027 208 : std::vector<dof_id_type> ordered_node_list;
1028 208 : std::vector<dof_id_type> ordered_elem_list;
1029 208 : MooseMeshUtils::makeOrderedNodeList(
1030 : node_assm, elem_id_list, ordered_node_list, ordered_elem_list);
1031 :
1032 : // For all the elements sharing the same off-boundary node, we need to know how many separated
1033 : // subdomains are involved
1034 : // If there are extra element ids defined on the mesh, they also want to retain their boundaries
1035 : // Only triangle elements that share a side can be merged
1036 208 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
1037 208 : std::vector<std::tuple<subdomain_id_type, std::vector<dof_id_type>, unsigned int>> blocks_info;
1038 624 : for (const auto & elem_id : ordered_elem_list)
1039 : {
1040 416 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
1041 : // Record all the element extra integers of the original quad element
1042 416 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
1043 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
1044 416 : if (!blocks_info.empty())
1045 : {
1046 400 : if (mesh.elem_ptr(elem_id)->subdomain_id() == std::get<0>(blocks_info.back()) &&
1047 192 : exist_extra_ids == std::get<1>(blocks_info.back()))
1048 : {
1049 192 : std::get<2>(blocks_info.back())++;
1050 192 : continue;
1051 : }
1052 : }
1053 224 : blocks_info.push_back(
1054 448 : std::make_tuple(mesh.elem_ptr(elem_id)->subdomain_id(), exist_extra_ids, 1));
1055 416 : }
1056 : // For each separated subdomain / set of extra ids, we try to improve the boundary elements
1057 208 : unsigned int side_counter = 0;
1058 432 : for (const auto & block_info : blocks_info)
1059 : {
1060 224 : const auto node_1 = mesh.node_ptr(ordered_node_list[side_counter]);
1061 : // we do not need to subtract 1 for node_2
1062 224 : const auto node_2 = mesh.node_ptr(ordered_node_list[side_counter + std::get<2>(block_info)]);
1063 224 : const auto node_0 = mesh.node_ptr(tri_group.first);
1064 224 : const Point v1 = *node_1 - *node_0;
1065 224 : const Point v2 = *node_2 - *node_0;
1066 224 : const Real angle = std::acos(v1 * v2 / v1.norm() / v2.norm()) / M_PI * 180.0;
1067 0 : const std::vector<dof_id_type> block_elems(ordered_elem_list.begin() + side_counter,
1068 224 : ordered_elem_list.begin() + side_counter +
1069 448 : std::get<2>(block_info));
1070 : // We assume that there are no sidesets defined inside a subdomain
1071 : // For the first TRI3 element, we want to check if its side defined by node_0 and node_1 is
1072 : // defined in any sidesets
1073 : unsigned short side_id_0;
1074 : unsigned short side_id_t;
1075 : bool is_inverse_0;
1076 : bool is_inverse_t;
1077 224 : elemSideLocator(mesh,
1078 224 : block_elems.front(),
1079 224 : tri_group.first,
1080 224 : ordered_node_list[side_counter],
1081 : side_id_0,
1082 : is_inverse_0);
1083 224 : elemSideLocator(mesh,
1084 224 : block_elems.back(),
1085 224 : ordered_node_list[side_counter + std::get<2>(block_info)],
1086 224 : tri_group.first,
1087 : side_id_t,
1088 : is_inverse_t);
1089 : // Collect boundary information of the identified sides
1090 224 : std::vector<boundary_id_type> side_0_boundary_ids;
1091 448 : boundary_info.boundary_ids(
1092 224 : mesh.elem_ptr(block_elems.front()), side_id_0, side_0_boundary_ids);
1093 224 : std::vector<boundary_id_type> side_t_boundary_ids;
1094 224 : boundary_info.boundary_ids(mesh.elem_ptr(block_elems.back()), side_id_t, side_t_boundary_ids);
1095 :
1096 : // Ideally we want this angle to be 60 degrees
1097 : // In reality, we want one TRI3 element if the angle is less than 90 degrees;
1098 : // we want two TRI3 elements if the angle is greater than 90 degrees and less than 135
1099 : // degrees; we want three TRI3 elements if the angle is greater than 135 degrees and less than
1100 : // 180 degrees.
1101 224 : if (angle < 90.0)
1102 : {
1103 152 : if (std::get<2>(block_info) > 1)
1104 : {
1105 128 : makeImprovedTriElement(mesh,
1106 64 : tri_group.first,
1107 64 : ordered_node_list[side_counter],
1108 64 : ordered_node_list[side_counter + std::get<2>(block_info)],
1109 64 : std::get<0>(block_info),
1110 64 : std::get<1>(block_info),
1111 : {boundary_to_improve},
1112 : side_0_boundary_ids,
1113 : side_t_boundary_ids);
1114 64 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1115 : }
1116 : }
1117 72 : else if (angle < 135.0)
1118 : {
1119 : // We can just add the middle node because there's nothing on the other side
1120 56 : const auto node_m = mesh.add_point((*node_1 + *node_2) / 2.0);
1121 112 : makeImprovedTriElement(mesh,
1122 56 : tri_group.first,
1123 56 : ordered_node_list[side_counter],
1124 : node_m->id(),
1125 56 : std::get<0>(block_info),
1126 56 : std::get<1>(block_info),
1127 : {boundary_to_improve},
1128 : side_0_boundary_ids,
1129 112 : std::vector<boundary_id_type>());
1130 168 : makeImprovedTriElement(mesh,
1131 56 : tri_group.first,
1132 : node_m->id(),
1133 56 : ordered_node_list[side_counter + std::get<2>(block_info)],
1134 56 : std::get<0>(block_info),
1135 56 : std::get<1>(block_info),
1136 : {boundary_to_improve},
1137 112 : std::vector<boundary_id_type>(),
1138 : side_t_boundary_ids);
1139 56 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1140 : }
1141 : else
1142 : {
1143 16 : const auto node_m1 = mesh.add_point((*node_1 * 2.0 + *node_2) / 3.0);
1144 16 : const auto node_m2 = mesh.add_point((*node_1 + *node_2 * 2.0) / 3.0);
1145 32 : makeImprovedTriElement(mesh,
1146 16 : tri_group.first,
1147 16 : ordered_node_list[side_counter],
1148 : node_m1->id(),
1149 16 : std::get<0>(block_info),
1150 16 : std::get<1>(block_info),
1151 : {boundary_to_improve},
1152 : side_0_boundary_ids,
1153 32 : std::vector<boundary_id_type>());
1154 48 : makeImprovedTriElement(mesh,
1155 16 : tri_group.first,
1156 : node_m1->id(),
1157 : node_m2->id(),
1158 16 : std::get<0>(block_info),
1159 16 : std::get<1>(block_info),
1160 : {boundary_to_improve},
1161 32 : std::vector<boundary_id_type>(),
1162 32 : std::vector<boundary_id_type>());
1163 48 : makeImprovedTriElement(mesh,
1164 16 : tri_group.first,
1165 : node_m2->id(),
1166 16 : ordered_node_list[side_counter + std::get<2>(block_info)],
1167 16 : std::get<0>(block_info),
1168 16 : std::get<1>(block_info),
1169 : {boundary_to_improve},
1170 32 : std::vector<boundary_id_type>(),
1171 : side_t_boundary_ids);
1172 16 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1173 : }
1174 224 : side_counter += std::get<2>(block_info);
1175 224 : }
1176 : // TODO: Need to check if the new element is inverted?
1177 208 : }
1178 : // Delete the original elements
1179 336 : for (const auto & elem_to_remove : elems_to_remove)
1180 328 : mesh.delete_elem(mesh.elem_ptr(elem_to_remove));
1181 8 : mesh.contract();
1182 8 : }
1183 :
1184 : void
1185 224 : makeImprovedTriElement(ReplicatedMesh & mesh,
1186 : const dof_id_type node_id_0,
1187 : const dof_id_type node_id_1,
1188 : const dof_id_type node_id_2,
1189 : const subdomain_id_type subdomain_id,
1190 : const std::vector<dof_id_type> & extra_elem_ids,
1191 : const std::vector<boundary_id_type> & boundary_ids_for_side_1,
1192 : const std::vector<boundary_id_type> & boundary_ids_for_side_0,
1193 : const std::vector<boundary_id_type> & boundary_ids_for_side_2)
1194 : {
1195 224 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1196 224 : Elem * elem_Tri3_new = mesh.add_elem(new Tri3);
1197 224 : elem_Tri3_new->set_node(0, mesh.node_ptr(node_id_0));
1198 224 : elem_Tri3_new->set_node(1, mesh.node_ptr(node_id_1));
1199 224 : elem_Tri3_new->set_node(2, mesh.node_ptr(node_id_2));
1200 264 : for (const auto & boundary_id_for_side_0 : boundary_ids_for_side_0)
1201 40 : boundary_info.add_side(elem_Tri3_new, 0, boundary_id_for_side_0);
1202 448 : for (const auto & boundary_id_for_side_1 : boundary_ids_for_side_1)
1203 224 : boundary_info.add_side(elem_Tri3_new, 1, boundary_id_for_side_1);
1204 224 : for (const auto & boundary_id_for_side_2 : boundary_ids_for_side_2)
1205 0 : boundary_info.add_side(elem_Tri3_new, 2, boundary_id_for_side_2);
1206 224 : elem_Tri3_new->subdomain_id() = subdomain_id;
1207 : // Retain element extra integers
1208 224 : for (unsigned int j = 0; j < extra_elem_ids.size(); j++)
1209 : {
1210 0 : elem_Tri3_new->set_extra_integer(j, extra_elem_ids[j]);
1211 : }
1212 224 : }
1213 :
1214 : bool
1215 448 : elemSideLocator(ReplicatedMesh & mesh,
1216 : const dof_id_type elem_id,
1217 : const dof_id_type node_id_0,
1218 : const dof_id_type node_id_1,
1219 : unsigned short & side_id,
1220 : bool & is_inverse)
1221 : {
1222 448 : Elem * elem = mesh.elem_ptr(elem_id);
1223 896 : for (unsigned short i = 0; i < elem->n_sides(); i++)
1224 : {
1225 2088 : if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_0 &&
1226 1192 : elem->side_ptr(i)->node_ptr(1)->id() == node_id_1)
1227 : {
1228 144 : side_id = i;
1229 144 : is_inverse = false;
1230 144 : return true;
1231 : }
1232 1880 : else if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_1 &&
1233 1128 : elem->side_ptr(i)->node_ptr(1)->id() == node_id_0)
1234 : {
1235 304 : side_id = i;
1236 304 : is_inverse = true;
1237 304 : return true;
1238 : }
1239 : }
1240 0 : return false;
1241 : }
1242 : }
|
