seems to work

src/coreComponents/mesh/PerforationData.cpp

@@ -239,17 +239,31 @@
     // compute the bounding box of the element
     real64 boxDims[ 3 ];
     computationalGeometry::getBoundingBox( ei,
-                                          subRegion.nodeList(),
-                                          nodeManager.referencePosition(),
-                                          boxDims );
+                                           subRegion.nodeList(),
+                                           nodeManager.referencePosition(),
+                                           boxDims );
 
     // dx and dz from bounding box
     dx = boxDims[ 0 ];
     dy = boxDims[ 1 ];
+    dz = boxDims[ 2 ];
 
-    // dz is computed as vol / (dx * dy)
-    dz  = subRegion.getElementVolume()[ei];
-    dz /= dx * dy;
+    if( dx < 1e-10 )
+    {
+      dx  = subRegion.getElementVolume()[ei];
+      dx /= dy * dz;
+    }
+    else if( dy < 1e-10 )
+    {
+      dy  = subRegion.getElementVolume()[ei];
+      dy /= dx * dz;
+    }
+    else
+    {
+      // dz is computed as vol / (dx * dy)
+      dz  = subRegion.getElementVolume()[ei];
+      dz /= dx * dy;
+    }
   } );
 }
 

src/coreComponents/mesh/utilities/ComputationalGeometry.hpp

@@ -434,104 +434,104 @@
 }
 
 template< typename POLYGON_TYPE, typename POINT_TYPE >
 bool isPointInPolygon2d( POLYGON_TYPE const & polygon, integer n, POINT_TYPE const & point )
 {
   integer count = 0;
 
   for( integer i = 0; i < n; i++ )
   {
     auto const & p1 = polygon[i];
     auto const & p2 = polygon[(i + 1) % n];
 
     if((point[1] > std::min( p1[1], p2[1] )) &&
        (point[1] <= std::max( p1[1], p2[1] )) &&
        (point[0] <= std::max( p1[0], p2[0] )))
     {
       real64 const xIntersect = (point[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0];
       if( std::abs( p1[0] - p2[0] ) < 1e-10 || point[0] <= xIntersect )
       {
         count++;
       }
     }
   }
 
   return count % 2 == 1;
 }
 
 template< typename POLYGON_TYPE, typename POINT_TYPE >
 bool isPointInPolygon3d( POLYGON_TYPE const & polygon, integer const n, POINT_TYPE const & point )
 {
   // Check if the point lies in the plane of the polygon
   auto const & p0 = polygon[0];
   POINT_TYPE normal = {0, 0, 0};
   for( integer i = 1; i < n - 1; i++ )
   {
     auto const & p1 = polygon[i];
     auto const & p2 = polygon[i + 1];
     normal[0] += (p1[1] - p0[1]) * (p2[2] - p0[2]) - (p1[2] - p0[2]) * (p2[1] - p0[1]);
     normal[1] += (p1[2] - p0[2]) * (p2[0] - p0[0]) - (p1[0] - p0[0]) * (p2[2] - p0[2]);
     normal[2] += (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p1[1] - p0[1]) * (p2[0] - p0[0]);
   }
 
   real64 d = -(normal[0] * p0[0] + normal[1] * p0[1] + normal[2] * p0[2]);
   real64 dist = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2] + d;
 
   if( std::abs( dist ) > 1e-6 )
   {
     return false;
   }
 
   // Determine the dominant component of the normal vector
   int dominantIndex = 0;
   if( std::abs( normal[1] ) > std::abs( normal[0] ))
   {
     dominantIndex = 1;
   }
   if( std::abs( normal[2] ) > std::abs( normal[dominantIndex] ))
   {
     dominantIndex = 2;
   }
 
   // Project the polygon and the point onto a 2D plane
   POLYGON_TYPE projectedPolygon( n );
   POINT_TYPE projectedPoint;
   for( integer i = 0; i < n; i++ )
   {
     projectedPolygon[i][0] = polygon[i][1];
     projectedPolygon[i][1] = polygon[i][2];
   }
   if( dominantIndex == 0 )  // X is dominant, project onto YZ plane
   {
     for( int i = 0; i < n; i++ )
     {
       projectedPolygon[i][0] = polygon[i][1];
       projectedPolygon[i][1] = polygon[i][2];
     }
     projectedPoint[0] = point[1];
     projectedPoint[1] = point[2];
   }
   else if( dominantIndex == 1 )  // Y is dominant, project onto XZ plane
   {
     for( int i = 0; i < n; i++ )
     {
       projectedPolygon[i][0] = polygon[i][0];
       projectedPolygon[i][1] = polygon[i][2];
     }
     projectedPoint[0] = point[0];
     projectedPoint[1] = point[2];
   }
   else  // Z is dominant, project onto XY plane
   {
     for( int i = 0; i < n; i++ )
     {
       projectedPolygon[i][0] = polygon[i][0];
       projectedPolygon[i][1] = polygon[i][1];
     }
     projectedPoint[0] = point[0];
     projectedPoint[1] = point[1];
   }
 
   return isPointInPolygon2d( projectedPolygon, n, projectedPoint );
 }
 
 /**
@@ -946,8 +946,6 @@
                      arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & pointCoordinates,
                      VEC_TYPE && boxDims )
 {
-  std::cout << "getBoundingBox " << elemIndex << std::endl;
-
   // This holds the min coordinates of the set in each direction
   R1Tensor minCoords = { LvArray::NumericLimits< real64 >::max,
                          LvArray::NumericLimits< real64 >::max,
@@ -957,7 +955,6 @@
   LvArray::tensorOps::fill< 3 >( boxDims, LvArray::NumericLimits< real64 >::lowest );
 
   // loop over all the vertices of the element to get the min and max coords
-  std::cout << pointIndices[elemIndex].size() << std::endl;
   for( localIndex a = 0; a < pointIndices[elemIndex].size(); ++a )
   {
     localIndex const id = pointIndices( elemIndex, a );