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 "MeshRepairGenerator.h"
11 : #include "CastUniquePointer.h"
12 : #include "MooseMeshUtils.h"
13 :
14 : #include "libmesh/mesh_tools.h"
15 : #include "libmesh/mesh_modification.h"
16 : #include "libmesh/face_c0polygon.h"
17 :
18 : registerMooseObject("MooseApp", MeshRepairGenerator);
19 :
20 : InputParameters
21 3212 : MeshRepairGenerator::validParams()
22 : {
23 :
24 3212 : InputParameters params = MeshGenerator::validParams();
25 6424 : params.addClassDescription(
26 : "Mesh generator to perform various improvement / fixing operations on an input mesh");
27 12848 : params.addRequiredParam<MeshGeneratorName>("input",
28 : "Name of the mesh generator providing the mesh");
29 :
30 12848 : params.addParam<bool>("fix_node_overlap", false, "Whether to merge overlapping nodes");
31 9636 : params.addParam<Real>(
32 6424 : "node_overlap_tol", 1e-8, "Absolute tolerance for merging overlapping nodes");
33 :
34 9636 : params.addParam<bool>(
35 6424 : "fix_elements_orientation", false, "Whether to flip elements with negative volumes");
36 :
37 9636 : params.addParam<bool>("separate_blocks_by_element_types",
38 6424 : false,
39 : "Create new blocks if multiple element types are present in a block");
40 :
41 9636 : params.addParam<bool>("merge_boundary_ids_with_same_name",
42 6424 : false,
43 : "Merge boundaries if they have the same name but different boundary IDs");
44 :
45 9636 : params.addParam<bool>(
46 : "renumber_contiguously",
47 6424 : false,
48 : "Whether to renumber the elements of the mesh so the numbering is contiguous");
49 :
50 6424 : params.addParam<bool>(
51 6424 : "split_nonconvex_polygons", false, "Split non-convex polygons to form convex ones");
52 :
53 3212 : return params;
54 0 : }
55 :
56 74 : MeshRepairGenerator::MeshRepairGenerator(const InputParameters & parameters)
57 : : MeshGenerator(parameters),
58 74 : _input(getMesh("input")),
59 148 : _fix_overlapping_nodes(getParam<bool>("fix_node_overlap")),
60 148 : _node_overlap_tol(getParam<Real>("node_overlap_tol")),
61 148 : _fix_element_orientation(getParam<bool>("fix_elements_orientation")),
62 148 : _elem_type_separation(getParam<bool>("separate_blocks_by_element_types")),
63 148 : _boundary_id_merge(getParam<bool>("merge_boundary_ids_with_same_name")),
64 222 : _split_nonconvex_polygons(getParam<bool>("split_nonconvex_polygons"))
65 : {
66 45 : if (!_fix_overlapping_nodes && !_fix_element_orientation && !_elem_type_separation &&
67 149 : !_boundary_id_merge && !getParam<bool>("renumber_contiguously") && !_split_nonconvex_polygons)
68 0 : mooseError("No specific item to fix. Are any of the parameters misspelled?");
69 74 : }
70 :
71 : std::unique_ptr<MeshBase>
72 67 : MeshRepairGenerator::generate()
73 : {
74 67 : std::unique_ptr<MeshBase> mesh = std::move(_input);
75 :
76 : // We're trying to repair a potentially broken mesh; we'll just
77 : // start with a full prepare rather than trying to be efficient and
78 : // risking missing something.
79 67 : mesh->prepare_for_use();
80 :
81 : // Blanket ban on distributed. This can be relaxed for some operations if needed
82 67 : if (!mesh->is_serial())
83 0 : mooseError("MeshRepairGenerator requires a serial mesh. The mesh should not be distributed.");
84 :
85 67 : if (_fix_overlapping_nodes)
86 25 : fixOverlappingNodes(mesh);
87 :
88 : // Flip orientation of elements to keep positive volumes
89 67 : if (_fix_element_orientation)
90 9 : MeshTools::Modification::orient_elements(*mesh);
91 :
92 : // Disambiguate any block that has elements of multiple types
93 67 : if (_elem_type_separation)
94 9 : separateSubdomainsByElementType(mesh);
95 :
96 : // Assign a single boundary ID to boundaries that have the same name
97 67 : if (_boundary_id_merge)
98 9 : MooseMeshUtils::mergeBoundaryIDsWithSameName(*mesh);
99 :
100 : // Renumber the mesh despite any mesh flag
101 201 : if (getParam<bool>("renumber_contiguously"))
102 : {
103 8 : const auto prev_status = mesh->allow_renumbering();
104 8 : mesh->allow_renumbering(true);
105 8 : mesh->renumber_nodes_and_elements();
106 8 : mesh->allow_renumbering(prev_status);
107 : }
108 :
109 : // Split non-convex polygons to make convex ones
110 67 : if (_split_nonconvex_polygons)
111 7 : splitNonConvexPolygons(mesh);
112 :
113 67 : mesh->unset_is_prepared();
114 67 : return mesh;
115 0 : }
116 :
117 : void
118 25 : MeshRepairGenerator::fixOverlappingNodes(std::unique_ptr<MeshBase> & mesh) const
119 : {
120 25 : unsigned int num_fixed_nodes = 0;
121 25 : auto pl = mesh->sub_point_locator();
122 25 : pl->set_close_to_point_tol(_node_overlap_tol);
123 :
124 25 : std::unordered_set<dof_id_type> nodes_removed;
125 : // loop on nodes
126 517 : for (auto & node : mesh->node_ptr_range())
127 : {
128 : // this node has already been removed
129 246 : if (nodes_removed.count(node->id()))
130 66 : continue;
131 :
132 : // find all the elements around this node
133 180 : std::set<const Elem *> elements;
134 180 : (*pl)(*node, elements);
135 :
136 426 : for (auto & elem : elements)
137 : {
138 246 : bool found = false;
139 1103 : for (auto & elem_node : elem->node_ref_range())
140 : {
141 1037 : if (node->id() == elem_node.id())
142 : {
143 180 : found = true;
144 180 : break;
145 : }
146 : }
147 246 : if (!found)
148 : {
149 440 : for (auto & elem_node : elem->node_ref_range())
150 : {
151 374 : if (elem_node.id() == node->id())
152 0 : continue;
153 374 : const Real tol = _node_overlap_tol;
154 : // Compares the coordinates
155 374 : const auto x_node = (*node)(0);
156 374 : const auto x_elem_node = elem_node(0);
157 374 : const auto y_node = (*node)(1);
158 374 : const auto y_elem_node = elem_node(1);
159 374 : const auto z_node = (*node)(2);
160 374 : const auto z_elem_node = elem_node(2);
161 :
162 374 : if (MooseUtils::absoluteFuzzyEqual(x_node, x_elem_node, tol) &&
163 472 : MooseUtils::absoluteFuzzyEqual(y_node, y_elem_node, tol) &&
164 98 : MooseUtils::absoluteFuzzyEqual(z_node, z_elem_node, tol))
165 : {
166 : // Merging two nodes from the same element is almost never a good idea
167 66 : if (elem->get_node_index(node) != libMesh::invalid_uint)
168 : {
169 0 : _console << "Two overlapping nodes in element " << elem->id() << " right by "
170 0 : << elem->vertex_average() << ".\n They will not be stitched" << std::endl;
171 0 : continue;
172 : }
173 :
174 : // Coordinates are the same but it's not the same node
175 : // Replace the node in the element
176 66 : const_cast<Elem *>(elem)->set_node(elem->get_node_index(&elem_node), node);
177 66 : nodes_removed.insert(elem_node.id());
178 :
179 66 : num_fixed_nodes++;
180 66 : if (num_fixed_nodes < 10)
181 66 : _console << "Stitching nodes " << *node << " and " << elem_node
182 66 : << std::endl;
183 0 : else if (num_fixed_nodes == 10)
184 0 : _console << "Node stitching will now proceed silently." << std::endl;
185 : }
186 : }
187 : }
188 : }
189 205 : }
190 25 : _console << "Number of overlapping nodes which got merged: " << num_fixed_nodes << std::endl;
191 25 : if (mesh->allow_renumbering())
192 9 : mesh->renumber_nodes_and_elements();
193 : else
194 : {
195 16 : mesh->remove_orphaned_nodes();
196 16 : mesh->update_parallel_id_counts();
197 : }
198 25 : }
199 :
200 : void
201 9 : MeshRepairGenerator::separateSubdomainsByElementType(std::unique_ptr<MeshBase> & mesh) const
202 : {
203 9 : std::set<subdomain_id_type> ids;
204 9 : mesh->subdomain_ids(ids);
205 : // loop on sub-domain
206 18 : for (const auto id : ids)
207 : {
208 : // Gather all the element types and blocks
209 : // ElemType defines an enum for geometric element types
210 9 : std::set<ElemType> types;
211 : // loop on elements within this sub-domain
212 27 : for (auto & elem : mesh->active_subdomain_elements_ptr_range(id))
213 27 : types.insert(elem->type());
214 :
215 : // This call must be performed on all processes
216 9 : auto next_block_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
217 :
218 9 : if (types.size() > 1)
219 : {
220 9 : subdomain_id_type i = 0;
221 27 : for (const auto type : types)
222 : {
223 18 : auto new_id = next_block_id + i++;
224 : // Create blocks when a block has multiple element types
225 18 : mesh->subdomain_name(new_id) = mesh->subdomain_name(id) + "_" + Moose::stringify(type);
226 :
227 : // Re-assign elements to the new blocks
228 72 : for (auto elem : mesh->active_subdomain_elements_ptr_range(id))
229 27 : if (elem->type() == type)
230 36 : elem->subdomain_id() = new_id;
231 : }
232 : }
233 9 : }
234 9 : }
235 :
236 : void
237 7 : MeshRepairGenerator::splitNonConvexPolygons(std::unique_ptr<MeshBase> & mesh) const
238 : {
239 : // Counters to keep track of what happened in the splitting
240 7 : unsigned int num_polygons = 0;
241 7 : unsigned int num_nonconvex = 0;
242 7 : unsigned int num_triangulated = 0;
243 7 : std::vector<Elem *> elems_to_delete;
244 7 : std::vector<std::unique_ptr<libMesh::C0Polygon>> elements_to_add_to_mesh;
245 :
246 : // Missing parallel packing for polys
247 7 : if (!mesh->is_serial())
248 0 : mooseError("MeshRepairGenerator requires a serial mesh for this operation. "
249 : "The mesh should not be distributed.");
250 :
251 35 : for (auto elem : mesh->element_ptr_range())
252 : {
253 : // We are not splitting non-convex quads. Maybe later.
254 28 : if (elem->type() == libMesh::C0POLYGON)
255 : {
256 28 : num_polygons++;
257 :
258 : // Lambda function to determine if an element is convex
259 154 : auto is_convex = [](const Elem * elem,
260 : std::vector<std::pair<int, Real>> & obtuse_angle_nodes) -> bool
261 : {
262 : mooseAssert(elem, "Should have an elem");
263 : mooseAssert(obtuse_angle_nodes.empty(), "Should not be empty");
264 :
265 : // Find the normal to the element plane
266 154 : const auto elem_centroid = elem->vertex_average();
267 154 : const auto n_nodes = elem->n_nodes();
268 154 : Point plane_normal;
269 : // Two nodes could be aligned with the centroid in a non-convex polygon
270 154 : for (const auto i : make_range(n_nodes))
271 : {
272 154 : const auto v1 = elem_centroid - elem->point(i);
273 154 : const auto v2 = elem_centroid - elem->point((i + 1) % n_nodes);
274 154 : if (!MooseUtils::absoluteFuzzyEqual((v1.cross(v2)).norm_sq(), 0))
275 : {
276 154 : plane_normal = v1.cross(v2).unit();
277 154 : break;
278 : }
279 : }
280 :
281 : // Look for angles > 180 deg with the cross product
282 154 : bool convex = true;
283 756 : for (const auto i : make_range(n_nodes))
284 : {
285 602 : const auto n1 = elem->point(i);
286 602 : const auto n2 = elem->point((i + 1) % n_nodes);
287 602 : const auto n3 = elem->point((i + 2) % n_nodes);
288 602 : const auto top_dir = (n2 - n1).cross(n3 - n2);
289 :
290 602 : if (plane_normal * top_dir <= 0)
291 : {
292 105 : convex = false;
293 105 : if (!MooseUtils::absoluteFuzzyEqual(plane_normal * top_dir, 0))
294 105 : obtuse_angle_nodes.push_back(
295 210 : std::make_pair<int, Real>(i, plane_normal * top_dir / top_dir.norm()));
296 : // n1, n2 and n3 are likely aligned
297 : else
298 0 : obtuse_angle_nodes.push_back(std::make_pair<int, Real>(i, -1));
299 : }
300 : }
301 154 : return convex;
302 : };
303 :
304 28 : std::vector<std::pair<int, Real>> parent_obtuse_angles;
305 28 : const auto convex = is_convex(elem, parent_obtuse_angles);
306 :
307 28 : if (!convex)
308 28 : num_nonconvex++;
309 : // Move to next one if it is convex, no need to fix it
310 : else
311 0 : continue;
312 :
313 : // Lambda function to split a polygon into two polygons at the specified nodes
314 : // TODO: we could try to split a polygon into more than two polygons
315 63 : auto cut_polygon = [&mesh](const Elem * elem, unsigned int node_i_cut1, unsigned node_i_cut2)
316 : -> std::pair<std::unique_ptr<libMesh::C0Polygon>, std::unique_ptr<libMesh::C0Polygon>>
317 : {
318 63 : const auto cut1 = std::min(node_i_cut1, node_i_cut2);
319 63 : const auto cut2 = std::max(node_i_cut1, node_i_cut2);
320 :
321 : // Count the number of sides on each side of the cut
322 63 : const auto ns1 = cut2 - cut1 + 1;
323 63 : const auto ns2 = elem->n_nodes() - ns1 + 2;
324 : mooseAssert(ns1 >= 3, "Should cut at least a triangle");
325 : mooseAssert(ns2 >= 3, "Should cut at least a triangle");
326 :
327 : // Create the children and add their nodes
328 63 : auto child1 = std::make_unique<libMesh::C0Polygon>(ns1);
329 280 : for (const auto i_n : make_range(cut1, cut2 + 1))
330 217 : child1->set_node(i_n - cut1, mesh->node_ptr(elem->node_ptr(i_n)->id()));
331 :
332 63 : auto child2 = std::make_unique<libMesh::C0Polygon>(ns2);
333 280 : for (const auto i_n : make_range(cut2, elem->n_nodes() + cut1 + 1))
334 434 : child2->set_node((i_n - cut2) % child2->n_nodes(),
335 217 : mesh->node_ptr(elem->node_ptr(i_n % elem->n_nodes())->id()));
336 :
337 126 : return std::make_pair(std::move(child1), std::move(child2));
338 63 : };
339 :
340 : // If not convex, try to split it.
341 : // Let's first try to split it recursively at every large angle
342 : // NOTE: These vectors of unique pointer keep the memory ownership of the elements
343 : // during the loop attempting to fix the polygon with the heuristic below
344 28 : bool failed = false;
345 28 : std::vector<std::unique_ptr<libMesh::C0Polygon>> elements_to_add_to_mesh_temp;
346 28 : std::vector<std::pair<std::unique_ptr<Elem>, std::vector<std::pair<int, Real>>>> elems_to_cut;
347 :
348 : // Create a clone of the first element to simplify logic
349 : std::unique_ptr<libMesh::C0Polygon> base_elem =
350 28 : std::make_unique<libMesh::C0Polygon>(elem->n_nodes());
351 196 : for (const auto i_n : make_range(elem->n_nodes()))
352 168 : base_elem->set_node(i_n, elem->node_ptr(i_n));
353 28 : elems_to_cut.push_back(std::make_pair(std::move(base_elem), parent_obtuse_angles));
354 :
355 91 : while (!failed && !elems_to_cut.empty())
356 : {
357 63 : auto & [current_elem, large_angles] = elems_to_cut.back();
358 :
359 : // Split the element on one side at the node with the worst obtuse angle
360 : // TODO: try all of them instead to see if a single cut can fix the element
361 : auto worst_angle_it =
362 63 : std::min_element(large_angles.begin(),
363 : large_angles.end(),
364 42 : [](std::pair<int, Real> lhs, std::pair<int, Real> rhs) -> bool
365 42 : { return lhs.second < rhs.second; });
366 63 : unsigned int worst_angle_i = (*worst_angle_it).first;
367 63 : const auto n_nodes = current_elem->n_nodes();
368 :
369 : // Split the element on the other side at roughly the opposite node
370 : // TODO: we could try all of the 'other' nodes here too to see if a single cut can fix the
371 : // element NOTE: if trying multiple cuts (not implemented), we should remove the opposite
372 : // from the large_angles
373 63 : const auto opposite = ((worst_angle_i + n_nodes / 2) % n_nodes);
374 :
375 63 : auto [child1, child2] = cut_polygon(current_elem.get(), worst_angle_i, opposite);
376 :
377 : // Check the two fragments
378 63 : std::vector<std::pair<int, Real>> child1_large_angles;
379 63 : bool is_convex1 = is_convex(child1.get(), child1_large_angles);
380 63 : std::vector<std::pair<int, Real>> child2_large_angles;
381 63 : bool is_convex2 = is_convex(child2.get(), child2_large_angles);
382 :
383 : // Can we keep going?
384 63 : if (child1->n_nodes() < 4 && !is_convex1)
385 0 : failed = true;
386 63 : if (child2->n_nodes() < 4 && !is_convex2)
387 0 : failed = true;
388 :
389 : // We managed to cut the elements into a convex part, but it's flipped
390 : // so the element is likely "outside" of the starting element
391 63 : if (is_convex1 && child1->volume() <= 0)
392 0 : failed = true;
393 63 : if (is_convex2 && child2->volume() <= 0)
394 0 : failed = true;
395 :
396 : // Current element is done
397 : // NOTE: if trying multiple cuts (not implemented), we should not erase yet, we should
398 : // instead only remove the cut we tried from the 'large_angles' vector used to pick cuts
399 63 : large_angles.erase(worst_angle_it);
400 63 : elems_to_cut.pop_back();
401 :
402 : // Add non convex elements to the list of elements to cut
403 63 : if (!is_convex1)
404 14 : elems_to_cut.push_back(std::make_pair(std::move(child1), child1_large_angles));
405 63 : if (!is_convex2)
406 21 : elems_to_cut.push_back(std::make_pair(std::move(child2), child2_large_angles));
407 :
408 : // Add convex elements to the mesh
409 : // NOTE: if trying multiple cuts, we should not do this yet
410 63 : if (is_convex1)
411 49 : elements_to_add_to_mesh_temp.push_back(std::move(child1));
412 63 : if (is_convex2)
413 42 : elements_to_add_to_mesh_temp.push_back(std::move(child2));
414 63 : }
415 28 : bool fixed_it = !failed;
416 :
417 : // Keep all the cut elements, to be added at the end to avoid invalidating iterators
418 28 : if (fixed_it)
419 119 : for (auto & elem_uptr : elements_to_add_to_mesh_temp)
420 : {
421 91 : elem_uptr->inherit_data_from(*elem);
422 91 : elements_to_add_to_mesh.push_back(std::move(elem_uptr));
423 : }
424 :
425 : // Heuristic did not work, just use a triangulation
426 : // NOTE: This is not the triangulation of the polygon. It is simple though
427 28 : if (!fixed_it)
428 : {
429 0 : const auto centroid_node = mesh->add_point(elem->true_centroid());
430 0 : num_triangulated++;
431 0 : const auto * poly = dynamic_cast<libMesh::C0Polygon *>(elem);
432 0 : const auto n_sides = poly->n_sides();
433 0 : for (const auto i_tr : make_range(n_sides))
434 : {
435 0 : auto new_elem = std::make_unique<libMesh::C0Polygon>(3);
436 0 : new_elem->set_node(0, mesh->node_ptr(elem->node_ptr(i_tr)->id()));
437 0 : new_elem->set_node(1, centroid_node);
438 0 : new_elem->set_node(2, mesh->node_ptr(elem->node_ptr((i_tr + 1) % n_sides)->id()));
439 :
440 : // Check for degenerate case
441 0 : if (MooseUtils::absoluteFuzzyEqual(
442 0 : (*(Point *)centroid_node - *(Point *)elem->node_ptr(i_tr))
443 0 : .cross(*(Point *)centroid_node -
444 0 : *(Point *)elem->node_ptr((i_tr + 1) % n_sides))
445 0 : .norm_sq(),
446 0 : 0))
447 0 : mooseError("Manual triangulation failed as two consecutive nodes are aligned with the "
448 : "centroid");
449 :
450 0 : new_elem->inherit_data_from(*elem);
451 0 : elements_to_add_to_mesh.push_back(std::move(new_elem));
452 0 : }
453 : }
454 : // Element got fixed, can be deleted after (can't while using the range)
455 28 : elems_to_delete.push_back(elem);
456 28 : }
457 7 : }
458 : // Delete the original element
459 35 : for (auto elem : elems_to_delete)
460 28 : mesh->delete_elem(elem);
461 : // Add the new ones
462 98 : for (auto & new_elem_ptr : elements_to_add_to_mesh)
463 91 : mesh->add_elem(std::move(new_elem_ptr));
464 :
465 7 : _console << "Number of non-convex polygons which got split into convex polygons: "
466 7 : << num_nonconvex << std::endl;
467 7 : if (num_triangulated)
468 0 : _console << "Number of non-convex polygons split using a triangulation: " << num_triangulated
469 0 : << ", using heuristic: " << num_nonconvex - num_triangulated << std::endl;
470 7 : if (!num_polygons)
471 0 : mooseWarning("No C0 polygons in mesh: the polyon convexity fix did nothing");
472 7 : }
|
