From 5df5e4f4689a13efdf3a91ffdfed1aa73ad59533 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Tue, 19 Sep 2023 10:34:55 -0600 Subject: [PATCH] Add check on side normal orientation flip internally refs #25278 --- .../meshgenerators/MeshDiagnosticsGenerator.h | 8 +- framework/src/mesh/MooseMesh.C | 2 +- .../meshgenerators/MeshDiagnosticsGenerator.C | 74 ++++++++++++++++++- 3 files changed, 77 insertions(+), 7 deletions(-) diff --git a/framework/include/meshgenerators/MeshDiagnosticsGenerator.h b/framework/include/meshgenerators/MeshDiagnosticsGenerator.h index 49f02cb64c14..3a496f422d19 100644 --- a/framework/include/meshgenerators/MeshDiagnosticsGenerator.h +++ b/framework/include/meshgenerators/MeshDiagnosticsGenerator.h @@ -50,11 +50,11 @@ class MeshDiagnosticsGenerator : public MeshGenerator * Utility routine to output the final diagnostics level in the desired mode * @param msg the message to output * @param log_level the log level to output the message at - * @param may_error if set to false, prevents erroring from the log, despite the log level - * may_error is used to avoid erroring when the log is requested but there are no issues so it - * should just say "0 problems" with an info message + * @param problem_detected if set to false, prevents erroring from the log, despite the log level + * problem_detected is used to avoid erroring when the log is requested but there are no issues so + * it should just say "0 problems" with an info message */ - void diagnosticsLog(std::string msg, const MooseEnum & log_level, bool may_error) const; + void diagnosticsLog(std::string msg, const MooseEnum & log_level, bool problem_detected) const; /** * Utility routine to re-order a vector of nodes so that they form a valid quad diff --git a/framework/src/mesh/MooseMesh.C b/framework/src/mesh/MooseMesh.C index bd22fb49fd11..966918c05da6 100644 --- a/framework/src/mesh/MooseMesh.C +++ b/framework/src/mesh/MooseMesh.C @@ -3428,7 +3428,7 @@ MooseMesh::buildFiniteVolumeInfo() const boundary_ids.clear(); // We initialize the weights/other information in faceInfo. If the neighbor does not exist - // or is remote (so when we are on some sort of mesh boundary), we initiualize the ghost + // or is remote (so when we are on some sort of mesh boundary), we initialize the ghost // cell and use it to compute the weights corresponding to the faceInfo. if (!neighbor || neighbor == remote_elem) fi.computeBoundaryCoefficients(); diff --git a/framework/src/meshgenerators/MeshDiagnosticsGenerator.C b/framework/src/meshgenerators/MeshDiagnosticsGenerator.C index f02992478f28..ee187361b922 100644 --- a/framework/src/meshgenerators/MeshDiagnosticsGenerator.C +++ b/framework/src/meshgenerators/MeshDiagnosticsGenerator.C @@ -191,6 +191,76 @@ MeshDiagnosticsGenerator::checkSidesetsOrientation(const std::unique_ptr(side_tuples[index]) != bid) + continue; + const auto elem_ptr = mesh->elem_ptr(std::get<0>(side_tuples[index])); + const auto side_id = std::get<1>(side_tuples[index]); + + // Get side normal + const std::unique_ptr face = elem_ptr->build_side_ptr(side_id); + std::unique_ptr fe(FEBase::build(elem_ptr->dim(), FEType(elem_ptr->default_order()))); + QGauss qface(elem_ptr->dim() - 1, CONSTANT); + fe->attach_quadrature_rule(&qface); + const std::vector & normals = fe->get_normals(); + fe->reinit(elem_ptr, side_id); + mooseAssert(normals.size() == 1, "We expected only one normal here"); + const auto side_normal = normals[0]; + + // Compare to the sideset normals of neighbor sides in that sideset + for (const auto neighbor : elem_ptr->neighbor_ptr_range()) + if (neighbor) + for (const auto neigh_side_index : neighbor->side_index_range()) + { + // Check that the neighbor side is also in the sideset being examined + if (!boundary_info.has_boundary_id(neighbor, neigh_side_index, bid)) + continue; + + // We re-init everything for the neighbor in case it's a different dimension + std::unique_ptr fe_neighbor( + FEBase::build(neighbor->dim(), FEType(neighbor->default_order()))); + QGauss qface(neighbor->dim() - 1, CONSTANT); + fe_neighbor->attach_quadrature_rule(&qface); + const std::vector & neigh_normals = fe_neighbor->get_normals(); + fe_neighbor->reinit(neighbor, neigh_side_index); + mooseAssert(neigh_normals.size() == 1, "We expected only one normal here"); + const auto neigh_side_normal = neigh_normals[0]; + + // Check the angle by computing the dot product + if (neigh_side_normal * side_normal <= 0) + { + num_normals_flipping++; + steepest_side_angles = + std::min(std::acos(neigh_side_normal * side_normal), steepest_side_angles); + if (num_normals_flipping <= _num_outputs) + _console << "Side normals changed by more than pi/2 for sideset " + << sideset_full_name << " between side " << side_id << " of element " + << elem_ptr->id() << " and side " << neigh_side_index + << " of neighbor element " << neighbor->id() << std::endl; + else if (num_normals_flipping == _num_outputs + 1) + _console << "Maximum output reached for sideset normal flipping check. Silencing " + "output from now on" + << std::endl; + } + } + } + + if (num_normals_flipping) + message = "Sideset " + sideset_full_name + + " has two neighboring sides with a very large angle. Largest angle detected: " + + std::to_string(steepest_side_angles) + " rad."; + else + message = "Sideset " + sideset_full_name + + " does not appear to have side-to-neighbor-side orientation flips. All neighbor " + "sides normal differ by less than pi/2"; + + diagnosticsLog(message, _check_sidesets_orientation, num_normals_flipping); } } @@ -1213,11 +1283,11 @@ MeshDiagnosticsGenerator::checkLocalJacobians(const std::unique_ptr & void MeshDiagnosticsGenerator::diagnosticsLog(std::string msg, const MooseEnum & log_level, - bool may_error) const + bool problem_detected) const { mooseAssert(log_level != "NO_CHECK", "We should not be outputting logs if the check had been disabled"); - if (log_level == "INFO" || !may_error) + if (log_level == "INFO" || !problem_detected) mooseInfoRepeated(msg); else if (log_level == "WARNING") mooseWarning(msg);