Add check on side normal orientation flip internally

refs idaholab#25278

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

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();

framework/src/meshgenerators/MeshDiagnosticsGenerator.C

@@ -191,6 +191,76 @@ MeshDiagnosticsGenerator::checkSidesetsOrientation(const std::unique_ptr<MeshBas
                 " is consistently oriented with regards to the blocks it neighbors";
 
     diagnosticsLog(message, _check_sidesets_orientation, flipped_pairs.size());
+
+    // Now check that there is no sideset radically flipping from one side to another
+    // We'll consider pi / 2 to be the most steep angle we'll pass
+    unsigned int num_normals_flipping = 0;
+    Real steepest_side_angles = 1;
+    for (const auto index : index_range(side_tuples))
+    {
+      if (std::get<2>(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<const Elem> face = elem_ptr->build_side_ptr(side_id);
+      std::unique_ptr<FEBase> 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<Point> & 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<FEBase> 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<Point> & 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<MeshBase> &
 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);