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 27 : 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 27 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
39 27 : auto bdry_side_list = boundary_info.build_side_list();
40 : // Only select the boundaries_to_conform
41 27 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> slc_bdry_side_list;
42 2731 : for (const auto i : index_range(bdry_side_list))
43 4476 : if (std::get<2>(bdry_side_list[i]) == external_boundary_id ||
44 1772 : std::find(other_boundaries_to_conform.begin(),
45 : other_boundaries_to_conform.end(),
46 6248 : std::get<2>(bdry_side_list[i])) != other_boundaries_to_conform.end())
47 1700 : 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 27 : std::vector<dof_id_type> crossed_elems_to_remove;
52 6539 : 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 6512 : unsigned short removal_side_count = 0;
57 32560 : for (const auto i : make_range((*elem_it)->n_vertices()))
58 : {
59 : // First check if the vertex is on the XY-Plane
60 26048 : 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 26048 : if (lineSideDeterminator((*elem_it)->point(i)(0),
64 26048 : (*elem_it)->point(i)(1),
65 26048 : bdry_pars[0],
66 26048 : bdry_pars[1],
67 26048 : bdry_pars[2],
68 : side_to_remove))
69 16212 : removal_side_count++;
70 : }
71 6512 : if (removal_side_count == (*elem_it)->n_vertices())
72 : {
73 3698 : (*elem_it)->subdomain_id() = block_id_to_remove;
74 3698 : continue;
75 : }
76 : // Check the average of the vertices of the element
77 2814 : if (lineSideDeterminator((*elem_it)->vertex_average()(0),
78 5628 : (*elem_it)->vertex_average()(1),
79 2814 : bdry_pars[0],
80 2814 : bdry_pars[1],
81 2814 : bdry_pars[2],
82 : side_to_remove))
83 422 : crossed_elems_to_remove.push_back((*elem_it)->id());
84 27 : }
85 : // Check each crossed element to see if removing it would lead to boundary moving
86 449 : for (const auto & elem_id : crossed_elems_to_remove)
87 : {
88 422 : bool remove_flag = true;
89 1758 : for (const auto i : make_range(mesh.elem_ptr(elem_id)->n_sides()))
90 : {
91 1448 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i) != nullptr)
92 2155 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id() != block_id_to_remove &&
93 769 : std::find(crossed_elems_to_remove.begin(),
94 : crossed_elems_to_remove.end(),
95 1538 : mesh.elem_ptr(elem_id)->neighbor_ptr(i)->id()) ==
96 2155 : crossed_elems_to_remove.end())
97 : {
98 552 : if (mesh.elem_ptr(elem_id)->subdomain_id() !=
99 552 : mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id())
100 : {
101 112 : remove_flag = false;
102 112 : break;
103 : }
104 : }
105 : }
106 422 : if (remove_flag)
107 310 : 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 27 : std::vector<dof_id_type> node_list;
114 2531 : for (auto elem_it = mesh.active_subdomain_set_elements_begin(subdomain_ids_set);
115 2531 : elem_it != mesh.active_subdomain_set_elements_end(subdomain_ids_set);
116 : elem_it++)
117 : {
118 12520 : for (const auto i : make_range((*elem_it)->n_sides()))
119 : {
120 10016 : if ((*elem_it)->neighbor_ptr(i) != nullptr)
121 9528 : if ((*elem_it)->neighbor_ptr(i)->subdomain_id() == block_id_to_remove)
122 : {
123 726 : node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(0)->id());
124 726 : node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(1)->id());
125 726 : boundary_info.add_side(*elem_it, i, trimming_section_boundary_id);
126 726 : 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 27 : }
131 : // Remove duplicate nodes
132 27 : const auto unique_it = std::unique(node_list.begin(), node_list.end());
133 54 : 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 27 : std::vector<bool> node_list_flag(node_list.size(), false);
137 27 : std::vector<Point> node_list_point(node_list.size(), Point(0.0, 0.0, 0.0));
138 : // Loop over all the selected sides
139 1727 : for (const auto i : index_range(slc_bdry_side_list))
140 : {
141 : // Get the two node ids of the side
142 1700 : dof_id_type side_id_0 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
143 1700 : ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
144 1700 : ->node_ptr(0)
145 1700 : ->id();
146 1700 : dof_id_type side_id_1 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
147 1700 : ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
148 1700 : ->node_ptr(1)
149 1700 : ->id();
150 : // True means the selected bdry node is in the node list of the trimming interface
151 : bool side_id_0_in =
152 1700 : !(std::find(node_list.begin(), node_list.end(), side_id_0) == node_list.end());
153 : bool side_id_1_in =
154 1700 : !(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 1700 : bool side_node_0_remove = lineSideDeterminator((*mesh.node_ptr(side_id_0))(0),
158 1700 : (*mesh.node_ptr(side_id_0))(1),
159 1700 : bdry_pars[0],
160 1700 : bdry_pars[1],
161 1700 : bdry_pars[2],
162 : side_to_remove);
163 1700 : bool side_node_1_remove = lineSideDeterminator((*mesh.node_ptr(side_id_1))(0),
164 1700 : (*mesh.node_ptr(side_id_1))(1),
165 1700 : bdry_pars[0],
166 1700 : bdry_pars[1],
167 1700 : bdry_pars[2],
168 : side_to_remove);
169 : // If both nodes of that side are involved in the trimming interface
170 1700 : 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 1692 : 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 67 : node_list_flag[std::distance(
181 67 : node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_0))] = true;
182 67 : const Point p0 = *mesh.node_ptr(side_id_0);
183 67 : const Point p1 = *mesh.node_ptr(side_id_1);
184 :
185 67 : node_list_point[std::distance(node_list.begin(),
186 67 : std::find(node_list.begin(), node_list.end(), side_id_0))] =
187 67 : twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
188 67 : }
189 : // If node 1 is on the trimming interface, and the side is cut by the trimming line
190 1625 : 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 1215 : for (const auto i : index_range(node_list))
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 1188 : if (node_list_flag[i])
211 115 : *(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 1073 : const Real x0 = (*(mesh.node_ptr(node_list[i])))(0);
217 1073 : const Real y0 = (*(mesh.node_ptr(node_list[i])))(1);
218 1073 : (*(mesh.node_ptr(node_list[i])))(0) =
219 1073 : (bdry_pars[1] * (bdry_pars[1] * x0 - bdry_pars[0] * y0) - bdry_pars[0] * bdry_pars[2]) /
220 1073 : (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
221 1073 : (*(mesh.node_ptr(node_list[i])))(1) =
222 1073 : (bdry_pars[0] * (-bdry_pars[1] * x0 + bdry_pars[0] * y0) - bdry_pars[1] * bdry_pars[2]) /
223 1073 : (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
224 : }
225 : }
226 :
227 : // Delete the block
228 4035 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
229 4035 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
230 : elem_it++)
231 4035 : mesh.delete_elem(*elem_it);
232 27 : mesh.contract();
233 27 : mesh.find_neighbors();
234 : // Delete zero volume elements
235 27 : std::vector<dof_id_type> zero_elems;
236 2531 : for (auto elem_it = mesh.elements_begin(); elem_it != mesh.elements_end(); elem_it++)
237 : {
238 2504 : if (MooseUtils::absoluteFuzzyEqual((*elem_it)->volume(), 0.0))
239 : {
240 80 : for (const auto i : make_range((*elem_it)->n_sides()))
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 27 : }
255 43 : for (const auto & zero_elem : zero_elems)
256 16 : mesh.delete_elem(mesh.elem_ptr(zero_elem));
257 27 : mesh.contract();
258 : // As we modified the side_list, it is safer to clear the node_list
259 27 : boundary_info.clear_boundary_node_ids();
260 27 : mesh.prepare_for_use();
261 27 : }
262 :
263 : bool
264 51264 : pointOnLine(const Real px,
265 : const Real py,
266 : const Real param_1,
267 : const Real param_2,
268 : const Real param_3,
269 : const Real dis_tol)
270 : {
271 51264 : return std::abs(px * param_1 + py * param_2 + param_3) <= dis_tol;
272 : }
273 :
274 : bool
275 82526 : lineSideDeterminator(const Real px,
276 : const Real py,
277 : const Real param_1,
278 : const Real param_2,
279 : const Real param_3,
280 : const bool direction_param,
281 : const Real dis_tol)
282 : {
283 82526 : const Real tmp = px * param_1 + py * param_2 + param_3;
284 82526 : return direction_param ? tmp >= dis_tol : tmp <= dis_tol;
285 : }
286 :
287 : Point
288 1763 : twoLineIntersection(const Real param_11,
289 : const Real param_12,
290 : const Real param_13,
291 : const Real param_21,
292 : const Real param_22,
293 : const Real param_23)
294 : {
295 : return Point(
296 1763 : (param_12 * param_23 - param_22 * param_13) / (param_11 * param_22 - param_21 * param_12),
297 1763 : (param_13 * param_21 - param_23 * param_11) / (param_11 * param_22 - param_21 * param_12),
298 1763 : 0.0);
299 : }
300 :
301 : Point
302 1763 : twoPointandLineIntersection(const Point & pt1,
303 : const Point & pt2,
304 : const Real param_1,
305 : const Real param_2,
306 : const Real param_3)
307 : {
308 5289 : return twoLineIntersection(param_1,
309 : param_2,
310 : param_3,
311 1763 : pt2(1) - pt1(1),
312 1763 : pt1(0) - pt2(0),
313 3526 : pt2(0) * pt1(1) - pt1(0) * pt2(1));
314 : }
315 :
316 : bool
317 27 : quasiTriElementsFixer(ReplicatedMesh & mesh,
318 : const std::set<subdomain_id_type> & subdomain_ids_set,
319 : const subdomain_id_type tri_elem_subdomain_shift,
320 : const SubdomainName tri_elem_subdomain_name_suffix)
321 : {
322 27 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
323 : // Define the subdomain id shift for the new TRI3 element subdomain(s)
324 27 : const subdomain_id_type max_subdomain_id(*subdomain_ids_set.rbegin());
325 27 : const subdomain_id_type tri_subdomain_id_shift =
326 27 : 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 27 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
332 27 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
333 27 : 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 2515 : 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 2488 : const auto elem_angles = vertex_angles(*elem);
341 2488 : const auto elem_distances = vertex_distances(*elem);
342 : // Type 1
343 2488 : if (MooseUtils::absoluteFuzzyEqual(elem_angles.front().first, M_PI, 0.001))
344 : {
345 202 : bad_elems_rec.push_back(std::make_tuple(elem, elem_angles.front().second, false, true));
346 202 : continue;
347 : }
348 : // Type 2
349 2286 : if (MooseUtils::absoluteFuzzyEqual(elem_distances.front().first, 0.0))
350 : {
351 8 : bad_elems_rec.push_back(std::make_tuple(elem, elem_distances.front().second, false, false));
352 : }
353 2717 : }
354 27 : std::set<subdomain_id_type> new_subdomain_ids;
355 : // Loop over all the identified degenerate QUAD elements
356 237 : for (const auto & bad_elem : bad_elems_rec)
357 : {
358 210 : std::vector<boundary_id_type> elem_bdry_container_0;
359 210 : std::vector<boundary_id_type> elem_bdry_container_1;
360 210 : std::vector<boundary_id_type> elem_bdry_container_2;
361 :
362 210 : Elem * elem_0 = std::get<0>(bad_elem);
363 210 : 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 202 : Elem * elem_1 = elem_0->neighbor_ptr(std::get<1>(bad_elem));
370 202 : Elem * elem_2 = elem_0->neighbor_ptr((std::get<1>(bad_elem) - 1) % elem_0->n_vertices());
371 202 : if ((elem_1 != nullptr || elem_2 != nullptr))
372 0 : throw MooseException("The input mesh has degenerate quad element before trimming.");
373 : }
374 210 : mesh.get_boundary_info().boundary_ids(elem_0, std::get<1>(bad_elem), elem_bdry_container_0);
375 210 : mesh.get_boundary_info().boundary_ids(
376 210 : elem_0, (std::get<1>(bad_elem) + 1) % elem_0->n_vertices(), elem_bdry_container_1);
377 210 : mesh.get_boundary_info().boundary_ids(
378 210 : elem_0, (std::get<1>(bad_elem) + 2) % elem_0->n_vertices(), elem_bdry_container_2);
379 210 : mesh.get_boundary_info().boundary_ids(
380 210 : 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 210 : 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 210 : auto pt0 = elem_0->node_ptr((std::get<1>(bad_elem) + 1) % elem_0->n_vertices());
386 210 : auto pt1 = elem_0->node_ptr((std::get<1>(bad_elem) + 2) % elem_0->n_vertices());
387 210 : 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 210 : for (const auto j : make_range(n_elem_extra_ids))
390 0 : exist_extra_ids[j] = elem_0->get_extra_integer(j);
391 : // Delete the degenerate QUAD element
392 210 : mesh.delete_elem(elem_0);
393 : // Create the new TRI element
394 210 : Elem * elem_Tri3 = mesh.add_elem(new Tri3);
395 210 : elem_Tri3->set_node(0, pt0);
396 210 : elem_Tri3->set_node(1, pt1);
397 210 : elem_Tri3->set_node(2, pt2);
398 : // Retain the boundary information
399 420 : for (auto bdry_id : elem_bdry_container_0)
400 210 : boundary_info.add_side(elem_Tri3, 2, bdry_id);
401 250 : for (auto bdry_id : elem_bdry_container_1)
402 40 : boundary_info.add_side(elem_Tri3, 0, bdry_id);
403 277 : for (auto bdry_id : elem_bdry_container_2)
404 67 : boundary_info.add_side(elem_Tri3, 1, bdry_id);
405 : // Assign subdomain id for the TRI element by shifting its original subdomain id
406 210 : elem_Tri3->subdomain_id() = elem_block_id + tri_subdomain_id_shift;
407 210 : new_subdomain_ids.emplace(elem_block_id + tri_subdomain_id_shift);
408 : // Retain element extra integers
409 210 : for (const auto j : make_range(n_elem_extra_ids))
410 0 : elem_Tri3->set_extra_integer(j, exist_extra_ids[j]);
411 210 : }
412 : // Assign subdomain names for the new TRI elements
413 115 : for (auto & nid : new_subdomain_ids)
414 : {
415 91 : const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
416 91 : if (MooseMeshUtils::getSubdomainID(
417 182 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
418 182 : : old_name) +
419 364 : "_" + tri_elem_subdomain_name_suffix,
420 91 : mesh) != Moose::INVALID_BLOCK_ID)
421 3 : throw MooseException("The new subdomain name already exists in the mesh.");
422 88 : mesh.subdomain_name(nid) =
423 176 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
424 176 : : old_name) +
425 264 : "_" + tri_elem_subdomain_name_suffix;
426 88 : mooseWarning("Degenerate QUAD elements have been converted into TRI elements with a new "
427 88 : "subdomain name: " +
428 176 : mesh.subdomain_name(nid) + ".");
429 91 : }
430 48 : return bad_elems_rec.size();
431 33 : }
432 :
433 : std::vector<std::pair<Real, unsigned int>>
434 2488 : vertex_angles(const Elem & elem)
435 : {
436 2488 : std::vector<std::pair<Real, unsigned int>> angles;
437 2488 : const unsigned int n_vertices = elem.n_vertices();
438 :
439 12440 : for (const auto i : make_range(n_vertices))
440 : {
441 9952 : Point v1 = (*elem.node_ptr((i - 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
442 9952 : Point v2 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
443 9952 : Real tmp = v1 * v2 / v1.norm() / v2.norm();
444 9952 : if (tmp > 1.0)
445 0 : tmp = 1.0;
446 9952 : else if (tmp < -1.0)
447 43 : tmp = -1.0;
448 9952 : angles.push_back(std::make_pair(acos(tmp), i));
449 : }
450 2488 : std::sort(angles.begin(), angles.end(), std::greater<>());
451 2488 : return angles;
452 0 : }
453 :
454 : std::vector<std::pair<Real, unsigned int>>
455 2488 : vertex_distances(const Elem & elem)
456 : {
457 2488 : std::vector<std::pair<Real, unsigned int>> distances;
458 2488 : const unsigned int n_vertices = elem.n_vertices();
459 :
460 12440 : for (const auto i : make_range(n_vertices))
461 : {
462 9952 : Point v1 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
463 9952 : distances.push_back(std::make_pair(v1.norm(), i));
464 : }
465 2488 : std::sort(distances.begin(), distances.end());
466 2488 : return distances;
467 0 : }
468 :
469 : void
470 800 : triElemSplitter(ReplicatedMesh & mesh,
471 : const dof_id_type elem_id,
472 : const unsigned short node_shift,
473 : const dof_id_type nid_3,
474 : const dof_id_type nid_4,
475 : const subdomain_id_type single_elem_side_id,
476 : const subdomain_id_type double_elem_side_id)
477 : {
478 800 : const auto elem_old = mesh.elem_ptr(elem_id);
479 800 : const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
480 800 : const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
481 800 : const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
482 :
483 : const bool m1_side_flag =
484 800 : MooseUtils::absoluteFuzzyEqual((*(mesh.node_ptr(nid_3)) - *(mesh.node_ptr(nid_0)))
485 800 : .cross(*(mesh.node_ptr(nid_1)) - *(mesh.node_ptr(nid_0)))
486 800 : .norm(),
487 1600 : 0.0);
488 800 : const dof_id_type nid_m1 = m1_side_flag ? nid_3 : nid_4;
489 800 : const dof_id_type nid_m2 = m1_side_flag ? nid_4 : nid_3;
490 : // Build boundary information of the mesh
491 800 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
492 800 : auto bdry_side_list = boundary_info.build_side_list();
493 : // Create a list of sidesets involving the element to be split
494 800 : std::vector<std::vector<boundary_id_type>> elem_side_list;
495 800 : elem_side_list.resize(3);
496 115680 : for (const auto i : index_range(bdry_side_list))
497 : {
498 114880 : if (std::get<0>(bdry_side_list[i]) == elem_id)
499 : {
500 352 : elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
501 176 : std::get<2>(bdry_side_list[i]));
502 : }
503 : }
504 :
505 800 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
506 800 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
507 : // Record all the element extra integers of the original element
508 800 : for (const auto j : make_range(n_elem_extra_ids))
509 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
510 :
511 800 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
512 800 : elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
513 800 : elem_Tri3_0->set_node(1, mesh.node_ptr(nid_m1));
514 800 : elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m2));
515 800 : elem_Tri3_0->subdomain_id() = single_elem_side_id;
516 800 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
517 800 : elem_Tri3_1->set_node(0, mesh.node_ptr(nid_1));
518 800 : elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m2));
519 800 : elem_Tri3_1->set_node(2, mesh.node_ptr(nid_m1));
520 800 : elem_Tri3_1->subdomain_id() = double_elem_side_id;
521 800 : Elem * elem_Tri3_2 = mesh.add_elem(new Tri3);
522 800 : elem_Tri3_2->set_node(0, mesh.node_ptr(nid_2));
523 800 : elem_Tri3_2->set_node(1, mesh.node_ptr(nid_m2));
524 800 : elem_Tri3_2->set_node(2, mesh.node_ptr(nid_1));
525 800 : elem_Tri3_2->subdomain_id() = double_elem_side_id;
526 : // Retain element extra integers
527 800 : for (const auto j : make_range(n_elem_extra_ids))
528 : {
529 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
530 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
531 0 : elem_Tri3_2->set_extra_integer(j, exist_extra_ids[j]);
532 : }
533 :
534 : // Add sideset information to the new elements
535 848 : for (const auto & side_info_0 : elem_side_list[0])
536 : {
537 48 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
538 48 : boundary_info.add_side(elem_Tri3_1, 2, side_info_0);
539 : }
540 848 : for (const auto & side_info_1 : elem_side_list[1])
541 48 : boundary_info.add_side(elem_Tri3_2, 2, side_info_1);
542 880 : for (const auto & side_info_2 : elem_side_list[2])
543 : {
544 80 : boundary_info.add_side(elem_Tri3_0, 2, side_info_2);
545 80 : boundary_info.add_side(elem_Tri3_2, 0, side_info_2);
546 : }
547 800 : }
548 :
549 : void
550 96 : triElemSplitter(ReplicatedMesh & mesh,
551 : const dof_id_type elem_id,
552 : const unsigned short node_shift,
553 : const dof_id_type nid_m,
554 : const subdomain_id_type first_elem_side_id,
555 : const subdomain_id_type second_elem_side_id)
556 : {
557 96 : const auto elem_old = mesh.elem_ptr(elem_id);
558 96 : const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
559 96 : const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
560 96 : const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
561 : // Build boundary information of the mesh
562 96 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
563 96 : auto bdry_side_list = boundary_info.build_side_list();
564 : // Create a list of sidesets involving the element to be split
565 96 : std::vector<std::vector<boundary_id_type>> elem_side_list;
566 96 : elem_side_list.resize(3);
567 10560 : for (const auto i : index_range(bdry_side_list))
568 : {
569 10464 : if (std::get<0>(bdry_side_list[i]) == elem_id)
570 : {
571 96 : elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
572 48 : std::get<2>(bdry_side_list[i]));
573 : }
574 : }
575 :
576 96 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
577 96 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
578 : // Record all the element extra integers of the original element
579 96 : for (const auto j : make_range(n_elem_extra_ids))
580 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
581 :
582 96 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
583 96 : elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
584 96 : elem_Tri3_0->set_node(1, mesh.node_ptr(nid_1));
585 96 : elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m));
586 96 : elem_Tri3_0->subdomain_id() = first_elem_side_id;
587 96 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
588 96 : elem_Tri3_1->set_node(0, mesh.node_ptr(nid_0));
589 96 : elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m));
590 96 : elem_Tri3_1->set_node(2, mesh.node_ptr(nid_2));
591 96 : elem_Tri3_1->subdomain_id() = second_elem_side_id;
592 : // Retain element extra integers
593 96 : for (const auto j : make_range(n_elem_extra_ids))
594 : {
595 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
596 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
597 : }
598 :
599 : // Add sideset information to the new elements
600 112 : for (const auto & side_info_0 : elem_side_list[0])
601 16 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
602 96 : for (const auto & side_info_1 : elem_side_list[1])
603 : {
604 0 : boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
605 0 : boundary_info.add_side(elem_Tri3_1, 1, side_info_1);
606 : }
607 128 : for (const auto & side_info_2 : elem_side_list[2])
608 32 : boundary_info.add_side(elem_Tri3_1, 2, side_info_2);
609 96 : }
610 :
611 : void
612 672 : quadElemSplitter(ReplicatedMesh & mesh,
613 : const dof_id_type elem_id,
614 : const subdomain_id_type tri_elem_subdomain_shift)
615 : {
616 : // Build boundary information of the mesh
617 672 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
618 672 : auto bdry_side_list = boundary_info.build_side_list();
619 : // Create a list of sidesets involving the element to be split
620 672 : std::vector<std::vector<boundary_id_type>> elem_side_list;
621 672 : elem_side_list.resize(4);
622 81056 : for (const auto i : index_range(bdry_side_list))
623 : {
624 80384 : if (std::get<0>(bdry_side_list[i]) == elem_id)
625 : {
626 326 : elem_side_list[std::get<1>(bdry_side_list[i])].push_back(std::get<2>(bdry_side_list[i]));
627 : }
628 : }
629 :
630 672 : auto node_0 = mesh.elem_ptr(elem_id)->node_ptr(0);
631 672 : auto node_1 = mesh.elem_ptr(elem_id)->node_ptr(1);
632 672 : auto node_2 = mesh.elem_ptr(elem_id)->node_ptr(2);
633 672 : auto node_3 = mesh.elem_ptr(elem_id)->node_ptr(3);
634 :
635 672 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
636 672 : 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 672 : for (const auto j : make_range(n_elem_extra_ids))
639 0 : 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 672 : if (std::abs((*node_1 - *node_0).cross(*node_3 - *node_0).norm() -
644 672 : (*node_1 - *node_2).cross(*node_3 - *node_2).norm()) >
645 672 : std::abs((*node_0 - *node_1).cross(*node_2 - *node_1).norm() -
646 672 : (*node_0 - *node_3).cross(*node_2 - *node_3).norm()))
647 : {
648 196 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
649 196 : elem_Tri3_0->set_node(0, node_0);
650 196 : elem_Tri3_0->set_node(1, node_1);
651 196 : elem_Tri3_0->set_node(2, node_2);
652 196 : elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
653 196 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
654 196 : elem_Tri3_1->set_node(0, node_0);
655 196 : elem_Tri3_1->set_node(1, node_2);
656 196 : elem_Tri3_1->set_node(2, node_3);
657 196 : elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
658 : // Retain element extra integers
659 196 : for (const auto j : make_range(n_elem_extra_ids))
660 : {
661 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
662 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
663 : }
664 :
665 : // Add sideset information to the new elements
666 196 : for (const auto & side_info_0 : elem_side_list[0])
667 0 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
668 292 : for (const auto & side_info_1 : elem_side_list[1])
669 96 : boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
670 196 : for (const auto & side_info_2 : elem_side_list[2])
671 0 : boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
672 196 : for (const auto & side_info_3 : elem_side_list[3])
673 0 : boundary_info.add_side(elem_Tri3_1, 2, side_info_3);
674 : }
675 : else
676 : {
677 476 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
678 476 : elem_Tri3_0->set_node(0, node_0);
679 476 : elem_Tri3_0->set_node(1, node_1);
680 476 : elem_Tri3_0->set_node(2, node_3);
681 476 : elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
682 476 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
683 476 : elem_Tri3_1->set_node(0, node_1);
684 476 : elem_Tri3_1->set_node(1, node_2);
685 476 : elem_Tri3_1->set_node(2, node_3);
686 476 : elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
687 : // Retain element extra integers
688 476 : for (const auto j : make_range(n_elem_extra_ids))
689 : {
690 0 : elem_Tri3_0->set_extra_integer(j, exist_extra_ids[j]);
691 0 : elem_Tri3_1->set_extra_integer(j, exist_extra_ids[j]);
692 : }
693 :
694 : // Add sideset information to the new elements
695 492 : for (const auto & side_info_0 : elem_side_list[0])
696 16 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
697 604 : for (const auto & side_info_1 : elem_side_list[1])
698 128 : boundary_info.add_side(elem_Tri3_1, 0, side_info_1);
699 514 : for (const auto & side_info_2 : elem_side_list[2])
700 38 : boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
701 524 : for (const auto & side_info_3 : elem_side_list[3])
702 48 : boundary_info.add_side(elem_Tri3_0, 2, side_info_3);
703 : }
704 672 : }
705 :
706 : void
707 43 : quadToTriOnLine(ReplicatedMesh & mesh,
708 : const std::vector<Real> & cut_line_params,
709 : const dof_id_type tri_subdomain_id_shift,
710 : const SubdomainName tri_elem_subdomain_name_suffix)
711 : {
712 : // Preprocess: find all the quad elements that are across the cutting line
713 43 : std::vector<dof_id_type> cross_elems_quad;
714 43 : std::set<subdomain_id_type> new_subdomain_ids;
715 6811 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
716 : elem_it++)
717 : {
718 6768 : if ((*elem_it)->n_vertices() == 4)
719 : {
720 6768 : std::vector<unsigned short> node_side_rec;
721 33840 : for (const auto i : make_range(4))
722 : {
723 27072 : const Point v_point = (*elem_it)->point(i);
724 27072 : if (!pointOnLine(
725 27072 : v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2]))
726 26616 : node_side_rec.push_back(lineSideDeterminator(v_point(0),
727 26616 : v_point(1),
728 26616 : cut_line_params[0],
729 26616 : cut_line_params[1],
730 26616 : 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 6768 : const auto num_nodes = std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0);
740 6768 : if (num_nodes != (int)node_side_rec.size() && num_nodes > 0)
741 : {
742 672 : cross_elems_quad.push_back((*elem_it)->id());
743 672 : new_subdomain_ids.emplace((*elem_it)->subdomain_id() + tri_subdomain_id_shift);
744 : }
745 6768 : }
746 43 : }
747 : // Then convert these quad elements into tri elements
748 715 : for (const auto & cross_elem_quad : cross_elems_quad)
749 : {
750 672 : quadElemSplitter(mesh, cross_elem_quad, tri_subdomain_id_shift);
751 672 : mesh.delete_elem(mesh.elem_ptr(cross_elem_quad));
752 : }
753 139 : for (auto & nid : new_subdomain_ids)
754 : {
755 99 : const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
756 99 : if (MooseMeshUtils::getSubdomainID(
757 198 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
758 198 : : old_name) +
759 396 : "_" + tri_elem_subdomain_name_suffix,
760 99 : mesh) != Moose::INVALID_BLOCK_ID)
761 3 : throw MooseException("The new subdomain name already exists in the mesh.");
762 96 : mesh.subdomain_name(nid) =
763 192 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
764 192 : : old_name) +
765 288 : "_" + tri_elem_subdomain_name_suffix;
766 96 : mooseWarning("QUAD elements have been converted into TRI elements with a new "
767 96 : "subdomain name: " +
768 192 : mesh.subdomain_name(nid) + ".");
769 99 : }
770 40 : mesh.contract();
771 46 : }
772 :
773 : void
774 40 : lineRemoverCutElemTri(ReplicatedMesh & mesh,
775 : const std::vector<Real> & cut_line_params,
776 : const subdomain_id_type block_id_to_remove,
777 : const boundary_id_type new_boundary_id)
778 : {
779 : // Find all the elements that are across the cutting line
780 40 : std::vector<dof_id_type> cross_elems;
781 : // A vector for element specific information
782 40 : std::vector<std::vector<std::pair<dof_id_type, dof_id_type>>> node_pairs_vec;
783 : // A set for unique pairs
784 40 : std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs_unique_vec;
785 6376 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
786 : elem_it++)
787 : {
788 6336 : const auto n_vertices = (*elem_it)->n_vertices();
789 6336 : unsigned int n_points_on_line = 0;
790 6336 : std::vector<unsigned short> node_side_rec(n_vertices, 0);
791 30528 : for (const auto i : make_range(n_vertices))
792 : {
793 : // First check if the vertex is in the XY Plane
794 24192 : if (!MooseUtils::absoluteFuzzyEqual((*elem_it)->point(i)(2), 0.0))
795 0 : mooseError("MooseMeshXYCuttingUtils::lineRemoverCutElemTri() only works for 2D meshes in "
796 : "XY Plane.");
797 24192 : const Point v_point = (*elem_it)->point(i);
798 24192 : if (pointOnLine(
799 24192 : v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2]))
800 544 : ++n_points_on_line;
801 : else
802 23648 : node_side_rec[i] = lineSideDeterminator(v_point(0),
803 23648 : v_point(1),
804 23648 : cut_line_params[0],
805 23648 : cut_line_params[1],
806 23648 : 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 6336 : const unsigned int num_nodes = std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0);
816 6336 : if (num_nodes == node_side_rec.size() - n_points_on_line)
817 : {
818 3552 : (*elem_it)->subdomain_id() = block_id_to_remove;
819 : }
820 2784 : else if (num_nodes > 0)
821 : {
822 896 : if ((*elem_it)->n_vertices() != 3 || (*elem_it)->n_nodes() != 3)
823 0 : mooseError("The element across the cutting line is not TRI3, which is not supported.");
824 896 : cross_elems.push_back((*elem_it)->id());
825 : // Then we need to check pairs of nodes that are on the different side
826 896 : std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs;
827 3584 : for (const auto i : index_range(node_side_rec))
828 : {
829 : // first node on removal side and second node on retaining side
830 2688 : if (node_side_rec[i] > 0 && node_side_rec[(i + 1) % node_side_rec.size()] == 0)
831 : {
832 : // Removal side first
833 896 : node_pairs.push_back(
834 896 : std::make_pair((*elem_it)->node_ptr(i)->id(),
835 896 : (*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id()));
836 896 : node_pairs_unique_vec.push_back(node_pairs.back());
837 : }
838 : // first node on retaining side and second node on removal side
839 1792 : else if (node_side_rec[i] == 0 && node_side_rec[(i + 1) % node_side_rec.size()] > 0)
840 : {
841 : // Removal side first
842 896 : node_pairs.push_back(
843 896 : std::make_pair((*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id(),
844 896 : (*elem_it)->node_ptr(i)->id()));
845 896 : node_pairs_unique_vec.push_back(node_pairs.back());
846 : }
847 : }
848 896 : node_pairs_vec.push_back(node_pairs);
849 896 : }
850 6376 : }
851 40 : auto vec_ip = std::unique(node_pairs_unique_vec.begin(), node_pairs_unique_vec.end());
852 80 : 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 40 : std::vector<Node *> nodes_on_line;
856 : // whether the on-line node is overlapped with the node pairs or a brand new node
857 40 : std::vector<unsigned short> nodes_on_line_overlap;
858 1688 : for (const auto & node_pair : node_pairs_unique_vec)
859 : {
860 1648 : const Point pt1 = *mesh.node_ptr(node_pair.first);
861 1648 : const Point pt2 = *mesh.node_ptr(node_pair.second);
862 1648 : const Point pt_line = twoPointandLineIntersection(
863 1648 : pt1, pt2, cut_line_params[0], cut_line_params[1], cut_line_params[2]);
864 1648 : if ((pt_line - pt1).norm() < libMesh::TOLERANCE)
865 : {
866 0 : nodes_on_line.push_back(mesh.node_ptr(node_pair.first));
867 0 : nodes_on_line_overlap.push_back(1);
868 : }
869 1648 : else if ((pt_line - pt2).norm() < libMesh::TOLERANCE)
870 : {
871 96 : nodes_on_line.push_back(mesh.node_ptr(node_pair.second));
872 96 : nodes_on_line_overlap.push_back(2);
873 : }
874 : else
875 : {
876 1552 : nodes_on_line.push_back(mesh.add_point(pt_line));
877 1552 : nodes_on_line_overlap.push_back(0);
878 : }
879 : }
880 :
881 : // make new elements
882 936 : for (const auto i : index_range(cross_elems))
883 : {
884 : // Only TRI elements are involved after preprocessing
885 896 : auto cross_elem = mesh.elem_ptr(cross_elems[i]);
886 896 : auto node_0 = cross_elem->node_ptr(0);
887 896 : auto node_1 = cross_elem->node_ptr(1);
888 896 : auto node_2 = cross_elem->node_ptr(2);
889 1792 : const std::vector<dof_id_type> tri_nodes = {node_0->id(), node_1->id(), node_2->id()};
890 :
891 896 : 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 896 : node_pairs_vec[i][0]));
895 896 : 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 896 : node_pairs_vec[i][1]));
899 896 : auto node_3 = nodes_on_line[online_node_index_1];
900 896 : auto node_4 = nodes_on_line[online_node_index_2];
901 896 : const auto node_3_overlap_flag = nodes_on_line_overlap[online_node_index_1];
902 896 : const auto node_4_overlap_flag = nodes_on_line_overlap[online_node_index_2];
903 : // Most common case, no overlapped nodes
904 896 : 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 800 : 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 800 : common_node_side ? block_id_to_remove : cross_elem->subdomain_id();
910 : const subdomain_id_type block_id_to_assign_2 =
911 800 : 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 800 : common_node_side ? node_pairs_vec[i][0].first : node_pairs_vec[i][0].second;
916 :
917 1600 : triElemSplitter(mesh,
918 : cross_elem->id(),
919 1600 : 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 800 : mesh.delete_elem(cross_elem);
926 800 : }
927 : // both node_3 and node_4 are overlapped
928 96 : 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 0 : cross_elem->subdomain_id() = lineSideDeterminator(cross_elem->vertex_average()(0),
933 0 : cross_elem->vertex_average()(1),
934 0 : cut_line_params[0],
935 0 : cut_line_params[1],
936 0 : cut_line_params[2],
937 : true)
938 0 : ? block_id_to_remove
939 0 : : cross_elem->subdomain_id();
940 : }
941 : // node_3 or node_4 is overlapped
942 : else
943 : {
944 96 : const auto node_3_finder = std::distance(
945 96 : tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_3->id()));
946 96 : const auto node_4_finder = std::distance(
947 96 : 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 96 : const dof_id_type node_id = node_3_finder < node_4_finder ? node_4->id() : node_3->id();
951 96 : const auto node_finder = std::min(node_3_finder, node_4_finder);
952 :
953 144 : triElemSplitter(
954 : mesh,
955 : cross_elem->id(),
956 : node_finder,
957 : node_id,
958 96 : tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
959 48 : ? block_id_to_remove
960 48 : : cross_elem->subdomain_id(),
961 96 : tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
962 48 : ? cross_elem->subdomain_id()
963 : : block_id_to_remove);
964 96 : mesh.delete_elem(cross_elem);
965 : }
966 896 : }
967 40 : 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.
971 40 : mesh.find_neighbors();
972 40 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
973 8072 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
974 : elem_it++)
975 : {
976 8032 : if ((*elem_it)->subdomain_id() != block_id_to_remove)
977 : {
978 14448 : for (const auto j : make_range((*elem_it)->n_sides()))
979 : {
980 11280 : if ((*elem_it)->neighbor_ptr(j) != nullptr)
981 10704 : if ((*elem_it)->neighbor_ptr(j)->subdomain_id() == block_id_to_remove)
982 960 : boundary_info.add_side(*elem_it, j, new_boundary_id);
983 : }
984 : }
985 40 : }
986 :
987 : // Delete the block to remove
988 4904 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
989 4904 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
990 : elem_it++)
991 4904 : mesh.delete_elem(*elem_it);
992 40 : mesh.contract();
993 40 : }
994 :
995 : void
996 43 : lineRemoverCutElem(ReplicatedMesh & mesh,
997 : const std::vector<Real> & cut_line_params,
998 : const dof_id_type tri_subdomain_id_shift,
999 : const SubdomainName tri_elem_subdomain_name_suffix,
1000 : const subdomain_id_type block_id_to_remove,
1001 : const boundary_id_type new_boundary_id,
1002 : const bool improve_boundary_tri_elems)
1003 : {
1004 : // Convert any quad elements crossed by the line into tri elements
1005 43 : 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 40 : lineRemoverCutElemTri(mesh, cut_line_params, block_id_to_remove, new_boundary_id);
1009 :
1010 40 : if (improve_boundary_tri_elems)
1011 8 : boundaryTriElemImprover(mesh, new_boundary_id);
1012 40 : }
1013 :
1014 : void
1015 8 : boundaryTriElemImprover(ReplicatedMesh & mesh, const boundary_id_type boundary_to_improve)
1016 : {
1017 8 : if (!MooseMeshUtils::hasBoundaryID(mesh, boundary_to_improve))
1018 0 : mooseError(
1019 : "MooseMeshXYCuttingUtils::boundaryTriElemImprover(): The boundary_to_improve provided "
1020 : "does not exist in the given mesh.");
1021 8 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1022 8 : 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 8 : tri3_elem_info;
1031 904 : for (const auto & side : side_list)
1032 : {
1033 896 : if (std::get<2>(side) == boundary_to_improve)
1034 : {
1035 416 : Elem * elem = mesh.elem_ptr(std::get<0>(side));
1036 416 : if (elem->type() == TRI3)
1037 : {
1038 416 : const auto key_node_id = elem->node_id((std::get<1>(side) + 2) % 3);
1039 416 : const auto value_elem_id = elem->id();
1040 416 : const auto value_node_id_1 = elem->node_id(std::get<1>(side));
1041 416 : const auto value_node_id_2 = elem->node_id((std::get<1>(side) + 1) % 3);
1042 832 : tri3_elem_info[key_node_id].push_back(
1043 832 : 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 8 : 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 216 : 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 208 : std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
1055 208 : std::vector<dof_id_type> elem_id_list;
1056 624 : for (const auto & tri : tri_group.second)
1057 : {
1058 416 : node_assm.push_back(std::make_pair(std::get<1>(tri), std::get<2>(tri)));
1059 416 : elem_id_list.push_back(std::get<0>(tri));
1060 : }
1061 208 : std::vector<dof_id_type> ordered_node_list;
1062 208 : std::vector<dof_id_type> ordered_elem_list;
1063 208 : MooseMeshUtils::makeOrderedNodeList(
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 208 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
1071 208 : std::vector<std::tuple<subdomain_id_type, std::vector<dof_id_type>, unsigned int>> blocks_info;
1072 624 : for (const auto & elem_id : ordered_elem_list)
1073 : {
1074 416 : 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 416 : for (const auto j : make_range(n_elem_extra_ids))
1077 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
1078 416 : if (!blocks_info.empty())
1079 : {
1080 400 : if (mesh.elem_ptr(elem_id)->subdomain_id() == std::get<0>(blocks_info.back()) &&
1081 192 : exist_extra_ids == std::get<1>(blocks_info.back()))
1082 : {
1083 192 : std::get<2>(blocks_info.back())++;
1084 192 : continue;
1085 : }
1086 : }
1087 224 : blocks_info.push_back(
1088 448 : std::make_tuple(mesh.elem_ptr(elem_id)->subdomain_id(), exist_extra_ids, 1));
1089 416 : }
1090 : // For each separated subdomain / set of extra ids, we try to improve the boundary elements
1091 208 : unsigned int side_counter = 0;
1092 432 : for (const auto & block_info : blocks_info)
1093 : {
1094 224 : const auto node_1 = mesh.node_ptr(ordered_node_list[side_counter]);
1095 : // we do not need to subtract 1 for node_2
1096 224 : const auto node_2 = mesh.node_ptr(ordered_node_list[side_counter + std::get<2>(block_info)]);
1097 224 : const auto node_0 = mesh.node_ptr(tri_group.first);
1098 224 : const Point v1 = *node_1 - *node_0;
1099 224 : const Point v2 = *node_2 - *node_0;
1100 224 : const Real angle = std::acos(v1 * v2 / v1.norm() / v2.norm()) / M_PI * 180.0;
1101 0 : const std::vector<dof_id_type> block_elems(ordered_elem_list.begin() + side_counter,
1102 448 : ordered_elem_list.begin() + side_counter +
1103 448 : 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 224 : elemSideLocator(mesh,
1112 224 : block_elems.front(),
1113 224 : tri_group.first,
1114 224 : ordered_node_list[side_counter],
1115 : side_id_0,
1116 : is_inverse_0);
1117 224 : elemSideLocator(mesh,
1118 224 : block_elems.back(),
1119 224 : ordered_node_list[side_counter + std::get<2>(block_info)],
1120 224 : tri_group.first,
1121 : side_id_t,
1122 : is_inverse_t);
1123 : // Collect boundary information of the identified sides
1124 224 : std::vector<boundary_id_type> side_0_boundary_ids;
1125 448 : boundary_info.boundary_ids(
1126 224 : mesh.elem_ptr(block_elems.front()), side_id_0, side_0_boundary_ids);
1127 224 : std::vector<boundary_id_type> side_t_boundary_ids;
1128 224 : 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 224 : if (angle < 90.0)
1136 : {
1137 152 : if (std::get<2>(block_info) > 1)
1138 : {
1139 192 : makeImprovedTriElement(mesh,
1140 64 : tri_group.first,
1141 64 : ordered_node_list[side_counter],
1142 64 : ordered_node_list[side_counter + std::get<2>(block_info)],
1143 64 : std::get<0>(block_info),
1144 64 : std::get<1>(block_info),
1145 : {boundary_to_improve},
1146 : side_0_boundary_ids,
1147 : side_t_boundary_ids);
1148 64 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1149 : }
1150 : }
1151 72 : else if (angle < 135.0)
1152 : {
1153 : // We can just add the middle node because there's nothing on the other side
1154 56 : const auto node_m = mesh.add_point((*node_1 + *node_2) / 2.0);
1155 168 : makeImprovedTriElement(mesh,
1156 56 : tri_group.first,
1157 56 : ordered_node_list[side_counter],
1158 : node_m->id(),
1159 56 : std::get<0>(block_info),
1160 56 : std::get<1>(block_info),
1161 : {boundary_to_improve},
1162 : side_0_boundary_ids,
1163 112 : std::vector<boundary_id_type>());
1164 224 : makeImprovedTriElement(mesh,
1165 56 : tri_group.first,
1166 : node_m->id(),
1167 56 : ordered_node_list[side_counter + std::get<2>(block_info)],
1168 56 : std::get<0>(block_info),
1169 56 : std::get<1>(block_info),
1170 : {boundary_to_improve},
1171 112 : std::vector<boundary_id_type>(),
1172 : side_t_boundary_ids);
1173 56 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1174 : }
1175 : else
1176 : {
1177 16 : const auto node_m1 = mesh.add_point((*node_1 * 2.0 + *node_2) / 3.0);
1178 16 : const auto node_m2 = mesh.add_point((*node_1 + *node_2 * 2.0) / 3.0);
1179 48 : makeImprovedTriElement(mesh,
1180 16 : tri_group.first,
1181 16 : ordered_node_list[side_counter],
1182 : node_m1->id(),
1183 16 : std::get<0>(block_info),
1184 16 : std::get<1>(block_info),
1185 : {boundary_to_improve},
1186 : side_0_boundary_ids,
1187 32 : std::vector<boundary_id_type>());
1188 64 : makeImprovedTriElement(mesh,
1189 16 : tri_group.first,
1190 : node_m1->id(),
1191 : node_m2->id(),
1192 16 : std::get<0>(block_info),
1193 16 : std::get<1>(block_info),
1194 : {boundary_to_improve},
1195 32 : std::vector<boundary_id_type>(),
1196 32 : std::vector<boundary_id_type>());
1197 64 : makeImprovedTriElement(mesh,
1198 16 : tri_group.first,
1199 : node_m2->id(),
1200 16 : ordered_node_list[side_counter + std::get<2>(block_info)],
1201 16 : std::get<0>(block_info),
1202 16 : std::get<1>(block_info),
1203 : {boundary_to_improve},
1204 32 : std::vector<boundary_id_type>(),
1205 : side_t_boundary_ids);
1206 16 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1207 : }
1208 224 : side_counter += std::get<2>(block_info);
1209 224 : }
1210 : // TODO: Need to check if the new element is inverted?
1211 208 : }
1212 : // Delete the original elements
1213 336 : for (const auto & elem_to_remove : elems_to_remove)
1214 328 : mesh.delete_elem(mesh.elem_ptr(elem_to_remove));
1215 8 : mesh.contract();
1216 8 : }
1217 :
1218 : void
1219 224 : makeImprovedTriElement(ReplicatedMesh & mesh,
1220 : const dof_id_type node_id_0,
1221 : const dof_id_type node_id_1,
1222 : const dof_id_type node_id_2,
1223 : const subdomain_id_type subdomain_id,
1224 : const std::vector<dof_id_type> & extra_elem_ids,
1225 : const std::vector<boundary_id_type> & boundary_ids_for_side_1,
1226 : const std::vector<boundary_id_type> & boundary_ids_for_side_0,
1227 : const std::vector<boundary_id_type> & boundary_ids_for_side_2)
1228 : {
1229 224 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1230 224 : Elem * elem_Tri3_new = mesh.add_elem(new Tri3);
1231 224 : elem_Tri3_new->set_node(0, mesh.node_ptr(node_id_0));
1232 224 : elem_Tri3_new->set_node(1, mesh.node_ptr(node_id_1));
1233 224 : elem_Tri3_new->set_node(2, mesh.node_ptr(node_id_2));
1234 264 : for (const auto & boundary_id_for_side_0 : boundary_ids_for_side_0)
1235 40 : boundary_info.add_side(elem_Tri3_new, 0, boundary_id_for_side_0);
1236 448 : for (const auto & boundary_id_for_side_1 : boundary_ids_for_side_1)
1237 224 : boundary_info.add_side(elem_Tri3_new, 1, boundary_id_for_side_1);
1238 224 : for (const auto & boundary_id_for_side_2 : boundary_ids_for_side_2)
1239 0 : boundary_info.add_side(elem_Tri3_new, 2, boundary_id_for_side_2);
1240 224 : elem_Tri3_new->subdomain_id() = subdomain_id;
1241 : // Retain element extra integers
1242 224 : for (const auto j : index_range(extra_elem_ids))
1243 : {
1244 0 : elem_Tri3_new->set_extra_integer(j, extra_elem_ids[j]);
1245 : }
1246 224 : }
1247 :
1248 : bool
1249 448 : elemSideLocator(ReplicatedMesh & mesh,
1250 : const dof_id_type elem_id,
1251 : const dof_id_type node_id_0,
1252 : const dof_id_type node_id_1,
1253 : unsigned short & side_id,
1254 : bool & is_inverse)
1255 : {
1256 448 : Elem * elem = mesh.elem_ptr(elem_id);
1257 896 : for (unsigned short i = 0; i < elem->n_sides(); i++)
1258 : {
1259 2088 : if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_0 &&
1260 1192 : elem->side_ptr(i)->node_ptr(1)->id() == node_id_1)
1261 : {
1262 144 : side_id = i;
1263 144 : is_inverse = false;
1264 144 : return true;
1265 : }
1266 1880 : else if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_1 &&
1267 1128 : elem->side_ptr(i)->node_ptr(1)->id() == node_id_0)
1268 : {
1269 304 : side_id = i;
1270 304 : is_inverse = true;
1271 304 : return true;
1272 : }
1273 : }
1274 0 : return false;
1275 : }
1276 : }
|
