Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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 "Boundary2DDelaunayGenerator.h"
11 :
12 : #include "MooseMeshUtils.h"
13 : #include "CastUniquePointer.h"
14 :
15 : #include "libmesh/boundary_info.h"
16 : #include "libmesh/poly2tri_triangulator.h"
17 : #include "libmesh/mesh_triangle_holes.h"
18 : #include "libmesh/mesh_modification.h"
19 : #include "libmesh/parallel_algebra.h"
20 :
21 : registerMooseObject("MooseApp", Boundary2DDelaunayGenerator);
22 :
23 : InputParameters
24 15155 : Boundary2DDelaunayGenerator::validParams()
25 : {
26 15155 : InputParameters params = SurfaceDelaunayGeneratorBase::validParams();
27 15155 : params += FunctionParserUtils<false>::validParams();
28 :
29 30310 : params.addClassDescription(
30 : "Mesh generator that convert a 2D surface given as one or a few boundaries of a 3D mesh into "
31 : "a 2D mesh using Delaunay triangulation.");
32 60620 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
33 60620 : params.addRequiredParam<std::vector<BoundaryName>>("boundary_names", "The boundaries to be used");
34 45465 : params.addParam<std::vector<std::vector<BoundaryName>>>(
35 : "hole_boundary_names",
36 30310 : std::vector<std::vector<BoundaryName>>(),
37 : "The optional boundaries to be used as the holes in the mesh during triangulation. Note that "
38 : "this is a vector of vectors, which allows each hole to be defined as a combination of "
39 : "multiple boundaries.");
40 :
41 60620 : params.addParam<std::string>(
42 : "level_set",
43 : "Level set used to achieve more accurate reverse projection compared to interpolation.");
44 45465 : params.addParam<unsigned int>(
45 : "max_level_set_correction_iterations",
46 30310 : 3,
47 : "Maximum number of iterations to correct the nodes based on the level set function.");
48 75775 : params.addRangeCheckedParam<Real>(
49 : "max_angle_deviation",
50 30310 : 60.0,
51 : "max_angle_deviation>0 & max_angle_deviation<90",
52 : "Maximum angle deviation from the global average normal vector in the input mesh.");
53 :
54 45465 : params.addParam<BoundaryName>(
55 : "output_external_boundary_name",
56 : "",
57 : "The optional name of the external boundary of the mesh to generate. If not provided, the "
58 : "external boundary will be have a trivial id of 0.");
59 :
60 15155 : return params;
61 0 : }
62 :
63 205 : Boundary2DDelaunayGenerator::Boundary2DDelaunayGenerator(const InputParameters & parameters)
64 : : SurfaceDelaunayGeneratorBase(parameters),
65 : FunctionParserUtils<false>(parameters),
66 205 : _input(getMesh("input")),
67 410 : _boundary_names(getParam<std::vector<BoundaryName>>("boundary_names")),
68 410 : _hole_boundary_names(getParam<std::vector<std::vector<BoundaryName>>>("hole_boundary_names")),
69 205 : _max_level_set_correction_iterations(
70 410 : getParam<unsigned int>("max_level_set_correction_iterations")),
71 410 : _max_angle_deviation(getParam<Real>("max_angle_deviation")),
72 820 : _output_external_boundary_name(getParam<BoundaryName>("output_external_boundary_name"))
73 : {
74 615 : if (isParamValid("level_set"))
75 : {
76 49 : _func_level_set = std::make_shared<SymFunction>();
77 : // set FParser internal feature flags
78 49 : setParserFeatureFlags(_func_level_set);
79 147 : if (isParamValid("constant_names") && isParamValid("constant_expressions"))
80 0 : addFParserConstants(_func_level_set,
81 : getParam<std::vector<std::string>>("constant_names"),
82 : getParam<std::vector<std::string>>("constant_expressions"));
83 245 : if (_func_level_set->Parse(getParam<std::string>("level_set"), "x,y,z") >= 0)
84 0 : mooseError("Invalid function f(x,y,z)\n",
85 0 : _func_level_set,
86 : "\nin CutMeshByLevelSetGenerator ",
87 0 : name(),
88 : ".\n",
89 0 : _func_level_set->ErrorMsg());
90 :
91 49 : _func_params.resize(3);
92 : }
93 205 : }
94 :
95 : std::unique_ptr<MeshBase>
96 115 : Boundary2DDelaunayGenerator::generate()
97 : {
98 115 : std::unique_ptr<MeshBase> mesh_3d = std::move(_input);
99 :
100 : // Generate a new 2D block based on the sidesets
101 115 : const auto new_block_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh_3d);
102 : try
103 : {
104 242 : MooseMeshUtils::createSubdomainFromSidesets(
105 349 : mesh_3d, _boundary_names, new_block_id, SubdomainName(), type());
106 : }
107 4 : catch (MooseException & e)
108 : {
109 8 : if (((std::string)e.what()).compare(0, 12, "The sideset ") == 0)
110 8 : paramError("boundary_names", e.what());
111 : else
112 0 : mooseError(e.what());
113 0 : }
114 :
115 : // If holes are provided, we need to create new blocks for them too
116 111 : std::vector<subdomain_id_type> hole_block_ids;
117 128 : for (const auto & hole : _hole_boundary_names)
118 : {
119 21 : hole_block_ids.push_back(MooseMeshUtils::getNextFreeSubdomainID(*mesh_3d));
120 : try
121 : {
122 71 : MooseMeshUtils::createSubdomainFromSidesets(
123 67 : mesh_3d, hole, hole_block_ids.back(), SubdomainName(), type());
124 : }
125 4 : catch (MooseException & e)
126 : {
127 8 : if (((std::string)e.what()).compare(0, 12, "The sideset ") == 0)
128 8 : paramError("hole_boundary_names", e.what());
129 : else
130 0 : mooseError(e.what());
131 0 : }
132 : }
133 :
134 : // Create a 2D mesh form the 2D block
135 107 : auto mesh_2d = buildMeshBaseObject();
136 321 : MooseMeshUtils::convertBlockToMesh(mesh_3d, mesh_2d, {std::to_string(new_block_id)});
137 : // If holes are provided, we need to create a 2D mesh for each hole
138 107 : std::vector<std::unique_ptr<MeshBase>> hole_meshes_2d;
139 124 : for (const auto & hole_block_id : hole_block_ids)
140 : {
141 17 : hole_meshes_2d.push_back(buildMeshBaseObject());
142 68 : MooseMeshUtils::convertBlockToMesh(
143 34 : mesh_3d, hole_meshes_2d.back(), {std::to_string(hole_block_id)});
144 : }
145 :
146 : // We do not need the 3D mesh anymore
147 107 : mesh_3d->clear();
148 :
149 206 : return General2DDelaunay(mesh_2d, hole_meshes_2d);
150 223 : }
151 :
152 : Point
153 8029 : Boundary2DDelaunayGenerator::elemNormal(const Elem & elem)
154 : {
155 : mooseAssert(elem.n_vertices() == 3 || elem.n_vertices() == 4, "unsupported element type.");
156 : // Only the first three vertices are used to calculate the normal vector
157 8029 : const Point & p0 = *elem.node_ptr(0);
158 8029 : const Point & p1 = *elem.node_ptr(1);
159 8029 : const Point & p2 = *elem.node_ptr(2);
160 :
161 8029 : if (elem.n_vertices() == 4)
162 : {
163 0 : const Point & p3 = *elem.node_ptr(3);
164 0 : return ((p2 - p0).cross(p3 - p1)).unit();
165 : }
166 :
167 8029 : return ((p2 - p1).cross(p0 - p1)).unit();
168 : }
169 :
170 : Point
171 103 : Boundary2DDelaunayGenerator::meshNormal2D(const MeshBase & mesh)
172 : {
173 103 : Point mesh_norm = Point(0.0, 0.0, 0.0);
174 103 : Real mesh_area = 0.0;
175 :
176 : // Check all the elements' normal vectors
177 2787 : for (const auto & elem : mesh.active_local_element_ptr_range())
178 : {
179 2684 : const Real elem_area = elem->volume();
180 2684 : mesh_norm += elemNormal(*elem) * elem_area;
181 2684 : mesh_area += elem_area;
182 103 : }
183 103 : mesh.comm().sum(mesh_norm);
184 103 : mesh.comm().sum(mesh_area);
185 103 : mesh_norm /= mesh_area;
186 103 : return mesh_norm.unit();
187 : }
188 :
189 : Real
190 103 : Boundary2DDelaunayGenerator::meshNormalDeviation2D(const MeshBase & mesh, const Point & global_norm)
191 : {
192 103 : Real max_deviation(0.0);
193 : // Check all the elements' deviation from the global normal vector
194 2787 : for (const auto & elem : mesh.active_local_element_ptr_range())
195 : {
196 2684 : const Real elem_deviation = std::acos(global_norm * elemNormal(*elem)) / M_PI * 180.0;
197 2684 : max_deviation = std::max(max_deviation, elem_deviation);
198 103 : }
199 103 : mesh.comm().max(max_deviation);
200 103 : return max_deviation;
201 : }
202 :
203 : Real
204 12139 : Boundary2DDelaunayGenerator::levelSetEvaluator(const Point & point)
205 : {
206 48556 : return evaluate(_func_level_set, std::vector<Real>({point(0), point(1), point(2), 0}));
207 : }
208 :
209 : void
210 1320 : Boundary2DDelaunayGenerator::levelSetCorrection(Node & node)
211 : {
212 : // Based on the given level set, we try to move the node in its normal direction
213 1320 : const Real diff = libMesh::TOLERANCE * 10.0; // A small value to perturb the node
214 1320 : const Real original_eval = levelSetEvaluator(node);
215 1320 : const Real xp_eval = levelSetEvaluator(node + Point(diff, 0.0, 0.0));
216 1320 : const Real yp_eval = levelSetEvaluator(node + Point(0.0, diff, 0.0));
217 1320 : const Real zp_eval = levelSetEvaluator(node + Point(0.0, 0.0, diff));
218 1320 : const Real xm_eval = levelSetEvaluator(node - Point(diff, 0.0, 0.0));
219 1320 : const Real ym_eval = levelSetEvaluator(node - Point(0.0, diff, 0.0));
220 1320 : const Real zm_eval = levelSetEvaluator(node - Point(0.0, 0.0, diff));
221 1320 : const Point grad = Point((xp_eval - xm_eval) / (2.0 * diff),
222 1320 : (yp_eval - ym_eval) / (2.0 * diff),
223 1320 : (zp_eval - zm_eval) / (2.0 * diff));
224 1320 : const Real xyz_diff = -original_eval / grad.contract(grad);
225 1320 : node(0) += xyz_diff * grad(0);
226 1320 : node(1) += xyz_diff * grad(1);
227 1320 : node(2) += xyz_diff * grad(2);
228 1320 : }
229 :
230 : std::unique_ptr<MeshBase>
231 107 : Boundary2DDelaunayGenerator::General2DDelaunay(
232 : std::unique_ptr<MeshBase> & mesh_2d, std::vector<std::unique_ptr<MeshBase>> & hole_meshes_2d)
233 : {
234 : // If a level set is provided, we need to check if the nodes in the original 2D mesh match the
235 : // level set
236 107 : if (_func_level_set)
237 : {
238 949 : for (const auto & node : mesh_2d->node_ptr_range())
239 : {
240 904 : if (std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE)
241 : {
242 8 : paramError("level_set",
243 : "The level set function does not match the nodes in the given boundary of the "
244 : "input mesh.");
245 : }
246 45 : }
247 : }
248 :
249 : // Easier to work with a TRI3 mesh
250 : // all_tri() also prepares the mesh for use
251 103 : mesh_2d->prepare_for_use();
252 103 : MeshTools::Modification::all_tri(*mesh_2d);
253 103 : const auto mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(*mesh_2d);
254 3499 : for (const auto & elem : mesh_2d->active_element_ptr_range())
255 13584 : for (const auto & i_side : elem->side_index_range())
256 10188 : if (elem->neighbor_ptr(i_side) == nullptr)
257 1675 : mesh_2d->get_boundary_info().add_side(elem, i_side, mesh_2d_ext_bdry);
258 :
259 : // Create a clone of the 2D mesh to be used for the 1D mesh generation
260 103 : auto mesh_2d_dummy = dynamic_pointer_cast<MeshBase>(mesh_2d->clone());
261 : // Generate a new 1D block based on the external boundary
262 103 : const auto new_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*mesh_2d_dummy);
263 :
264 412 : MooseMeshUtils::createSubdomainFromSidesets(
265 309 : mesh_2d_dummy, {std::to_string(mesh_2d_ext_bdry)}, new_block_id_1d, SubdomainName(), type());
266 :
267 : // Create a 1D mesh form the 1D block
268 103 : auto mesh_1d = buildMeshBaseObject();
269 309 : MooseMeshUtils::convertBlockToMesh(mesh_2d_dummy, mesh_1d, {std::to_string(new_block_id_1d)});
270 103 : mesh_2d_dummy->clear();
271 :
272 : // If we have holes, we need to create a 1D mesh for each hole
273 103 : std::vector<std::unique_ptr<MeshBase>> hole_meshes_1d;
274 116 : for (auto & hole_mesh_2d : hole_meshes_2d)
275 : {
276 : // As we do not need these holes for reverse projection, we do not need to convert them to TRI3
277 : // meshes, but we still need to create a 1D mesh for each hole
278 13 : hole_mesh_2d->find_neighbors();
279 13 : const auto hole_mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(*hole_mesh_2d);
280 65 : for (const auto & elem : hole_mesh_2d->active_element_ptr_range())
281 260 : for (const auto & i_side : elem->side_index_range())
282 208 : if (elem->neighbor_ptr(i_side) == nullptr)
283 117 : hole_mesh_2d->get_boundary_info().add_side(elem, i_side, mesh_2d_ext_bdry);
284 13 : const auto new_hole_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*hole_mesh_2d);
285 52 : MooseMeshUtils::createSubdomainFromSidesets(hole_mesh_2d,
286 13 : {std::to_string(hole_mesh_2d_ext_bdry)},
287 : new_hole_block_id_1d,
288 26 : SubdomainName(),
289 13 : type());
290 : // Create a 1D mesh form the 1D block
291 13 : hole_meshes_1d.push_back(buildMeshBaseObject());
292 52 : MooseMeshUtils::convertBlockToMesh(
293 26 : hole_mesh_2d, hole_meshes_1d.back(), {std::to_string(new_hole_block_id_1d)});
294 13 : hole_mesh_2d->clear();
295 : }
296 :
297 : // Find centroid of the 2D mesh
298 103 : const Point centroid = MooseMeshUtils::meshCentroidCalculator(*mesh_2d);
299 : // calculate an average normal vector of the 2D mesh
300 103 : const Point mesh_norm = meshNormal2D(*mesh_2d);
301 : // Check the deviation of the mesh normal vector from the global average normal vector
302 103 : if (meshNormalDeviation2D(*mesh_2d, mesh_norm) > _max_angle_deviation)
303 4 : paramError("boundary_names",
304 : "The normal vector of some elements in the 2D mesh deviates too much from the "
305 4 : "global average normal vector. The maximum deviation is " +
306 8 : std::to_string(_max_angle_deviation) +
307 : ". Consider dividing the boundary into several parts to "
308 : "reduce the angle deviation.");
309 :
310 : // Move both 2d and 1d meshes to the centroid of the 2D mesh
311 99 : MeshTools::Modification::translate(*mesh_1d, -centroid(0), -centroid(1), -centroid(2));
312 99 : MeshTools::Modification::translate(*mesh_2d, -centroid(0), -centroid(1), -centroid(2));
313 : // Alsop need to translate the 1D hole meshes if applicable
314 108 : for (auto & hole_mesh_1d : hole_meshes_1d)
315 9 : MeshTools::Modification::translate(*hole_mesh_1d, -centroid(0), -centroid(1), -centroid(2));
316 :
317 : // Calculate the Euler angles to rotate the meshes so that the 2D mesh is close to the XY plane
318 : // (i.e., the normal vector of the 2D mesh is aligned with the Z axis)
319 99 : const Real theta = std::acos(mesh_norm(2)) / M_PI * 180.0;
320 : const Real phi =
321 99 : (MooseUtils::absoluteFuzzyLessThan(mesh_norm(2), 1.0) ? std::atan2(mesh_norm(1), mesh_norm(0))
322 99 : : 0.0) /
323 99 : M_PI * 180.0;
324 99 : MeshTools::Modification::rotate(*mesh_1d, 90.0 - phi, theta, 0.0);
325 99 : MeshTools::Modification::rotate(*mesh_2d, 90.0 - phi, theta, 0.0);
326 : // Also rotate the 1D hole meshes if applicable
327 108 : for (auto & hole_mesh_1d : hole_meshes_1d)
328 9 : MeshTools::Modification::rotate(*hole_mesh_1d, 90.0 - phi, theta, 0.0);
329 :
330 : // Clone the 2D mesh to be used for reverse projection later
331 99 : auto mesh_2d_xyz = dynamic_pointer_cast<MeshBase>(mesh_2d->clone());
332 :
333 : // Project the 2D mesh to the XY plane so that XYDelaunay can be used
334 2538 : for (const auto & node : mesh_2d->node_ptr_range())
335 2538 : (*node)(2) = 0;
336 : // Project the 1D mesh to the XY plane as well
337 1575 : for (const auto & node : mesh_1d->node_ptr_range())
338 1575 : (*node)(2) = 0;
339 : // Also project the 1D hole meshes to the XY plane if applicable
340 108 : for (auto & hole_mesh_1d : hole_meshes_1d)
341 : {
342 81 : for (const auto & node : hole_mesh_1d->node_ptr_range())
343 81 : (*node)(2) = 0;
344 : }
345 :
346 99 : std::vector<libMesh::TriangulatorInterface::MeshedHole> meshed_holes;
347 99 : std::vector<libMesh::TriangulatorInterface::Hole *> triangulator_hole_ptrs(hole_meshes_1d.size());
348 99 : meshed_holes.reserve(hole_meshes_1d.size());
349 108 : for (auto hole_i : index_range(hole_meshes_1d))
350 : {
351 9 : hole_meshes_1d[hole_i]->prepare_for_use();
352 9 : meshed_holes.emplace_back(*hole_meshes_1d[hole_i]);
353 9 : triangulator_hole_ptrs[hole_i] = &meshed_holes.back();
354 : }
355 :
356 : // Finally, triangulation
357 : std::unique_ptr<UnstructuredMesh> mesh =
358 99 : dynamic_pointer_cast<UnstructuredMesh>(std::move(mesh_1d));
359 :
360 99 : Poly2TriTriangulator poly2tri(*mesh);
361 99 : poly2tri.triangulation_type() = TriangulatorInterface::PSLG;
362 :
363 99 : poly2tri.set_interpolate_boundary_points(0);
364 99 : poly2tri.set_refine_boundary_allowed(false);
365 99 : poly2tri.set_verify_hole_boundaries(false);
366 99 : poly2tri.desired_area() = 0;
367 99 : poly2tri.minimum_angle() = 0; // Not yet supported
368 99 : poly2tri.smooth_after_generating() = false;
369 99 : if (!triangulator_hole_ptrs.empty())
370 9 : poly2tri.attach_hole_list(&triangulator_hole_ptrs);
371 : // Future TODO: correct the area function based on the local normal vector
372 99 : if (_use_auto_area_func)
373 90 : poly2tri.set_auto_area_function(this->comm(),
374 90 : _auto_area_function_num_points,
375 90 : _auto_area_function_power,
376 90 : _auto_area_func_default_size,
377 90 : _auto_area_func_default_size_dist);
378 99 : poly2tri.triangulate();
379 :
380 : // Reverse the projection based on the original 2D mesh
381 2760 : for (const auto & node : mesh->node_ptr_range())
382 : {
383 2661 : bool node_mod = false;
384 : // Try to find the element in mesh_2d that contains the new node
385 47083 : for (const auto & elem : mesh_2d->active_element_ptr_range())
386 : {
387 47083 : if (elem->contains_point(Point((*node)(0), (*node)(1), 0.0)))
388 : {
389 : // Element id
390 2661 : const auto elem_id = elem->id();
391 : // element in xyz_in_xyz
392 2661 : const Elem & elem_xyz = *mesh_2d_xyz->elem_ptr(elem_id);
393 :
394 2661 : const Point elem_normal = elemNormal(elem_xyz);
395 2661 : const Point & elem_p = *mesh_2d_xyz->elem_ptr(elem_id)->node_ptr(0);
396 :
397 : // if the x and y values of the node is the same as the elem_p's first node, we can just
398 : // move it to that node's position
399 3561 : if (MooseUtils::absoluteFuzzyEqual((*node)(0), elem_p(0)) &&
400 900 : MooseUtils::absoluteFuzzyEqual((*node)(1), elem_p(1)))
401 : {
402 189 : (*node)(2) = elem_p(2);
403 189 : node_mod = true;
404 189 : break;
405 : }
406 : // Otherwise, we need to find a position inside the 2D element
407 : // It has the same x and y coordinates as the node in the projected mesh;
408 2472 : (*node)(2) = elem_p(2) - (((*node)(0) - elem_p(0)) * elem_normal(0) +
409 2472 : ((*node)(1) - elem_p(1)) * elem_normal(1)) /
410 2472 : elem_normal(2);
411 2472 : node_mod = true;
412 2472 : break;
413 : }
414 2661 : }
415 2661 : if (!node_mod)
416 0 : mooseError("Node not found in mesh_in_xy");
417 99 : }
418 :
419 : // Rotate the mesh back
420 99 : MeshTools::Modification::rotate(*mesh, 0.0, -theta, phi - 90.0);
421 : // Translate the mesh back
422 99 : MeshTools::Modification::translate(*mesh, centroid(0), centroid(1), centroid(2));
423 :
424 : // Correct the nodes based on the level set function
425 99 : if (_func_level_set)
426 : {
427 1160 : for (const auto & node : mesh->node_ptr_range())
428 : {
429 1115 : unsigned int iter_ct = 0;
430 4430 : while (iter_ct < _max_level_set_correction_iterations &&
431 1995 : std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE * libMesh::TOLERANCE)
432 : {
433 1320 : levelSetCorrection(*node);
434 1320 : ++iter_ct;
435 : }
436 45 : }
437 : }
438 :
439 : // Assign the external boundary name if provided
440 99 : if (!_output_external_boundary_name.empty())
441 : {
442 : const std::vector<BoundaryID> output_boundary_id =
443 54 : MooseMeshUtils::getBoundaryIDs(*mesh, {_output_external_boundary_name}, true);
444 :
445 18 : libMesh::MeshTools::Modification::change_boundary_id(*mesh, 0, output_boundary_id[0]);
446 18 : mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = _output_external_boundary_name;
447 18 : }
448 :
449 99 : mesh->set_isnt_prepared();
450 :
451 198 : return mesh;
452 349 : }
|
