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