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 "SurfaceSubdomainsDelaunayRemesher.h"
11 :
12 : #include "MooseMeshUtils.h"
13 : #include "CastUniquePointer.h"
14 :
15 : #include "libmesh/boundary_info.h"
16 : #include "libmesh/mesh_triangle_holes.h"
17 : #include "libmesh/mesh_modification.h"
18 : #include "libmesh/parallel_algebra.h"
19 : #include "libmesh/poly2tri_triangulator.h"
20 : #include "libmesh/unstructured_mesh.h"
21 :
22 : #include "libmesh/mesh_base.h"
23 : #include "libmesh/replicated_mesh.h"
24 :
25 : registerMooseObject("MooseApp", SurfaceSubdomainsDelaunayRemesher);
26 :
27 : InputParameters
28 3125 : SurfaceSubdomainsDelaunayRemesher::validParams()
29 : {
30 3125 : InputParameters params = SurfaceDelaunayGeneratorBase::validParams();
31 3125 : params += LevelSetMeshingHelper::validParams();
32 12500 : params.renameParameterGroup("Parsed expression advanced", "Level set shape parsed expression");
33 :
34 6250 : params.addClassDescription(
35 : "Mesh generator that re-meshes a 2D surface mesh given as one or more subdomains into "
36 : "a 2D surface mesh using Delaunay triangulation.");
37 :
38 : // Definition of the region to re-mesh
39 12500 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
40 12500 : params.addParam<std::vector<std::vector<SubdomainName>>>(
41 : "subdomain_names", {}, "The surface mesh subdomains to be re-meshed");
42 12500 : params.addParam<std::vector<SubdomainName>>(
43 : "exclude_subdomain_names", {}, "The surface mesh subdomains that should not be re-meshed");
44 :
45 : // Surface re-meshing parameters
46 18750 : params.addRangeCheckedParam<Real>(
47 : "max_edge_length",
48 : "max_edge_length>0",
49 : "Maximum length of an edge in the 1D meshes around each subdomain. Only a single value is "
50 : "currently allowed as the boundary refinement must be consistent between all subdomains.");
51 : // NOTE: we could make this a per-included-subdomain parameter, we would just need to identify
52 : // interfaces between each subdomain and use the minimum (or average or maximum, just needs to
53 : // be consistent) edge length on those interfaces.
54 12500 : params.addParam<std::vector<unsigned int>>(
55 : "interpolate_boundaries",
56 : {1},
57 : "Ratio of points to add to the boundaries. Can be set to a single value for all groups at "
58 : "once, or specified individually. A single value is recommended as the boundary refinement "
59 : "must be consistent between all subdomains.");
60 12500 : params.addParam<std::vector<Real>>(
61 : "desired_areas",
62 : {0},
63 : "Target element size when triangulating projection of the subdomain group. Can be set to a "
64 : "single value for all groups at once, or specified individually. Default of 0 means no "
65 : "constraint.");
66 12500 : params.addParamNamesToGroup("max_edge_length interpolate_boundaries desired_areas",
67 : "Delaunay triangulation");
68 :
69 : // Parameters for stitching meshes at the end
70 9375 : params.addParam<bool>("avoid_merging_subdomains",
71 6250 : false,
72 : "Whether to prevent merging subdomains by offsetting ids. The first mesh "
73 : "in the input will keep the same subdomains ids, the others will have "
74 : "offsets. All subdomain names will remain valid");
75 9375 : params.addParam<bool>("clear_stitching_boundaries",
76 6250 : true,
77 : "Whether to clear the boundaries between the subdomains being stitched "
78 : "together after they were re-meshed with triangles");
79 12500 : MooseEnum algorithm("BINARY EXHAUSTIVE", "BINARY");
80 12500 : params.addParam<MooseEnum>(
81 : "stitching_algorithm",
82 : algorithm,
83 : "Control the use of binary search for the nodes of the stitched surfaces.");
84 9375 : params.addParam<bool>(
85 6250 : "verbose_stitching", false, "Whether mesh stitching should have verbose output.");
86 9375 : params.addParamNamesToGroup(
87 : "avoid_merging_subdomains clear_stitching_boundaries stitching_algorithm verbose_stitching",
88 : "Stitching triangularized meshes");
89 :
90 6250 : return params;
91 3125 : }
92 :
93 32 : SurfaceSubdomainsDelaunayRemesher::SurfaceSubdomainsDelaunayRemesher(
94 32 : const InputParameters & parameters)
95 : : SurfaceDelaunayGeneratorBase(parameters),
96 : LevelSetMeshingHelper(parameters),
97 32 : _input(getMesh("input")),
98 64 : _subdomain_names(getParam<std::vector<std::vector<SubdomainName>>>("subdomain_names")),
99 64 : _interpolate_boundaries(getParam<std::vector<unsigned int>>("interpolate_boundaries")),
100 128 : _desired_areas(getParam<std::vector<Real>>("desired_areas"))
101 : {
102 112 : if (isParamSetByUser("subdomain_names") && isParamSetByUser("exclude_subdomain_names"))
103 0 : paramError("exclude_subdomain_names",
104 : "Excluding subdomain names is only to be set when 'subdomain_names' is not set");
105 32 : }
106 :
107 : std::unique_ptr<MeshBase>
108 32 : SurfaceSubdomainsDelaunayRemesher::generate()
109 : {
110 32 : std::unique_ptr<MeshBase> mesh_3d = std::move(_input);
111 :
112 : // Select subdomain names if it has not been selected by the user
113 : // We form 1 group per subdomain by default (easiest to mesh)
114 32 : if (_subdomain_names.empty() || _subdomain_names[0].empty())
115 : {
116 24 : if (_subdomain_names.size())
117 0 : _subdomain_names.erase(_subdomain_names.begin());
118 :
119 : // Get all subdomains from the mesh
120 24 : std::set<subdomain_id_type> sub_ids;
121 24 : mesh_3d->subdomain_ids(sub_ids);
122 :
123 : // Copy only the ones that are not excluded
124 48 : const auto & excluded = getParam<std::vector<SubdomainName>>("exclude_subdomain_names");
125 288 : for (const auto & sub_id : sub_ids)
126 : {
127 264 : const auto & sn = mesh_3d->subdomain_name(sub_id) == "" ? std::to_string(sub_id)
128 264 : : mesh_3d->subdomain_name(sub_id);
129 264 : if (excluded.empty() || std::find(excluded.begin(), excluded.end(), sn) == excluded.end())
130 768 : _subdomain_names.push_back({sn});
131 264 : }
132 24 : }
133 32 : _num_groups = _subdomain_names.size();
134 32 : if (_interpolate_boundaries.size() == 1)
135 32 : _interpolate_boundaries.resize(_num_groups, _interpolate_boundaries[0]);
136 32 : if (_desired_areas.size() == 1)
137 32 : _desired_areas.resize(_num_groups, _desired_areas[0]);
138 32 : if (_interpolate_boundaries.size() != _num_groups)
139 0 : paramError("interpolate_boundaries",
140 : "Should be the same size as 'subdomain_names' or of size 1");
141 32 : if (_desired_areas.size() != _num_groups)
142 0 : paramError("desired_areas", "Should be the same size as 'subdomain_names' or of size 1");
143 :
144 : // Re-create surface meshes from each of the subdomains
145 32 : std::vector<std::unique_ptr<UnstructuredMesh>> remeshed_2d;
146 :
147 304 : for (const auto i : make_range(_num_groups))
148 : {
149 272 : auto mesh_2d = buildReplicatedMesh();
150 272 : MooseMeshUtils::convertBlockToMesh(*mesh_3d, *mesh_2d, _subdomain_names[i]);
151 :
152 : // Mesh the subdomains by groups
153 : // TODO: holes for each subdomain are not currently supported
154 272 : std::vector<std::unique_ptr<UnstructuredMesh>> no_holes = {};
155 272 : auto new_mesh = General2DDelaunay(*mesh_2d, /*holes*/ no_holes, i);
156 :
157 : // Keep the remeshed subdomains 2D meshes
158 272 : remeshed_2d.push_back(std::move(new_mesh));
159 272 : }
160 :
161 : // TODO: can we gain efficiency by using different boundary ids for each stitch?
162 32 : const auto primary_bcid = 1;
163 : // Prevent deleting a legitimate ID if we start conserving boundary IDs from the input mesh
164 : // to the output mesh
165 32 : const auto starting_stitching_id = MooseMeshUtils::getNextFreeBoundaryID(*mesh_3d);
166 32 : auto paired_bcid = starting_stitching_id;
167 64 : const auto verbose_stitching = getParam<bool>("verbose_stitching");
168 64 : const auto use_binary_search = getParam<MooseEnum>("stitching_algorithm") == "BINARY";
169 :
170 : // Stitch all the parts to the first one
171 32 : std::unique_ptr<UnstructuredMesh> full_mesh;
172 32 : full_mesh = std::move(remeshed_2d[0]);
173 :
174 : // Build a subdomain map
175 32 : auto & main_subdomain_map = full_mesh->set_subdomain_name_map();
176 :
177 304 : for (auto remesh_i : index_range(remeshed_2d))
178 : {
179 272 : if (remesh_i == 0)
180 32 : continue;
181 :
182 240 : UnstructuredMesh & remeshed = dynamic_cast<UnstructuredMesh &>(*remeshed_2d[remesh_i]);
183 :
184 : // NOTE: we cannot use MeshedHole because the meshes are no longer in the XY plane
185 : // Potential optimization: instead of using the full mesh outer boundary for stitching
186 : // use what is done in the XYDelaunayGenerator and re-build only the boundary near the
187 : // mesh that we are stitching
188 : // Potential optimization: we could keep the external boundaries from the triangulation phase
189 :
190 : // Create the boundary outside the remeshed mesh
191 240 : bool rem_has_external_bdy = false;
192 240 : MooseMeshUtils::addExternalBoundary(remeshed, paired_bcid, rem_has_external_bdy);
193 : mooseAssert(rem_has_external_bdy, "Subdomain mesh should have an external boundary");
194 :
195 : // We need to re-create the boundary outside the parent mesh because we are clearing
196 : // it on every stitch
197 240 : bool base_has_external_bdy = false;
198 240 : MooseMeshUtils::addExternalBoundary(*full_mesh, primary_bcid, base_has_external_bdy);
199 240 : if (!base_has_external_bdy)
200 0 : mooseWarning(
201 : "Full mesh should have an external boundary. This could occur if the surface mesh is "
202 : "disconnected in several parts and a part's surface mesh has just been fully remeshed");
203 :
204 : // Retrieve subdomain name map from the mesh to be stitched and insert it into the main
205 : // subdomain map
206 240 : const auto & increment_subdomain_map = remeshed.get_subdomain_name_map();
207 240 : main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());
208 :
209 240 : full_mesh->stitch_meshes(remeshed,
210 : primary_bcid,
211 : paired_bcid,
212 : TOLERANCE,
213 240 : /*clear_stitched_bcids*/ getParam<bool>("clear_stitching_boundaries"),
214 : verbose_stitching,
215 : use_binary_search,
216 : /*enforce_all_nodes_match_on_boundaries*/ false,
217 : /*merge_boundary_nodes_all_or_nothing*/ false,
218 720 : /*remap_subdomain_ids*/ !getParam<bool>("avoid_merging_subdomains"));
219 :
220 : // If we want to keep track of the stitching boundaries
221 240 : paired_bcid++;
222 : // TODO: when implementing hole meshes, we might want to also stitch the hole meshes
223 : }
224 :
225 : // We do not need the 3D mesh anymore
226 32 : mesh_3d->clear();
227 : // The stitching boundary should be removed by the stitcher, but this does not hurt
228 96 : if (getParam<bool>("clear_stitching_boundaries"))
229 272 : for (const auto bcid : make_range(starting_stitching_id, paired_bcid))
230 240 : full_mesh->get_boundary_info().remove_id(bcid);
231 : // Mesh is not prepared after stitching (TODO: fix this)
232 32 : full_mesh->unset_is_prepared();
233 :
234 64 : return full_mesh;
235 288 : }
236 :
237 : std::unique_ptr<UnstructuredMesh>
238 272 : SurfaceSubdomainsDelaunayRemesher::General2DDelaunay(
239 : UnstructuredMesh & mesh_2d,
240 : std::vector<std::unique_ptr<UnstructuredMesh>> & hole_meshes_2d,
241 : unsigned int group_i)
242 : {
243 272 : if (_verbose)
244 : {
245 272 : if (!mesh_2d.preparation().has_cached_elem_data)
246 272 : mesh_2d.cache_elem_data();
247 272 : _console << "Re-meshing mesh\n " << mesh_2d << std::endl;
248 1360 : _console << "with subdomains " << Moose::stringify(mesh_2d.get_subdomain_name_map())
249 272 : << std::endl;
250 272 : if (hole_meshes_2d.size())
251 0 : _console << "With " << hole_meshes_2d.size() << " holes:" << std::endl;
252 272 : for (const auto & hole_m : hole_meshes_2d)
253 : {
254 0 : if (!hole_m->preparation().has_cached_elem_data)
255 0 : hole_m->cache_elem_data();
256 0 : _console << "Hole subdomains " << Moose::stringify(hole_m->get_subdomain_name_map())
257 0 : << std::endl;
258 0 : _console << *hole_m << std::endl;
259 : }
260 : }
261 : // If a level set is provided, we need to check if the nodes in the original 2D mesh match the
262 : // level set
263 272 : if (_func_level_set)
264 : {
265 1872 : for (const auto & node : mesh_2d.node_ptr_range())
266 : {
267 1616 : if (std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE)
268 : {
269 0 : paramError("level_set",
270 : "The level set function does not match the nodes in the given boundary of the "
271 0 : "input mesh. Level set evaluates at: " +
272 0 : std::to_string(levelSetEvaluator(*node)) + " at node: " + node->get_info());
273 : }
274 256 : }
275 : }
276 :
277 : // Create the external boundary of the 2D mesh
278 : // Easier to work with a TRI3 mesh
279 : // all_tri() also prepares the mesh for use
280 272 : mesh_2d.prepare_for_use();
281 272 : bool has_external_side = false;
282 272 : MeshTools::Modification::all_tri(mesh_2d);
283 272 : const auto mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(mesh_2d);
284 1456 : for (const auto elem : mesh_2d.active_element_ptr_range())
285 4736 : for (const auto i_side : elem->side_index_range())
286 3552 : if (elem->neighbor_ptr(i_side) == nullptr)
287 : {
288 1728 : mesh_2d.get_boundary_info().add_side(elem, i_side, mesh_2d_ext_bdry);
289 1728 : has_external_side = true;
290 272 : }
291 :
292 : // Create a clone of the 2D mesh to be used for the 1D mesh generation
293 272 : auto mesh_2d_dummy = mesh_2d.clone();
294 : // Generate a new 1D block based on the external boundary
295 272 : const auto new_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*mesh_2d_dummy);
296 :
297 272 : if (has_external_side)
298 1360 : MooseMeshUtils::createSubdomainFromSidesets(*mesh_2d_dummy,
299 272 : {std::to_string(mesh_2d_ext_bdry)},
300 : new_block_id_1d,
301 544 : SubdomainName(),
302 272 : type());
303 :
304 : // Create a 1D mesh from the 1D block
305 272 : auto mesh_1d = buildMeshBaseObject();
306 272 : if (has_external_side)
307 816 : MooseMeshUtils::convertBlockToMesh(*mesh_2d_dummy, *mesh_1d, {std::to_string(new_block_id_1d)});
308 272 : mesh_2d_dummy->clear();
309 :
310 : // Add more nodes in the 1D mesh
311 : // We need to do this BEFORE the mesh is projected to avoid distortion in the node distances
312 : // The distortion is not the same for each subdomain, so we would end up with a non-conformal mesh
313 816 : if (isParamValid("max_edge_length"))
314 : {
315 : // We need to extract the points in the order of the boundary. This utility will build the
316 : // boundary as a loop
317 : // Add the nodeset from the sideset
318 : // NOTE: this is redoing the 1D mesh from the boundary
319 176 : mesh_2d.get_boundary_info().build_node_list_from_side_list({mesh_2d_ext_bdry});
320 88 : const auto boundary_1d = MooseMeshUtils::buildLoopBoundaryOf2DMesh(mesh_2d, mesh_2d_ext_bdry);
321 : // Extract the points
322 88 : std::vector<Point> mesh_1d_points;
323 88 : mesh_1d_points.reserve(boundary_1d->n_nodes());
324 648 : for (const auto & node : boundary_1d->node_ptr_range())
325 648 : mesh_1d_points.push_back(*node);
326 :
327 : // Re-create the 1D mesh with a polyline with a maximum edge length
328 88 : mesh_1d->clear();
329 : // we add a tiny offset to avoid roundoff differences
330 176 : MooseMeshUtils::buildPolyLineMesh(
331 352 : *mesh_1d, mesh_1d_points, true, "", "", getParam<Real>("max_edge_length") + 1e-8);
332 88 : }
333 :
334 : // Find centroid of the 2D mesh
335 272 : const Point centroid = MooseMeshUtils::meshCentroidCalculator(mesh_2d);
336 : // calculate an average normal vector of the 2D mesh
337 272 : const Point mesh_norm = meshNormal2D(mesh_2d);
338 : // Check the deviation of the mesh normal vector from the global average normal vector
339 272 : if (meshNormalDeviation2D(mesh_2d, mesh_norm) > _max_angle_deviation)
340 0 : paramError("subdomain_names",
341 : "The normal vector of some elements in the 2D mesh deviates too much from the "
342 0 : "global average normal vector. The maximum deviation found / allowed is " +
343 0 : std::to_string(meshNormalDeviation2D(mesh_2d, mesh_norm)) + " / " +
344 0 : std::to_string(_max_angle_deviation) +
345 : ". Consider dividing the boundary into several parts to "
346 : "reduce the angle deviation.");
347 :
348 : // Move both 2d and 1d meshes to the centroid of the 2D mesh
349 272 : MeshTools::Modification::translate(*mesh_1d, -centroid(0), -centroid(1), -centroid(2));
350 272 : MeshTools::Modification::translate(mesh_2d, -centroid(0), -centroid(1), -centroid(2));
351 :
352 : // Calculate the Euler angles to rotate the meshes so that the 2D mesh is close to the XY plane
353 : // (i.e., the normal vector of the 2D mesh is aligned with the Z axis)
354 272 : const Real theta = std::acos(mesh_norm(2)) / M_PI * 180.0;
355 : const Real phi =
356 272 : (MooseUtils::absoluteFuzzyLessThan(mesh_norm(2), 1.0) ? std::atan2(mesh_norm(1), mesh_norm(0))
357 272 : : 0.0) /
358 272 : M_PI * 180.0;
359 272 : MeshTools::Modification::rotate(*mesh_1d, 90.0 - phi, theta, 0.0);
360 272 : MeshTools::Modification::rotate(mesh_2d, 90.0 - phi, theta, 0.0);
361 :
362 : // Clone the 2D mesh to be used for reverse projection later
363 272 : auto mesh_2d_xyz = dynamic_pointer_cast<MeshBase>(mesh_2d.clone());
364 :
365 : // Project the 2D mesh to the XY plane so that XYDelaunay can be used
366 2000 : for (const auto & node : mesh_2d.node_ptr_range())
367 2000 : (*node)(2) = 0;
368 : // Project the 1D mesh to the XY plane as well
369 2560 : for (const auto & node : mesh_1d->node_ptr_range())
370 2560 : (*node)(2) = 0;
371 :
372 : // Finally, triangulation
373 : std::unique_ptr<UnstructuredMesh> mesh =
374 272 : dynamic_pointer_cast<UnstructuredMesh>(std::move(mesh_1d));
375 :
376 272 : Poly2TriTriangulator poly2tri(*mesh);
377 272 : poly2tri.triangulation_type() = TriangulatorInterface::PSLG;
378 :
379 272 : poly2tri.set_interpolate_boundary_points(_interpolate_boundaries[group_i]);
380 272 : poly2tri.set_verify_hole_boundaries(false);
381 272 : poly2tri.desired_area() = _desired_areas[group_i];
382 272 : poly2tri.minimum_angle() = 0; // Not yet supported
383 272 : poly2tri.smooth_after_generating() = false;
384 : // Future TODO: correct the area function based on the local normal vector
385 272 : if (_use_auto_area_func)
386 0 : poly2tri.set_auto_area_function(this->comm(),
387 0 : _auto_area_function_num_points,
388 0 : _auto_area_function_power,
389 0 : _auto_area_func_default_size,
390 0 : _auto_area_func_default_size_dist);
391 272 : poly2tri.triangulate();
392 :
393 : // Reverse the projection based on the original 2D mesh
394 5800 : for (const auto & node : mesh->node_ptr_range())
395 : {
396 5528 : bool node_mod = false;
397 : // Try to find the element in mesh_2d that contains the new node
398 14944 : for (const auto & elem : mesh_2d.active_element_ptr_range())
399 : {
400 14944 : if (elem->contains_point(Point((*node)(0), (*node)(1), 0.0)))
401 : {
402 : // Element id
403 5528 : const auto elem_id = elem->id();
404 : // element in xyz_in_xyz
405 5528 : const Elem & elem_xyz = *mesh_2d_xyz->elem_ptr(elem_id);
406 :
407 5528 : const Point elem_normal = elemNormal(elem_xyz);
408 5528 : const auto & elem_p = mesh_2d_xyz->elem_ptr(elem_id)->point(0);
409 :
410 : // if the x and y values of the node is the same as the elem_p's first node, we can just
411 : // move it to that node's position
412 5912 : if (MooseUtils::absoluteFuzzyEqual((*node)(0), elem_p(0)) &&
413 384 : MooseUtils::absoluteFuzzyEqual((*node)(1), elem_p(1)))
414 : {
415 296 : (*node)(2) = elem_p(2);
416 296 : node_mod = true;
417 296 : break;
418 : }
419 : // Otherwise, we need to find a position inside the 2D element
420 : // It has the same x and y coordinates as the node in the projected mesh;
421 5232 : (*node)(2) = elem_p(2) - (((*node)(0) - elem_p(0)) * elem_normal(0) +
422 5232 : ((*node)(1) - elem_p(1)) * elem_normal(1)) /
423 5232 : elem_normal(2);
424 5232 : node_mod = true;
425 5232 : break;
426 : }
427 5528 : }
428 5528 : if (!node_mod)
429 0 : mooseError("Node not found in mesh_in_xy");
430 272 : }
431 :
432 : // Rotate the mesh back
433 272 : MeshTools::Modification::rotate(*mesh, 0.0, -theta, phi - 90.0);
434 : // Translate the mesh back
435 272 : MeshTools::Modification::translate(*mesh, centroid(0), centroid(1), centroid(2));
436 :
437 : // Correct the nodes based on the level set function
438 272 : if (_func_level_set)
439 : {
440 5560 : for (const auto & node : mesh->node_ptr_range())
441 : {
442 5304 : unsigned int iter_ct = 0;
443 29048 : while (iter_ct < _max_level_set_correction_iterations &&
444 12680 : std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE * libMesh::TOLERANCE)
445 : {
446 11064 : levelSetCorrection(*node);
447 11064 : ++iter_ct;
448 : }
449 256 : }
450 : }
451 : // Give the old subdomain to all elements
452 272 : const auto common_id = mesh_2d.get_mesh_subdomains().begin();
453 5648 : for (auto & elem : mesh->active_element_ptr_range())
454 5648 : elem->subdomain_id() = *common_id;
455 :
456 : // Pass the subdomain names
457 272 : mesh->set_subdomain_name_map() = mesh_2d.get_subdomain_name_map();
458 : // Remove the boundaries, they get re-added (with a known ID) when stitching the pieces together
459 272 : mesh->get_boundary_info().remove_id(mesh_2d_ext_bdry);
460 : // Elements have changed, neighbors, ids etc
461 272 : mesh->unset_is_prepared();
462 :
463 544 : return mesh;
464 816 : }
|
