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 6334 : CutMeshByLevelSetGeneratorBase::validParams()
26 : {
27 6334 : InputParameters params = MeshGenerator::validParams();
28 6334 : params += FunctionParserUtils<false>::validParams();
29 :
30 25336 : params.addRequiredParam<MeshGeneratorName>("input", "The input mesh that needs to be trimmed.");
31 :
32 19002 : params.addParam<bool>(
33 : "generate_transition_layer",
34 12668 : 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 25336 : 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 19002 : params.addParam<BoundaryName>(
44 12668 : "cut_face_name", BoundaryName(), "The boundary name of the face generated by the cut.");
45 25336 : 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 25336 : 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 6334 : 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 6334 : return params;
63 0 : }
64 :
65 106 : CutMeshByLevelSetGeneratorBase::CutMeshByLevelSetGeneratorBase(const InputParameters & parameters)
66 : : MeshGenerator(parameters),
67 : FunctionParserUtils<false>(parameters),
68 106 : _input_name(getParam<MeshGeneratorName>("input")),
69 212 : _generate_transition_layer(getParam<bool>("generate_transition_layer")),
70 212 : _cut_face_name(getParam<BoundaryName>("cut_face_name")),
71 212 : _converted_tet_element_subdomain_name_suffix(
72 : getParam<SubdomainName>("converted_tet_element_subdomain_name_suffix")),
73 212 : _converted_pyramid_element_subdomain_name_suffix(
74 : getParam<SubdomainName>("converted_pyramid_element_subdomain_name_suffix")),
75 212 : _input(getMeshByName(_input_name))
76 : {
77 320 : _cut_face_id = isParamValid("cut_face_id") ? getParam<boundary_id_type>("cut_face_id") : -1;
78 106 : }
79 :
80 : std::unique_ptr<MeshBase>
81 106 : CutMeshByLevelSetGeneratorBase::generate()
82 : {
83 : // We're querying elem dim caches from our input mesh
84 106 : if (!_input->preparation().has_cached_elem_data)
85 103 : _input->cache_elem_data();
86 :
87 106 : auto replicated_mesh_ptr = dynamic_cast<ReplicatedMesh *>(_input.get());
88 106 : if (!replicated_mesh_ptr)
89 6 : paramError("input", "Input is not a replicated mesh, which is required");
90 203 : if (*(replicated_mesh_ptr->elem_dimensions().begin()) != 3 ||
91 203 : *(replicated_mesh_ptr->elem_dimensions().rbegin()) != 3)
92 6 : paramError(
93 : "input",
94 : "Only 3D meshes are supported. Use XYMeshLineCutter for 2D meshes if applicable. Mixed "
95 : "dimensional meshes are not supported at the moment.");
96 :
97 100 : ReplicatedMesh & mesh = *replicated_mesh_ptr;
98 :
99 100 : if (!_cut_face_name.empty())
100 : {
101 51 : if (MooseMeshUtils::hasBoundaryName(mesh, _cut_face_name))
102 : {
103 : const boundary_id_type exist_cut_face_id =
104 11 : MooseMeshUtils::getBoundaryID(_cut_face_name, mesh);
105 11 : if (_cut_face_id != -1 && _cut_face_id != exist_cut_face_id)
106 6 : paramError("cut_face_id",
107 : "The cut face boundary name and id are both provided, but they are inconsistent "
108 : "with an existing boundary name which has a different id.");
109 : else
110 8 : _cut_face_id = exist_cut_face_id;
111 : }
112 : else
113 : {
114 40 : if (_cut_face_id == -1)
115 8 : _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh);
116 40 : mesh.get_boundary_info().sideset_name(_cut_face_id) = _cut_face_name;
117 : }
118 : }
119 49 : else if (_cut_face_id == -1)
120 : {
121 41 : _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh);
122 : }
123 :
124 : // Subdomain ID for new utility blocks must be new
125 : // Find sid shifts
126 97 : const auto sid_shift_base = MooseMeshUtils::getNextFreeSubdomainID(mesh);
127 97 : const auto block_id_to_remove = sid_shift_base * 3;
128 :
129 : // For the boolean value in the pair, true means the element is crossed by the cutting plane
130 : // false means the element is on the remaining side
131 97 : std::vector<std::pair<dof_id_type, bool>> cross_and_remained_elems_pre_convert;
132 :
133 : // collect the subdomain ids of the original mesh
134 97 : std::set<subdomain_id_type> original_subdomain_ids;
135 14497 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
136 : elem_it++)
137 : {
138 14403 : original_subdomain_ids.emplace((*elem_it)->subdomain_id());
139 14403 : const unsigned int & n_vertices = (*elem_it)->n_vertices();
140 14403 : unsigned short elem_vertices_counter(0);
141 121499 : for (unsigned int i = 0; i < n_vertices; i++)
142 : {
143 : // We define elem_vertices_counter in this way so that those elements with one face on the
144 : // plane are also processed to have the cut face boundary assigned.
145 107096 : if (pointLevelSetRelation((*(*elem_it)->node_ptr(i))) !=
146 : PointLevelSetRelationIndex::level_set_in_side)
147 64763 : elem_vertices_counter++;
148 : }
149 14403 : if (elem_vertices_counter == n_vertices)
150 5152 : (*elem_it)->subdomain_id() = block_id_to_remove;
151 : else
152 : {
153 : // Check if any elements to be processed are not first order
154 9251 : if ((*elem_it)->default_order() != Order::FIRST)
155 3 : mooseError("Only first order elements are supported for cutting.");
156 : // If we are generating a transition layer, we only need to record the crossed elements
157 9248 : if (_generate_transition_layer && elem_vertices_counter > 0)
158 : {
159 1614 : cross_and_remained_elems_pre_convert.push_back(std::make_pair((*elem_it)->id(), true));
160 : // Since we are converting these elements to TET4, and we would keep the original type of
161 : // some elements in these blocks
162 3228 : if ((*elem_it)->type() != TET4)
163 1614 : (*elem_it)->subdomain_id() += sid_shift_base;
164 : }
165 7634 : else if (!_generate_transition_layer)
166 6912 : cross_and_remained_elems_pre_convert.push_back(
167 13824 : std::make_pair((*elem_it)->id(), elem_vertices_counter > 0));
168 : }
169 94 : }
170 :
171 : // Then we need to identify and make the transition layer
172 94 : std::vector<dof_id_type> transition_elems_list;
173 94 : std::vector<std::vector<unsigned int>> transition_elems_sides_list;
174 94 : if (_generate_transition_layer)
175 : {
176 30 : if (!mesh.is_prepared())
177 30 : mesh.find_neighbors();
178 : // First, we need to identify the retained elements that share the boundary with the crossed
179 : // elements.
180 1644 : for (const auto & elem_info : cross_and_remained_elems_pre_convert)
181 : {
182 11298 : for (const auto & i_side : make_range(mesh.elem_ptr(elem_info.first)->n_sides()))
183 : {
184 : // No need to work on a TRI3 side
185 9684 : if (mesh.elem_ptr(elem_info.first)->side_ptr(i_side)->type() == QUAD4)
186 : {
187 9684 : if (mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side) != nullptr &&
188 8172 : mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side)->subdomain_id() !=
189 17856 : block_id_to_remove &&
190 6968 : std::find(cross_and_remained_elems_pre_convert.begin(),
191 : cross_and_remained_elems_pre_convert.end(),
192 6968 : std::make_pair(mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side)->id(),
193 23620 : true)) == cross_and_remained_elems_pre_convert.end())
194 : {
195 : const auto & neighbor_elem_id =
196 1012 : mesh.elem_ptr(elem_info.first)->neighbor_ptr(i_side)->id();
197 : const auto & neighbor_elem_side_index =
198 1012 : mesh.elem_ptr(elem_info.first)
199 : ->neighbor_ptr(i_side)
200 1012 : ->which_neighbor_am_i(mesh.elem_ptr(elem_info.first));
201 1012 : const auto & id_found = std::find(
202 1012 : transition_elems_list.begin(), transition_elems_list.end(), neighbor_elem_id);
203 1012 : if (id_found == transition_elems_list.end())
204 : {
205 460 : transition_elems_list.push_back(neighbor_elem_id);
206 460 : transition_elems_sides_list.push_back(
207 920 : std::vector<unsigned int>({neighbor_elem_side_index}));
208 : }
209 : else
210 : {
211 : // If the element is already in the list, we just add the side index
212 552 : const auto index = std::distance(transition_elems_list.begin(), id_found);
213 552 : transition_elems_sides_list[index].push_back(neighbor_elem_side_index);
214 : }
215 : }
216 : }
217 : }
218 : }
219 : }
220 :
221 94 : std::vector<dof_id_type> converted_elems_ids_to_cut;
222 : // Then convert these non-TET4 elements into TET4 elements
223 94 : MooseMeshElementConversionUtils::convert3DMeshToAllTet4(mesh,
224 : cross_and_remained_elems_pre_convert,
225 : converted_elems_ids_to_cut,
226 : block_id_to_remove,
227 : false);
228 :
229 : // Make the transition layer if applicable
230 94 : auto & sideset_map = mesh.get_boundary_info().get_sideset_map();
231 554 : for (const auto & i_elem : index_range(transition_elems_list))
232 : {
233 460 : const auto & elem_id = transition_elems_list[i_elem];
234 460 : const auto & side_indices = transition_elems_sides_list[i_elem];
235 : // Find the involved sidesets of the element so that we can retain them
236 460 : std::vector<std::vector<boundary_id_type>> elem_side_info(mesh.elem_ptr(elem_id)->n_sides());
237 460 : auto side_range = sideset_map.equal_range(mesh.elem_ptr(elem_id));
238 798 : for (auto i = side_range.first; i != side_range.second; ++i)
239 338 : elem_side_info[i->second.first].push_back(i->second.second);
240 :
241 460 : MooseMeshElementConversionUtils::convertElem(
242 : mesh, elem_id, side_indices, elem_side_info, sid_shift_base);
243 460 : }
244 :
245 94 : std::vector<const Node *> new_on_plane_nodes;
246 : // We build the sideset information now, as we only need the information of the elements before
247 : // cutting
248 94 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
249 94 : const auto bdry_side_list = boundary_info.build_side_list();
250 : // Cut the TET4 Elements
251 35538 : for (const auto & converted_elems_id_to_cut : converted_elems_ids_to_cut)
252 : {
253 35444 : tet4ElemCutter(
254 : mesh, bdry_side_list, converted_elems_id_to_cut, block_id_to_remove, new_on_plane_nodes);
255 : }
256 :
257 : // If a transition layer is generated, we need to add some subdomain names
258 94 : if (_generate_transition_layer)
259 : {
260 : try
261 : {
262 30 : MooseMeshElementConversionUtils::assignConvertedElementsSubdomainNameSuffix(
263 : mesh,
264 : original_subdomain_ids,
265 : sid_shift_base,
266 30 : _converted_tet_element_subdomain_name_suffix,
267 30 : _converted_pyramid_element_subdomain_name_suffix);
268 : }
269 6 : catch (const std::exception & e)
270 : {
271 12 : if (((std::string)e.what()).compare(26, 4, "TET4") == 0)
272 6 : paramError("converted_tet_element_subdomain_name_suffix", e.what());
273 : else // if (((std::string)e.what()).compare(26, 7, "PYRAMID5") == 0)
274 6 : paramError("converted_pyramid_element_subdomain_name_suffix", e.what());
275 0 : }
276 : }
277 :
278 : // delete the original elements that were converted
279 488 : for (const auto & elem_id : transition_elems_list)
280 400 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
281 : // Delete the block to remove
282 43424 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
283 43424 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
284 : elem_it++)
285 43424 : mesh.delete_elem(*elem_it);
286 :
287 88 : mesh.contract();
288 88 : mesh.unset_is_prepared();
289 176 : return std::move(_input);
290 88 : }
291 :
292 : CutMeshByLevelSetGeneratorBase::PointLevelSetRelationIndex
293 248872 : CutMeshByLevelSetGeneratorBase::pointLevelSetRelation(const Point & point)
294 : {
295 248872 : const Real level_set_value = levelSetEvaluator(point);
296 248872 : if (MooseUtils::absoluteFuzzyLessThan(level_set_value, 0.0))
297 103725 : return PointLevelSetRelationIndex::level_set_in_side;
298 145147 : else if (MooseUtils::absoluteFuzzyGreaterThan(level_set_value, 0.0))
299 138284 : return PointLevelSetRelationIndex::level_set_out_side;
300 : else
301 6863 : return PointLevelSetRelationIndex::on_level_set;
302 : }
303 :
304 : Point
305 83812 : CutMeshByLevelSetGeneratorBase::pointPairLevelSetInterception(const Point & point1,
306 : const Point & point2)
307 : {
308 83812 : Real dist1 = levelSetEvaluator(point1);
309 83812 : Real dist2 = levelSetEvaluator(point2);
310 :
311 83812 : if (MooseUtils::absoluteFuzzyEqual(dist1, 0.0) || MooseUtils::absoluteFuzzyEqual(dist2, 0.0))
312 0 : mooseError("At least one of the two points are on the plane.");
313 83812 : if (MooseUtils::absoluteFuzzyGreaterThan(dist1 * dist2, 0.0))
314 0 : mooseError("The two points are on the same side of the plane.");
315 :
316 83812 : Point p1(point1);
317 83812 : Point p2(point2);
318 83812 : Real dist = abs(dist1) + abs(dist2);
319 83812 : Point mid_point;
320 :
321 : // Bisection method to find midpoint
322 3267928 : while (MooseUtils::absoluteFuzzyGreaterThan(dist, 0.0))
323 : {
324 3184116 : mid_point = 0.5 * (p1 + p2);
325 3184116 : const Real dist_mid = levelSetEvaluator(mid_point);
326 : // We do not need Fuzzy here as it will be covered by the while loop
327 3184116 : if (dist_mid == 0.0)
328 1984 : dist = 0.0;
329 3182132 : else if (dist_mid * dist1 < 0.0)
330 : {
331 1575828 : p2 = mid_point;
332 1575828 : dist2 = levelSetEvaluator(p2);
333 1575828 : dist = abs(dist1) + abs(dist2);
334 : }
335 : else
336 : {
337 1606304 : p1 = mid_point;
338 1606304 : dist1 = levelSetEvaluator(p1);
339 1606304 : dist = abs(dist1) + abs(dist2);
340 : }
341 : }
342 167624 : return mid_point;
343 : }
344 :
345 : const Node *
346 83812 : CutMeshByLevelSetGeneratorBase::nonDuplicateNodeCreator(
347 : ReplicatedMesh & mesh,
348 : std::vector<const Node *> & new_on_plane_nodes,
349 : const Point & new_point) const
350 : {
351 32054322 : for (const auto & new_on_plane_node : new_on_plane_nodes)
352 : {
353 32036990 : if (MooseUtils::absoluteFuzzyEqual((*new_on_plane_node - new_point).norm(), 0.0))
354 66480 : return new_on_plane_node;
355 : }
356 17332 : new_on_plane_nodes.push_back(mesh.add_point(new_point));
357 17332 : return new_on_plane_nodes.back();
358 : }
359 :
360 : void
361 35444 : CutMeshByLevelSetGeneratorBase::tet4ElemCutter(
362 : ReplicatedMesh & mesh,
363 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
364 : const dof_id_type elem_id,
365 : const subdomain_id_type & block_id_to_remove,
366 : std::vector<const Node *> & new_on_plane_nodes)
367 : {
368 : // Retrieve boundary information for the mesh
369 35444 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
370 : // Create a list of sidesets involving the element to be split
371 : // It might be complex to assign the boundary id to the new elements
372 : // In TET4, none of the four faces have the same normal vector
373 : // So we are using the normal vector to identify the faces that
374 : // need to be assigned the same boundary id
375 35444 : std::vector<Point> elem_side_normal_list;
376 35444 : elem_side_normal_list.resize(4);
377 177220 : for (const auto i : make_range(4))
378 : {
379 141776 : auto elem_side = mesh.elem_ptr(elem_id)->side_ptr(i);
380 283552 : elem_side_normal_list[i] = (*elem_side->node_ptr(1) - *elem_side->node_ptr(0))
381 141776 : .cross(*elem_side->node_ptr(2) - *elem_side->node_ptr(1))
382 141776 : .unit();
383 141776 : }
384 :
385 35444 : std::vector<std::vector<boundary_id_type>> elem_side_list;
386 35444 : MooseMeshElementConversionUtils::elementBoundaryInfoCollector(
387 : bdry_side_list, elem_id, 4, elem_side_list);
388 :
389 70888 : std::vector<PointLevelSetRelationIndex> node_plane_relation(4);
390 35444 : std::vector<const Node *> tet4_nodes(4);
391 35444 : std::vector<const Node *> tet4_nodes_on_plane;
392 35444 : std::vector<const Node *> tet4_nodes_outside_plane;
393 35444 : std::vector<const Node *> tet4_nodes_inside_plane;
394 : // Sort tetrahedral nodes based on their positioning wrt the plane
395 177220 : for (unsigned int i = 0; i < 4; i++)
396 : {
397 141776 : tet4_nodes[i] = mesh.elem_ptr(elem_id)->node_ptr(i);
398 141776 : node_plane_relation[i] = pointLevelSetRelation(*tet4_nodes[i]);
399 141776 : if (node_plane_relation[i] == PointLevelSetRelationIndex::on_level_set)
400 3740 : tet4_nodes_on_plane.push_back(tet4_nodes[i]);
401 138036 : else if (node_plane_relation[i] == PointLevelSetRelationIndex::level_set_out_side)
402 76644 : tet4_nodes_outside_plane.push_back(tet4_nodes[i]);
403 : else
404 61392 : tet4_nodes_inside_plane.push_back(tet4_nodes[i]);
405 : }
406 35444 : std::vector<Elem *> elems_tet4;
407 35444 : bool new_elements_created(false);
408 : // No cutting needed if no nodes are outside the plane
409 35444 : if (tet4_nodes_outside_plane.size() == 0)
410 : {
411 4312 : if (tet4_nodes_on_plane.size() == 3)
412 : {
413 : // Record the element for future cross section boundary assignment
414 320 : elems_tet4.push_back(mesh.elem_ptr(elem_id));
415 : }
416 : }
417 : // Remove the element if all the nodes are outside the plane
418 31132 : else if (tet4_nodes_inside_plane.size() == 0)
419 : {
420 5354 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
421 5354 : if (tet4_nodes_on_plane.size() == 3)
422 : {
423 : // I think the neighboring element will be handled,
424 : // So we do not need to assign the cross section boundary here
425 : }
426 : }
427 : // As we have nodes on both sides, six different scenarios are possible
428 : else
429 : {
430 25778 : new_elements_created = true;
431 25778 : if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 3)
432 : {
433 10820 : std::vector<const Node *> new_plane_nodes;
434 : // A smaller TET4 element is created, this solution is unique
435 43280 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
436 : {
437 32460 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
438 : mesh,
439 : new_on_plane_nodes,
440 64920 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0])));
441 : }
442 10820 : auto new_elem_tet4 = std::make_unique<Tet4>();
443 10820 : new_elem_tet4->set_node(0, const_cast<Node *>(tet4_nodes_inside_plane[0]));
444 10820 : new_elem_tet4->set_node(1, const_cast<Node *>(new_plane_nodes[0]));
445 10820 : new_elem_tet4->set_node(2, const_cast<Node *>(new_plane_nodes[1]));
446 10820 : new_elem_tet4->set_node(3, const_cast<Node *>(new_plane_nodes[2]));
447 10820 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
448 10820 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
449 10820 : }
450 14958 : else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 2)
451 : {
452 7566 : std::vector<const Node *> new_plane_nodes;
453 : // 3 smaller TET3 elements are created
454 22698 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
455 : {
456 45396 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
457 : {
458 30264 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
459 : mesh,
460 : new_on_plane_nodes,
461 60528 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_node_inside_plane)));
462 : }
463 : }
464 7566 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[1],
465 7566 : new_plane_nodes[3],
466 7566 : new_plane_nodes[1],
467 7566 : tet4_nodes_inside_plane[0],
468 7566 : new_plane_nodes[2],
469 15132 : new_plane_nodes[0]};
470 7566 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
471 7566 : std::vector<std::vector<const Node *>> optimized_node_list;
472 7566 : MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer(
473 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
474 :
475 30264 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
476 : {
477 22698 : auto new_elem_tet4 = std::make_unique<Tet4>();
478 22698 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
479 22698 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
480 22698 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
481 22698 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
482 22698 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
483 22698 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
484 22698 : }
485 7566 : }
486 7392 : else if (tet4_nodes_inside_plane.size() == 3 && tet4_nodes_outside_plane.size() == 1)
487 : {
488 6560 : std::vector<const Node *> new_plane_nodes;
489 : // 3 smaller Tet4 elements are created
490 26240 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
491 : {
492 19680 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
493 : mesh,
494 : new_on_plane_nodes,
495 39360 : pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0])));
496 : }
497 6560 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[0],
498 6560 : tet4_nodes_inside_plane[1],
499 6560 : tet4_nodes_inside_plane[2],
500 6560 : new_plane_nodes[0],
501 6560 : new_plane_nodes[1],
502 13120 : new_plane_nodes[2]};
503 6560 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
504 6560 : std::vector<std::vector<const Node *>> optimized_node_list;
505 6560 : MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer(
506 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
507 :
508 26240 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
509 : {
510 19680 : auto new_elem_tet4 = std::make_unique<Tet4>();
511 19680 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
512 19680 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
513 19680 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
514 19680 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
515 19680 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
516 19680 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
517 19680 : }
518 6560 : }
519 832 : else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 1)
520 : {
521 256 : auto new_plane_node = nonDuplicateNodeCreator(
522 : mesh,
523 : new_on_plane_nodes,
524 256 : pointPairLevelSetInterception(*tet4_nodes_inside_plane[0], *tet4_nodes_outside_plane[0]));
525 : // A smaller Tet4 is created, this solution is unique
526 256 : auto new_elem_tet4 = std::make_unique<Tet4>();
527 256 : new_elem_tet4->set_node(0, const_cast<Node *>(new_plane_node));
528 256 : new_elem_tet4->set_node(1, const_cast<Node *>(tet4_nodes_on_plane[0]));
529 256 : new_elem_tet4->set_node(2, const_cast<Node *>(tet4_nodes_on_plane[1]));
530 256 : new_elem_tet4->set_node(3, const_cast<Node *>(tet4_nodes_inside_plane[0]));
531 256 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
532 256 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
533 256 : }
534 576 : else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 2)
535 : {
536 284 : std::vector<const Node *> new_plane_nodes;
537 : // A smaller Tet4 element is created, this solution is unique
538 852 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
539 : {
540 568 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
541 : mesh,
542 : new_on_plane_nodes,
543 1136 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0])));
544 : }
545 284 : auto new_elem_tet4 = std::make_unique<Tet4>();
546 284 : new_elem_tet4->set_node(0, const_cast<Node *>(new_plane_nodes[0]));
547 284 : new_elem_tet4->set_node(1, const_cast<Node *>(new_plane_nodes[1]));
548 284 : new_elem_tet4->set_node(2, const_cast<Node *>(tet4_nodes_on_plane[0]));
549 284 : new_elem_tet4->set_node(3, const_cast<Node *>(tet4_nodes_inside_plane[0]));
550 284 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
551 284 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
552 284 : }
553 292 : else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 1)
554 : {
555 292 : std::vector<const Node *> new_plane_nodes;
556 : // 2 smaller TET4 elements are created
557 876 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
558 : {
559 584 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
560 : mesh,
561 : new_on_plane_nodes,
562 1168 : pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0])));
563 : }
564 292 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[0],
565 292 : tet4_nodes_inside_plane[1],
566 292 : new_plane_nodes[1],
567 292 : new_plane_nodes[0],
568 584 : tet4_nodes_on_plane[0]};
569 292 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
570 292 : std::vector<std::vector<const Node *>> optimized_node_list;
571 292 : MooseMeshElementConversionUtils::pyramidNodesToTetNodesDeterminer(
572 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
573 :
574 876 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
575 : {
576 584 : auto new_elem_tet4 = std::make_unique<Tet4>();
577 584 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
578 584 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
579 584 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
580 584 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
581 584 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
582 584 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
583 584 : }
584 292 : }
585 : else
586 0 : mooseError("Unexpected scenario.");
587 :
588 25778 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
589 : }
590 :
591 90086 : for (auto & elem_tet4 : elems_tet4)
592 : {
593 54642 : if (new_elements_created)
594 : {
595 54322 : if (elem_tet4->volume() < 0.0)
596 : {
597 22354 : Node * temp = elem_tet4->node_ptr(0);
598 22354 : elem_tet4->set_node(0, elem_tet4->node_ptr(1));
599 22354 : elem_tet4->set_node(1, temp);
600 : }
601 : }
602 : // Find the boundary id of the new element
603 273210 : for (unsigned int i = 0; i < 4; i++)
604 : {
605 218568 : const Point & side_pt_0 = *elem_tet4->side_ptr(i)->node_ptr(0);
606 218568 : const Point & side_pt_1 = *elem_tet4->side_ptr(i)->node_ptr(1);
607 218568 : const Point & side_pt_2 = *elem_tet4->side_ptr(i)->node_ptr(2);
608 :
609 218568 : const Point side_normal = (side_pt_1 - side_pt_0).cross(side_pt_2 - side_pt_1).unit();
610 1092840 : for (unsigned int j = 0; j < 4; j++)
611 : {
612 874272 : if (new_elements_created)
613 : {
614 869152 : if (MooseUtils::absoluteFuzzyEqual(side_normal * elem_side_normal_list[j], 1.0))
615 : {
616 129184 : for (const auto & side_info : elem_side_list[j])
617 : {
618 2328 : boundary_info.add_side(elem_tet4, i, side_info);
619 : }
620 : }
621 : }
622 : }
623 437136 : if (MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_0), 0.0) &&
624 343272 : MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_1), 0.0) &&
625 343272 : MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_2), 0.0))
626 : {
627 :
628 33664 : boundary_info.add_side(elem_tet4, i, _cut_face_id);
629 : }
630 : }
631 : }
632 35444 : }
633 :
634 : Real
635 7285962 : CutMeshByLevelSetGeneratorBase::levelSetEvaluator(const Point & point)
636 : {
637 7285962 : _func_params[0] = point(0);
638 7285962 : _func_params[1] = point(1);
639 7285962 : _func_params[2] = point(2);
640 21857886 : return evaluate(_func_level_set);
641 : }
|
