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 "MeshDiagnosticsGenerator.h"
11 : #include "MooseMeshUtils.h"
12 : #include "CastUniquePointer.h"
13 : #include "MeshCoarseningUtils.h"
14 : #include "MeshBaseDiagnosticsUtils.h"
15 :
16 : #include "libmesh/mesh_tools.h"
17 : #include "libmesh/mesh_refinement.h"
18 : #include "libmesh/fe.h"
19 : #include "libmesh/quadrature_gauss.h"
20 : #include "libmesh/face_tri3.h"
21 : #include "libmesh/cell_tet4.h"
22 : #include "libmesh/face_quad4.h"
23 : #include "libmesh/cell_hex8.h"
24 : #include "libmesh/string_to_enum.h"
25 : #include "libmesh/enum_point_locator_type.h"
26 :
27 : // C++
28 : #include <cstring>
29 :
30 : registerMooseObject("MooseApp", MeshDiagnosticsGenerator);
31 :
32 : InputParameters
33 3940 : MeshDiagnosticsGenerator::validParams()
34 : {
35 :
36 3940 : InputParameters params = MeshGenerator::validParams();
37 :
38 15760 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to diagnose");
39 7880 : params.addClassDescription("Runs a series of diagnostics on the mesh to detect potential issues "
40 : "such as unsupported features");
41 :
42 : // Options for the output level
43 15760 : MooseEnum chk_option("NO_CHECK INFO WARNING ERROR", "NO_CHECK");
44 :
45 15760 : params.addParam<MooseEnum>(
46 : "examine_sidesets_orientation",
47 : chk_option,
48 : "whether to check that sidesets are consistently oriented using neighbor subdomains. If a "
49 : "sideset is inconsistently oriented within a subdomain, this will not be detected");
50 15760 : params.addParam<MooseEnum>(
51 : "check_for_watertight_sidesets",
52 : chk_option,
53 : "whether to check for external sides that are not assigned to any sidesets");
54 15760 : params.addParam<MooseEnum>(
55 : "check_for_watertight_nodesets",
56 : chk_option,
57 : "whether to check for external nodes that are not assigned to any nodeset");
58 15760 : params.addParam<std::vector<BoundaryName>>(
59 : "boundaries_to_check",
60 : {},
61 : "Names boundaries that should form a watertight envelope around the mesh. Defaults to all "
62 : "the boundaries combined.");
63 15760 : params.addParam<MooseEnum>(
64 : "examine_element_volumes", chk_option, "whether to examine volume of the elements");
65 15760 : params.addParam<Real>("minimum_element_volumes", 1e-16, "minimum size for element volume");
66 15760 : params.addParam<Real>("maximum_element_volumes", 1e16, "Maximum size for element volume");
67 :
68 15760 : params.addParam<MooseEnum>("examine_element_types",
69 : chk_option,
70 : "whether to look for multiple element types in the same sub-domain");
71 15760 : params.addParam<MooseEnum>(
72 : "examine_element_overlap", chk_option, "whether to find overlapping elements");
73 15760 : params.addParam<MooseEnum>(
74 : "examine_nonplanar_sides", chk_option, "whether to check element sides are planar");
75 15760 : params.addParam<MooseEnum>("examine_non_conformality",
76 : chk_option,
77 : "whether to examine the conformality of elements in the mesh");
78 15760 : params.addParam<MooseEnum>("examine_non_matching_edges",
79 : chk_option,
80 : "Whether to check if there are any intersecting edges");
81 15760 : params.addParam<Real>("intersection_tol", TOLERANCE, "tolerence for intersecting edges");
82 15760 : params.addParam<Real>("nonconformal_tol", TOLERANCE, "tolerance for element non-conformality");
83 15760 : params.addParam<MooseEnum>(
84 : "search_for_adaptivity_nonconformality",
85 : chk_option,
86 : "whether to check for non-conformality arising from adaptive mesh refinement");
87 15760 : params.addParam<MooseEnum>("check_local_jacobian",
88 : chk_option,
89 : "whether to check the local Jacobian for bad (non-positive) values");
90 15760 : params.addParam<MooseEnum>(
91 : "check_polygons", chk_option, "Whether to check that all C0 polygons are convex");
92 7880 : params.addParam<unsigned int>(
93 : "log_length_limit",
94 7880 : 10,
95 : "How many problematic element/nodes/sides/etc are explicitly reported on by each check");
96 7880 : return params;
97 3940 : }
98 :
99 441 : MeshDiagnosticsGenerator::MeshDiagnosticsGenerator(const InputParameters & parameters)
100 : : MeshGenerator(parameters),
101 441 : _input(getMesh("input")),
102 882 : _check_sidesets_orientation(getParam<MooseEnum>("examine_sidesets_orientation")),
103 882 : _check_watertight_sidesets(getParam<MooseEnum>("check_for_watertight_sidesets")),
104 882 : _check_watertight_nodesets(getParam<MooseEnum>("check_for_watertight_nodesets")),
105 882 : _watertight_boundary_names(getParam<std::vector<BoundaryName>>("boundaries_to_check")),
106 882 : _check_element_volumes(getParam<MooseEnum>("examine_element_volumes")),
107 882 : _min_volume(getParam<Real>("minimum_element_volumes")),
108 882 : _max_volume(getParam<Real>("maximum_element_volumes")),
109 882 : _check_element_types(getParam<MooseEnum>("examine_element_types")),
110 882 : _check_element_overlap(getParam<MooseEnum>("examine_element_overlap")),
111 882 : _check_non_planar_sides(getParam<MooseEnum>("examine_nonplanar_sides")),
112 882 : _check_non_conformal_mesh(getParam<MooseEnum>("examine_non_conformality")),
113 882 : _non_conformality_tol(getParam<Real>("nonconformal_tol")),
114 882 : _check_non_matching_edges(getParam<MooseEnum>("examine_non_matching_edges")),
115 882 : _non_matching_edge_tol(getParam<Real>("intersection_tol")),
116 882 : _check_adaptivity_non_conformality(
117 : getParam<MooseEnum>("search_for_adaptivity_nonconformality")),
118 882 : _check_local_jacobian(getParam<MooseEnum>("check_local_jacobian")),
119 882 : _check_polygons(getParam<MooseEnum>("check_polygons")),
120 1764 : _num_outputs(getParam<unsigned int>("log_length_limit"))
121 : {
122 : // Check that no secondary parameters have been passed with the main check disabled
123 1323 : if ((isParamSetByUser("minimum_element_volumes") ||
124 1764 : isParamSetByUser("maximum_element_volumes")) &&
125 6 : _check_element_volumes == "NO_CHECK")
126 0 : paramError("examine_element_volumes",
127 : "You must set this parameter to true to trigger element size checks");
128 1323 : if (isParamSetByUser("nonconformal_tol") && _check_non_conformal_mesh == "NO_CHECK")
129 0 : paramError("examine_non_conformality",
130 : "You must set this parameter to true to trigger mesh conformality check");
131 828 : if (_check_sidesets_orientation == "NO_CHECK" && _check_watertight_sidesets == "NO_CHECK" &&
132 373 : _check_watertight_nodesets == "NO_CHECK" && _check_element_volumes == "NO_CHECK" &&
133 323 : _check_element_types == "NO_CHECK" && _check_element_overlap == "NO_CHECK" &&
134 130 : _check_non_planar_sides == "NO_CHECK" && _check_non_conformal_mesh == "NO_CHECK" &&
135 109 : _check_adaptivity_non_conformality == "NO_CHECK" && _check_local_jacobian == "NO_CHECK" &&
136 828 : _check_non_matching_edges == "NO_CHECK" && _check_polygons == "NO_CHECK")
137 3 : mooseError("You need to turn on at least one diagnostic. Did you misspell a parameter?");
138 438 : }
139 :
140 : std::unique_ptr<MeshBase>
141 305 : MeshDiagnosticsGenerator::generate()
142 : {
143 305 : std::unique_ptr<MeshBase> mesh = std::move(_input);
144 :
145 : // Most of the checks assume we have the full mesh
146 305 : if (!mesh->is_serial())
147 0 : mooseError("Only serialized meshes are supported");
148 :
149 : // We prepare for use at the beginning to facilitate diagnosis
150 : // This deliberately does not trust the mesh to know whether it's already prepared or not
151 305 : mesh->prepare_for_use();
152 :
153 : // check that specified boundary is valid, convert BoundaryNames to BoundaryIDs, and sort
154 305 : for (const auto & boundary_name : _watertight_boundary_names)
155 : {
156 0 : if (!MooseMeshUtils::hasBoundaryName(*mesh, boundary_name))
157 0 : mooseError("User specified boundary_to_check \'", boundary_name, "\' does not exist");
158 : }
159 305 : _watertight_boundaries = MooseMeshUtils::getBoundaryIDs(*mesh, _watertight_boundary_names, false);
160 305 : std::sort(_watertight_boundaries.begin(), _watertight_boundaries.end());
161 :
162 305 : if (_check_sidesets_orientation != "NO_CHECK")
163 54 : checkSidesetsOrientation(mesh);
164 :
165 305 : if (_check_watertight_sidesets != "NO_CHECK")
166 14 : checkWatertightSidesets(mesh);
167 :
168 305 : if (_check_watertight_nodesets != "NO_CHECK")
169 14 : checkWatertightNodesets(mesh);
170 :
171 305 : if (_check_element_volumes != "NO_CHECK")
172 46 : checkElementVolumes(mesh);
173 :
174 299 : if (_check_element_types != "NO_CHECK")
175 43 : checkElementTypes(mesh);
176 :
177 299 : if (_check_element_overlap != "NO_CHECK")
178 127 : checkElementOverlap(mesh);
179 :
180 299 : if (_check_non_planar_sides != "NO_CHECK")
181 39 : checkNonPlanarSides(mesh);
182 :
183 299 : if (_check_non_conformal_mesh != "NO_CHECK")
184 134 : checkNonConformalMesh(mesh);
185 :
186 299 : if (_check_adaptivity_non_conformality != "NO_CHECK")
187 198 : checkNonConformalMeshFromAdaptivity(mesh);
188 :
189 284 : if (_check_local_jacobian != "NO_CHECK")
190 47 : checkLocalJacobians(mesh);
191 :
192 284 : if (_check_non_matching_edges != "NO_CHECK")
193 15 : checkNonMatchingEdges(mesh);
194 :
195 284 : if (_check_polygons != "NO_CHECK")
196 14 : checkPolygons(mesh);
197 :
198 568 : return dynamic_pointer_cast<MeshBase>(mesh);
199 284 : }
200 :
201 : void
202 54 : MeshDiagnosticsGenerator::checkSidesetsOrientation(const std::unique_ptr<MeshBase> & mesh) const
203 : {
204 54 : auto & boundary_info = mesh->get_boundary_info();
205 54 : auto side_tuples = boundary_info.build_side_list();
206 :
207 148 : for (const auto bid : boundary_info.get_boundary_ids())
208 : {
209 : // This check only looks at subdomains on both sides of the sideset
210 : // it wont pick up if the sideset is changing orientations while inside of a subdomain
211 94 : std::set<std::pair<subdomain_id_type, subdomain_id_type>> block_neighbors;
212 1414 : for (const auto index : index_range(side_tuples))
213 : {
214 1320 : if (std::get<2>(side_tuples[index]) != bid)
215 864 : continue;
216 456 : const auto elem_ptr = mesh->elem_ptr(std::get<0>(side_tuples[index]));
217 456 : if (elem_ptr->neighbor_ptr(std::get<1>(side_tuples[index])))
218 56 : block_neighbors.insert(std::make_pair(
219 56 : elem_ptr->subdomain_id(),
220 56 : elem_ptr->neighbor_ptr(std::get<1>(side_tuples[index]))->subdomain_id()));
221 : }
222 :
223 : // Check that there is no flipped pair
224 94 : std::set<std::pair<subdomain_id_type, subdomain_id_type>> flipped_pairs;
225 122 : for (const auto & block_pair_1 : block_neighbors)
226 84 : for (const auto & block_pair_2 : block_neighbors)
227 56 : if (block_pair_1 != block_pair_2)
228 28 : if (block_pair_1.first == block_pair_2.second &&
229 28 : block_pair_1.second == block_pair_2.first)
230 28 : flipped_pairs.insert(block_pair_1);
231 :
232 94 : std::string message;
233 : const std::string sideset_full_name =
234 94 : boundary_info.sideset_name(bid) + " (" + std::to_string(bid) + ")";
235 94 : if (!flipped_pairs.empty())
236 : {
237 14 : std::string block_pairs_string = "";
238 42 : for (const auto & pair : flipped_pairs)
239 : block_pairs_string +=
240 56 : " [" + mesh->subdomain_name(pair.first) + " (" + std::to_string(pair.first) + "), " +
241 84 : mesh->subdomain_name(pair.second) + " (" + std::to_string(pair.second) + ")]";
242 28 : message = "Inconsistent orientation of sideset " + sideset_full_name +
243 14 : " with regards to subdomain pairs" + block_pairs_string;
244 14 : }
245 : else
246 160 : message = "Sideset " + sideset_full_name +
247 80 : " is consistently oriented with regards to the blocks it neighbors";
248 :
249 94 : diagnosticsLog(message, _check_sidesets_orientation, flipped_pairs.size());
250 :
251 : // Now check that there is no sideset radically flipping from one side's normal to another
252 : // side next to it, in the same sideset
253 : // We'll consider pi / 2 to be the most steep angle we'll pass
254 94 : unsigned int num_normals_flipping = 0;
255 94 : Real steepest_side_angles = 0;
256 1414 : for (const auto & [elem_id, side_id, side_bid] : side_tuples)
257 : {
258 1320 : if (side_bid != bid)
259 864 : continue;
260 456 : const auto & elem_ptr = mesh->elem_ptr(elem_id);
261 :
262 : // Get side normal
263 456 : const std::unique_ptr<const Elem> face = elem_ptr->build_side_ptr(side_id);
264 : std::unique_ptr<libMesh::FEBase> fe(
265 456 : libMesh::FEBase::build(elem_ptr->dim(), libMesh::FEType(elem_ptr->default_order())));
266 456 : libMesh::QGauss qface(elem_ptr->dim() - 1, CONSTANT);
267 456 : fe->attach_quadrature_rule(&qface);
268 456 : const auto & normals = fe->get_normals();
269 456 : fe->reinit(elem_ptr, side_id);
270 : mooseAssert(normals.size() == 1, "We expected only one normal here");
271 456 : const auto & side_normal = normals[0];
272 :
273 : // Compare to the sideset normals of neighbor sides in that sideset
274 2568 : for (const auto neighbor : elem_ptr->neighbor_ptr_range())
275 2112 : if (neighbor)
276 8080 : for (const auto neigh_side_index : neighbor->side_index_range())
277 : {
278 : // Check that the neighbor side is also in the sideset being examined
279 6768 : if (!boundary_info.has_boundary_id(neighbor, neigh_side_index, bid))
280 5872 : continue;
281 :
282 : // We re-init everything for the neighbor in case it's a different dimension
283 : std::unique_ptr<libMesh::FEBase> fe_neighbor(libMesh::FEBase::build(
284 896 : neighbor->dim(), libMesh::FEType(neighbor->default_order())));
285 896 : libMesh::QGauss qface(neighbor->dim() - 1, CONSTANT);
286 896 : fe_neighbor->attach_quadrature_rule(&qface);
287 896 : const auto & neigh_normals = fe_neighbor->get_normals();
288 896 : fe_neighbor->reinit(neighbor, neigh_side_index);
289 : mooseAssert(neigh_normals.size() == 1, "We expected only one normal here");
290 896 : const auto & neigh_side_normal = neigh_normals[0];
291 :
292 : // Check the angle by computing the dot product
293 896 : if (neigh_side_normal * side_normal <= 0)
294 : {
295 56 : num_normals_flipping++;
296 56 : steepest_side_angles =
297 56 : std::max(std::acos(neigh_side_normal * side_normal), steepest_side_angles);
298 56 : if (num_normals_flipping <= _num_outputs)
299 56 : _console << "Side normals changed by more than pi/2 for sideset "
300 56 : << sideset_full_name << " between side " << side_id << " of element "
301 56 : << elem_ptr->id() << " and side " << neigh_side_index
302 56 : << " of neighbor element " << neighbor->id() << std::endl;
303 0 : else if (num_normals_flipping == _num_outputs + 1)
304 : _console << "Maximum output reached for sideset normal flipping check. Silencing "
305 0 : "output from now on"
306 0 : << std::endl;
307 : }
308 896 : }
309 456 : }
310 :
311 94 : if (num_normals_flipping)
312 28 : message = "Sideset " + sideset_full_name +
313 28 : " has two neighboring sides with a very large angle. Largest angle detected: " +
314 56 : std::to_string(steepest_side_angles) + " rad (" +
315 42 : std::to_string(steepest_side_angles * 180 / libMesh::pi) + " degrees).";
316 : else
317 160 : message = "Sideset " + sideset_full_name +
318 : " does not appear to have side-to-neighbor-side orientation flips. All neighbor "
319 80 : "sides normal differ by less than pi/2";
320 :
321 94 : diagnosticsLog(message, _check_sidesets_orientation, num_normals_flipping);
322 94 : }
323 54 : }
324 :
325 : void
326 14 : MeshDiagnosticsGenerator::checkWatertightSidesets(const std::unique_ptr<MeshBase> & mesh) const
327 : {
328 : /*
329 : Algorithm Overview:
330 : 1) Loop through all elements
331 : 2) For each element loop through all its sides
332 : 3) If it has no neighbors it's an external side
333 : 4) If external check if it's part of a sideset
334 : */
335 14 : if (mesh->mesh_dimension() < 2)
336 0 : mooseError("The sideset check only works for 2D and 3D meshes");
337 14 : auto & boundary_info = mesh->get_boundary_info();
338 14 : boundary_info.build_side_list();
339 14 : const auto sideset_map = boundary_info.get_sideset_map();
340 14 : unsigned int num_faces_without_sideset = 0;
341 :
342 182 : for (const auto elem : mesh->active_element_ptr_range())
343 : {
344 532 : for (auto i : elem->side_index_range())
345 : {
346 : // Check if side is external
347 448 : if (elem->neighbor_ptr(i) == nullptr)
348 : {
349 : // If external get the boundary ids associated with this side
350 224 : std::vector<boundary_id_type> boundary_ids;
351 224 : auto side_range = sideset_map.equal_range(elem);
352 728 : for (const auto & itr : as_range(side_range))
353 504 : if (itr.second.first == i)
354 182 : boundary_ids.push_back(i);
355 : // get intersection of boundary_ids and _watertight_boundaries
356 : std::vector<boundary_id_type> intersections =
357 224 : findBoundaryOverlap(_watertight_boundaries, boundary_ids);
358 :
359 224 : bool no_specified_ids = boundary_ids.empty();
360 224 : bool specified_ids = !_watertight_boundaries.empty() && intersections.empty();
361 224 : std::string message;
362 224 : if (mesh->mesh_dimension() == 3)
363 336 : message = "Element " + std::to_string(elem->id()) +
364 168 : " contains an external face which has not been assigned to ";
365 : else
366 112 : message = "Element " + std::to_string(elem->id()) +
367 56 : " contains an external edge which has not been assigned to ";
368 224 : if (no_specified_ids)
369 42 : message = message + "a sideset";
370 182 : else if (specified_ids)
371 0 : message = message + "one of the specified sidesets";
372 224 : if ((no_specified_ids || specified_ids) && num_faces_without_sideset < _num_outputs)
373 : {
374 42 : _console << message << std::endl;
375 42 : num_faces_without_sideset++;
376 : }
377 224 : }
378 : }
379 14 : }
380 14 : std::string message;
381 14 : if (mesh->mesh_dimension() == 3)
382 14 : message = "Number of external element faces that have not been assigned to a sideset: " +
383 21 : std::to_string(num_faces_without_sideset);
384 : else
385 14 : message = "Number of external element edges that have not been assigned to a sideset: " +
386 21 : std::to_string(num_faces_without_sideset);
387 14 : diagnosticsLog(message, _check_watertight_sidesets, num_faces_without_sideset);
388 14 : }
389 :
390 : void
391 14 : MeshDiagnosticsGenerator::checkWatertightNodesets(const std::unique_ptr<MeshBase> & mesh) const
392 : {
393 : /*
394 : Diagnostic Overview:
395 : 1) Mesh precheck
396 : 2) Loop through all elements
397 : 3) Loop through all sides of that element
398 : 4) If side is external loop through its nodes
399 : 5) If node is not associated with any nodeset add to list
400 : 6) Print out node id
401 : */
402 14 : if (mesh->mesh_dimension() < 2)
403 0 : mooseError("The nodeset check only works for 2D and 3D meshes");
404 14 : auto & boundary_info = mesh->get_boundary_info();
405 14 : unsigned int num_nodes_without_nodeset = 0;
406 14 : std::set<dof_id_type> checked_nodes_id;
407 :
408 182 : for (const auto elem : mesh->active_element_ptr_range())
409 : {
410 532 : for (const auto i : elem->side_index_range())
411 : {
412 : // Check if side is external
413 448 : if (elem->neighbor_ptr(i) == nullptr)
414 : {
415 : // Side is external, now check nodes
416 224 : auto side = elem->side_ptr(i);
417 224 : const auto & node_list = side->get_nodes();
418 1008 : for (unsigned int j = 0; j < side->n_nodes(); j++)
419 : {
420 784 : const auto node = node_list[j];
421 784 : if (checked_nodes_id.count(node->id()))
422 28 : continue;
423 : // get vector of node's boundaries (in most cases it will only have one)
424 756 : std::vector<boundary_id_type> boundary_ids;
425 756 : boundary_info.boundary_ids(node, boundary_ids);
426 : std::vector<boundary_id_type> intersection =
427 756 : findBoundaryOverlap(_watertight_boundaries, boundary_ids);
428 :
429 756 : bool no_specified_ids = boundary_info.n_boundary_ids(node) == 0;
430 756 : bool specified_ids = !_watertight_boundaries.empty() && intersection.empty();
431 : std::string message =
432 1512 : "Node " + std::to_string(node->id()) +
433 756 : " is on an external boundary of the mesh, but has not been assigned to ";
434 756 : if (no_specified_ids)
435 14 : message = message + "a nodeset";
436 742 : else if (specified_ids)
437 0 : message = message + "one of the specified nodesets";
438 756 : if ((no_specified_ids || specified_ids) && num_nodes_without_nodeset < _num_outputs)
439 : {
440 14 : checked_nodes_id.insert(node->id());
441 14 : num_nodes_without_nodeset++;
442 14 : _console << message << std::endl;
443 : }
444 756 : }
445 224 : }
446 : }
447 14 : }
448 14 : std::string message;
449 28 : message = "Number of external nodes that have not been assigned to a nodeset: " +
450 42 : std::to_string(num_nodes_without_nodeset);
451 14 : diagnosticsLog(message, _check_watertight_nodesets, num_nodes_without_nodeset);
452 14 : }
453 :
454 : std::vector<boundary_id_type>
455 980 : MeshDiagnosticsGenerator::findBoundaryOverlap(
456 : const std::vector<boundary_id_type> & watertight_boundaries,
457 : std::vector<boundary_id_type> & boundary_ids) const
458 : {
459 : // Only the boundary_ids vector is sorted here. watertight_boundaries has to be sorted beforehand
460 : // Returns their intersection (elements that they share)
461 980 : std::sort(boundary_ids.begin(), boundary_ids.end());
462 980 : std::vector<boundary_id_type> boundary_intersection;
463 980 : std::set_intersection(watertight_boundaries.begin(),
464 : watertight_boundaries.end(),
465 : boundary_ids.begin(),
466 : boundary_ids.end(),
467 : std::back_inserter(boundary_intersection));
468 980 : return boundary_intersection;
469 0 : }
470 :
471 : void
472 46 : MeshDiagnosticsGenerator::checkElementVolumes(const std::unique_ptr<MeshBase> & mesh) const
473 : {
474 46 : unsigned int num_tiny_elems = 0;
475 46 : unsigned int num_negative_elems = 0;
476 46 : unsigned int num_big_elems = 0;
477 : // loop elements within the mesh (assumes replicated)
478 12346 : for (auto & elem : mesh->active_element_ptr_range())
479 : {
480 6150 : Real vol = elem->volume();
481 :
482 6150 : if (vol <= _min_volume)
483 : {
484 3 : if (num_tiny_elems < _num_outputs)
485 3 : _console << "Element with volume below threshold detected : \n"
486 3 : << "id " << elem->id() << " near point " << elem->vertex_average() << std::endl;
487 0 : else if (num_tiny_elems == _num_outputs)
488 0 : _console << "Maximum output reached, log is silenced" << std::endl;
489 3 : num_tiny_elems++;
490 : }
491 6150 : if (vol < 0)
492 : {
493 0 : if (num_negative_elems < _num_outputs)
494 0 : _console << "Element with negative volume detected : \n"
495 0 : << "id " << elem->id() << " near point " << elem->vertex_average() << std::endl;
496 0 : else if (num_negative_elems == _num_outputs)
497 0 : _console << "Maximum output reached, log is silenced" << std::endl;
498 0 : num_negative_elems++;
499 : }
500 6150 : if (vol >= _max_volume)
501 : {
502 3 : if (num_big_elems < _num_outputs)
503 3 : _console << "Element with volume above threshold detected : \n"
504 3 : << elem->get_info() << std::endl;
505 0 : else if (num_big_elems == _num_outputs)
506 0 : _console << "Maximum output reached, log is silenced" << std::endl;
507 3 : num_big_elems++;
508 : }
509 46 : }
510 46 : diagnosticsLog("Number of elements below prescribed minimum volume : " +
511 43 : std::to_string(num_tiny_elems),
512 46 : _check_element_volumes,
513 : num_tiny_elems);
514 43 : diagnosticsLog("Number of elements with negative volume : " + std::to_string(num_negative_elems),
515 43 : _check_element_volumes,
516 : num_negative_elems);
517 43 : diagnosticsLog("Number of elements above prescribed maximum volume : " +
518 40 : std::to_string(num_big_elems),
519 43 : _check_element_volumes,
520 : num_big_elems);
521 40 : }
522 :
523 : void
524 43 : MeshDiagnosticsGenerator::checkElementTypes(const std::unique_ptr<MeshBase> & mesh) const
525 : {
526 43 : std::set<subdomain_id_type> ids;
527 43 : mesh->subdomain_ids(ids);
528 : // loop on sub-domain
529 334 : for (auto & id : ids)
530 : {
531 : // ElemType defines an enum for geometric element types
532 291 : std::set<ElemType> types;
533 : // loop on elements within this sub-domain
534 6441 : for (auto & elem : mesh->active_subdomain_elements_ptr_range(id))
535 6441 : types.insert(elem->type());
536 :
537 291 : std::string elem_type_names = "";
538 585 : for (auto & elem_type : types)
539 294 : elem_type_names += " " + Moose::stringify(elem_type);
540 :
541 582 : _console << "Element type in subdomain " + mesh->subdomain_name(id) + " (" +
542 1164 : std::to_string(id) + ") :" + elem_type_names
543 291 : << std::endl;
544 291 : if (types.size() > 1)
545 3 : diagnosticsLog("Two different element types in subdomain " + std::to_string(id),
546 3 : _check_element_types,
547 : true);
548 291 : }
549 43 : }
550 :
551 : void
552 127 : MeshDiagnosticsGenerator::checkElementOverlap(const std::unique_ptr<MeshBase> & mesh) const
553 : {
554 : {
555 127 : unsigned int num_elem_overlaps = 0;
556 127 : auto pl = mesh->sub_point_locator();
557 : // loop on nodes, assumed replicated mesh
558 25247 : for (auto & node : mesh->node_ptr_range())
559 : {
560 : // find all the elements around this node
561 12560 : std::set<const Elem *> elements;
562 12560 : (*pl)(*node, elements);
563 :
564 70109 : for (auto & elem : elements)
565 : {
566 57549 : if (!elem->contains_point(*node))
567 0 : continue;
568 :
569 : // not overlapping inside the element if part of its nodes
570 57549 : bool found = false;
571 212107 : for (auto & elem_node : elem->node_ref_range())
572 212100 : if (*node == elem_node)
573 : {
574 57542 : found = true;
575 57542 : break;
576 : }
577 : // not overlapping inside the element if right on its side
578 57549 : bool on_a_side = false;
579 347905 : for (const auto & elem_side_index : elem->side_index_range())
580 290356 : if (elem->side_ptr(elem_side_index)->contains_point(*node, _non_conformality_tol))
581 153228 : on_a_side = true;
582 57549 : if (!found && !on_a_side)
583 : {
584 7 : num_elem_overlaps++;
585 7 : if (num_elem_overlaps < _num_outputs)
586 7 : _console << "Element overlap detected at : " << *node << std::endl;
587 0 : else if (num_elem_overlaps == _num_outputs)
588 0 : _console << "Maximum output reached, log is silenced" << std::endl;
589 : }
590 : }
591 12687 : }
592 :
593 127 : diagnosticsLog("Number of elements overlapping (node-based heuristics): " +
594 127 : Moose::stringify(num_elem_overlaps),
595 127 : _check_element_overlap,
596 : num_elem_overlaps);
597 127 : num_elem_overlaps = 0;
598 :
599 : // loop on all elements in mesh: assumes a replicated mesh
600 22043 : for (auto & elem : mesh->active_element_ptr_range())
601 : {
602 : // find all the elements around the centroid of this element
603 10958 : std::set<const Elem *> overlaps;
604 10958 : (*pl)(elem->vertex_average(), overlaps);
605 :
606 10958 : if (overlaps.size() > 1)
607 : {
608 7 : num_elem_overlaps++;
609 7 : if (num_elem_overlaps < _num_outputs)
610 7 : _console << "Element overlap detected with element : " << elem->id() << " near point "
611 7 : << elem->vertex_average() << std::endl;
612 0 : else if (num_elem_overlaps == _num_outputs)
613 0 : _console << "Maximum output reached, log is silenced" << std::endl;
614 : }
615 11085 : }
616 127 : diagnosticsLog("Number of elements overlapping (centroid-based heuristics): " +
617 127 : Moose::stringify(num_elem_overlaps),
618 127 : _check_element_overlap,
619 : num_elem_overlaps);
620 127 : }
621 127 : }
622 :
623 : void
624 39 : MeshDiagnosticsGenerator::checkNonPlanarSides(const std::unique_ptr<MeshBase> & mesh) const
625 : {
626 39 : unsigned int sides_non_planar = 0;
627 : // loop on all elements in mesh: assumes a replicated mesh
628 10805 : for (auto & elem : mesh->active_element_ptr_range())
629 : {
630 21553 : for (auto i : make_range(elem->n_sides()))
631 : {
632 16170 : auto side = elem->side_ptr(i);
633 16170 : std::vector<const Point *> nodes;
634 48594 : for (auto & node : side->node_ref_range())
635 32424 : nodes.emplace_back(&node);
636 :
637 16170 : if (nodes.size() <= 3)
638 16128 : continue;
639 : // First vector of the base
640 42 : const RealVectorValue v1 = *nodes[0] - *nodes[1];
641 :
642 : // Find another node so that we can form a basis. It should just be node 0, 1, 2
643 : // to form two independent vectors, but degenerate elements can make them aligned
644 42 : bool aligned = true;
645 42 : unsigned int third_node_index = 2;
646 42 : RealVectorValue v2;
647 84 : while (aligned && third_node_index < nodes.size())
648 : {
649 42 : v2 = *nodes[0] - *nodes[third_node_index++];
650 42 : aligned = MooseUtils::absoluteFuzzyEqual(v1 * v2 - v1.norm() * v2.norm(), 0);
651 : }
652 :
653 : // Degenerate element, could not find a third node that is not aligned
654 42 : if (aligned)
655 0 : continue;
656 :
657 42 : bool found_non_planar = false;
658 :
659 84 : for (auto node_offset : make_range(nodes.size() - 3))
660 : {
661 42 : RealVectorValue v3 = *nodes[0] - *nodes[node_offset + 3];
662 42 : bool planar = MooseUtils::absoluteFuzzyEqual(v2.cross(v1) * v3, 0);
663 42 : if (!planar)
664 7 : found_non_planar = true;
665 : }
666 :
667 42 : if (found_non_planar)
668 : {
669 7 : sides_non_planar++;
670 7 : if (sides_non_planar < _num_outputs)
671 7 : _console << "Nonplanar side detected for side " << i
672 7 : << " of element :" << elem->get_info() << std::endl;
673 0 : else if (sides_non_planar == _num_outputs)
674 0 : _console << "Maximum output reached, log is silenced" << std::endl;
675 : }
676 32298 : }
677 39 : }
678 39 : diagnosticsLog("Number of non-planar element sides detected: " +
679 39 : Moose::stringify(sides_non_planar),
680 39 : _check_non_planar_sides,
681 : sides_non_planar);
682 39 : }
683 :
684 : void
685 134 : MeshDiagnosticsGenerator::checkNonConformalMesh(const std::unique_ptr<MeshBase> & mesh) const
686 : {
687 134 : unsigned int num_nonconformal_nodes = 0;
688 134 : MeshBaseDiagnosticsUtils::checkNonConformalMesh(
689 134 : mesh, _console, _num_outputs, _non_conformality_tol, num_nonconformal_nodes);
690 134 : diagnosticsLog("Number of non-conformal nodes: " + Moose::stringify(num_nonconformal_nodes),
691 134 : _check_non_conformal_mesh,
692 : num_nonconformal_nodes);
693 134 : }
694 :
695 : void
696 198 : MeshDiagnosticsGenerator::checkNonConformalMeshFromAdaptivity(
697 : const std::unique_ptr<MeshBase> & mesh) const
698 : {
699 198 : unsigned int num_likely_AMR_created_nonconformality = 0;
700 198 : auto pl = mesh->sub_point_locator();
701 198 : pl->set_close_to_point_tol(_non_conformality_tol);
702 :
703 : // We have to make a copy because adding the new parent element to the mesh
704 : // will modify the mesh for the analysis of the next nodes
705 : // Make a copy of the mesh, add this element
706 198 : auto mesh_copy = mesh->clone();
707 198 : libMesh::MeshRefinement mesh_refiner(*mesh_copy);
708 :
709 : // loop on nodes, assumes a replicated mesh
710 38786 : for (auto & node : mesh->node_ptr_range())
711 : {
712 : // find all the elements around this node
713 19294 : std::set<const Elem *> elements_around;
714 19294 : (*pl)(*node, elements_around);
715 :
716 : // Keep track of the refined elements and the coarse element
717 19294 : std::set<const Elem *> fine_elements;
718 19294 : std::set<const Elem *> coarse_elements;
719 :
720 : // loop through the set of elements near this node
721 105102 : for (auto elem : elements_around)
722 : {
723 : // If the node is not part of this element's nodes, it is a
724 : // case of non-conformality
725 85808 : bool node_on_elem = false;
726 :
727 85808 : if (elem->get_node_index(node) != libMesh::invalid_uint)
728 : {
729 84628 : node_on_elem = true;
730 : // non-vertex nodes are not cause for the kind of non-conformality we are looking for
731 84628 : if (!elem->is_vertex(elem->get_node_index(node)))
732 10196 : continue;
733 : }
734 :
735 : // Keep track of all the elements this node is a part of. They are potentially the
736 : // 'fine' (refined) elements next to a coarser element
737 75612 : if (node_on_elem)
738 74432 : fine_elements.insert(elem);
739 : // Else, the node is not part of the element considered, so if the element had been part
740 : // of an AMR-created non-conformality, this element is on the coarse side
741 75612 : if (!node_on_elem)
742 1180 : coarse_elements.insert(elem);
743 : }
744 :
745 : // all the elements around contained the node as one of their nodes
746 : // if the coarse and refined sides are not stitched together, this check can fail,
747 : // as nodes that are physically near one element are not part of it because of the lack of
748 : // stitching (overlapping nodes)
749 19294 : if (fine_elements.size() == elements_around.size())
750 14905 : continue;
751 :
752 4389 : if (fine_elements.empty())
753 3730 : continue;
754 :
755 : // Depending on the type of element, we already know the number of elements we expect
756 : // to be part of this set of likely refined candidates for a given non-conformal node to
757 : // examine. We can only decide if it was born out of AMR if it's the center node of the face
758 : // of a coarse element near refined elements
759 659 : const auto elem_type = (*fine_elements.begin())->type();
760 806 : if ((elem_type == QUAD4 || elem_type == QUAD8 || elem_type == QUAD9) &&
761 147 : fine_elements.size() != 2)
762 112 : continue;
763 932 : else if ((elem_type == HEX8 || elem_type == HEX20 || elem_type == HEX27) &&
764 385 : fine_elements.size() != 4)
765 308 : continue;
766 273 : else if ((elem_type == TRI3 || elem_type == TRI6 || elem_type == TRI7) &&
767 34 : fine_elements.size() != 3)
768 0 : continue;
769 332 : else if ((elem_type == TET4 || elem_type == TET10 || elem_type == TET14) &&
770 93 : (fine_elements.size() % 2 != 0))
771 0 : continue;
772 :
773 : // only one coarse element in front of refined elements except for tets. Whatever we're
774 : // looking at is not the interface between coarse and refined elements
775 : // Tets are split on their edges (rather than the middle of a face) so there could be any
776 : // number of coarse elements in front of the node non-conformality created by refinement
777 239 : if (elem_type != TET4 && elem_type != TET10 && elem_type != TET14 && coarse_elements.size() > 1)
778 0 : continue;
779 :
780 : // There exists non-conformality, the node should have been a node of all the elements
781 : // that are close enough to the node, and it is not
782 :
783 : // Nodes of the tentative parent element
784 239 : std::vector<const Node *> tentative_coarse_nodes;
785 :
786 : // For quads and hexes, there is one (quad) or four (hexes) sides that are tied to this node
787 : // at the non-conformal interface between the refined elements and a coarse element
788 239 : if (elem_type == QUAD4 || elem_type == QUAD8 || elem_type == QUAD9 || elem_type == HEX8 ||
789 134 : elem_type == HEX20 || elem_type == HEX27)
790 : {
791 112 : const auto elem = *fine_elements.begin();
792 :
793 : // Find which sides (of the elements) the node considered is part of
794 112 : std::vector<Elem *> node_on_sides;
795 112 : unsigned int side_inside_parent = std::numeric_limits<unsigned int>::max();
796 368 : for (auto i : make_range(elem->n_sides()))
797 : {
798 368 : const auto side = elem->side_ptr(i);
799 368 : std::vector<const Node *> other_nodes_on_side;
800 368 : bool node_on_side = false;
801 1674 : for (const auto & elem_node : side->node_ref_range())
802 : {
803 1306 : if (*node == elem_node)
804 203 : node_on_side = true;
805 : else
806 1103 : other_nodes_on_side.emplace_back(&elem_node);
807 : }
808 : // node is on the side, but is it the side that goes away from the coarse element?
809 368 : if (node_on_side)
810 : {
811 : // if all the other nodes on this side are in one of the other potentially refined
812 : // elements, it's one of the side(s) (4 sides in a 3D hex for example) inside the
813 : // parent
814 203 : bool all_side_nodes_are_shared = true;
815 694 : for (const auto & other_node : other_nodes_on_side)
816 : {
817 491 : bool shared_with_a_fine_elem = false;
818 2337 : for (const auto & other_elem : fine_elements)
819 3201 : if (other_elem != elem &&
820 1355 : other_elem->get_node_index(other_node) != libMesh::invalid_uint)
821 554 : shared_with_a_fine_elem = true;
822 :
823 491 : if (!shared_with_a_fine_elem)
824 91 : all_side_nodes_are_shared = false;
825 : }
826 203 : if (all_side_nodes_are_shared)
827 : {
828 112 : side_inside_parent = i;
829 : // We stop examining sides, it does not matter which side we pick inside the parent
830 112 : break;
831 : }
832 : }
833 480 : }
834 112 : if (side_inside_parent == std::numeric_limits<unsigned int>::max())
835 0 : continue;
836 :
837 : // Gather the other potential elements in the refined element:
838 : // they are point neighbors of the node that is shared between all the elements flagged
839 : // for the non-conformality
840 : // Find shared node
841 112 : const auto interior_side = elem->side_ptr(side_inside_parent);
842 112 : const Node * interior_node = nullptr;
843 292 : for (const auto & other_node : interior_side->node_ref_range())
844 : {
845 292 : if (other_node == *node)
846 31 : continue;
847 261 : bool in_all_node_neighbor_elements = true;
848 1235 : for (auto other_elem : fine_elements)
849 : {
850 974 : if (other_elem->get_node_index(&other_node) == libMesh::invalid_uint)
851 298 : in_all_node_neighbor_elements = false;
852 : }
853 261 : if (in_all_node_neighbor_elements)
854 : {
855 112 : interior_node = &other_node;
856 112 : break;
857 : }
858 : }
859 : // Did not find interior node. Probably not AMR
860 112 : if (!interior_node)
861 0 : continue;
862 :
863 : // Add point neighbors of interior node to list of potentially refined elements
864 112 : std::set<const Elem *> all_elements;
865 112 : elem->find_point_neighbors(*interior_node, all_elements);
866 :
867 112 : if (elem_type == QUAD4 || elem_type == QUAD8 || elem_type == QUAD9)
868 : {
869 35 : const auto success = MeshCoarseningUtils::getFineElementsFromInteriorNode(
870 : *interior_node, *node, *elem, tentative_coarse_nodes, fine_elements);
871 35 : if (!success)
872 0 : continue;
873 35 : }
874 : // For hexes we first look at the fine-neighbors of the non-conformality
875 : // then the fine elements neighbors of the center 'node' of the potential parent
876 : else
877 : {
878 : // Get the coarse neighbor side to be able to recognize nodes that should become part of
879 : // the coarse parent
880 77 : const auto & coarse_elem = *coarse_elements.begin();
881 77 : unsigned short coarse_side_i = 0;
882 385 : for (const auto & coarse_side_index : coarse_elem->side_index_range())
883 : {
884 385 : const auto coarse_side_ptr = coarse_elem->side_ptr(coarse_side_index);
885 : // The side of interest is the side that contains the non-conformality
886 385 : if (!coarse_side_ptr->close_to_point(*node, 10 * _non_conformality_tol))
887 308 : continue;
888 : else
889 : {
890 77 : coarse_side_i = coarse_side_index;
891 77 : break;
892 : }
893 385 : }
894 77 : const auto coarse_side = coarse_elem->side_ptr(coarse_side_i);
895 :
896 : // We did not find the side of the coarse neighbor near the refined elements
897 : // Try again at another node
898 77 : if (!coarse_side)
899 0 : continue;
900 :
901 : // We cant directly use the coarse neighbor nodes
902 : // - The user might be passing a disjoint mesh
903 : // - There could two levels of refinement separating the coarse neighbor and its refined
904 : // counterparts
905 : // We use the fine element nodes
906 77 : unsigned int i = 0;
907 77 : tentative_coarse_nodes.resize(4);
908 385 : for (const auto & elem_1 : fine_elements)
909 3108 : for (const auto & coarse_node : elem_1->node_ref_range())
910 : {
911 2800 : bool node_shared = false;
912 14000 : for (const auto & elem_2 : fine_elements)
913 : {
914 11200 : if (elem_2 != elem_1)
915 8400 : if (elem_2->get_node_index(&coarse_node) != libMesh::invalid_uint)
916 3332 : node_shared = true;
917 : }
918 : // A node for the coarse parent will appear in only one fine neighbor (not shared)
919 : // and will lay on the side of the coarse neighbor
920 : // We only care about the coarse neighbor vertex nodes
921 3164 : if (!node_shared && coarse_side->close_to_point(coarse_node, _non_conformality_tol) &&
922 364 : elem_1->is_vertex(elem_1->get_node_index(&coarse_node)))
923 308 : tentative_coarse_nodes[i++] = &coarse_node;
924 : mooseAssert(i <= 5, "We went too far in this index");
925 : }
926 :
927 : // We did not find 4 coarse nodes. Mesh might be disjoint and the coarse element does not
928 : // contain the fine elements nodes we found
929 77 : if (i != 4)
930 0 : continue;
931 :
932 : // Need to order these nodes to form a valid quad / base of an hex
933 : // We go around the axis formed by the node and the interior node
934 77 : Point axis = *interior_node - *node;
935 77 : const auto start_circle = elem->vertex_average();
936 77 : MeshCoarseningUtils::reorderNodes(
937 : tentative_coarse_nodes, *interior_node, start_circle, axis);
938 77 : tentative_coarse_nodes.resize(8);
939 :
940 : // Use the neighbors of the fine elements that contain these nodes to get the vertex
941 : // nodes
942 385 : for (const auto & elem : fine_elements)
943 : {
944 : // Find the index of the coarse node for the starting element
945 308 : unsigned int node_index = 0;
946 770 : for (const auto & coarse_node : tentative_coarse_nodes)
947 : {
948 770 : if (elem->get_node_index(coarse_node) != libMesh::invalid_uint)
949 308 : break;
950 462 : node_index++;
951 : }
952 :
953 : // Get the neighbor element that is part of the fine elements to coarsen together
954 2156 : for (const auto & neighbor : elem->neighbor_ptr_range())
955 1848 : if (all_elements.count(neighbor) && !fine_elements.count(neighbor))
956 : {
957 : // Find the coarse node for the neighbor
958 308 : const Node * coarse_elem_node = nullptr;
959 1386 : for (const auto & fine_node : neighbor->node_ref_range())
960 : {
961 1386 : if (!neighbor->is_vertex(neighbor->get_node_index(&fine_node)))
962 0 : continue;
963 1386 : bool node_shared = false;
964 12474 : for (const auto & elem_2 : all_elements)
965 20790 : if (elem_2 != neighbor &&
966 9702 : elem_2->get_node_index(&fine_node) != libMesh::invalid_uint)
967 2926 : node_shared = true;
968 1386 : if (!node_shared)
969 : {
970 308 : coarse_elem_node = &fine_node;
971 308 : break;
972 : }
973 : }
974 : // Insert the coarse node at the right place
975 308 : tentative_coarse_nodes[node_index + 4] = coarse_elem_node;
976 : mooseAssert(node_index + 4 < tentative_coarse_nodes.size(), "Indexed too far");
977 : mooseAssert(coarse_elem_node, "Did not find last coarse element node");
978 : }
979 : }
980 77 : }
981 :
982 : // No need to separate fine elements near the non-conformal node and away from it
983 112 : fine_elements = all_elements;
984 224 : }
985 : // For TRI elements, we use the fine triangle element at the center of the potential
986 : // coarse triangle element
987 127 : else if (elem_type == TRI3 || elem_type == TRI6 || elem_type == TRI7)
988 : {
989 : // Find the center element
990 : // It's the only element that shares a side with both of the other elements near the node
991 : // considered
992 34 : const Elem * center_elem = nullptr;
993 136 : for (const auto refined_elem_1 : fine_elements)
994 : {
995 102 : unsigned int num_neighbors = 0;
996 408 : for (const auto refined_elem_2 : fine_elements)
997 : {
998 306 : if (refined_elem_1 == refined_elem_2)
999 102 : continue;
1000 204 : if (refined_elem_1->has_neighbor(refined_elem_2))
1001 136 : num_neighbors++;
1002 : }
1003 102 : if (num_neighbors >= 2)
1004 34 : center_elem = refined_elem_1;
1005 : }
1006 : // Did not find the center fine element, probably not AMR
1007 34 : if (!center_elem)
1008 0 : continue;
1009 : // Now get the tentative coarse element nodes
1010 136 : for (const auto refined_elem : fine_elements)
1011 : {
1012 102 : if (refined_elem == center_elem)
1013 34 : continue;
1014 314 : for (const auto & other_node : refined_elem->node_ref_range())
1015 344 : if (center_elem->get_node_index(&other_node) == libMesh::invalid_uint &&
1016 98 : refined_elem->is_vertex(refined_elem->get_node_index(&other_node)))
1017 68 : tentative_coarse_nodes.push_back(&other_node);
1018 : }
1019 :
1020 : // Get the final tentative new coarse element node, on the other side of the center
1021 : // element from the non-conformality
1022 34 : unsigned int center_side_opposite_node = std::numeric_limits<unsigned int>::max();
1023 136 : for (auto side_index : center_elem->side_index_range())
1024 102 : if (center_elem->side_ptr(side_index)->get_node_index(node) == libMesh::invalid_uint)
1025 34 : center_side_opposite_node = side_index;
1026 : const auto neighbor_on_other_side_of_opposite_center_side =
1027 34 : center_elem->neighbor_ptr(center_side_opposite_node);
1028 :
1029 : // Element is on a boundary, cannot form a coarse element
1030 34 : if (!neighbor_on_other_side_of_opposite_center_side)
1031 0 : continue;
1032 :
1033 34 : fine_elements.insert(neighbor_on_other_side_of_opposite_center_side);
1034 157 : for (const auto & tri_node : neighbor_on_other_side_of_opposite_center_side->node_ref_range())
1035 123 : if (neighbor_on_other_side_of_opposite_center_side->is_vertex(
1036 348 : neighbor_on_other_side_of_opposite_center_side->get_node_index(&tri_node)) &&
1037 225 : center_elem->side_ptr(center_side_opposite_node)->get_node_index(&tri_node) ==
1038 : libMesh::invalid_uint)
1039 34 : tentative_coarse_nodes.push_back(&tri_node);
1040 :
1041 : mooseAssert(center_side_opposite_node != std::numeric_limits<unsigned int>::max(),
1042 : "Did not find the side opposite the non-conformality");
1043 : mooseAssert(tentative_coarse_nodes.size() == 3,
1044 : "We are forming a coarsened triangle element");
1045 34 : }
1046 : // For TET elements, it's very different because the non-conformality does not happen inside
1047 : // of a face, but on an edge of one or more coarse elements
1048 93 : else if (elem_type == TET4 || elem_type == TET10 || elem_type == TET14)
1049 : {
1050 : // There are 4 tets on the tips of the coarsened tet and 4 tets inside
1051 : // let's identify all of them
1052 93 : std::set<const Elem *> tips_tets;
1053 93 : std::set<const Elem *> inside_tets;
1054 :
1055 : // pick a coarse element and work with its fine neighbors
1056 93 : const Elem * coarse_elem = nullptr;
1057 93 : std::set<const Elem *> fine_tets;
1058 113 : for (auto & coarse_one : coarse_elements)
1059 : {
1060 1285 : for (const auto & elem : fine_elements)
1061 : // for two levels of refinement across, this is not working
1062 : // we would need a "has_face_embedded_in_this_other_ones_face" routine
1063 1172 : if (elem->has_neighbor(coarse_one))
1064 279 : fine_tets.insert(elem);
1065 :
1066 113 : if (fine_tets.size())
1067 : {
1068 93 : coarse_elem = coarse_one;
1069 93 : break;
1070 : }
1071 : }
1072 : // There's no coarse element neighbor to a group of finer tets, not AMR
1073 93 : if (!coarse_elem)
1074 0 : continue;
1075 :
1076 : // There is one last point neighbor of the node that is sandwiched between two neighbors
1077 554 : for (const auto & elem : fine_elements)
1078 : {
1079 554 : int num_face_neighbors = 0;
1080 2216 : for (const auto & tet : fine_tets)
1081 1662 : if (tet->has_neighbor(elem))
1082 369 : num_face_neighbors++;
1083 554 : if (num_face_neighbors == 2)
1084 : {
1085 93 : fine_tets.insert(elem);
1086 93 : break;
1087 : }
1088 : }
1089 :
1090 : // There should be two other nodes with non-conformality near this coarse element
1091 : // Find both, as they will be nodes of the rest of the elements to add to the potential
1092 : // fine tet list. They are shared by two of the fine tets we have already found
1093 93 : std::set<const Node *> other_nodes;
1094 465 : for (const auto & tet_1 : fine_tets)
1095 : {
1096 3396 : for (const auto & node_1 : tet_1->node_ref_range())
1097 : {
1098 3024 : if (&node_1 == node)
1099 372 : continue;
1100 2652 : if (!tet_1->is_vertex(tet_1->get_node_index(&node_1)))
1101 1536 : continue;
1102 5580 : for (const auto & tet_2 : fine_tets)
1103 : {
1104 4464 : if (tet_2 == tet_1)
1105 1116 : continue;
1106 3348 : if (tet_2->get_node_index(&node_1) != libMesh::invalid_uint)
1107 : // check that it's near the coarse element as well
1108 1352 : if (coarse_elem->close_to_point(node_1, 10 * _non_conformality_tol))
1109 744 : other_nodes.insert(&node_1);
1110 : }
1111 : }
1112 : }
1113 : mooseAssert(other_nodes.size() == 2,
1114 : "Should find only two extra non-conformal nodes near the coarse element");
1115 :
1116 : // Now we can go towards this tip element next to two non-conformalities
1117 506 : for (const auto & tet_1 : fine_tets)
1118 : {
1119 2065 : for (const auto & neighbor : tet_1->neighbor_ptr_range())
1120 1652 : if (neighbor->get_node_index(*other_nodes.begin()) != libMesh::invalid_uint &&
1121 719 : neighbor->is_vertex(neighbor->get_node_index(*other_nodes.begin())) &&
1122 2698 : neighbor->get_node_index(*other_nodes.rbegin()) != libMesh::invalid_uint &&
1123 1979 : neighbor->is_vertex(neighbor->get_node_index(*other_nodes.rbegin())))
1124 327 : fine_tets.insert(neighbor);
1125 : }
1126 : // Now that the element next to the time is in the fine_tets, we can get the tip
1127 599 : for (const auto & tet_1 : fine_tets)
1128 : {
1129 2530 : for (const auto & neighbor : tet_1->neighbor_ptr_range())
1130 2024 : if (neighbor->get_node_index(*other_nodes.begin()) != libMesh::invalid_uint &&
1131 923 : neighbor->is_vertex(neighbor->get_node_index(*other_nodes.begin())) &&
1132 3385 : neighbor->get_node_index(*other_nodes.rbegin()) != libMesh::invalid_uint &&
1133 2462 : neighbor->is_vertex(neighbor->get_node_index(*other_nodes.rbegin())))
1134 438 : fine_tets.insert(neighbor);
1135 : }
1136 :
1137 : // Get the sandwiched tets between the tets we already found
1138 744 : for (const auto & tet_1 : fine_tets)
1139 3255 : for (const auto & neighbor : tet_1->neighbor_ptr_range())
1140 19246 : for (const auto & tet_2 : fine_tets)
1141 16642 : if (tet_1 != tet_2 && tet_2->has_neighbor(neighbor) && neighbor != coarse_elem)
1142 1640 : fine_tets.insert(neighbor);
1143 :
1144 : // tips tests are the only ones to have a node that is shared by no other tet in the group
1145 744 : for (const auto & tet_1 : fine_tets)
1146 : {
1147 651 : unsigned int unshared_nodes = 0;
1148 5943 : for (const auto & other_node : tet_1->node_ref_range())
1149 : {
1150 5292 : if (!tet_1->is_vertex(tet_1->get_node_index(&other_node)))
1151 2688 : continue;
1152 2604 : bool node_shared = false;
1153 20832 : for (const auto & tet_2 : fine_tets)
1154 18228 : if (tet_2 != tet_1 && tet_2->get_node_index(&other_node) != libMesh::invalid_uint)
1155 7998 : node_shared = true;
1156 2604 : if (!node_shared)
1157 279 : unshared_nodes++;
1158 : }
1159 651 : if (unshared_nodes == 1)
1160 279 : tips_tets.insert(tet_1);
1161 372 : else if (unshared_nodes == 0)
1162 372 : inside_tets.insert(tet_1);
1163 : else
1164 0 : mooseError("Did not expect a tet to have two unshared vertex nodes here");
1165 : }
1166 :
1167 : // Finally grab the last tip of the tentative coarse tet. It shares:
1168 : // - 3 nodes with the other tips, only one with each
1169 : // - 1 face with only one tet of the fine tet group
1170 : // - it has a node that no other fine tet shares (the tip node)
1171 465 : for (const auto & tet : inside_tets)
1172 : {
1173 1860 : for (const auto & neighbor : tet->neighbor_ptr_range())
1174 : {
1175 : // Check that it shares a face with no other potential fine tet
1176 1488 : bool shared_with_another_tet = false;
1177 11904 : for (const auto & tet_2 : fine_tets)
1178 : {
1179 10416 : if (tet_2 == tet)
1180 1488 : continue;
1181 8928 : if (tet_2->has_neighbor(neighbor))
1182 1581 : shared_with_another_tet = true;
1183 : }
1184 1488 : if (shared_with_another_tet)
1185 837 : continue;
1186 :
1187 : // Used to count the nodes shared with tip tets. Can only be 1 per tip tet
1188 651 : std::vector<const Node *> tip_nodes_shared;
1189 651 : unsigned int unshared_nodes = 0;
1190 5943 : for (const auto & other_node : neighbor->node_ref_range())
1191 : {
1192 5292 : if (!neighbor->is_vertex(neighbor->get_node_index(&other_node)))
1193 2688 : continue;
1194 :
1195 : // Check for being a node-neighbor of the 3 other tip tets
1196 10788 : for (const auto & tip_tet : tips_tets)
1197 : {
1198 8184 : if (neighbor == tip_tet)
1199 1116 : continue;
1200 :
1201 : // we could break here but we want to check that no other tip shares that node
1202 7068 : if (tip_tet->get_node_index(&other_node) != libMesh::invalid_uint)
1203 2139 : tip_nodes_shared.push_back(&other_node);
1204 : }
1205 : // Check for having a node shared with no other tet
1206 2604 : bool node_shared = false;
1207 20832 : for (const auto & tet_2 : fine_tets)
1208 18228 : if (tet_2 != neighbor && tet_2->get_node_index(&other_node) != libMesh::invalid_uint)
1209 7161 : node_shared = true;
1210 2604 : if (!node_shared)
1211 651 : unshared_nodes++;
1212 : }
1213 651 : if (tip_nodes_shared.size() == 3 && unshared_nodes == 1)
1214 93 : tips_tets.insert(neighbor);
1215 651 : }
1216 : }
1217 :
1218 : // append the missing fine tets (inside the coarse element, away from the node considered)
1219 : // into the fine elements set for the check on "did it refine the tentative coarse tet
1220 : // onto the same fine tets"
1221 93 : fine_elements.clear();
1222 465 : for (const auto & elem : tips_tets)
1223 372 : fine_elements.insert(elem);
1224 465 : for (const auto & elem : inside_tets)
1225 372 : fine_elements.insert(elem);
1226 :
1227 : // get the vertex of the coarse element from the tip tets
1228 465 : for (const auto & tip : tips_tets)
1229 : {
1230 930 : for (const auto & node : tip->node_ref_range())
1231 : {
1232 930 : bool outside = true;
1233 :
1234 930 : const auto id = tip->get_node_index(&node);
1235 930 : if (!tip->is_vertex(id))
1236 0 : continue;
1237 4650 : for (const auto & tet : inside_tets)
1238 3720 : if (tet->get_node_index(&node) != libMesh::invalid_uint)
1239 1488 : outside = false;
1240 930 : if (outside)
1241 : {
1242 372 : tentative_coarse_nodes.push_back(&node);
1243 : // only one tip node per tip tet
1244 372 : break;
1245 : }
1246 : }
1247 : }
1248 :
1249 93 : std::sort(tentative_coarse_nodes.begin(), tentative_coarse_nodes.end());
1250 186 : tentative_coarse_nodes.erase(
1251 93 : std::unique(tentative_coarse_nodes.begin(), tentative_coarse_nodes.end()),
1252 93 : tentative_coarse_nodes.end());
1253 :
1254 : // The group of fine elements ended up having less or more than 4 tips, so it's clearly
1255 : // not forming a coarse tetrahedral
1256 93 : if (tentative_coarse_nodes.size() != 4)
1257 0 : continue;
1258 186 : }
1259 : else
1260 : {
1261 0 : mooseInfo("Unsupported element type ",
1262 : elem_type,
1263 : ". Skipping detection for this node and all future nodes near only this "
1264 : "element type");
1265 0 : continue;
1266 : }
1267 :
1268 : // Check the fine element types: if not all the same then it's not uniform AMR
1269 1875 : for (auto elem : fine_elements)
1270 1636 : if (elem->type() != elem_type)
1271 0 : continue;
1272 :
1273 : // Check the number of coarse element nodes gathered
1274 1469 : for (const auto & check_node : tentative_coarse_nodes)
1275 1230 : if (check_node == nullptr)
1276 0 : continue;
1277 :
1278 : // Form a parent, of a low order type as we only have the extreme vertex nodes
1279 239 : std::unique_ptr<Elem> parent = Elem::build(Elem::first_order_equivalent_type(elem_type));
1280 239 : auto parent_ptr = mesh_copy->add_elem(parent.release());
1281 :
1282 : // Set the nodes to the coarse element
1283 1469 : for (auto i : index_range(tentative_coarse_nodes))
1284 1230 : parent_ptr->set_node(i, mesh_copy->node_ptr(tentative_coarse_nodes[i]->id()));
1285 :
1286 : // Refine this parent
1287 239 : parent_ptr->set_refinement_flag(Elem::REFINE);
1288 239 : parent_ptr->refine(mesh_refiner);
1289 239 : const auto num_children = parent_ptr->n_children();
1290 :
1291 : // Compare with the original set of elements
1292 : // We already know the child share the exterior node. If they share the same vertex
1293 : // average as the group of unrefined elements we will call this good enough for now
1294 : // For tetrahedral elements we cannot rely on the children all matching as the choice in
1295 : // the diagonal selection can be made differently. We'll just say 4 matching children is
1296 : // good enough for the heuristic
1297 239 : unsigned int num_children_match = 0;
1298 1875 : for (const auto & child : parent_ptr->child_ref_range())
1299 : {
1300 7120 : for (const auto & potential_children : fine_elements)
1301 6996 : if (MooseUtils::absoluteFuzzyEqual(child.vertex_average()(0),
1302 6996 : potential_children->vertex_average()(0),
1303 6996 : _non_conformality_tol) &&
1304 3346 : MooseUtils::absoluteFuzzyEqual(child.vertex_average()(1),
1305 3346 : potential_children->vertex_average()(1),
1306 10342 : _non_conformality_tol) &&
1307 1940 : MooseUtils::absoluteFuzzyEqual(child.vertex_average()(2),
1308 8936 : potential_children->vertex_average()(2),
1309 1940 : _non_conformality_tol))
1310 : {
1311 1512 : num_children_match++;
1312 1512 : break;
1313 : }
1314 : }
1315 :
1316 239 : if (num_children_match == num_children ||
1317 31 : ((elem_type == TET4 || elem_type == TET10 || elem_type == TET14) &&
1318 : num_children_match == 4))
1319 : {
1320 239 : num_likely_AMR_created_nonconformality++;
1321 239 : if (num_likely_AMR_created_nonconformality < _num_outputs)
1322 : {
1323 221 : _console << "Detected non-conformality likely created by AMR near" << *node
1324 221 : << Moose::stringify(elem_type)
1325 221 : << " elements that could be merged into a coarse element:" << std::endl;
1326 1713 : for (const auto & elem : fine_elements)
1327 1492 : _console << elem->id() << " ";
1328 221 : _console << std::endl;
1329 : }
1330 18 : else if (num_likely_AMR_created_nonconformality == _num_outputs)
1331 3 : _console << "Maximum log output reached, silencing output" << std::endl;
1332 : }
1333 57602 : }
1334 :
1335 198 : diagnosticsLog(
1336 381 : "Number of non-conformal nodes likely due to mesh refinement detected by heuristic: " +
1337 183 : Moose::stringify(num_likely_AMR_created_nonconformality),
1338 198 : _check_adaptivity_non_conformality,
1339 : num_likely_AMR_created_nonconformality);
1340 183 : pl->unset_close_to_point_tol();
1341 183 : }
1342 :
1343 : void
1344 47 : MeshDiagnosticsGenerator::checkLocalJacobians(const std::unique_ptr<MeshBase> & mesh) const
1345 : {
1346 47 : unsigned int num_bad_elem_qp_jacobians = 0;
1347 : // Get a high-ish order quadrature
1348 47 : auto qrule_dimension = mesh->mesh_dimension();
1349 47 : libMesh::QGauss qrule(qrule_dimension, FIFTH);
1350 :
1351 : // Use a constant monomial
1352 47 : const libMesh::FEType fe_type(CONSTANT, libMesh::MONOMIAL);
1353 :
1354 : // Initialize a basic constant monomial shape function everywhere
1355 47 : std::unique_ptr<libMesh::FEBase> fe_elem;
1356 47 : if (mesh->mesh_dimension() == 1)
1357 0 : fe_elem = std::make_unique<libMesh::FEMonomial<1>>(fe_type);
1358 47 : if (mesh->mesh_dimension() == 2)
1359 39 : fe_elem = std::make_unique<libMesh::FEMonomial<2>>(fe_type);
1360 : else
1361 8 : fe_elem = std::make_unique<libMesh::FEMonomial<3>>(fe_type);
1362 :
1363 47 : fe_elem->get_JxW();
1364 47 : fe_elem->attach_quadrature_rule(&qrule);
1365 :
1366 : // Check elements (assumes serialized mesh)
1367 6198 : for (const auto & elem : mesh->element_ptr_range())
1368 : {
1369 : // Handle mixed-dimensional meshes
1370 6151 : if (qrule_dimension != elem->dim())
1371 : {
1372 : // Re-initialize a quadrature
1373 0 : qrule_dimension = elem->dim();
1374 0 : qrule = libMesh::QGauss(qrule_dimension, FIFTH);
1375 :
1376 : // Re-initialize a monomial FE
1377 0 : if (elem->dim() == 1)
1378 0 : fe_elem = std::make_unique<libMesh::FEMonomial<1>>(fe_type);
1379 0 : if (elem->dim() == 2)
1380 0 : fe_elem = std::make_unique<libMesh::FEMonomial<2>>(fe_type);
1381 : else
1382 0 : fe_elem = std::make_unique<libMesh::FEMonomial<3>>(fe_type);
1383 :
1384 0 : fe_elem->get_JxW();
1385 0 : fe_elem->attach_quadrature_rule(&qrule);
1386 : }
1387 :
1388 : try
1389 : {
1390 6151 : fe_elem->reinit(elem);
1391 : }
1392 7 : catch (std::exception & e)
1393 : {
1394 7 : if (!strstr(e.what(), "Jacobian"))
1395 0 : throw;
1396 :
1397 7 : num_bad_elem_qp_jacobians++;
1398 7 : if (num_bad_elem_qp_jacobians < _num_outputs)
1399 7 : _console << "Bad Jacobian found in element " << elem->id() << " near point "
1400 7 : << elem->vertex_average() << std::endl;
1401 0 : else if (num_bad_elem_qp_jacobians == _num_outputs)
1402 0 : _console << "Maximum log output reached, silencing output" << std::endl;
1403 7 : }
1404 47 : }
1405 47 : diagnosticsLog("Number of elements with a bad Jacobian: " +
1406 47 : Moose::stringify(num_bad_elem_qp_jacobians),
1407 47 : _check_local_jacobian,
1408 : num_bad_elem_qp_jacobians);
1409 :
1410 47 : unsigned int num_bad_side_qp_jacobians = 0;
1411 : // Get a high-ish order side quadrature
1412 47 : auto qrule_side_dimension = mesh->mesh_dimension() - 1;
1413 47 : libMesh::QGauss qrule_side(qrule_side_dimension, FIFTH);
1414 :
1415 : // Use the side quadrature now
1416 47 : fe_elem->attach_quadrature_rule(&qrule_side);
1417 :
1418 : // Check element sides
1419 6198 : for (const auto & elem : mesh->element_ptr_range())
1420 : {
1421 : // Handle mixed-dimensional meshes
1422 6151 : if (int(qrule_side_dimension) != elem->dim() - 1)
1423 : {
1424 0 : qrule_side_dimension = elem->dim() - 1;
1425 0 : qrule_side = libMesh::QGauss(qrule_side_dimension, FIFTH);
1426 :
1427 : // Re-initialize a side FE
1428 0 : if (elem->dim() == 1)
1429 0 : fe_elem = std::make_unique<libMesh::FEMonomial<1>>(fe_type);
1430 0 : if (elem->dim() == 2)
1431 0 : fe_elem = std::make_unique<libMesh::FEMonomial<2>>(fe_type);
1432 : else
1433 0 : fe_elem = std::make_unique<libMesh::FEMonomial<3>>(fe_type);
1434 :
1435 0 : fe_elem->get_JxW();
1436 0 : fe_elem->attach_quadrature_rule(&qrule_side);
1437 : }
1438 :
1439 26915 : for (const auto & side : elem->side_index_range())
1440 : {
1441 : try
1442 : {
1443 20764 : fe_elem->reinit(elem, side);
1444 : }
1445 28 : catch (std::exception & e)
1446 : {
1447 : // In 2D dbg/devel modes libMesh could hit
1448 : // libmesh_assert_not_equal_to on a side reinit
1449 28 : if (!strstr(e.what(), "Jacobian") && !strstr(e.what(), "det != 0"))
1450 0 : throw;
1451 :
1452 28 : num_bad_side_qp_jacobians++;
1453 28 : if (num_bad_side_qp_jacobians < _num_outputs)
1454 28 : _console << "Bad Jacobian found in side " << side << " of element" << elem->id()
1455 28 : << " near point " << elem->vertex_average() << std::endl;
1456 0 : else if (num_bad_side_qp_jacobians == _num_outputs)
1457 0 : _console << "Maximum log output reached, silencing output" << std::endl;
1458 28 : }
1459 : }
1460 47 : }
1461 47 : diagnosticsLog("Number of element sides with bad Jacobians: " +
1462 47 : Moose::stringify(num_bad_side_qp_jacobians),
1463 47 : _check_local_jacobian,
1464 : num_bad_side_qp_jacobians);
1465 47 : }
1466 :
1467 : void
1468 15 : MeshDiagnosticsGenerator::checkNonMatchingEdges(const std::unique_ptr<MeshBase> & mesh) const
1469 : {
1470 : /*Algorithm Overview
1471 : 1)Prechecks
1472 : a)This algorithm only works for 3D so check for that first
1473 : 2)Loop
1474 : a)Loop through every element
1475 : b)For each element get the edges associated with it
1476 : c)For each edge check overlap with any edges of nearby elements
1477 : d)Have check to make sure the same pair of edges are not being tested twice for overlap
1478 : 3)Overlap check
1479 : a)Shortest line that connects both lines is perpendicular to both lines
1480 : b)A good overview of the math for finding intersecting lines can be found
1481 : here->paulbourke.net/geometry/pointlineplane/
1482 : */
1483 15 : if (mesh->mesh_dimension() != 3)
1484 : {
1485 0 : mooseWarning("The edge intersection algorithm only works with 3D meshes. "
1486 : "'examine_non_matching_edges' is skipped");
1487 0 : return;
1488 : }
1489 15 : if (!mesh->is_serial())
1490 0 : mooseError("Only serialized/replicated meshes are supported");
1491 15 : unsigned int num_intersecting_edges = 0;
1492 :
1493 : // Create map of element to bounding box to avoing reinitializing the same bounding box multiple
1494 : // times
1495 15 : std::unordered_map<Elem *, BoundingBox> bounding_box_map;
1496 1455 : for (const auto elem : mesh->active_element_ptr_range())
1497 : {
1498 1440 : const auto boundingBox = elem->loose_bounding_box();
1499 1440 : bounding_box_map.insert({elem, boundingBox});
1500 15 : }
1501 :
1502 15 : std::unique_ptr<PointLocatorBase> point_locator = mesh->sub_point_locator();
1503 15 : std::set<std::array<dof_id_type, 4>> overlapping_edges_nodes;
1504 1455 : for (const auto elem : mesh->active_element_ptr_range())
1505 : {
1506 : // loop through elem's nodes and find nearby elements with it
1507 1440 : std::set<const Elem *> candidate_elems;
1508 1440 : std::set<const Elem *> nearby_elems;
1509 10272 : for (unsigned int i = 0; i < elem->n_nodes(); i++)
1510 : {
1511 8832 : (*point_locator)(elem->point(i), candidate_elems);
1512 8832 : nearby_elems.insert(candidate_elems.begin(), candidate_elems.end());
1513 : }
1514 1440 : std::vector<std::unique_ptr<const Elem>> elem_edges(elem->n_edges());
1515 14688 : for (auto i : elem->edge_index_range())
1516 13248 : elem_edges[i] = elem->build_edge_ptr(i);
1517 37352 : for (const auto other_elem : nearby_elems)
1518 : {
1519 : // If they're the same element then there's no need to check for overlap
1520 35912 : if (elem->id() >= other_elem->id())
1521 18676 : continue;
1522 :
1523 17236 : std::vector<std::unique_ptr<const Elem>> other_edges(other_elem->n_edges());
1524 156364 : for (auto j : other_elem->edge_index_range())
1525 139128 : other_edges[j] = other_elem->build_edge_ptr(j);
1526 156364 : for (auto & edge : elem_edges)
1527 : {
1528 1402440 : for (auto & other_edge : other_edges)
1529 : {
1530 : // Get nodes from edges
1531 1263312 : const Node * n1 = edge->get_nodes()[0];
1532 1263312 : const Node * n2 = edge->get_nodes()[1];
1533 1263312 : const Node * n3 = other_edge->get_nodes()[0];
1534 1263312 : const Node * n4 = other_edge->get_nodes()[1];
1535 :
1536 : // Create array<dof_id_type, 4> to check against set
1537 1263312 : std::array<dof_id_type, 4> node_id_array = {n1->id(), n2->id(), n3->id(), n4->id()};
1538 1263312 : std::sort(node_id_array.begin(), node_id_array.end());
1539 :
1540 : // Check if the edges have already been added to our check_edges list
1541 1263312 : if (overlapping_edges_nodes.count(node_id_array))
1542 : {
1543 305 : continue;
1544 : }
1545 :
1546 : // Check element/edge type
1547 1263007 : if (edge->type() != EDGE2)
1548 : {
1549 0 : std::string element_message = "Edge of type " + Utility::enum_to_string(edge->type()) +
1550 0 : " was found in cell " + std::to_string(elem->id()) +
1551 0 : " which is of type " +
1552 0 : Utility::enum_to_string(elem->type()) + '\n' +
1553 : "The edge intersection check only works for EDGE2 "
1554 0 : "elements.\nThis message will not be output again";
1555 0 : mooseDoOnce(_console << element_message << std::endl);
1556 0 : continue;
1557 0 : }
1558 1263007 : if (other_edge->type() != EDGE2)
1559 0 : continue;
1560 :
1561 : // Now compare edge with other_edge
1562 1263007 : Point intersection_coords;
1563 1263007 : bool overlap = MeshBaseDiagnosticsUtils::checkFirstOrderEdgeOverlap(
1564 1263007 : *edge, *other_edge, intersection_coords, _non_matching_edge_tol);
1565 1263007 : if (overlap)
1566 : {
1567 : // Add the nodes that make up the 2 edges to the vector overlapping_edges_nodes
1568 14 : overlapping_edges_nodes.insert(node_id_array);
1569 14 : num_intersecting_edges += 2;
1570 14 : if (num_intersecting_edges < _num_outputs)
1571 : {
1572 : // Print error message
1573 14 : std::string elem_id = std::to_string(elem->id());
1574 14 : std::string other_elem_id = std::to_string(other_elem->id());
1575 14 : std::string x_coord = std::to_string(intersection_coords(0));
1576 14 : std::string y_coord = std::to_string(intersection_coords(1));
1577 14 : std::string z_coord = std::to_string(intersection_coords(2));
1578 28 : std::string message = "Intersecting edges found between elements " + elem_id +
1579 28 : " and " + other_elem_id + " near point (" + x_coord + ", " +
1580 14 : y_coord + ", " + z_coord + ")";
1581 14 : _console << message << std::endl;
1582 14 : }
1583 : }
1584 : }
1585 : }
1586 17236 : }
1587 1455 : }
1588 15 : diagnosticsLog("Number of intersecting element edges: " +
1589 15 : Moose::stringify(num_intersecting_edges),
1590 15 : _check_non_matching_edges,
1591 : num_intersecting_edges);
1592 15 : }
1593 :
1594 : void
1595 14 : MeshDiagnosticsGenerator::checkPolygons(const std::unique_ptr<MeshBase> & mesh) const
1596 : {
1597 14 : unsigned int num_polygons = 0;
1598 14 : unsigned int num_nonconvex = 0;
1599 14 : unsigned int num_nonplanar = 0;
1600 14 : unsigned int num_flat_consecutive_sides = 0;
1601 :
1602 252 : for (const auto & elem : mesh->element_ptr_range())
1603 119 : if (elem->type() == libMesh::C0POLYGON)
1604 : {
1605 119 : num_polygons++;
1606 119 : const auto n_nodes = elem->n_nodes();
1607 119 : Point base_top_dir(0, 0, 0);
1608 119 : bool nonconvex = false;
1609 119 : bool nonplanar = false;
1610 581 : for (const auto & i : make_range(n_nodes))
1611 : {
1612 462 : const auto n1 = elem->point(i);
1613 462 : const auto n2 = elem->point((i + 1) % n_nodes);
1614 462 : const auto n3 = elem->point((i + 2) % n_nodes);
1615 : // can't be const with unit
1616 462 : Point top_dir = (n2 - n1).cross(n3 - n2);
1617 :
1618 462 : if (top_dir.norm_sq() > 0 && base_top_dir.norm() == 0)
1619 : {
1620 119 : base_top_dir = top_dir.unit();
1621 119 : continue;
1622 : }
1623 343 : if (base_top_dir * top_dir < 0)
1624 70 : nonconvex = true;
1625 343 : if (top_dir.norm_sq() > 0)
1626 343 : top_dir = top_dir.unit();
1627 : else
1628 0 : num_flat_consecutive_sides++;
1629 413 : if (!MooseUtils::absoluteFuzzyEqual((top_dir - base_top_dir).norm_sq(), 0, TOLERANCE) &&
1630 413 : !MooseUtils::absoluteFuzzyEqual((top_dir + base_top_dir).norm_sq(), 0, TOLERANCE))
1631 0 : nonplanar = true;
1632 : }
1633 :
1634 119 : if (nonconvex)
1635 : {
1636 28 : num_nonconvex++;
1637 28 : if (num_nonconvex < _num_outputs)
1638 28 : _console << "Non convex C0 polygon detected:" << elem->get_info() << std::endl;
1639 0 : else if (num_nonconvex == _num_outputs)
1640 0 : _console << "Ouptut limit reached for non-convex polygons" << std::endl;
1641 : }
1642 119 : if (nonplanar)
1643 : {
1644 0 : num_nonplanar++;
1645 0 : if (num_nonconvex < _num_outputs)
1646 0 : _console << "Non planar C0 polygon detected:" << elem->get_info() << std::endl;
1647 0 : else if (num_nonconvex == _num_outputs)
1648 0 : _console << "Ouptut limit reached for non-planar polygons" << std::endl;
1649 : }
1650 14 : }
1651 :
1652 14 : if (!num_polygons)
1653 0 : mooseWarning("No C0 polygons in geometry: polyon check did nothing");
1654 : else
1655 : {
1656 14 : diagnosticsLog("Number of non convex polygons: " + Moose::stringify(num_nonconvex),
1657 14 : _check_polygons,
1658 : num_nonconvex);
1659 14 : diagnosticsLog("Number of non planar polygons: " + Moose::stringify(num_nonplanar),
1660 14 : _check_polygons,
1661 : num_nonplanar);
1662 14 : diagnosticsLog("Number of colinear consecutive sides of polygons: " +
1663 14 : Moose::stringify(num_flat_consecutive_sides),
1664 14 : _check_polygons,
1665 : num_flat_consecutive_sides);
1666 : }
1667 14 : }
1668 :
1669 : void
1670 1127 : MeshDiagnosticsGenerator::diagnosticsLog(std::string msg,
1671 : const MooseEnum & log_level,
1672 : bool problem_detected) const
1673 : {
1674 : mooseAssert(log_level != "NO_CHECK",
1675 : "We should not be outputting logs if the check had been disabled");
1676 1127 : if (log_level == "INFO" || !problem_detected)
1677 1106 : mooseInfoRepeated(msg);
1678 21 : else if (log_level == "WARNING")
1679 3 : mooseWarning(msg);
1680 18 : else if (log_level == "ERROR")
1681 18 : mooseError(msg);
1682 : else
1683 0 : mooseError("Should not reach here");
1684 1106 : }
|
