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 : #include "CutMeshByLevelSetGeneratorBase.h"
11 : #include "MooseMeshElementConversionUtils.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/cell_tet4.h"
20 :
21 : // C++ includes
22 : #include <cmath>
23 :
24 : InputParameters
25 28776 : CutMeshByLevelSetGeneratorBase::validParams()
26 : {
27 28776 : InputParameters params = MeshGenerator::validParams();
28 28776 : params += FunctionParserUtils<false>::validParams();
29 :
30 115104 : params.addRequiredParam<MeshGeneratorName>("input", "The input mesh that needs to be trimmed.");
31 :
32 86328 : params.addParam<bool>(
33 : "generate_transition_layer",
34 57552 : false,
35 : "Whether to generate a transition layer for the cut mesh. "
36 : "If false, the entire input mesh will be converted to TET4 elements to facilitate the "
37 : "cutting; if true, only the elements near the cut face will be converted with a transition "
38 : "layer to maintain compatibility with the original mesh.");
39 :
40 115104 : params.addParam<boundary_id_type>("cut_face_id",
41 : "The boundary id of the face generated by the cut. An "
42 : "id will be automatically assigned if not provided.");
43 86328 : params.addParam<BoundaryName>(
44 57552 : "cut_face_name", BoundaryName(), "The boundary name of the face generated by the cut.");
45 115104 : params.addParam<SubdomainName>("converted_tet_element_subdomain_name_suffix",
46 : "to_tet",
47 : "The suffix to be added to the original subdomain name for the "
48 : "subdomains containing the elements converted to TET4. This is "
49 : "only applicable when transition layer is generated.");
50 115104 : params.addParam<SubdomainName>(
51 : "converted_pyramid_element_subdomain_name_suffix",
52 : "to_pyramid",
53 : "The suffix to be added to the original subdomain name for the subdomains containing the "
54 : "elements converted to PYRAMID5. This is only applicable when transition layer is "
55 : "generated.");
56 :
57 28776 : params.addClassDescription(
58 : "This CutMeshByLevelSetGeneratorBase object is designed to be the base class of mesh "
59 : "generator that cuts a 3D mesh based on an analytic level set function. The level set "
60 : "function could be provided explicitly or indirectly.");
61 :
62 28776 : return params;
63 0 : }
64 :
65 123 : CutMeshByLevelSetGeneratorBase::CutMeshByLevelSetGeneratorBase(const InputParameters & parameters)
66 : : MeshGenerator(parameters),
67 : FunctionParserUtils<false>(parameters),
68 123 : _input_name(getParam<MeshGeneratorName>("input")),
69 246 : _generate_transition_layer(getParam<bool>("generate_transition_layer")),
70 246 : _cut_face_name(getParam<BoundaryName>("cut_face_name")),
71 246 : _converted_tet_element_subdomain_name_suffix(
72 : getParam<SubdomainName>("converted_tet_element_subdomain_name_suffix")),
73 246 : _converted_pyramid_element_subdomain_name_suffix(
74 : getParam<SubdomainName>("converted_pyramid_element_subdomain_name_suffix")),
75 246 : _input(getMeshByName(_input_name))
76 : {
77 370 : _cut_face_id = isParamValid("cut_face_id") ? getParam<boundary_id_type>("cut_face_id") : -1;
78 123 : }
79 :
80 : std::unique_ptr<MeshBase>
81 123 : CutMeshByLevelSetGeneratorBase::generate()
82 : {
83 123 : auto replicated_mesh_ptr = dynamic_cast<ReplicatedMesh *>(_input.get());
84 123 : if (!replicated_mesh_ptr)
85 8 : paramError("input", "Input is not a replicated mesh, which is required");
86 234 : if (*(replicated_mesh_ptr->elem_dimensions().begin()) != 3 ||
87 234 : *(replicated_mesh_ptr->elem_dimensions().rbegin()) != 3)
88 12 : paramError(
89 : "input",
90 : "Only 3D meshes are supported. Use XYMeshLineCutter for 2D meshes if applicable. Mixed "
91 : "dimensional meshes are not supported at the moment.");
92 :
93 115 : ReplicatedMesh & mesh = *replicated_mesh_ptr;
94 :
95 115 : if (!_cut_face_name.empty())
96 : {
97 58 : if (MooseMeshUtils::hasBoundaryName(mesh, _cut_face_name))
98 : {
99 : const boundary_id_type exist_cut_face_id =
100 13 : MooseMeshUtils::getBoundaryID(_cut_face_name, mesh);
101 13 : if (_cut_face_id != -1 && _cut_face_id != exist_cut_face_id)
102 8 : paramError("cut_face_id",
103 : "The cut face boundary name and id are both provided, but they are inconsistent "
104 : "with an existing boundary name which has a different id.");
105 : else
106 9 : _cut_face_id = exist_cut_face_id;
107 : }
108 : else
109 : {
110 45 : if (_cut_face_id == -1)
111 9 : _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh);
112 45 : mesh.get_boundary_info().sideset_name(_cut_face_id) = _cut_face_name;
113 : }
114 : }
115 57 : else if (_cut_face_id == -1)
116 : {
117 48 : _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh);
118 : }
119 :
120 : // Subdomain ID for new utility blocks must be new
121 : // Find sid shifts
122 111 : const auto sid_shift_base = MooseMeshUtils::getNextFreeSubdomainID(mesh);
123 111 : const auto block_id_to_remove = sid_shift_base * 3;
124 :
125 : // For the boolean value in the pair, true means the element is crossed by the cutting plane
126 : // false means the element is on the remaining side
127 111 : std::vector<std::pair<dof_id_type, bool>> cross_and_remained_elems_pre_convert;
128 :
129 : // collect the subdomain ids of the original mesh
130 111 : std::set<subdomain_id_type> original_subdomain_ids;
131 16391 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
132 : elem_it++)
133 : {
134 16284 : original_subdomain_ids.emplace((*elem_it)->subdomain_id());
135 16284 : const unsigned int & n_vertices = (*elem_it)->n_vertices();
136 16284 : unsigned short elem_vertices_counter(0);
137 137412 : for (unsigned int i = 0; i < n_vertices; i++)
138 : {
139 : // We define elem_vertices_counter in this way so that those elements with one face on the
140 : // plane are also processed to have the cut face boundary assigned.
141 121128 : if (pointLevelSetRelation((*(*elem_it)->node_ptr(i))) !=
142 : PointLevelSetRelationIndex::level_set_in_side)
143 73184 : elem_vertices_counter++;
144 : }
145 16284 : if (elem_vertices_counter == n_vertices)
146 5816 : (*elem_it)->subdomain_id() = block_id_to_remove;
147 : else
148 : {
149 : // Check if any elements to be processed are not first order
150 10468 : if ((*elem_it)->default_order() != Order::FIRST)
151 4 : mooseError("Only first order elements are supported for cutting.");
152 : // If we are generating a transition layer, we only need to record the crossed elements
153 10464 : if (_generate_transition_layer && elem_vertices_counter > 0)
154 : {
155 1857 : cross_and_remained_elems_pre_convert.push_back(std::make_pair((*elem_it)->id(), true));
156 : // Since we are converting these elements to TET4, and we would keep the original type of
157 : // some elements in these blocks
158 3714 : if ((*elem_it)->type() != TET4)
159 1857 : (*elem_it)->subdomain_id() += sid_shift_base;
160 : }
161 8607 : else if (!_generate_transition_layer)
162 7776 : cross_and_remained_elems_pre_convert.push_back(
163 15552 : std::make_pair((*elem_it)->id(), elem_vertices_counter > 0));
164 : }
165 107 : }
166 :
167 : // Then we need to identify and make the transition layer
168 107 : std::vector<dof_id_type> transition_elems_list;
169 107 : std::vector<std::vector<unsigned int>> transition_elems_sides_list;
170 107 : if (_generate_transition_layer)
171 : {
172 35 : if (!mesh.is_prepared())
173 35 : mesh.find_neighbors();
174 : // First, we need to identify the retained elements that share the boundary with the crossed
175 : // elements.
176 1892 : for (const auto & elem_info : cross_and_remained_elems_pre_convert)
177 : {
178 12999 : for (const auto & i_side : make_range(mesh.elem_ptr(elem_info.first)->n_sides()))
179 : {
180 : // No need to work on a TRI3 side
181 11142 : if (mesh.elem_ptr(elem_info.first)->side_ptr(i_side)->type() == QUAD4)
182 : {
183 11142 : if (mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side) != nullptr &&
184 9386 : mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side)->subdomain_id() !=
185 20528 : block_id_to_remove &&
186 8004 : std::find(cross_and_remained_elems_pre_convert.begin(),
187 : cross_and_remained_elems_pre_convert.end(),
188 8004 : std::make_pair(mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side)->id(),
189 27150 : true)) == cross_and_remained_elems_pre_convert.end())
190 : {
191 : const auto & neighbor_elem_id =
192 1166 : mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side)->id();
193 : const auto & neighbor_elem_side_index =
194 1166 : mesh.elem_ptr(elem_info.first)
195 : ->neighbor_ptr(i_side)
196 1166 : ->which_neighbor_am_i(mesh.elem_ptr(elem_info.first));
197 1166 : const auto & id_found = std::find(
198 1166 : transition_elems_list.begin(), transition_elems_list.end(), neighbor_elem_id);
199 1166 : if (id_found == transition_elems_list.end())
200 : {
201 530 : transition_elems_list.push_back(neighbor_elem_id);
202 530 : transition_elems_sides_list.push_back(
203 1060 : std::vector<unsigned int>({neighbor_elem_side_index}));
204 : }
205 : else
206 : {
207 : // If the element is already in the list, we just add the side index
208 636 : const auto index = std::distance(transition_elems_list.begin(), id_found);
209 636 : transition_elems_sides_list[index].push_back(neighbor_elem_side_index);
210 : }
211 : }
212 : }
213 : }
214 : }
215 : }
216 :
217 107 : std::vector<dof_id_type> converted_elems_ids_to_cut;
218 : // Then convert these non-TET4 elements into TET4 elements
219 107 : MooseMeshElementConversionUtils::convert3DMeshToAllTet4(mesh,
220 : cross_and_remained_elems_pre_convert,
221 : converted_elems_ids_to_cut,
222 : block_id_to_remove,
223 : false);
224 :
225 : // Make the transition layer if applicable
226 107 : auto & sideset_map = mesh.get_boundary_info().get_sideset_map();
227 637 : for (const auto & i_elem : index_range(transition_elems_list))
228 : {
229 530 : const auto & elem_id = transition_elems_list[i_elem];
230 530 : const auto & side_indices = transition_elems_sides_list[i_elem];
231 : // Find the involved sidesets of the element so that we can retain them
232 530 : std::vector<std::vector<boundary_id_type>> elem_side_info(mesh.elem_ptr(elem_id)->n_sides());
233 530 : auto side_range = sideset_map.equal_range(mesh.elem_ptr(elem_id));
234 929 : for (auto i = side_range.first; i != side_range.second; ++i)
235 399 : elem_side_info[i->second.first].push_back(i->second.second);
236 :
237 530 : MooseMeshElementConversionUtils::convertElem(
238 : mesh, elem_id, side_indices, elem_side_info, sid_shift_base);
239 530 : }
240 :
241 107 : std::vector<const Node *> new_on_plane_nodes;
242 : // We build the sideset information now, as we only need the information of the elements before
243 : // cutting
244 107 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
245 107 : const auto bdry_side_list = boundary_info.build_side_list();
246 : // Cut the TET4 Elements
247 40229 : for (const auto & converted_elems_id_to_cut : converted_elems_ids_to_cut)
248 : {
249 40122 : tet4ElemCutter(
250 : mesh, bdry_side_list, converted_elems_id_to_cut, block_id_to_remove, new_on_plane_nodes);
251 : }
252 :
253 : // If a transition layer is generated, we need to add some subdomain names
254 107 : if (_generate_transition_layer)
255 : {
256 : try
257 : {
258 35 : MooseMeshElementConversionUtils::assignConvertedElementsSubdomainNameSuffix(
259 : mesh,
260 : original_subdomain_ids,
261 : sid_shift_base,
262 35 : _converted_tet_element_subdomain_name_suffix,
263 35 : _converted_pyramid_element_subdomain_name_suffix);
264 : }
265 8 : catch (const std::exception & e)
266 : {
267 16 : if (((std::string)e.what()).compare(26, 4, "TET4") == 0)
268 8 : paramError("converted_tet_element_subdomain_name_suffix", e.what());
269 : else // if (((std::string)e.what()).compare(26, 7, "PYRAMID5") == 0)
270 8 : paramError("converted_pyramid_element_subdomain_name_suffix", e.what());
271 0 : }
272 : }
273 :
274 : // delete the original elements that were converted
275 549 : for (const auto & elem_id : transition_elems_list)
276 450 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
277 : // Delete the block to remove
278 48852 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
279 48852 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
280 : elem_it++)
281 48852 : mesh.delete_elem(*elem_it);
282 :
283 99 : mesh.contract();
284 99 : mesh.set_isnt_prepared();
285 198 : return std::move(_input);
286 99 : }
287 :
288 : CutMeshByLevelSetGeneratorBase::PointLevelSetRelationIndex
289 281616 : CutMeshByLevelSetGeneratorBase::pointLevelSetRelation(const Point & point)
290 : {
291 281616 : const Real level_set_value = levelSetEvaluator(point);
292 281616 : if (MooseUtils::absoluteFuzzyLessThan(level_set_value, 0.0))
293 117520 : return PointLevelSetRelationIndex::level_set_in_side;
294 164096 : else if (MooseUtils::absoluteFuzzyGreaterThan(level_set_value, 0.0))
295 156342 : return PointLevelSetRelationIndex::level_set_out_side;
296 : else
297 7754 : return PointLevelSetRelationIndex::on_level_set;
298 : }
299 :
300 : Point
301 95001 : CutMeshByLevelSetGeneratorBase::pointPairLevelSetInterception(const Point & point1,
302 : const Point & point2)
303 : {
304 95001 : Real dist1 = levelSetEvaluator(point1);
305 95001 : Real dist2 = levelSetEvaluator(point2);
306 :
307 95001 : if (MooseUtils::absoluteFuzzyEqual(dist1, 0.0) || MooseUtils::absoluteFuzzyEqual(dist2, 0.0))
308 0 : mooseError("At least one of the two points are on the plane.");
309 95001 : if (MooseUtils::absoluteFuzzyGreaterThan(dist1 * dist2, 0.0))
310 0 : mooseError("The two points are on the same side of the plane.");
311 :
312 95001 : Point p1(point1);
313 95001 : Point p2(point2);
314 95001 : Real dist = abs(dist1) + abs(dist2);
315 95001 : Point mid_point;
316 :
317 : // Bisection method to find midpoint
318 3704444 : while (MooseUtils::absoluteFuzzyGreaterThan(dist, 0.0))
319 : {
320 3609443 : mid_point = 0.5 * (p1 + p2);
321 3609443 : const Real dist_mid = levelSetEvaluator(mid_point);
322 : // We do not need Fuzzy here as it will be covered by the while loop
323 3609443 : if (dist_mid == 0.0)
324 2232 : dist = 0.0;
325 3607211 : else if (dist_mid * dist1 < 0.0)
326 : {
327 1785974 : p2 = mid_point;
328 1785974 : dist2 = levelSetEvaluator(p2);
329 1785974 : dist = abs(dist1) + abs(dist2);
330 : }
331 : else
332 : {
333 1821237 : p1 = mid_point;
334 1821237 : dist1 = levelSetEvaluator(p1);
335 1821237 : dist = abs(dist1) + abs(dist2);
336 : }
337 : }
338 190002 : return mid_point;
339 : }
340 :
341 : const Node *
342 95001 : CutMeshByLevelSetGeneratorBase::nonDuplicateNodeCreator(
343 : ReplicatedMesh & mesh,
344 : std::vector<const Node *> & new_on_plane_nodes,
345 : const Point & new_point) const
346 : {
347 36110131 : for (const auto & new_on_plane_node : new_on_plane_nodes)
348 : {
349 36090460 : if (MooseUtils::absoluteFuzzyEqual((*new_on_plane_node - new_point).norm(), 0.0))
350 75330 : return new_on_plane_node;
351 : }
352 19671 : new_on_plane_nodes.push_back(mesh.add_point(new_point));
353 19671 : return new_on_plane_nodes.back();
354 : }
355 :
356 : void
357 40122 : CutMeshByLevelSetGeneratorBase::tet4ElemCutter(
358 : ReplicatedMesh & mesh,
359 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
360 : const dof_id_type elem_id,
361 : const subdomain_id_type & block_id_to_remove,
362 : std::vector<const Node *> & new_on_plane_nodes)
363 : {
364 : // Retrieve boundary information for the mesh
365 40122 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
366 : // Create a list of sidesets involving the element to be split
367 : // It might be complex to assign the boundary id to the new elements
368 : // In TET4, none of the four faces have the same normal vector
369 : // So we are using the normal vector to identify the faces that
370 : // need to be assigned the same boundary id
371 40122 : std::vector<Point> elem_side_normal_list;
372 40122 : elem_side_normal_list.resize(4);
373 200610 : for (const auto i : make_range(4))
374 : {
375 160488 : auto elem_side = mesh.elem_ptr(elem_id)->side_ptr(i);
376 320976 : elem_side_normal_list[i] = (*elem_side->node_ptr(1) - *elem_side->node_ptr(0))
377 160488 : .cross(*elem_side->node_ptr(2) - *elem_side->node_ptr(1))
378 160488 : .unit();
379 160488 : }
380 :
381 40122 : std::vector<std::vector<boundary_id_type>> elem_side_list;
382 40122 : MooseMeshElementConversionUtils::elementBoundaryInfoCollector(
383 : bdry_side_list, elem_id, 4, elem_side_list);
384 :
385 80244 : std::vector<PointLevelSetRelationIndex> node_plane_relation(4);
386 40122 : std::vector<const Node *> tet4_nodes(4);
387 40122 : std::vector<const Node *> tet4_nodes_on_plane;
388 40122 : std::vector<const Node *> tet4_nodes_outside_plane;
389 40122 : std::vector<const Node *> tet4_nodes_inside_plane;
390 : // Sort tetrahedral nodes based on their positioning wrt the plane
391 200610 : for (unsigned int i = 0; i < 4; i++)
392 : {
393 160488 : tet4_nodes[i] = mesh.elem_ptr(elem_id)->node_ptr(i);
394 160488 : node_plane_relation[i] = pointLevelSetRelation(*tet4_nodes[i]);
395 160488 : if (node_plane_relation[i] == PointLevelSetRelationIndex::on_level_set)
396 4230 : tet4_nodes_on_plane.push_back(tet4_nodes[i]);
397 156258 : else if (node_plane_relation[i] == PointLevelSetRelationIndex::level_set_out_side)
398 86682 : tet4_nodes_outside_plane.push_back(tet4_nodes[i]);
399 : else
400 69576 : tet4_nodes_inside_plane.push_back(tet4_nodes[i]);
401 : }
402 40122 : std::vector<Elem *> elems_tet4;
403 40122 : bool new_elements_created(false);
404 : // No cutting needed if no nodes are outside the plane
405 40122 : if (tet4_nodes_outside_plane.size() == 0)
406 : {
407 4876 : if (tet4_nodes_on_plane.size() == 3)
408 : {
409 : // Record the element for future cross section boundary assignment
410 360 : elems_tet4.push_back(mesh.elem_ptr(elem_id));
411 : }
412 : }
413 : // Remove the element if all the nodes are outside the plane
414 35246 : else if (tet4_nodes_inside_plane.size() == 0)
415 : {
416 6027 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
417 6027 : if (tet4_nodes_on_plane.size() == 3)
418 : {
419 : // I think the neighboring element will be handled,
420 : // So we do not need to assign the cross section boundary here
421 : }
422 : }
423 : // As we have nodes on both sides, six different scenarios are possible
424 : else
425 : {
426 29219 : new_elements_created = true;
427 29219 : if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 3)
428 : {
429 12245 : std::vector<const Node *> new_plane_nodes;
430 : // A smaller TET4 element is created, this solution is unique
431 48980 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
432 : {
433 36735 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
434 : mesh,
435 : new_on_plane_nodes,
436 73470 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0])));
437 : }
438 12245 : auto new_elem_tet4 = std::make_unique<Tet4>();
439 12245 : new_elem_tet4->set_node(0, const_cast<Node *>(tet4_nodes_inside_plane[0]));
440 12245 : new_elem_tet4->set_node(1, const_cast<Node *>(new_plane_nodes[0]));
441 12245 : new_elem_tet4->set_node(2, const_cast<Node *>(new_plane_nodes[1]));
442 12245 : new_elem_tet4->set_node(3, const_cast<Node *>(new_plane_nodes[2]));
443 12245 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
444 12245 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
445 12245 : }
446 16974 : else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 2)
447 : {
448 8583 : std::vector<const Node *> new_plane_nodes;
449 : // 3 smaller TET3 elements are created
450 25749 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
451 : {
452 51498 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
453 : {
454 34332 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
455 : mesh,
456 : new_on_plane_nodes,
457 68664 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_node_inside_plane)));
458 : }
459 : }
460 8583 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[1],
461 8583 : new_plane_nodes[3],
462 8583 : new_plane_nodes[1],
463 8583 : tet4_nodes_inside_plane[0],
464 8583 : new_plane_nodes[2],
465 17166 : new_plane_nodes[0]};
466 8583 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
467 8583 : std::vector<std::vector<const Node *>> optimized_node_list;
468 8583 : MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer(
469 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
470 :
471 34332 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
472 : {
473 25749 : auto new_elem_tet4 = std::make_unique<Tet4>();
474 25749 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
475 25749 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
476 25749 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
477 25749 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
478 25749 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
479 25749 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
480 25749 : }
481 8583 : }
482 8391 : else if (tet4_nodes_inside_plane.size() == 3 && tet4_nodes_outside_plane.size() == 1)
483 : {
484 7440 : std::vector<const Node *> new_plane_nodes;
485 : // 3 smaller Tet4 elements are created
486 29760 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
487 : {
488 22320 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
489 : mesh,
490 : new_on_plane_nodes,
491 44640 : pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0])));
492 : }
493 7440 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[0],
494 7440 : tet4_nodes_inside_plane[1],
495 7440 : tet4_nodes_inside_plane[2],
496 7440 : new_plane_nodes[0],
497 7440 : new_plane_nodes[1],
498 14880 : new_plane_nodes[2]};
499 7440 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
500 7440 : std::vector<std::vector<const Node *>> optimized_node_list;
501 7440 : MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer(
502 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
503 :
504 29760 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
505 : {
506 22320 : auto new_elem_tet4 = std::make_unique<Tet4>();
507 22320 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
508 22320 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
509 22320 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
510 22320 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
511 22320 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
512 22320 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
513 22320 : }
514 7440 : }
515 951 : else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 1)
516 : {
517 288 : auto new_plane_node = nonDuplicateNodeCreator(
518 : mesh,
519 : new_on_plane_nodes,
520 288 : pointPairLevelSetInterception(*tet4_nodes_inside_plane[0], *tet4_nodes_outside_plane[0]));
521 : // A smaller Tet4 is created, this solution is unique
522 288 : auto new_elem_tet4 = std::make_unique<Tet4>();
523 288 : new_elem_tet4->set_node(0, const_cast<Node *>(new_plane_node));
524 288 : new_elem_tet4->set_node(1, const_cast<Node *>(tet4_nodes_on_plane[0]));
525 288 : new_elem_tet4->set_node(2, const_cast<Node *>(tet4_nodes_on_plane[1]));
526 288 : new_elem_tet4->set_node(3, const_cast<Node *>(tet4_nodes_inside_plane[0]));
527 288 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
528 288 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
529 288 : }
530 663 : else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 2)
531 : {
532 327 : std::vector<const Node *> new_plane_nodes;
533 : // A smaller Tet4 element is created, this solution is unique
534 981 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
535 : {
536 654 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
537 : mesh,
538 : new_on_plane_nodes,
539 1308 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0])));
540 : }
541 327 : auto new_elem_tet4 = std::make_unique<Tet4>();
542 327 : new_elem_tet4->set_node(0, const_cast<Node *>(new_plane_nodes[0]));
543 327 : new_elem_tet4->set_node(1, const_cast<Node *>(new_plane_nodes[1]));
544 327 : new_elem_tet4->set_node(2, const_cast<Node *>(tet4_nodes_on_plane[0]));
545 327 : new_elem_tet4->set_node(3, const_cast<Node *>(tet4_nodes_inside_plane[0]));
546 327 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
547 327 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
548 327 : }
549 336 : else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 1)
550 : {
551 336 : std::vector<const Node *> new_plane_nodes;
552 : // 2 smaller TET4 elements are created
553 1008 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
554 : {
555 672 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
556 : mesh,
557 : new_on_plane_nodes,
558 1344 : pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0])));
559 : }
560 336 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[0],
561 336 : tet4_nodes_inside_plane[1],
562 336 : new_plane_nodes[1],
563 336 : new_plane_nodes[0],
564 672 : tet4_nodes_on_plane[0]};
565 336 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
566 336 : std::vector<std::vector<const Node *>> optimized_node_list;
567 336 : MooseMeshElementConversionUtils::pyramidNodesToTetNodesDeterminer(
568 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
569 :
570 1008 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
571 : {
572 672 : auto new_elem_tet4 = std::make_unique<Tet4>();
573 672 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
574 672 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
575 672 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
576 672 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
577 672 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
578 672 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
579 672 : }
580 336 : }
581 : else
582 0 : mooseError("Unexpected scenario.");
583 :
584 29219 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
585 : }
586 :
587 102083 : for (auto & elem_tet4 : elems_tet4)
588 : {
589 61961 : if (new_elements_created)
590 : {
591 61601 : if (elem_tet4->volume() < 0.0)
592 : {
593 25277 : Node * temp = elem_tet4->node_ptr(0);
594 25277 : elem_tet4->set_node(0, elem_tet4->node_ptr(1));
595 25277 : elem_tet4->set_node(1, temp);
596 : }
597 : }
598 : // Find the boundary id of the new element
599 309805 : for (unsigned int i = 0; i < 4; i++)
600 : {
601 247844 : const Point & side_pt_0 = *elem_tet4->side_ptr(i)->node_ptr(0);
602 247844 : const Point & side_pt_1 = *elem_tet4->side_ptr(i)->node_ptr(1);
603 247844 : const Point & side_pt_2 = *elem_tet4->side_ptr(i)->node_ptr(2);
604 :
605 247844 : const Point side_normal = (side_pt_1 - side_pt_0).cross(side_pt_2 - side_pt_1).unit();
606 1239220 : for (unsigned int j = 0; j < 4; j++)
607 : {
608 991376 : if (new_elements_created)
609 : {
610 985616 : if (MooseUtils::absoluteFuzzyEqual(side_normal * elem_side_normal_list[j], 1.0))
611 : {
612 146552 : for (const auto & side_info : elem_side_list[j])
613 : {
614 2714 : boundary_info.add_side(elem_tet4, i, side_info);
615 : }
616 : }
617 : }
618 : }
619 495688 : if (MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_0), 0.0) &&
620 389231 : MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_1), 0.0) &&
621 389231 : MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_2), 0.0))
622 : {
623 :
624 38162 : boundary_info.add_side(elem_tet4, i, _cut_face_id);
625 : }
626 : }
627 : }
628 40122 : }
629 :
630 : Real
631 8258821 : CutMeshByLevelSetGeneratorBase::levelSetEvaluator(const Point & point)
632 : {
633 8258821 : _func_params[0] = point(0);
634 8258821 : _func_params[1] = point(1);
635 8258821 : _func_params[2] = point(2);
636 24776463 : return evaluate(_func_level_set);
637 : }
|
