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 28690 : CutMeshByLevelSetGeneratorBase::validParams()
26 : {
27 28690 : InputParameters params = MeshGenerator::validParams();
28 28690 : params += FunctionParserUtils<false>::validParams();
29 :
30 28690 : params.addRequiredParam<MeshGeneratorName>("input", "The input mesh that needs to be trimmed.");
31 :
32 28690 : params.addParam<boundary_id_type>("cut_face_id",
33 : "The boundary id of the face generated by the cut. An "
34 : "id will be automatically assigned if not provided.");
35 86070 : params.addParam<BoundaryName>(
36 57380 : "cut_face_name", BoundaryName(), "The boundary name of the face generated by the cut.");
37 :
38 28690 : params.addClassDescription(
39 : "This CutMeshByLevelSetGeneratorBase object is designed to be the base class of mesh "
40 : "generator that cuts a 3D mesh based on an analytic level set function. The level set "
41 : "function could be provided explicitly or indirectly.");
42 :
43 28690 : return params;
44 0 : }
45 :
46 80 : CutMeshByLevelSetGeneratorBase::CutMeshByLevelSetGeneratorBase(const InputParameters & parameters)
47 : : MeshGenerator(parameters),
48 : FunctionParserUtils<false>(parameters),
49 80 : _input_name(getParam<MeshGeneratorName>("input")),
50 80 : _cut_face_name(getParam<BoundaryName>("cut_face_name")),
51 160 : _input(getMeshByName(_input_name))
52 : {
53 80 : _cut_face_id = isParamValid("cut_face_id") ? getParam<boundary_id_type>("cut_face_id") : -1;
54 80 : }
55 :
56 : std::unique_ptr<MeshBase>
57 80 : CutMeshByLevelSetGeneratorBase::generate()
58 : {
59 80 : auto replicated_mesh_ptr = dynamic_cast<ReplicatedMesh *>(_input.get());
60 80 : if (!replicated_mesh_ptr)
61 4 : paramError("input", "Input is not a replicated mesh, which is required");
62 148 : if (*(replicated_mesh_ptr->elem_dimensions().begin()) != 3 ||
63 148 : *(replicated_mesh_ptr->elem_dimensions().rbegin()) != 3)
64 4 : paramError(
65 : "input",
66 : "Only 3D meshes are supported. Use XYMeshLineCutter for 2D meshes if applicable. Mixed "
67 : "dimensional meshes are not supported at the moment.");
68 :
69 72 : ReplicatedMesh & mesh = *replicated_mesh_ptr;
70 :
71 72 : if (!_cut_face_name.empty())
72 : {
73 44 : if (MooseMeshUtils::hasBoundaryName(mesh, _cut_face_name))
74 : {
75 : const boundary_id_type exist_cut_face_id =
76 12 : MooseMeshUtils::getBoundaryID(_cut_face_name, mesh);
77 12 : if (_cut_face_id != -1 && _cut_face_id != exist_cut_face_id)
78 4 : paramError("cut_face_id",
79 : "The cut face boundary name and id are both provided, but they are inconsistent "
80 : "with an existing boundary name which has a different id.");
81 : else
82 8 : _cut_face_id = exist_cut_face_id;
83 : }
84 : else
85 : {
86 32 : if (_cut_face_id == -1)
87 8 : _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh);
88 32 : mesh.get_boundary_info().sideset_name(_cut_face_id) = _cut_face_name;
89 : }
90 : }
91 28 : else if (_cut_face_id == -1)
92 : {
93 20 : _cut_face_id = MooseMeshUtils::getNextFreeBoundaryID(mesh);
94 : }
95 :
96 : // Subdomain ID for new utility blocks must be new
97 68 : std::set<subdomain_id_type> subdomain_ids_set;
98 68 : mesh.subdomain_ids(subdomain_ids_set);
99 68 : const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin();
100 68 : const subdomain_id_type block_id_to_remove = max_subdomain_id + 1;
101 :
102 : // For the boolean value in the pair, true means the element is crossed by the cutting plane
103 : // false means the element is on the remaining side
104 68 : std::vector<std::pair<dof_id_type, bool>> cross_and_remained_elems_pre_convert;
105 :
106 11332 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
107 : elem_it++)
108 : {
109 11268 : const unsigned int & n_vertices = (*elem_it)->n_vertices();
110 11268 : unsigned short elem_vertices_counter(0);
111 93284 : for (unsigned int i = 0; i < n_vertices; i++)
112 : {
113 : // We define elem_vertices_counter in this way so that those elements with one face on the
114 : // plane are also processed to have the cut face boundary assigned.
115 82016 : if (pointLevelSetRelation((*(*elem_it)->node_ptr(i))) !=
116 : PointLevelSetRelationIndex::level_set_in_side)
117 50948 : elem_vertices_counter++;
118 : }
119 11268 : if (elem_vertices_counter == n_vertices)
120 4352 : (*elem_it)->subdomain_id() = block_id_to_remove;
121 : else
122 : {
123 : // Check if any elements to be processed are not first order
124 6916 : if ((*elem_it)->default_order() != Order::FIRST)
125 4 : mooseError("Only first order elements are supported for cutting.");
126 6912 : cross_and_remained_elems_pre_convert.push_back(
127 13824 : std::make_pair((*elem_it)->id(), elem_vertices_counter > 0));
128 : }
129 64 : }
130 :
131 64 : std::vector<dof_id_type> converted_elems_ids_to_cut;
132 : // Then convert these non-TET4 elements into TET4 elements
133 64 : MooseMeshElementConversionUtils::convert3DMeshToAllTet4(mesh,
134 : cross_and_remained_elems_pre_convert,
135 : converted_elems_ids_to_cut,
136 : block_id_to_remove,
137 : false);
138 :
139 64 : std::vector<const Node *> new_on_plane_nodes;
140 : // We build the sideset information now, as we only need the information of the elements before
141 : // cutting
142 64 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
143 64 : const auto bdry_side_list = boundary_info.build_side_list();
144 : // Cut the TET4 Elements
145 25824 : for (const auto & converted_elems_id_to_cut : converted_elems_ids_to_cut)
146 : {
147 25760 : tet4ElemCutter(
148 : mesh, bdry_side_list, converted_elems_id_to_cut, block_id_to_remove, new_on_plane_nodes);
149 : }
150 :
151 : // Delete the block to remove
152 33600 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
153 33600 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
154 : elem_it++)
155 33600 : mesh.delete_elem(*elem_it);
156 :
157 64 : mesh.contract();
158 64 : mesh.set_isnt_prepared();
159 128 : return std::move(_input);
160 64 : }
161 :
162 : CutMeshByLevelSetGeneratorBase::PointLevelSetRelationIndex
163 185056 : CutMeshByLevelSetGeneratorBase::pointLevelSetRelation(const Point & point)
164 : {
165 185056 : const Real level_set_value = levelSetEvaluator(point);
166 185056 : if (MooseUtils::absoluteFuzzyLessThan(level_set_value, 0.0))
167 75804 : return PointLevelSetRelationIndex::level_set_in_side;
168 109252 : else if (MooseUtils::absoluteFuzzyGreaterThan(level_set_value, 0.0))
169 105312 : return PointLevelSetRelationIndex::level_set_out_side;
170 : else
171 3940 : return PointLevelSetRelationIndex::on_level_set;
172 : }
173 :
174 : Point
175 60280 : CutMeshByLevelSetGeneratorBase::pointPairLevelSetInterception(const Point & point1,
176 : const Point & point2)
177 : {
178 60280 : Real dist1 = levelSetEvaluator(point1);
179 60280 : Real dist2 = levelSetEvaluator(point2);
180 :
181 60280 : if (MooseUtils::absoluteFuzzyEqual(dist1, 0.0) || MooseUtils::absoluteFuzzyEqual(dist2, 0.0))
182 0 : mooseError("At least one of the two points are on the plane.");
183 60280 : if (MooseUtils::absoluteFuzzyGreaterThan(dist1 * dist2, 0.0))
184 0 : mooseError("The two points are on the same side of the plane.");
185 :
186 60280 : Point p1(point1);
187 60280 : Point p2(point2);
188 60280 : Real dist = abs(dist1) + abs(dist2);
189 60280 : Point mid_point;
190 :
191 : // Bisection method to find midpoint
192 2324288 : while (MooseUtils::absoluteFuzzyGreaterThan(dist, 0.0))
193 : {
194 2264008 : mid_point = 0.5 * (p1 + p2);
195 2264008 : const Real dist_mid = levelSetEvaluator(mid_point);
196 : // We do not need Fuzzy here as it will be covered by the while loop
197 2264008 : if (dist_mid == 0.0)
198 1984 : dist = 0.0;
199 2262024 : else if (dist_mid * dist1 < 0.0)
200 : {
201 1121248 : p2 = mid_point;
202 1121248 : dist2 = levelSetEvaluator(p2);
203 1121248 : dist = abs(dist1) + abs(dist2);
204 : }
205 : else
206 : {
207 1140776 : p1 = mid_point;
208 1140776 : dist1 = levelSetEvaluator(p1);
209 1140776 : dist = abs(dist1) + abs(dist2);
210 : }
211 : }
212 120560 : return mid_point;
213 : }
214 :
215 : const Node *
216 60280 : CutMeshByLevelSetGeneratorBase::nonDuplicateNodeCreator(
217 : ReplicatedMesh & mesh,
218 : std::vector<const Node *> & new_on_plane_nodes,
219 : const Point & new_point) const
220 : {
221 28502976 : for (const auto & new_on_plane_node : new_on_plane_nodes)
222 : {
223 28490664 : if (MooseUtils::absoluteFuzzyEqual((*new_on_plane_node - new_point).norm(), 0.0))
224 47968 : return new_on_plane_node;
225 : }
226 12312 : new_on_plane_nodes.push_back(mesh.add_point(new_point));
227 12312 : return new_on_plane_nodes.back();
228 : }
229 :
230 : void
231 25760 : CutMeshByLevelSetGeneratorBase::tet4ElemCutter(
232 : ReplicatedMesh & mesh,
233 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
234 : const dof_id_type elem_id,
235 : const subdomain_id_type & block_id_to_remove,
236 : std::vector<const Node *> & new_on_plane_nodes)
237 : {
238 : // Retrieve boundary information for the mesh
239 25760 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
240 : // Create a list of sidesets involving the element to be split
241 : // It might be complex to assign the boundary id to the new elements
242 : // In TET4, none of the four faces have the same normal vector
243 : // So we are using the normal vector to identify the faces that
244 : // need to be assigned the same boundary id
245 25760 : std::vector<Point> elem_side_normal_list;
246 25760 : elem_side_normal_list.resize(4);
247 128800 : for (const auto i : make_range(4))
248 : {
249 103040 : auto elem_side = mesh.elem_ptr(elem_id)->side_ptr(i);
250 206080 : elem_side_normal_list[i] = (*elem_side->node_ptr(1) - *elem_side->node_ptr(0))
251 103040 : .cross(*elem_side->node_ptr(2) - *elem_side->node_ptr(1))
252 103040 : .unit();
253 103040 : }
254 :
255 25760 : std::vector<std::vector<boundary_id_type>> elem_side_list;
256 25760 : MooseMeshElementConversionUtils::elementBoundaryInfoCollector(
257 : bdry_side_list, elem_id, 4, elem_side_list);
258 :
259 25760 : std::vector<PointLevelSetRelationIndex> node_plane_relation(4);
260 25760 : std::vector<const Node *> tet4_nodes(4);
261 25760 : std::vector<const Node *> tet4_nodes_on_plane;
262 25760 : std::vector<const Node *> tet4_nodes_outside_plane;
263 25760 : std::vector<const Node *> tet4_nodes_inside_plane;
264 : // Sort tetrahedral nodes based on their positioning wrt the plane
265 128800 : for (unsigned int i = 0; i < 4; i++)
266 : {
267 103040 : tet4_nodes[i] = mesh.elem_ptr(elem_id)->node_ptr(i);
268 103040 : node_plane_relation[i] = pointLevelSetRelation(*tet4_nodes[i]);
269 103040 : if (node_plane_relation[i] == PointLevelSetRelationIndex::on_level_set)
270 1952 : tet4_nodes_on_plane.push_back(tet4_nodes[i]);
271 101088 : else if (node_plane_relation[i] == PointLevelSetRelationIndex::level_set_out_side)
272 56352 : tet4_nodes_outside_plane.push_back(tet4_nodes[i]);
273 : else
274 44736 : tet4_nodes_inside_plane.push_back(tet4_nodes[i]);
275 : }
276 25760 : std::vector<Elem *> elems_tet4;
277 25760 : bool new_elements_created(false);
278 : // No cutting needed if no nodes are outside the plane
279 25760 : if (tet4_nodes_outside_plane.size() == 0)
280 : {
281 2976 : if (tet4_nodes_on_plane.size() == 3)
282 : {
283 : // Record the element for future cross section boundary assignment
284 64 : elems_tet4.push_back(mesh.elem_ptr(elem_id));
285 : }
286 : }
287 : // Remove the element if all the nodes are outside the plane
288 22784 : else if (tet4_nodes_inside_plane.size() == 0)
289 : {
290 4160 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
291 4160 : if (tet4_nodes_on_plane.size() == 3)
292 : {
293 : // I think the neighboring element will be handled,
294 : // So we do not need to assign the cross section boundary here
295 : }
296 : }
297 : // As we have nodes on both sides, six different scenarios are possible
298 : else
299 : {
300 18624 : new_elements_created = true;
301 18624 : if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 3)
302 : {
303 7800 : std::vector<const Node *> new_plane_nodes;
304 : // A smaller TET4 element is created, this solution is unique
305 31200 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
306 : {
307 23400 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
308 : mesh,
309 : new_on_plane_nodes,
310 46800 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0])));
311 : }
312 7800 : auto new_elem_tet4 = std::make_unique<Tet4>();
313 7800 : new_elem_tet4->set_node(0, const_cast<Node *>(tet4_nodes_inside_plane[0]));
314 7800 : new_elem_tet4->set_node(1, const_cast<Node *>(new_plane_nodes[0]));
315 7800 : new_elem_tet4->set_node(2, const_cast<Node *>(new_plane_nodes[1]));
316 7800 : new_elem_tet4->set_node(3, const_cast<Node *>(new_plane_nodes[2]));
317 7800 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
318 7800 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
319 7800 : }
320 10824 : else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 2)
321 : {
322 5328 : std::vector<const Node *> new_plane_nodes;
323 : // 3 smaller TET3 elements are created
324 15984 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
325 : {
326 31968 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
327 : {
328 21312 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
329 : mesh,
330 : new_on_plane_nodes,
331 42624 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_node_inside_plane)));
332 : }
333 : }
334 5328 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[1],
335 5328 : new_plane_nodes[3],
336 5328 : new_plane_nodes[1],
337 5328 : tet4_nodes_inside_plane[0],
338 5328 : new_plane_nodes[2],
339 5328 : new_plane_nodes[0]};
340 5328 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
341 5328 : std::vector<std::vector<const Node *>> optimized_node_list;
342 5328 : MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer(
343 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
344 :
345 21312 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
346 : {
347 15984 : auto new_elem_tet4 = std::make_unique<Tet4>();
348 15984 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
349 15984 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
350 15984 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
351 15984 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
352 15984 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
353 15984 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
354 15984 : }
355 5328 : }
356 5496 : else if (tet4_nodes_inside_plane.size() == 3 && tet4_nodes_outside_plane.size() == 1)
357 : {
358 4832 : std::vector<const Node *> new_plane_nodes;
359 : // 3 smaller Tet4 elements are created
360 19328 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
361 : {
362 14496 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
363 : mesh,
364 : new_on_plane_nodes,
365 28992 : pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0])));
366 : }
367 4832 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[0],
368 4832 : tet4_nodes_inside_plane[1],
369 4832 : tet4_nodes_inside_plane[2],
370 4832 : new_plane_nodes[0],
371 4832 : new_plane_nodes[1],
372 4832 : new_plane_nodes[2]};
373 4832 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
374 4832 : std::vector<std::vector<const Node *>> optimized_node_list;
375 4832 : MooseMeshElementConversionUtils::prismNodesToTetNodesDeterminer(
376 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
377 :
378 19328 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
379 : {
380 14496 : auto new_elem_tet4 = std::make_unique<Tet4>();
381 14496 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
382 14496 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
383 14496 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
384 14496 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
385 14496 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
386 14496 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
387 14496 : }
388 4832 : }
389 664 : else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 1)
390 : {
391 256 : auto new_plane_node = nonDuplicateNodeCreator(
392 : mesh,
393 : new_on_plane_nodes,
394 256 : pointPairLevelSetInterception(*tet4_nodes_inside_plane[0], *tet4_nodes_outside_plane[0]));
395 : // A smaller Tet4 is created, this solution is unique
396 256 : auto new_elem_tet4 = std::make_unique<Tet4>();
397 256 : new_elem_tet4->set_node(0, const_cast<Node *>(new_plane_node));
398 256 : new_elem_tet4->set_node(1, const_cast<Node *>(tet4_nodes_on_plane[0]));
399 256 : new_elem_tet4->set_node(2, const_cast<Node *>(tet4_nodes_on_plane[1]));
400 256 : new_elem_tet4->set_node(3, const_cast<Node *>(tet4_nodes_inside_plane[0]));
401 256 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
402 256 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
403 256 : }
404 408 : else if (tet4_nodes_inside_plane.size() == 1 && tet4_nodes_outside_plane.size() == 2)
405 : {
406 200 : std::vector<const Node *> new_plane_nodes;
407 : // A smaller Tet4 element is created, this solution is unique
408 600 : for (const auto & tet4_node_outside_plane : tet4_nodes_outside_plane)
409 : {
410 400 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
411 : mesh,
412 : new_on_plane_nodes,
413 800 : pointPairLevelSetInterception(*tet4_node_outside_plane, *tet4_nodes_inside_plane[0])));
414 : }
415 200 : auto new_elem_tet4 = std::make_unique<Tet4>();
416 200 : new_elem_tet4->set_node(0, const_cast<Node *>(new_plane_nodes[0]));
417 200 : new_elem_tet4->set_node(1, const_cast<Node *>(new_plane_nodes[1]));
418 200 : new_elem_tet4->set_node(2, const_cast<Node *>(tet4_nodes_on_plane[0]));
419 200 : new_elem_tet4->set_node(3, const_cast<Node *>(tet4_nodes_inside_plane[0]));
420 200 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
421 200 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
422 200 : }
423 208 : else if (tet4_nodes_inside_plane.size() == 2 && tet4_nodes_outside_plane.size() == 1)
424 : {
425 208 : std::vector<const Node *> new_plane_nodes;
426 : // 2 smaller TET4 elements are created
427 624 : for (const auto & tet4_node_inside_plane : tet4_nodes_inside_plane)
428 : {
429 416 : new_plane_nodes.push_back(nonDuplicateNodeCreator(
430 : mesh,
431 : new_on_plane_nodes,
432 832 : pointPairLevelSetInterception(*tet4_node_inside_plane, *tet4_nodes_outside_plane[0])));
433 : }
434 208 : std::vector<const Node *> new_elems_nodes = {tet4_nodes_inside_plane[0],
435 208 : tet4_nodes_inside_plane[1],
436 208 : new_plane_nodes[1],
437 208 : new_plane_nodes[0],
438 208 : tet4_nodes_on_plane[0]};
439 208 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
440 208 : std::vector<std::vector<const Node *>> optimized_node_list;
441 208 : MooseMeshElementConversionUtils::pyramidNodesToTetNodesDeterminer(
442 : new_elems_nodes, rotated_tet_face_indices, optimized_node_list);
443 :
444 624 : for (unsigned int i = 0; i < optimized_node_list.size(); i++)
445 : {
446 416 : auto new_elem_tet4 = std::make_unique<Tet4>();
447 416 : new_elem_tet4->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
448 416 : new_elem_tet4->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
449 416 : new_elem_tet4->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
450 416 : new_elem_tet4->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
451 416 : new_elem_tet4->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
452 416 : elems_tet4.push_back(mesh.add_elem(std::move(new_elem_tet4)));
453 416 : }
454 208 : }
455 : else
456 0 : mooseError("Unexpected scenario.");
457 :
458 18624 : mesh.elem_ptr(elem_id)->subdomain_id() = block_id_to_remove;
459 : }
460 :
461 64976 : for (auto & elem_tet4 : elems_tet4)
462 : {
463 39216 : if (new_elements_created)
464 : {
465 39152 : if (elem_tet4->volume() < 0.0)
466 : {
467 16928 : Node * temp = elem_tet4->node_ptr(0);
468 16928 : elem_tet4->set_node(0, elem_tet4->node_ptr(1));
469 16928 : elem_tet4->set_node(1, temp);
470 : }
471 : }
472 : // Find the boundary id of the new element
473 196080 : for (unsigned int i = 0; i < 4; i++)
474 : {
475 156864 : const Point & side_pt_0 = *elem_tet4->side_ptr(i)->node_ptr(0);
476 156864 : const Point & side_pt_1 = *elem_tet4->side_ptr(i)->node_ptr(1);
477 156864 : const Point & side_pt_2 = *elem_tet4->side_ptr(i)->node_ptr(2);
478 :
479 156864 : const Point side_normal = (side_pt_1 - side_pt_0).cross(side_pt_2 - side_pt_1).unit();
480 784320 : for (unsigned int j = 0; j < 4; j++)
481 : {
482 627456 : if (new_elements_created)
483 : {
484 626432 : if (MooseUtils::absoluteFuzzyEqual(side_normal * elem_side_normal_list[j], 1.0))
485 : {
486 92864 : for (const auto & side_info : elem_side_list[j])
487 : {
488 1264 : boundary_info.add_side(elem_tet4, i, side_info);
489 : }
490 : }
491 : }
492 : }
493 313728 : if (MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_0), 0.0) &&
494 246696 : MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_1), 0.0) &&
495 246696 : MooseUtils::absoluteFuzzyEqual(levelSetEvaluator(side_pt_2), 0.0))
496 : {
497 :
498 24016 : boundary_info.add_side(elem_tet4, i, _cut_face_id);
499 : }
500 : }
501 : }
502 25760 : }
503 :
504 : Real
505 5193760 : CutMeshByLevelSetGeneratorBase::levelSetEvaluator(const Point & point)
506 : {
507 5193760 : _func_params[0] = point(0);
508 5193760 : _func_params[1] = point(1);
509 5193760 : _func_params[2] = point(2);
510 5193760 : return evaluate(_func_level_set);
511 : }
|
