Merge pull request #26852 from GiudGiud/PR_rt2

Misc ray tracing changes

modules/heat_transfer/src/userobjects/ViewFactorBase.C

@@ -278,7 +278,7 @@ ViewFactorBase::checkAndNormalizeViewFactor()
   Real max_sum_deviation_after_normalization = maxDevRowSum();
   Real max_reciprocity_deviation_after_normalization = maxDevReciprocity();
 
-  // print symmary
+  // print summary
   _console << std::endl;
   _console << COLOR_CYAN << "Summary of the view factor computation"
            << "\n"

modules/heat_transfer/src/userobjects/ViewFactorRayStudy.C

@@ -11,6 +11,7 @@
 
 // Local includes
 #include "ViewFactorRayBC.h"
+#include "GeometryUtils.h"
 
 // libMesh includes
 #include "libmesh/parallel_algebra.h"
@@ -265,6 +266,7 @@ ViewFactorRayStudy::generateRays()
   reserveRayBuffer(num_local_rays);
 
   Point direction;
+  unsigned int num_rays_skipped = 0;
 
   // loop through all starting points and spawn rays from each for each point and angle
   for (const auto & start_elem : _start_elems)
@@ -312,6 +314,59 @@ ViewFactorRayStudy::generateRays()
                                 : std::cos(_2d_aq_angles[l]) * _2d_aq_weights[l];
         const auto start_weight = start_elem._weights[start_i] * awf;
 
+        // Skip the ray if it exists the domain through the non-planar side it is starting from.
+        // We do not expect there are any neighbor elements to track it on if it exits the
+        // non-planar side.
+        bool intersection_found = false;
+        if (_is_3d && start_elem._start_elem &&
+            !start_elem._start_elem->neighbor_ptr(start_elem._incoming_side) &&
+            sideIsNonPlanar(start_elem._start_elem, start_elem._incoming_side))
+        {
+          // Find edge on side that is 'in front' of the future ray
+          Point intersection_point(std::numeric_limits<Real>::max(), -1, -1);
+          const auto side_elem = start_elem._start_elem->side_ptr(start_elem._incoming_side);
+          Point proj_dir;
+          for (const auto edge_i : side_elem->side_index_range())
+          {
+            const auto edge_1 = side_elem->side_ptr(edge_i);
+            // Project direction onto (start_point, node 1, node 2)
+            const auto d1 = *edge_1->node_ptr(0) - start_elem._points[start_i];
+            const auto d2 = *edge_1->node_ptr(1) - start_elem._points[start_i];
+            const auto d1_unit = d1.unit();
+            const auto d2_unit = d2.unit();
+            // If the starting point is aligned with the edge, it wont cross it
+            if (MooseUtils::absoluteFuzzyEqual(std::abs(d1_unit * d2_unit), 1))
+              continue;
+            const auto normal = (d1_unit.cross(d2_unit)).unit();
+
+            // One of the nodes must be in front of the start point following the direction
+            if (d1 * direction < 0 && d2 * direction < 0)
+              continue;
+
+            proj_dir = (direction - (direction * normal) * normal).unit();
+
+            // Only the side of interest will have the projected direction in between d1 and d2
+            if ((proj_dir * d2_unit > d1_unit * d2_unit) &&
+                (proj_dir * d1_unit > d1_unit * d2_unit))
+            {
+              const auto dist = geom_utils::distanceFromLine(
+                  start_elem._points[start_i], *edge_1->node_ptr(0), *edge_1->node_ptr(1));
+              // Ortho-normalize the base on the plane
+              intersection_point = start_elem._points[start_i] + dist * proj_dir;
+              intersection_found = true;
+              break;
+            }
+          }
+
+          // Skip the ray if it goes out of the element
+          const auto grazing_dir = (intersection_point - start_elem._points[start_i]).unit();
+          if (intersection_found && inward_normal * direction < inward_normal * grazing_dir)
+          {
+            num_rays_skipped++;
+            continue;
+          }
+        }
+
         // Acquire a Ray and fill with the starting information
         std::shared_ptr<Ray> ray = acquireRay();
         ray->setStart(
@@ -324,6 +379,10 @@ ViewFactorRayStudy::generateRays()
         moveRayToBuffer(ray);
       }
   }
+  if (num_rays_skipped)
+    mooseInfo(num_rays_skipped,
+              " rays were skipped as they exited the mesh at their starting point through "
+              "non-planar sides.");
 }
 
 void

modules/heat_transfer/test/tests/view_factors/tests

@@ -63,4 +63,14 @@
     design = 'RayTracingViewFactor.md'
     requirement = 'The system shall compute view factors for unobstructed, planar surfaces in three-dimensional meshes using ray tracing.'
   []
+
+  [ray3D_nonplanar]
+    type = RunApp
+    input = 'view_factor_3d_non_planar_face.i'
+    expect_out = '732 rays were skipped as they exited the mesh at their starting point through non-planar sides.'
+    design = 'RayTracingViewFactor.md'
+    requirement = 'The system shall be able to skip rays that exit the mesh when starting from non-planar faces in three-dimensional problems.'
+    recover = false
+    mesh_mode = REPLICATED
+  []
 []

modules/heat_transfer/test/tests/view_factors/view_factor_3d_non_planar_face.i

@@ -0,0 +1,107 @@
+[GlobalParams]
+  view_factor_object_name = rt_vf
+[]
+
+[Mesh]
+  [with_a_non_planar]
+    type = ElementGenerator
+    nodal_positions = '-1 -1 0
+                       1 -1 0
+                       1 1 0
+                       -1 1 0
+                       -1 -1 2
+                       1 -1 2
+                       1 1 2
+                       -1 1 1'
+    element_connectivity = '0 1 2 3 4 5 6 7'
+    elem_type = 'HEX8'
+  []
+  [sides]
+    type = AllSideSetsByNormalsGenerator
+    input = with_a_non_planar
+  []
+  [rename]
+    type = RenameBoundaryGenerator
+    input = sides
+    old_boundary = '1 2 3 4 5 6'
+    new_boundary = 'bottom front left back right top'
+  []
+[]
+
+[UserObjects]
+  [vf_study]
+    type = ViewFactorRayStudy
+    execute_on = INITIAL
+    boundary = 'bottom front left back right top'
+    face_order = CONSTANT
+    polar_quad_order = 3
+    azimuthal_quad_order = 200
+    face_type = GAUSS
+    warn_non_planar = false
+  []
+
+  [rt_vf]
+    type = RayTracingViewFactor
+    boundary = 'bottom front left back right top'
+    execute_on = INITIAL
+    ray_study_name = vf_study
+    normalize_view_factor = false
+  []
+[]
+
+[RayBCs]
+  [vf]
+    type = ViewFactorRayBC
+    boundary = 'left right front back bottom top'
+  []
+[]
+
+## For convenience, the "view_factor_object_name" for these
+## PPs are set in global params for switching between methods
+[Postprocessors]
+  [top_bottom]
+    type = ViewFactorPP
+    from_boundary = top
+    to_boundary = bottom
+  []
+  [top_left]
+    type = ViewFactorPP
+    from_boundary = top
+    to_boundary = left
+  []
+  [top_right]
+    type = ViewFactorPP
+    from_boundary = top
+    to_boundary = right
+  []
+  [top_front]
+    type = ViewFactorPP
+    from_boundary = top
+    to_boundary = front
+  []
+  [top_back]
+    type = ViewFactorPP
+    from_boundary = top
+    to_boundary = back
+  []
+  [sum]
+    type = ParsedPostprocessor
+    function = 'top_back + top_bottom + top_front + top_right + top_left'
+    pp_names = 'top_back top_bottom top_front top_right top_left'
+  []
+[]
+
+[Problem]
+  solve = false
+[]
+
+[Executioner]
+  type = Steady
+  [Quadrature] # higher order quadrature for unobstructed
+    order = SECOND
+  []
+[]
+
+[Outputs]
+  csv = true
+[]

modules/ray_tracing/include/raytracing/TraceRayTools.h

@@ -11,6 +11,7 @@
 
 // MOOSE includes
 #include "StaticallyAllocatedSet.h"
+#include "Conversion.h"
 
 // Local includes
 #include "DebugRay.h"
@@ -127,7 +128,7 @@ bool isWithinSegment(const Point & segment1,
  * @param neighbor_set The set to fill the neighbors into
  * @param untested_set Set for internal use
  * @param next_untested_set Set for internal use
- * @param active_neighbor_children Temprorary vector for use in the search
+ * @param active_neighbor_children Temporary vector for use in the search
  */
 void findPointNeighbors(
     const Elem * const elem,
@@ -162,8 +163,9 @@ void findEdgeNeighbors(
  * allocated set and accepts a functor for the rejection/acceptance of an element.
  *
  * Returns the active neighbors that fit the criteria of keep_functor.
+ * Does not the return the current element (elem)
  *
- * @param elem The element
+ * @param elem The element for which we are searching for the neighbors
  * @param neighbor_set The set to fill the neighbors into
  * @param untested_set Set for internal use
  * @param next_untested_set Set for internal use
@@ -186,16 +188,15 @@ findNeighbors(
   untested_set.clear();
   next_untested_set.clear();
 
-  neighbor_set.insert(elem);
   untested_set.insert(elem);
 
   while (!untested_set.empty())
   {
     // Loop over all the elements in the patch that haven't already been tested
-    for (const Elem * elem : untested_set)
-      for (auto current_neighbor : elem->neighbor_ptr_range())
-        if (current_neighbor &&
-            current_neighbor != remote_elem) // we have a real neighbor on elem side
+    for (const Elem * const test_elem : untested_set)
+      for (const auto * const current_neighbor : test_elem->neighbor_ptr_range())
+        if (current_neighbor && current_neighbor != remote_elem &&
+            current_neighbor != elem) // we have a real neighbor on elem side
         {
           if (current_neighbor->active()) // ... if it is active
           {
@@ -209,10 +210,11 @@ findNeighbors(
           else // ... the neighbor is *not* active,
           {    // ... so add *all* neighboring
                // active children that touch p
-            current_neighbor->active_family_tree_by_neighbor(active_neighbor_children, elem);
+            current_neighbor->active_family_tree_by_neighbor(active_neighbor_children, test_elem);
 
             for (const Elem * current_child : active_neighbor_children)
-              if (!neighbor_set.contains(current_child) && keep_functor(current_child))
+              if (!neighbor_set.contains(current_child) && keep_functor(current_child) &&
+                  current_child != elem)
               {
                 next_untested_set.insert(current_child);
                 neighbor_set.insert(current_child);
@@ -594,9 +596,11 @@ sideIntersectedByLine(const Elem * elem,
   if (segment_vertex != SEGMENT_VERTEX_NONE)
   {
     intersected_extrema.setVertex(T::side_nodes_map[side][segment_vertex]);
-    mooseAssert(intersected_extrema.vertexPoint(elem).absolute_fuzzy_equals(intersection_point,
-                                                                            TRACE_TOLERANCE),
-                "Doesn't intersect vertex");
+    mooseAssert(
+        intersected_extrema.vertexPoint(elem).absolute_fuzzy_equals(intersection_point,
+                                                                    3 * TRACE_TOLERANCE),
+        "Doesn't intersect vertex at: " + Moose::stringify(intersected_extrema.vertexPoint(elem)) +
+            " tentative intersection point: " + Moose::stringify(intersection_point));
   }
 
   return intersected;
@@ -860,7 +864,7 @@ sideIntersectedByLine(const Elem * elem,
  */
 bool isTraceableElem(const Elem * elem);
 /**
- * @return Whether or not the element is traceable with adapativity
+ * @return Whether or not the element is traceable with adaptivity
  */
 bool isAdaptivityTraceableElem(const Elem * elem);
 

modules/ray_tracing/src/raytracing/TraceRayTools.C

@@ -129,7 +129,7 @@ findPointNeighbors(
   // Helper for avoiding extraneous allocation when building side elements
   std::unique_ptr<const Elem> side_helper;
 
-  auto contains_point = [&point, &info, &side_helper](const Elem * const candidate)
+  auto contains_point = [&point, &info, &side_helper, &elem](const Elem * const candidate)
   {
     if (candidate->contains_point(point))
     {
@@ -141,7 +141,8 @@ findPointNeighbors(
           sides.push_back(s);
       }
 
-      if (!sides.empty())
+      // Dont add the local element
+      if (!sides.empty() && candidate != elem)
       {
         info.emplace_back(candidate, std::move(sides));
         return true;
@@ -167,7 +168,7 @@ findPointNeighbors(
   std::set<const Elem *> point_neighbors;
   elem->find_point_neighbors(point, point_neighbors);
   for (const auto & point_neighbor : point_neighbors)
-    if (!neighbor_set.contains(point_neighbor))
+    if (!neighbor_set.contains(point_neighbor) && point_neighbor != elem)
       mooseError("Missed a point neighbor");
 #endif
 }
@@ -241,7 +242,8 @@ findNodeNeighbors(
   elem->find_point_neighbors(*node, point_neighbors);
   for (const auto & point_neighbor : point_neighbors)
     for (const auto & neighbor_node : point_neighbor->node_ref_range())
-      if (node == &neighbor_node && !neighbor_set.contains(point_neighbor))
+      if (node == &neighbor_node && !neighbor_set.contains(point_neighbor) &&
+          point_neighbor != elem)
         mooseError("Missed a node neighbor");
 #endif
 }

modules/ray_tracing/src/userobjects/RayTracingStudy.C

@@ -527,6 +527,7 @@ RayTracingStudy::subdomainHMaxSetup()
         "tolerances used in computing ray intersections. This warning suggests that the\n"
         "approximate element size is not a good approximation. This is likely due to poor\n"
         "element aspect ratios.\n\n"
+        "This warning is only output for the first element affected.\n"
         "To disable this warning, set warn_subdomain_hmax = false.\n";
 
     for (const auto & elem : *_mesh.getActiveLocalElementRange())
@@ -537,13 +538,19 @@ RayTracingStudy::subdomainHMaxSetup()
 
       const auto hmax_rel = hmax / max_hmax;
       if (hmax_rel < 1.e-2 || hmax_rel > 1.e2)
-        mooseWarning(
-            warn_prefix, "Element hmax varies significantly from subdomain hmax.\n", warn_suffix);
+        mooseDoOnce(mooseWarning(warn_prefix,
+                                 "Element hmax varies significantly from subdomain hmax.\n",
+                                 warn_suffix,
+                                 "First element affected:\n",
+                                 Moose::stringify(*elem)););
 
       const auto h_rel = max_hmax / hmin;
       if (h_rel > 1.e2)
-        mooseWarning(
-            warn_prefix, "Element hmin varies significantly from subdomain hmax.\n", warn_suffix);
+        mooseDoOnce(mooseWarning(warn_prefix,
+                                 "Element hmin varies significantly from subdomain hmax.\n",
+                                 warn_suffix,
+                                 "First element affected:\n",
+                                 Moose::stringify(*elem)););
     }
   }
 }