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 "BoundingBoxIntersectionHelper.h"
      11             : 
      12             : // libMesh includes
      13             : #include "libmesh/bounding_box.h"
      14             : #include "libmesh/face_quad4.h"
      15             : #include "libmesh/cell_hex8.h"
      16             : 
      17             : #include "TraceRayTools.h"
      18             : #include "DebugRay.h"
      19             : 
      20             : using namespace libMesh;
      21             : 
      22        1535 : BoundingBoxIntersectionHelper::BoundingBoxIntersectionHelper(const BoundingBox & bbox,
      23        1535 :                                                              const unsigned int dim)
      24        1535 :   : _comm(), _mesh(std::make_unique<Mesh>(_comm, dim))
      25             : {
      26             :   // Add the nodes that represent our bounding box
      27        1535 :   _mesh->add_point(Point(bbox.min()(0), bbox.min()(1), bbox.min()(2)), 0, 0);
      28        1535 :   _mesh->add_point(Point(bbox.max()(0), bbox.min()(1), bbox.min()(2)), 1, 0);
      29        1535 :   if (dim > 1)
      30             :   {
      31        1493 :     _mesh->add_point(Point(bbox.max()(0), bbox.max()(1), bbox.min()(2)), 2, 0);
      32        1493 :     _mesh->add_point(Point(bbox.min()(0), bbox.max()(1), bbox.min()(2)), 3, 0);
      33             : 
      34        1493 :     if (dim == 3)
      35             :     {
      36        1036 :       _mesh->add_point(Point(bbox.min()(0), bbox.min()(1), bbox.max()(2)), 4, 0);
      37        1036 :       _mesh->add_point(Point(bbox.max()(0), bbox.min()(1), bbox.max()(2)), 5, 0);
      38        1036 :       _mesh->add_point(Point(bbox.max()(0), bbox.max()(1), bbox.max()(2)), 6, 0);
      39        1036 :       _mesh->add_point(Point(bbox.min()(0), bbox.max()(1), bbox.max()(2)), 7, 0);
      40             :     }
      41             :   }
      42             : 
      43             :   // Build the element and set the nodes accordingly
      44             :   Elem * elem;
      45        1535 :   if (dim == 1)
      46          42 :     elem = Elem::build(EDGE2).release();
      47        1493 :   else if (dim == 2)
      48         457 :     elem = Elem::build(QUAD4).release();
      49             :   else
      50        1036 :     elem = Elem::build(HEX8).release();
      51        1535 :   elem = _mesh->add_elem(elem);
      52        1535 :   elem->subdomain_id() = 0;
      53             :   elem->processor_id(0);
      54        1535 :   elem->set_id() = 0;
      55       11735 :   for (unsigned int n = 0; n < _mesh->n_nodes(); ++n)
      56       10200 :     elem->set_node(n, _mesh->node_ptr(n));
      57             : 
      58             :   // All done with the mesh
      59             :   _mesh->skip_partitioning(true);
      60        1535 :   _mesh->prepare_for_use();
      61             : 
      62             :   // This costs a little bit so why not cache it now
      63        1535 :   _hmax = elem->hmax();
      64        1535 : }
      65             : 
      66             : Point
      67     1251090 : BoundingBoxIntersectionHelper::intersection(const Point & point, const Point & direction) const
      68             : {
      69             :   mooseAssert(std::abs(1.0 - direction.norm()) < TOLERANCE, "Unnormalized direction");
      70             : 
      71     1251090 :   const Elem * elem = _mesh->elem_ptr(0);
      72             : 
      73             :   // We're going to check all possible intersections and keep the one that's furthest away
      74     1251090 :   Point best_intersection_point = RayTracingCommon::invalid_point;
      75             :   Real best_intersection_distance = -std::numeric_limits<Real>::max();
      76             : 
      77             :   // For use in the intersection routines
      78             :   Point intersection_point;
      79             :   ElemExtrema intersected_extrema;
      80             :   Real intersection_distance;
      81             :   bool intersected;
      82             : 
      83             :   // Check intersection with sides
      84     8720190 :   for (const auto s : elem->side_index_range())
      85             :   {
      86     7469100 :     intersection_point = RayTracingCommon::invalid_point;
      87             :     intersected_extrema.invalidate();
      88             : 
      89     7469100 :     if (elem->dim() == 3) // 3D: Intersect HEX8
      90     7394868 :       intersected = TraceRayTools::sideIntersectedByLine<Hex8>(elem,
      91             :                                                                point,
      92             :                                                                direction,
      93             :                                                                s,
      94             :                                                                _hmax,
      95             :                                                                intersection_point,
      96             :                                                                intersection_distance,
      97             :                                                                intersected_extrema,
      98     7394868 :                                                                _hmax
      99             : #ifdef DEBUG_RAY_INTERSECTIONS
     100             :                                                                ,
     101             :                                                                false
     102             : #endif
     103             :       );
     104       74232 :     else if (elem->dim() == 2) // 2D: Intersect QUAD4
     105             :     {
     106       74016 :       intersected = TraceRayTools::sideIntersectedByLine<Quad4>(elem,
     107             :                                                                 point,
     108             :                                                                 direction,
     109             :                                                                 s,
     110             :                                                                 _hmax,
     111             :                                                                 intersection_point,
     112             :                                                                 intersection_distance,
     113             :                                                                 intersected_extrema,
     114       74016 :                                                                 _hmax
     115             : #ifdef DEBUG_RAY_INTERSECTIONS
     116             :                                                                 ,
     117             :                                                                 false
     118             : #endif
     119             :       );
     120             :     }
     121             :     else // 1D: see if side point is between point and end far away
     122             :     {
     123             :       // Dummy end point outside of the element
     124             :       const Point dummy_end = point + 1.01 * direction * _hmax;
     125             : 
     126         216 :       intersected = TraceRayTools::isWithinSegment(point, dummy_end, elem->point(s));
     127         216 :       if (intersected)
     128             :       {
     129         216 :         intersection_distance = (elem->point(s) - point).norm();
     130         216 :         intersection_point = elem->point(s);
     131             :       }
     132             :     }
     133             : 
     134     7469100 :     if (intersected && intersection_distance > best_intersection_distance)
     135             :     {
     136             :       best_intersection_distance = intersection_distance;
     137     1260660 :       best_intersection_point = intersection_point;
     138             :     }
     139             :   }
     140             : 
     141     1251090 :   return best_intersection_point;
     142             : }