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