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 31 : 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 31 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
39 31 : auto bdry_side_list = boundary_info.build_side_list();
40 : // Only select the boundaries_to_conform
41 31 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> slc_bdry_side_list;
42 3103 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
43 5088 : if (std::get<2>(bdry_side_list[i]) == external_boundary_id ||
44 2016 : std::find(other_boundaries_to_conform.begin(),
45 : other_boundaries_to_conform.end(),
46 7104 : std::get<2>(bdry_side_list[i])) != other_boundaries_to_conform.end())
47 1920 : 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 31 : std::vector<dof_id_type> crossed_elems_to_remove;
52 7567 : 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 7536 : unsigned short removal_side_count = 0;
57 37680 : for (unsigned int i = 0; i < (*elem_it)->n_vertices(); i++)
58 : {
59 : // First check if the vertex is on the XY-Plane
60 30144 : 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 30144 : if (lineSideDeterminator((*elem_it)->point(i)(0),
64 30144 : (*elem_it)->point(i)(1),
65 30144 : bdry_pars[0],
66 30144 : bdry_pars[1],
67 30144 : bdry_pars[2],
68 : side_to_remove))
69 18656 : removal_side_count++;
70 : }
71 7536 : if (removal_side_count == (*elem_it)->n_vertices())
72 : {
73 4254 : (*elem_it)->subdomain_id() = block_id_to_remove;
74 4254 : continue;
75 : }
76 : // Check the average of the vertices of the element
77 3282 : if (lineSideDeterminator((*elem_it)->vertex_average()(0),
78 6564 : (*elem_it)->vertex_average()(1),
79 3282 : bdry_pars[0],
80 3282 : bdry_pars[1],
81 3282 : bdry_pars[2],
82 : side_to_remove))
83 486 : crossed_elems_to_remove.push_back((*elem_it)->id());
84 31 : }
85 : // Check each crossed element to see if removing it would lead to boundary moving
86 517 : for (const auto & elem_id : crossed_elems_to_remove)
87 : {
88 486 : bool remove_flag = true;
89 2034 : for (unsigned int i = 0; i < mesh.elem_ptr(elem_id)->n_sides(); i++)
90 : {
91 1674 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i) != nullptr)
92 2490 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id() != block_id_to_remove &&
93 887 : std::find(crossed_elems_to_remove.begin(),
94 : crossed_elems_to_remove.end(),
95 1774 : mesh.elem_ptr(elem_id)->neighbor_ptr(i)->id()) ==
96 2490 : crossed_elems_to_remove.end())
97 : {
98 636 : if (mesh.elem_ptr(elem_id)->subdomain_id() !=
99 636 : mesh.elem_ptr(elem_id)->neighbor_ptr(i)->subdomain_id())
100 : {
101 126 : remove_flag = false;
102 126 : break;
103 : }
104 : }
105 : }
106 486 : if (remove_flag)
107 360 : 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 31 : std::vector<dof_id_type> node_list;
114 2953 : for (auto elem_it = mesh.active_subdomain_set_elements_begin(subdomain_ids_set);
115 2953 : elem_it != mesh.active_subdomain_set_elements_end(subdomain_ids_set);
116 : elem_it++)
117 : {
118 14610 : for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
119 : {
120 11688 : if ((*elem_it)->neighbor_ptr(i) != nullptr)
121 11124 : if ((*elem_it)->neighbor_ptr(i)->subdomain_id() == block_id_to_remove)
122 : {
123 838 : node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(0)->id());
124 838 : node_list.push_back((*elem_it)->side_ptr(i)->node_ptr(1)->id());
125 838 : boundary_info.add_side(*elem_it, i, trimming_section_boundary_id);
126 838 : 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 31 : }
131 : // Remove duplicate nodes
132 31 : const auto unique_it = std::unique(node_list.begin(), node_list.end());
133 31 : 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 31 : std::vector<bool> node_list_flag(node_list.size(), false);
137 31 : std::vector<Point> node_list_point(node_list.size(), Point(0.0, 0.0, 0.0));
138 : // Loop over all the selected sides
139 1951 : for (unsigned int i = 0; i < slc_bdry_side_list.size(); i++)
140 : {
141 : // Get the two node ids of the side
142 1920 : dof_id_type side_id_0 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
143 1920 : ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
144 1920 : ->node_ptr(0)
145 1920 : ->id();
146 1920 : dof_id_type side_id_1 = mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i]))
147 1920 : ->side_ptr(std::get<1>(slc_bdry_side_list[i]))
148 1920 : ->node_ptr(1)
149 1920 : ->id();
150 : // True means the selected bdry node is in the node list of the trimming interface
151 : bool side_id_0_in =
152 1920 : !(std::find(node_list.begin(), node_list.end(), side_id_0) == node_list.end());
153 : bool side_id_1_in =
154 1920 : !(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 1920 : bool side_node_0_remove = lineSideDeterminator((*mesh.node_ptr(side_id_0))(0),
158 1920 : (*mesh.node_ptr(side_id_0))(1),
159 1920 : bdry_pars[0],
160 1920 : bdry_pars[1],
161 1920 : bdry_pars[2],
162 : side_to_remove);
163 1920 : bool side_node_1_remove = lineSideDeterminator((*mesh.node_ptr(side_id_1))(0),
164 1920 : (*mesh.node_ptr(side_id_1))(1),
165 1920 : bdry_pars[0],
166 1920 : bdry_pars[1],
167 1920 : bdry_pars[2],
168 : side_to_remove);
169 : // If both nodes of that side are involved in the trimming interface
170 1920 : 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 9 : boundary_info.remove_side(mesh.elem_ptr(std::get<0>(slc_bdry_side_list[i])),
174 9 : std::get<1>(slc_bdry_side_list[i]),
175 9 : 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 1911 : 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 76 : node_list_flag[std::distance(
181 76 : node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_0))] = true;
182 76 : const Point p0 = *mesh.node_ptr(side_id_0);
183 76 : const Point p1 = *mesh.node_ptr(side_id_1);
184 :
185 76 : node_list_point[std::distance(node_list.begin(),
186 76 : std::find(node_list.begin(), node_list.end(), side_id_0))] =
187 76 : twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
188 76 : }
189 : // If node 1 is on the trimming interface, and the side is cut by the trimming line
190 1835 : 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 54 : node_list_flag[std::distance(
194 54 : node_list.begin(), std::find(node_list.begin(), node_list.end(), side_id_1))] = true;
195 54 : const Point p0 = *mesh.node_ptr(side_id_0);
196 54 : const Point p1 = *mesh.node_ptr(side_id_1);
197 :
198 54 : node_list_point[std::distance(node_list.begin(),
199 54 : std::find(node_list.begin(), node_list.end(), side_id_1))] =
200 54 : twoPointandLineIntersection(p0, p1, bdry_pars[0], bdry_pars[1], bdry_pars[2]);
201 : }
202 : }
203 :
204 : // move nodes
205 1400 : 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 1369 : if (node_list_flag[i])
211 130 : *(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 1239 : const Real x0 = (*(mesh.node_ptr(node_list[i])))(0);
217 1239 : const Real y0 = (*(mesh.node_ptr(node_list[i])))(1);
218 1239 : (*(mesh.node_ptr(node_list[i])))(0) =
219 1239 : (bdry_pars[1] * (bdry_pars[1] * x0 - bdry_pars[0] * y0) - bdry_pars[0] * bdry_pars[2]) /
220 1239 : (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
221 1239 : (*(mesh.node_ptr(node_list[i])))(1) =
222 1239 : (bdry_pars[0] * (-bdry_pars[1] * x0 + bdry_pars[0] * y0) - bdry_pars[1] * bdry_pars[2]) /
223 1239 : (bdry_pars[0] * bdry_pars[0] + bdry_pars[1] * bdry_pars[1]);
224 : }
225 : }
226 :
227 : // Delete the block
228 4645 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
229 4645 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
230 : elem_it++)
231 4645 : mesh.delete_elem(*elem_it);
232 31 : mesh.contract();
233 31 : mesh.find_neighbors();
234 : // Delete zero volume elements
235 31 : std::vector<dof_id_type> zero_elems;
236 2953 : for (auto elem_it = mesh.elements_begin(); elem_it != mesh.elements_end(); elem_it++)
237 : {
238 2922 : if (MooseUtils::absoluteFuzzyEqual((*elem_it)->volume(), 0.0))
239 : {
240 90 : for (unsigned int i = 0; i < (*elem_it)->n_sides(); i++)
241 : {
242 72 : if ((*elem_it)->neighbor_ptr(i) != nullptr)
243 : {
244 36 : boundary_info.add_side((*elem_it)->neighbor_ptr(i),
245 36 : ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
246 : external_boundary_id);
247 36 : boundary_info.add_side((*elem_it)->neighbor_ptr(i),
248 36 : ((*elem_it)->neighbor_ptr(i))->which_neighbor_am_i(*elem_it),
249 : trimming_section_boundary_id);
250 : }
251 : }
252 18 : zero_elems.push_back((*elem_it)->id());
253 : }
254 31 : }
255 49 : for (const auto & zero_elem : zero_elems)
256 18 : mesh.delete_elem(mesh.elem_ptr(zero_elem));
257 31 : mesh.contract();
258 : // As we modified the side_list, it is safer to clear the node_list
259 31 : boundary_info.clear_boundary_node_ids();
260 31 : mesh.prepare_for_use();
261 31 : }
262 :
263 : bool
264 96273 : 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 96273 : const Real tmp = px * param_1 + py * param_2 + param_3;
273 96273 : return direction_param ? tmp >= dis_tol : tmp <= dis_tol;
274 : }
275 :
276 : Point
277 2488 : 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 2488 : (param_12 * param_23 - param_22 * param_13) / (param_11 * param_22 - param_21 * param_12),
286 2488 : (param_13 * param_21 - param_23 * param_11) / (param_11 * param_22 - param_21 * param_12),
287 2488 : 0.0);
288 : }
289 :
290 : Point
291 2488 : 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 7464 : return twoLineIntersection(param_1,
298 : param_2,
299 : param_3,
300 2488 : pt2(1) - pt1(1),
301 2488 : pt1(0) - pt2(0),
302 4976 : pt2(0) * pt1(1) - pt1(0) * pt2(1));
303 : }
304 :
305 : bool
306 31 : 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 31 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
312 : // Define the subdomain id shift for the new TRI3 element subdomain(s)
313 31 : const subdomain_id_type max_subdomain_id(*subdomain_ids_set.rbegin());
314 31 : const subdomain_id_type tri_subdomain_id_shift =
315 31 : 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 31 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
321 31 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
322 31 : 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 5839 : 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 2904 : const auto elem_angles = vertex_angles(*elem);
330 2904 : const auto elem_distances = vertex_distances(*elem);
331 : // Type 1
332 2904 : if (MooseUtils::absoluteFuzzyEqual(elem_angles.front().first, M_PI, 0.001))
333 : {
334 231 : bad_elems_rec.push_back(std::make_tuple(elem, elem_angles.front().second, false, true));
335 231 : continue;
336 : }
337 : // Type 2
338 2673 : if (MooseUtils::absoluteFuzzyEqual(elem_distances.front().first, 0.0))
339 : {
340 9 : bad_elems_rec.push_back(std::make_tuple(elem, elem_distances.front().second, false, false));
341 : }
342 3166 : }
343 31 : std::set<subdomain_id_type> new_subdomain_ids;
344 : // Loop over all the identified degenerate QUAD elements
345 271 : for (const auto & bad_elem : bad_elems_rec)
346 : {
347 240 : std::vector<boundary_id_type> elem_bdry_container_0;
348 240 : std::vector<boundary_id_type> elem_bdry_container_1;
349 240 : std::vector<boundary_id_type> elem_bdry_container_2;
350 :
351 240 : Elem * elem_0 = std::get<0>(bad_elem);
352 240 : 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 231 : Elem * elem_1 = elem_0->neighbor_ptr(std::get<1>(bad_elem));
359 231 : Elem * elem_2 = elem_0->neighbor_ptr((std::get<1>(bad_elem) - 1) % elem_0->n_vertices());
360 231 : if ((elem_1 != nullptr || elem_2 != nullptr))
361 0 : throw MooseException("The input mesh has degenerate quad element before trimming.");
362 : }
363 240 : mesh.get_boundary_info().boundary_ids(elem_0, std::get<1>(bad_elem), elem_bdry_container_0);
364 240 : mesh.get_boundary_info().boundary_ids(
365 240 : elem_0, (std::get<1>(bad_elem) + 1) % elem_0->n_vertices(), elem_bdry_container_1);
366 240 : mesh.get_boundary_info().boundary_ids(
367 240 : elem_0, (std::get<1>(bad_elem) + 2) % elem_0->n_vertices(), elem_bdry_container_2);
368 240 : mesh.get_boundary_info().boundary_ids(
369 240 : 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 240 : 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 240 : auto pt0 = elem_0->node_ptr((std::get<1>(bad_elem) + 1) % elem_0->n_vertices());
375 240 : auto pt1 = elem_0->node_ptr((std::get<1>(bad_elem) + 2) % elem_0->n_vertices());
376 240 : 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 240 : 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 240 : mesh.delete_elem(elem_0);
382 : // Create the new TRI element
383 240 : Elem * elem_Tri3 = mesh.add_elem(new Tri3);
384 240 : elem_Tri3->set_node(0, pt0);
385 240 : elem_Tri3->set_node(1, pt1);
386 240 : elem_Tri3->set_node(2, pt2);
387 : // Retain the boundary information
388 480 : for (auto bdry_id : elem_bdry_container_0)
389 240 : boundary_info.add_side(elem_Tri3, 2, bdry_id);
390 285 : for (auto bdry_id : elem_bdry_container_1)
391 45 : boundary_info.add_side(elem_Tri3, 0, bdry_id);
392 316 : for (auto bdry_id : elem_bdry_container_2)
393 76 : boundary_info.add_side(elem_Tri3, 1, bdry_id);
394 : // Assign subdomain id for the TRI element by shifting its original subdomain id
395 240 : elem_Tri3->subdomain_id() = elem_block_id + tri_subdomain_id_shift;
396 240 : new_subdomain_ids.emplace(elem_block_id + tri_subdomain_id_shift);
397 : // Retain element extra integers
398 240 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
399 0 : elem_Tri3->set_extra_integer(j, exist_extra_ids[j]);
400 240 : }
401 : // Assign subdomain names for the new TRI elements
402 130 : for (auto & nid : new_subdomain_ids)
403 : {
404 103 : const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
405 103 : if (MooseMeshUtils::getSubdomainID(
406 206 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
407 206 : : old_name) +
408 412 : "_" + tri_elem_subdomain_name_suffix,
409 103 : mesh) != Moose::INVALID_BLOCK_ID)
410 4 : throw MooseException("The new subdomain name already exists in the mesh.");
411 99 : mesh.subdomain_name(nid) =
412 198 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
413 198 : : old_name) +
414 297 : "_" + tri_elem_subdomain_name_suffix;
415 99 : mooseWarning("Degenerate QUAD elements have been converted into TRI elements with a new "
416 99 : "subdomain name: " +
417 198 : mesh.subdomain_name(nid) + ".");
418 103 : }
419 54 : return bad_elems_rec.size();
420 39 : }
421 :
422 : std::vector<std::pair<Real, unsigned int>>
423 2904 : vertex_angles(const Elem & elem)
424 : {
425 2904 : std::vector<std::pair<Real, unsigned int>> angles;
426 2904 : const unsigned int n_vertices = elem.n_vertices();
427 :
428 14520 : for (unsigned int i = 0; i < n_vertices; i++)
429 : {
430 11616 : Point v1 = (*elem.node_ptr((i - 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
431 11616 : Point v2 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
432 11616 : Real tmp = v1 * v2 / v1.norm() / v2.norm();
433 11616 : if (tmp > 1.0)
434 0 : tmp = 1.0;
435 11616 : else if (tmp < -1.0)
436 49 : tmp = -1.0;
437 11616 : angles.push_back(std::make_pair(acos(tmp), i));
438 : }
439 2904 : std::sort(angles.begin(), angles.end(), std::greater<>());
440 2904 : return angles;
441 0 : }
442 :
443 : std::vector<std::pair<Real, unsigned int>>
444 2904 : vertex_distances(const Elem & elem)
445 : {
446 2904 : std::vector<std::pair<Real, unsigned int>> distances;
447 2904 : const unsigned int n_vertices = elem.n_vertices();
448 :
449 14520 : for (unsigned int i = 0; i < n_vertices; i++)
450 : {
451 11616 : Point v1 = (*elem.node_ptr((i + 1) % n_vertices) - *elem.node_ptr(i % n_vertices));
452 11616 : distances.push_back(std::make_pair(v1.norm(), i));
453 : }
454 2904 : std::sort(distances.begin(), distances.end());
455 2904 : return distances;
456 0 : }
457 :
458 : void
459 900 : 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 900 : const auto elem_old = mesh.elem_ptr(elem_id);
468 900 : const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
469 900 : const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
470 900 : const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
471 :
472 : const bool m1_side_flag =
473 900 : MooseUtils::absoluteFuzzyEqual((*(mesh.node_ptr(nid_3)) - *(mesh.node_ptr(nid_0)))
474 900 : .cross(*(mesh.node_ptr(nid_1)) - *(mesh.node_ptr(nid_0)))
475 900 : .norm(),
476 1800 : 0.0);
477 900 : const dof_id_type nid_m1 = m1_side_flag ? nid_3 : nid_4;
478 900 : const dof_id_type nid_m2 = m1_side_flag ? nid_4 : nid_3;
479 : // Build boundary information of the mesh
480 900 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
481 900 : auto bdry_side_list = boundary_info.build_side_list();
482 : // Create a list of sidesets involving the element to be split
483 900 : std::vector<std::vector<boundary_id_type>> elem_side_list;
484 900 : elem_side_list.resize(3);
485 130140 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
486 : {
487 129240 : if (std::get<0>(bdry_side_list[i]) == elem_id)
488 : {
489 396 : elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
490 198 : std::get<2>(bdry_side_list[i]));
491 : }
492 : }
493 :
494 900 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
495 900 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
496 : // Record all the element extra integers of the original element
497 900 : 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 900 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
501 900 : elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
502 900 : elem_Tri3_0->set_node(1, mesh.node_ptr(nid_m1));
503 900 : elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m2));
504 900 : elem_Tri3_0->subdomain_id() = single_elem_side_id;
505 900 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
506 900 : elem_Tri3_1->set_node(0, mesh.node_ptr(nid_1));
507 900 : elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m2));
508 900 : elem_Tri3_1->set_node(2, mesh.node_ptr(nid_m1));
509 900 : elem_Tri3_1->subdomain_id() = double_elem_side_id;
510 900 : Elem * elem_Tri3_2 = mesh.add_elem(new Tri3);
511 900 : elem_Tri3_2->set_node(0, mesh.node_ptr(nid_2));
512 900 : elem_Tri3_2->set_node(1, mesh.node_ptr(nid_m2));
513 900 : elem_Tri3_2->set_node(2, mesh.node_ptr(nid_1));
514 900 : elem_Tri3_2->subdomain_id() = double_elem_side_id;
515 : // Retain element extra integers
516 900 : 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 954 : for (const auto & side_info_0 : elem_side_list[0])
525 : {
526 54 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
527 54 : boundary_info.add_side(elem_Tri3_1, 2, side_info_0);
528 : }
529 954 : for (const auto & side_info_1 : elem_side_list[1])
530 54 : boundary_info.add_side(elem_Tri3_2, 2, side_info_1);
531 990 : for (const auto & side_info_2 : elem_side_list[2])
532 : {
533 90 : boundary_info.add_side(elem_Tri3_0, 2, side_info_2);
534 90 : boundary_info.add_side(elem_Tri3_2, 0, side_info_2);
535 : }
536 900 : }
537 :
538 : void
539 108 : 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 108 : const auto elem_old = mesh.elem_ptr(elem_id);
547 108 : const dof_id_type nid_0 = elem_old->node_ptr(node_shift % 3)->id();
548 108 : const dof_id_type nid_1 = elem_old->node_ptr((1 + node_shift) % 3)->id();
549 108 : const dof_id_type nid_2 = elem_old->node_ptr((2 + node_shift) % 3)->id();
550 : // Build boundary information of the mesh
551 108 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
552 108 : auto bdry_side_list = boundary_info.build_side_list();
553 : // Create a list of sidesets involving the element to be split
554 108 : std::vector<std::vector<boundary_id_type>> elem_side_list;
555 108 : elem_side_list.resize(3);
556 11880 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
557 : {
558 11772 : if (std::get<0>(bdry_side_list[i]) == elem_id)
559 : {
560 108 : elem_side_list[(std::get<1>(bdry_side_list[i]) + 3 - node_shift) % 3].push_back(
561 54 : std::get<2>(bdry_side_list[i]));
562 : }
563 : }
564 :
565 108 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
566 108 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
567 : // Record all the element extra integers of the original element
568 108 : 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 108 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
572 108 : elem_Tri3_0->set_node(0, mesh.node_ptr(nid_0));
573 108 : elem_Tri3_0->set_node(1, mesh.node_ptr(nid_1));
574 108 : elem_Tri3_0->set_node(2, mesh.node_ptr(nid_m));
575 108 : elem_Tri3_0->subdomain_id() = first_elem_side_id;
576 108 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
577 108 : elem_Tri3_1->set_node(0, mesh.node_ptr(nid_0));
578 108 : elem_Tri3_1->set_node(1, mesh.node_ptr(nid_m));
579 108 : elem_Tri3_1->set_node(2, mesh.node_ptr(nid_2));
580 108 : elem_Tri3_1->subdomain_id() = second_elem_side_id;
581 : // Retain element extra integers
582 108 : 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 126 : for (const auto & side_info_0 : elem_side_list[0])
590 18 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
591 108 : 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 144 : for (const auto & side_info_2 : elem_side_list[2])
597 36 : boundary_info.add_side(elem_Tri3_1, 2, side_info_2);
598 108 : }
599 :
600 : void
601 901 : 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 901 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
607 901 : auto bdry_side_list = boundary_info.build_side_list();
608 : // Create a list of sidesets involving the element to be split
609 901 : std::vector<std::vector<boundary_id_type>> elem_side_list;
610 901 : elem_side_list.resize(4);
611 100453 : for (unsigned int i = 0; i < bdry_side_list.size(); i++)
612 : {
613 99552 : if (std::get<0>(bdry_side_list[i]) == elem_id)
614 : {
615 516 : elem_side_list[std::get<1>(bdry_side_list[i])].push_back(std::get<2>(bdry_side_list[i]));
616 : }
617 : }
618 :
619 901 : auto node_0 = mesh.elem_ptr(elem_id)->node_ptr(0);
620 901 : auto node_1 = mesh.elem_ptr(elem_id)->node_ptr(1);
621 901 : auto node_2 = mesh.elem_ptr(elem_id)->node_ptr(2);
622 901 : auto node_3 = mesh.elem_ptr(elem_id)->node_ptr(3);
623 :
624 901 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
625 901 : 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 901 : 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 901 : if (std::abs((*node_1 - *node_0).cross(*node_3 - *node_0).norm() -
633 901 : (*node_1 - *node_2).cross(*node_3 - *node_2).norm()) >
634 901 : std::abs((*node_0 - *node_1).cross(*node_2 - *node_1).norm() -
635 901 : (*node_0 - *node_3).cross(*node_2 - *node_3).norm()))
636 : {
637 254 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
638 254 : elem_Tri3_0->set_node(0, node_0);
639 254 : elem_Tri3_0->set_node(1, node_1);
640 254 : elem_Tri3_0->set_node(2, node_2);
641 254 : elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
642 254 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
643 254 : elem_Tri3_1->set_node(0, node_0);
644 254 : elem_Tri3_1->set_node(1, node_2);
645 254 : elem_Tri3_1->set_node(2, node_3);
646 254 : elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
647 : // Retain element extra integers
648 254 : 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 254 : for (const auto & side_info_0 : elem_side_list[0])
656 0 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
657 384 : for (const auto & side_info_1 : elem_side_list[1])
658 130 : boundary_info.add_side(elem_Tri3_0, 1, side_info_1);
659 254 : for (const auto & side_info_2 : elem_side_list[2])
660 0 : boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
661 254 : 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 647 : Elem * elem_Tri3_0 = mesh.add_elem(new Tri3);
667 647 : elem_Tri3_0->set_node(0, node_0);
668 647 : elem_Tri3_0->set_node(1, node_1);
669 647 : elem_Tri3_0->set_node(2, node_3);
670 647 : elem_Tri3_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
671 647 : Elem * elem_Tri3_1 = mesh.add_elem(new Tri3);
672 647 : elem_Tri3_1->set_node(0, node_1);
673 647 : elem_Tri3_1->set_node(1, node_2);
674 647 : elem_Tri3_1->set_node(2, node_3);
675 647 : elem_Tri3_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + tri_elem_subdomain_shift;
676 : // Retain element extra integers
677 647 : 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 683 : for (const auto & side_info_0 : elem_side_list[0])
685 36 : boundary_info.add_side(elem_Tri3_0, 0, side_info_0);
686 845 : for (const auto & side_info_1 : elem_side_list[1])
687 198 : boundary_info.add_side(elem_Tri3_1, 0, side_info_1);
688 709 : for (const auto & side_info_2 : elem_side_list[2])
689 62 : boundary_info.add_side(elem_Tri3_1, 1, side_info_2);
690 737 : for (const auto & side_info_3 : elem_side_list[3])
691 90 : boundary_info.add_side(elem_Tri3_0, 2, side_info_3);
692 : }
693 901 : }
694 :
695 : void
696 49 : 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 49 : std::vector<dof_id_type> cross_elems_quad;
703 49 : std::set<subdomain_id_type> new_subdomain_ids;
704 7873 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
705 : elem_it++)
706 : {
707 7824 : if ((*elem_it)->n_vertices() == 4)
708 : {
709 7824 : std::vector<unsigned short> node_side_rec;
710 39120 : for (unsigned int i = 0; i < 4; i++)
711 : {
712 31296 : const Point v_point = (*elem_it)->point(i);
713 31296 : node_side_rec.push_back(lineSideDeterminator(v_point(0),
714 31296 : v_point(1),
715 31296 : cut_line_params[0],
716 31296 : cut_line_params[1],
717 31296 : cut_line_params[2],
718 : true));
719 : }
720 11331 : if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) != 4 &&
721 3507 : std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) > 0)
722 : {
723 901 : cross_elems_quad.push_back((*elem_it)->id());
724 901 : new_subdomain_ids.emplace((*elem_it)->subdomain_id() + tri_subdomain_id_shift);
725 : }
726 7824 : }
727 49 : }
728 : // Then convert these quad elements into tri elements
729 950 : for (const auto & cross_elem_quad : cross_elems_quad)
730 : {
731 901 : quadElemSplitter(mesh, cross_elem_quad, tri_subdomain_id_shift);
732 901 : mesh.delete_elem(mesh.elem_ptr(cross_elem_quad));
733 : }
734 166 : for (auto & nid : new_subdomain_ids)
735 : {
736 121 : const SubdomainName old_name = mesh.subdomain_name(nid - tri_subdomain_id_shift);
737 121 : if (MooseMeshUtils::getSubdomainID(
738 242 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
739 242 : : old_name) +
740 484 : "_" + tri_elem_subdomain_name_suffix,
741 121 : mesh) != Moose::INVALID_BLOCK_ID)
742 4 : throw MooseException("The new subdomain name already exists in the mesh.");
743 117 : mesh.subdomain_name(nid) =
744 234 : (old_name.empty() ? (SubdomainName)(std::to_string(nid - tri_subdomain_id_shift))
745 234 : : old_name) +
746 351 : "_" + tri_elem_subdomain_name_suffix;
747 117 : mooseWarning("QUAD elements have been converted into TRI elements with a new "
748 117 : "subdomain name: " +
749 234 : mesh.subdomain_name(nid) + ".");
750 121 : }
751 45 : mesh.contract();
752 53 : }
753 :
754 : void
755 45 : 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 45 : std::vector<dof_id_type> cross_elems;
762 : // A vector for element specific information
763 45 : std::vector<std::vector<std::pair<dof_id_type, dof_id_type>>> node_pairs_vec;
764 : // A set for unique pairs
765 45 : std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs_unique_vec;
766 7290 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
767 : elem_it++)
768 : {
769 7245 : std::vector<unsigned short> node_side_rec;
770 7245 : const auto n_vertices = (*elem_it)->n_vertices();
771 7245 : node_side_rec.resize(n_vertices);
772 34695 : for (unsigned int i = 0; i < n_vertices; i++)
773 : {
774 : // First check if the vertex is in the XY Plane
775 27450 : 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 27450 : const Point v_point = (*elem_it)->point(i);
779 27450 : node_side_rec[i] = lineSideDeterminator(
780 27450 : v_point(0), v_point(1), cut_line_params[0], cut_line_params[1], cut_line_params[2], true);
781 : }
782 7245 : if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) == (int)node_side_rec.size())
783 : {
784 3852 : (*elem_it)->subdomain_id() = block_id_to_remove;
785 : }
786 3393 : else if (std::accumulate(node_side_rec.begin(), node_side_rec.end(), 0) > 0)
787 : {
788 1269 : 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 1269 : cross_elems.push_back((*elem_it)->id());
791 : // Then we need to check pairs of nodes that are on the different side
792 1269 : std::vector<std::pair<dof_id_type, dof_id_type>> node_pairs;
793 5076 : 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 3807 : if (node_side_rec[i] > 0 && node_side_rec[(i + 1) % node_side_rec.size()] == 0)
797 : {
798 : // Removal side first
799 1269 : node_pairs.push_back(
800 1269 : std::make_pair((*elem_it)->node_ptr(i)->id(),
801 1269 : (*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id()));
802 1269 : node_pairs_unique_vec.push_back(node_pairs.back());
803 : }
804 : // first node on retaining side and second node on removal side
805 2538 : else if (node_side_rec[i] == 0 && node_side_rec[(i + 1) % node_side_rec.size()] > 0)
806 : {
807 : // Removal side first
808 1269 : node_pairs.push_back(
809 1269 : std::make_pair((*elem_it)->node_ptr((i + 1) % node_side_rec.size())->id(),
810 1269 : (*elem_it)->node_ptr(i)->id()));
811 1269 : node_pairs_unique_vec.push_back(node_pairs.back());
812 : }
813 : }
814 1269 : node_pairs_vec.push_back(node_pairs);
815 1269 : }
816 7290 : }
817 45 : auto vec_ip = std::unique(node_pairs_unique_vec.begin(), node_pairs_unique_vec.end());
818 45 : 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 45 : std::vector<Node *> nodes_on_line;
822 : // whether the on-line node is overlapped with the node pairs or a brand new node
823 45 : std::vector<unsigned short> nodes_on_line_overlap;
824 2403 : for (const auto & node_pair : node_pairs_unique_vec)
825 : {
826 2358 : const Point pt1 = *mesh.node_ptr(node_pair.first);
827 2358 : const Point pt2 = *mesh.node_ptr(node_pair.second);
828 2358 : const Point pt_line = twoPointandLineIntersection(
829 2358 : pt1, pt2, cut_line_params[0], cut_line_params[1], cut_line_params[2]);
830 2358 : 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 2358 : else if ((pt_line - pt2).norm() < libMesh::TOLERANCE)
836 : {
837 612 : nodes_on_line.push_back(mesh.node_ptr(node_pair.second));
838 612 : nodes_on_line_overlap.push_back(2);
839 : }
840 : else
841 : {
842 1746 : nodes_on_line.push_back(mesh.add_point(pt_line));
843 1746 : nodes_on_line_overlap.push_back(0);
844 : }
845 : }
846 :
847 : // make new elements
848 1314 : for (unsigned int i = 0; i < cross_elems.size(); i++)
849 : {
850 : // Only TRI elements are involved after preprocessing
851 1269 : auto cross_elem = mesh.elem_ptr(cross_elems[i]);
852 1269 : auto node_0 = cross_elem->node_ptr(0);
853 1269 : auto node_1 = cross_elem->node_ptr(1);
854 1269 : auto node_2 = cross_elem->node_ptr(2);
855 1269 : const std::vector<dof_id_type> tri_nodes = {node_0->id(), node_1->id(), node_2->id()};
856 :
857 1269 : 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 1269 : node_pairs_vec[i][0]));
861 1269 : 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 1269 : node_pairs_vec[i][1]));
865 1269 : auto node_3 = nodes_on_line[online_node_index_1];
866 1269 : auto node_4 = nodes_on_line[online_node_index_2];
867 1269 : const auto node_3_overlap_flag = nodes_on_line_overlap[online_node_index_1];
868 1269 : const auto node_4_overlap_flag = nodes_on_line_overlap[online_node_index_2];
869 : // Most common case, no overlapped nodes
870 1269 : 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 900 : 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 900 : common_node_side ? block_id_to_remove : cross_elem->subdomain_id();
876 : const subdomain_id_type block_id_to_assign_2 =
877 900 : 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 900 : common_node_side ? node_pairs_vec[i][0].first : node_pairs_vec[i][0].second;
882 :
883 1800 : triElemSplitter(mesh,
884 : cross_elem->id(),
885 900 : 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 900 : mesh.delete_elem(cross_elem);
892 900 : }
893 : // both node_3 and node_4 are overlapped
894 369 : 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 522 : cross_elem->subdomain_id() = lineSideDeterminator(cross_elem->vertex_average()(0),
899 522 : cross_elem->vertex_average()(1),
900 261 : cut_line_params[0],
901 261 : cut_line_params[1],
902 261 : cut_line_params[2],
903 : true)
904 261 : ? block_id_to_remove
905 0 : : cross_elem->subdomain_id();
906 : }
907 : // node_3 or node_4 is overlapped
908 : else
909 : {
910 108 : const auto node_3_finder = std::distance(
911 108 : tri_nodes.begin(), std::find(tri_nodes.begin(), tri_nodes.end(), node_3->id()));
912 108 : const auto node_4_finder = std::distance(
913 108 : 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 108 : const dof_id_type node_id = node_3_finder < node_4_finder ? node_4->id() : node_3->id();
917 108 : const auto node_finder = std::min(node_3_finder, node_4_finder);
918 :
919 162 : triElemSplitter(
920 : mesh,
921 : cross_elem->id(),
922 : node_finder,
923 : node_id,
924 108 : tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
925 54 : ? block_id_to_remove
926 54 : : cross_elem->subdomain_id(),
927 108 : tri_nodes[(node_finder + 1) % 3] == node_pairs_vec[i][node_3_finder > node_4_finder].first
928 54 : ? cross_elem->subdomain_id()
929 : : block_id_to_remove);
930 108 : mesh.delete_elem(cross_elem);
931 : }
932 1269 : }
933 45 : 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 45 : mesh.find_neighbors();
938 45 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
939 9198 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
940 : elem_it++)
941 : {
942 9153 : if ((*elem_it)->subdomain_id() != block_id_to_remove)
943 : {
944 16254 : for (unsigned int j = 0; j < (*elem_it)->n_sides(); j++)
945 : {
946 12690 : if ((*elem_it)->neighbor_ptr(j) != nullptr)
947 12042 : if ((*elem_it)->neighbor_ptr(j)->subdomain_id() == block_id_to_remove)
948 1080 : boundary_info.add_side(*elem_it, j, new_boundary_id);
949 : }
950 : }
951 45 : }
952 :
953 : // Delete the block to remove
954 5634 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
955 5634 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
956 : elem_it++)
957 5634 : mesh.delete_elem(*elem_it);
958 45 : mesh.contract();
959 45 : }
960 :
961 : void
962 49 : 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 49 : 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 45 : lineRemoverCutElemTri(mesh, cut_line_params, block_id_to_remove, new_boundary_id);
975 :
976 45 : if (improve_boundary_tri_elems)
977 9 : boundaryTriElemImprover(mesh, new_boundary_id);
978 45 : }
979 :
980 : void
981 9 : boundaryTriElemImprover(ReplicatedMesh & mesh, const boundary_id_type boundary_to_improve)
982 : {
983 9 : 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 9 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
988 9 : 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 9 : tri3_elem_info;
997 1017 : for (const auto & side : side_list)
998 : {
999 1008 : if (std::get<2>(side) == boundary_to_improve)
1000 : {
1001 468 : Elem * elem = mesh.elem_ptr(std::get<0>(side));
1002 468 : if (elem->type() == TRI3)
1003 : {
1004 468 : const auto key_node_id = elem->node_id((std::get<1>(side) + 2) % 3);
1005 468 : const auto value_elem_id = elem->id();
1006 468 : const auto value_node_id_1 = elem->node_id(std::get<1>(side));
1007 468 : const auto value_node_id_2 = elem->node_id((std::get<1>(side) + 1) % 3);
1008 936 : tri3_elem_info[key_node_id].push_back(
1009 936 : 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 9 : 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 243 : 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 234 : std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
1021 234 : std::vector<dof_id_type> elem_id_list;
1022 702 : for (const auto & tri : tri_group.second)
1023 : {
1024 468 : node_assm.push_back(std::make_pair(std::get<1>(tri), std::get<2>(tri)));
1025 468 : elem_id_list.push_back(std::get<0>(tri));
1026 : }
1027 234 : std::vector<dof_id_type> ordered_node_list;
1028 234 : std::vector<dof_id_type> ordered_elem_list;
1029 234 : 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 234 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
1037 234 : std::vector<std::tuple<subdomain_id_type, std::vector<dof_id_type>, unsigned int>> blocks_info;
1038 702 : for (const auto & elem_id : ordered_elem_list)
1039 : {
1040 468 : 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 468 : 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 468 : if (!blocks_info.empty())
1045 : {
1046 450 : if (mesh.elem_ptr(elem_id)->subdomain_id() == std::get<0>(blocks_info.back()) &&
1047 216 : exist_extra_ids == std::get<1>(blocks_info.back()))
1048 : {
1049 216 : std::get<2>(blocks_info.back())++;
1050 216 : continue;
1051 : }
1052 : }
1053 252 : blocks_info.push_back(
1054 504 : std::make_tuple(mesh.elem_ptr(elem_id)->subdomain_id(), exist_extra_ids, 1));
1055 468 : }
1056 : // For each separated subdomain / set of extra ids, we try to improve the boundary elements
1057 234 : unsigned int side_counter = 0;
1058 486 : for (const auto & block_info : blocks_info)
1059 : {
1060 252 : const auto node_1 = mesh.node_ptr(ordered_node_list[side_counter]);
1061 : // we do not need to subtract 1 for node_2
1062 252 : const auto node_2 = mesh.node_ptr(ordered_node_list[side_counter + std::get<2>(block_info)]);
1063 252 : const auto node_0 = mesh.node_ptr(tri_group.first);
1064 252 : const Point v1 = *node_1 - *node_0;
1065 252 : const Point v2 = *node_2 - *node_0;
1066 252 : 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 252 : ordered_elem_list.begin() + side_counter +
1069 504 : 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 252 : elemSideLocator(mesh,
1078 252 : block_elems.front(),
1079 252 : tri_group.first,
1080 252 : ordered_node_list[side_counter],
1081 : side_id_0,
1082 : is_inverse_0);
1083 252 : elemSideLocator(mesh,
1084 252 : block_elems.back(),
1085 252 : ordered_node_list[side_counter + std::get<2>(block_info)],
1086 252 : tri_group.first,
1087 : side_id_t,
1088 : is_inverse_t);
1089 : // Collect boundary information of the identified sides
1090 252 : std::vector<boundary_id_type> side_0_boundary_ids;
1091 504 : boundary_info.boundary_ids(
1092 252 : mesh.elem_ptr(block_elems.front()), side_id_0, side_0_boundary_ids);
1093 252 : std::vector<boundary_id_type> side_t_boundary_ids;
1094 252 : 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 252 : if (angle < 90.0)
1102 : {
1103 171 : if (std::get<2>(block_info) > 1)
1104 : {
1105 144 : makeImprovedTriElement(mesh,
1106 72 : tri_group.first,
1107 72 : ordered_node_list[side_counter],
1108 72 : ordered_node_list[side_counter + std::get<2>(block_info)],
1109 72 : std::get<0>(block_info),
1110 72 : std::get<1>(block_info),
1111 : {boundary_to_improve},
1112 : side_0_boundary_ids,
1113 : side_t_boundary_ids);
1114 72 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1115 : }
1116 : }
1117 81 : else if (angle < 135.0)
1118 : {
1119 : // We can just add the middle node because there's nothing on the other side
1120 63 : const auto node_m = mesh.add_point((*node_1 + *node_2) / 2.0);
1121 126 : makeImprovedTriElement(mesh,
1122 63 : tri_group.first,
1123 63 : ordered_node_list[side_counter],
1124 : node_m->id(),
1125 63 : std::get<0>(block_info),
1126 63 : std::get<1>(block_info),
1127 : {boundary_to_improve},
1128 : side_0_boundary_ids,
1129 126 : std::vector<boundary_id_type>());
1130 189 : makeImprovedTriElement(mesh,
1131 63 : tri_group.first,
1132 : node_m->id(),
1133 63 : ordered_node_list[side_counter + std::get<2>(block_info)],
1134 63 : std::get<0>(block_info),
1135 63 : std::get<1>(block_info),
1136 : {boundary_to_improve},
1137 126 : std::vector<boundary_id_type>(),
1138 : side_t_boundary_ids);
1139 63 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1140 : }
1141 : else
1142 : {
1143 18 : const auto node_m1 = mesh.add_point((*node_1 * 2.0 + *node_2) / 3.0);
1144 18 : const auto node_m2 = mesh.add_point((*node_1 + *node_2 * 2.0) / 3.0);
1145 36 : makeImprovedTriElement(mesh,
1146 18 : tri_group.first,
1147 18 : ordered_node_list[side_counter],
1148 : node_m1->id(),
1149 18 : std::get<0>(block_info),
1150 18 : std::get<1>(block_info),
1151 : {boundary_to_improve},
1152 : side_0_boundary_ids,
1153 36 : std::vector<boundary_id_type>());
1154 54 : makeImprovedTriElement(mesh,
1155 18 : tri_group.first,
1156 : node_m1->id(),
1157 : node_m2->id(),
1158 18 : std::get<0>(block_info),
1159 18 : std::get<1>(block_info),
1160 : {boundary_to_improve},
1161 36 : std::vector<boundary_id_type>(),
1162 36 : std::vector<boundary_id_type>());
1163 54 : makeImprovedTriElement(mesh,
1164 18 : tri_group.first,
1165 : node_m2->id(),
1166 18 : ordered_node_list[side_counter + std::get<2>(block_info)],
1167 18 : std::get<0>(block_info),
1168 18 : std::get<1>(block_info),
1169 : {boundary_to_improve},
1170 36 : std::vector<boundary_id_type>(),
1171 : side_t_boundary_ids);
1172 18 : elems_to_remove.insert(elems_to_remove.end(), block_elems.begin(), block_elems.end());
1173 : }
1174 252 : side_counter += std::get<2>(block_info);
1175 252 : }
1176 : // TODO: Need to check if the new element is inverted?
1177 234 : }
1178 : // Delete the original elements
1179 378 : for (const auto & elem_to_remove : elems_to_remove)
1180 369 : mesh.delete_elem(mesh.elem_ptr(elem_to_remove));
1181 9 : mesh.contract();
1182 9 : }
1183 :
1184 : void
1185 252 : 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 252 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1196 252 : Elem * elem_Tri3_new = mesh.add_elem(new Tri3);
1197 252 : elem_Tri3_new->set_node(0, mesh.node_ptr(node_id_0));
1198 252 : elem_Tri3_new->set_node(1, mesh.node_ptr(node_id_1));
1199 252 : elem_Tri3_new->set_node(2, mesh.node_ptr(node_id_2));
1200 297 : for (const auto & boundary_id_for_side_0 : boundary_ids_for_side_0)
1201 45 : boundary_info.add_side(elem_Tri3_new, 0, boundary_id_for_side_0);
1202 504 : for (const auto & boundary_id_for_side_1 : boundary_ids_for_side_1)
1203 252 : boundary_info.add_side(elem_Tri3_new, 1, boundary_id_for_side_1);
1204 252 : 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 252 : elem_Tri3_new->subdomain_id() = subdomain_id;
1207 : // Retain element extra integers
1208 252 : 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 252 : }
1213 :
1214 : bool
1215 504 : 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 504 : Elem * elem = mesh.elem_ptr(elem_id);
1223 1008 : for (unsigned short i = 0; i < elem->n_sides(); i++)
1224 : {
1225 2349 : if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_0 &&
1226 1341 : elem->side_ptr(i)->node_ptr(1)->id() == node_id_1)
1227 : {
1228 162 : side_id = i;
1229 162 : is_inverse = false;
1230 162 : return true;
1231 : }
1232 2115 : else if (elem->side_ptr(i)->node_ptr(0)->id() == node_id_1 &&
1233 1269 : elem->side_ptr(i)->node_ptr(1)->id() == node_id_0)
1234 : {
1235 342 : side_id = i;
1236 342 : is_inverse = true;
1237 342 : return true;
1238 : }
1239 : }
1240 0 : return false;
1241 : }
1242 : }
|
