Merge pull request #5403 from OPM/feature-reservoir-surface-intersection-border

Feature reservoir surface intersection border

ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp

@@ -47,6 +47,9 @@
 #include "../cafHexInterpolator/cafHexInterpolator.h"
 #include "RivSectionFlattner.h"
 
+#include "RiaLogging.h"
+#include "clipper.hpp"
+
 cvf::ref<caf::DisplayCoordTransform> displayCoordTransform( const RimIntersection* intersection )
 {
     Rim3dView* rimView = nullptr;
@@ -177,8 +180,6 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays()
         std::array<cvf::Vec3d, 8> cellCorners;
         std::array<size_t, 8>     cornerIndices;
 
-        size_t startOfGeneratedTrianglesForNativeTriangles = outputTriangleVertices.size();
-
         for ( size_t ticIdx = 0; ticIdx < triIntersectedCellCandidates.size(); ++ticIdx )
         {
             size_t globalCellIdx = triIntersectedCellCandidates[ticIdx];
@@ -250,28 +251,6 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays()
                 }
             }
         }
-
-        // Add triangles for the part of the native triangle outside any gridcells
-
-        if ( startOfGeneratedTrianglesForNativeTriangles == outputTriangleVertices.size() )
-        {
-            // No triangles created, use the complete native triangle
-            outputTriangleVertices.push_back( cvf::Vec3f( p0 - displayModelOffset ) );
-            outputTriangleVertices.push_back( cvf::Vec3f( p1 - displayModelOffset ) );
-            outputTriangleVertices.push_back( cvf::Vec3f( p2 - displayModelOffset ) );
-
-            m_triangleToCellIdxMap.push_back( cvf::UNDEFINED_SIZE_T );
-
-            m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() );
-            m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() );
-            m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights() );
-        }
-        else
-        {
-            // Todo:
-            // Subtract the created triangles from the native triangle
-            // Add the remains
-        }
     }
 
     m_triangleVxes->assign( outputTriangleVertices );

ApplicationCode/ModelVisualization/Surfaces/RivSurfacePartMgr.cpp

@@ -111,7 +111,7 @@ void RivSurfacePartMgr::updateCellResultColor( size_t timeStepIndex )
                                                                                 m_intersectionFacesTextureCoords.p() );
     }
 
-    if ( m_nativeTrianglesPart.notNull() )
+    if ( false && m_nativeTrianglesPart.notNull() )
     {
         if ( !m_nativeVertexToCellIndexMap.size() )
         {
@@ -230,6 +230,8 @@ void RivSurfacePartMgr::appendIntersectionGeometryPartsToModel( cvf::ModelBasicL
         m_intersectionFaultGridLines->setTransform( scaleTransform );
         model->addPart( m_intersectionFaultGridLines.p() );
     }
+
+    appendNativeGeometryPartsToModel( model, scaleTransform );
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -241,9 +243,12 @@ void RivSurfacePartMgr::applySingleColor()
         caf::SurfaceEffectGenerator surfaceGen( cvf::Color4f( m_surfaceInView->surface()->color() ), caf::PO_1 );
         cvf::ref<cvf::Effect>       eff = surfaceGen.generateCachedEffect();
 
+        caf::SurfaceEffectGenerator surfaceGenBehind( cvf::Color4f( m_surfaceInView->surface()->color() ), caf::PO_2 );
+        cvf::ref<cvf::Effect>       effBehind = surfaceGenBehind.generateCachedEffect();
+
         if ( m_nativeTrianglesPart.notNull() )
         {
-            m_nativeTrianglesPart->setEffect( eff.p() );
+            m_nativeTrianglesPart->setEffect( effBehind.p() );
         }
 
         if ( m_intersectionFaces.notNull() )

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h

@@ -53,7 +53,7 @@ public:
 
     static std::vector<std::vector<cvf::Vec3d>>
         subtractPolygons( const std::vector<cvf::Vec3d>&              sourcePolygon,
-                          const std::vector<std::vector<cvf::Vec3d>>& polygonToSubtract );
+                          const std::vector<std::vector<cvf::Vec3d>>& polygonsToSubtract );
 
     enum ZInterpolationType
     {

ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -175,6 +175,31 @@ double GeometryTools::signedAreaPlanarPolygon( const cvf::Vec3d& planeNormal, co
     return signedArea;
 }
 
+//--------------------------------------------------------------------------------------------------
+/// This method below is more correct than the one above, both in naming and behaviour.
+/// Should be replaced, but is not done now to avoid possible sideeffects
+/// The difference is the sign of the area. The one below retuns correct sign according to the plane normal
+/// provided
+//--------------------------------------------------------------------------------------------------
+template <typename Vec3Type>
+double closestAxisSignedAreaPlanarPolygon( const cvf::Vec3d& planeNormal, const std::vector<Vec3Type>& polygon )
+{
+    int Z = cvf::GeometryTools::findClosestAxis( planeNormal );
+    int X = ( Z + 1 ) % 3;
+    int Y = ( Z + 2 ) % 3;
+
+    // Use Shoelace formula to calculate signed area.
+    // https://en.wikipedia.org/wiki/Shoelace_formula
+    double signedArea = 0.0;
+    for ( size_t i = 0; i < polygon.size(); ++i )
+    {
+        signedArea += ( polygon[( i + 1 ) % polygon.size()][X] + polygon[i][X] ) *
+                      ( polygon[( i + 1 ) % polygon.size()][Y] - polygon[i][Y] );
+    }
+
+    return ( planeNormal[Z] > 0 ) ? signedArea : -signedArea;
+}
+
 /*
    Determine the intersection point of two line segments
    From Paul Bourke, but modified to really handle coincident lines
@@ -880,12 +905,17 @@ bool EarClipTesselator::calculateTriangles( std::vector<size_t>* triangleIndices
         ++w; // v + 1;
         if ( w == m_polygonIndices.end() ) w = m_polygonIndices.begin(); // if (nv <= w) w = 0;
 
-        if ( isTriangleValid( u, v, w ) )
+        EarClipTesselator::TriangleStatus triStatus = calculateTriangleStatus( u, v, w );
+
+        if ( triStatus != INVALID_TRIANGLE )
         {
             // Indices of the vertices
-            triangleIndices->push_back( *u );
+            if ( triStatus == VALID_TRIANGLE )
-            triangleIndices->push_back( *v );
+            {
-            triangleIndices->push_back( *w );
+                triangleIndices->push_back( *u );
+                triangleIndices->push_back( *v );
+                triangleIndices->push_back( *w );
+            }
 
             // Remove v from remaining polygon
             m_polygonIndices.erase( v );
@@ -904,9 +934,9 @@ bool EarClipTesselator::calculateTriangles( std::vector<size_t>* triangleIndices
 /// Is this a valid triangle ? ( No points inside, and points not on a line. )
 //--------------------------------------------------------------------------------------------------
 
-bool EarClipTesselator::isTriangleValid( std::list<size_t>::const_iterator u,
+EarClipTesselator::TriangleStatus EarClipTesselator::calculateTriangleStatus( std::list<size_t>::const_iterator u,
-                                         std::list<size_t>::const_iterator v,
+                                                                              std::list<size_t>::const_iterator v,
-                                         std::list<size_t>::const_iterator w ) const
+                                                                              std::list<size_t>::const_iterator w ) const
 {
     CVF_ASSERT( m_X > -1 && m_Y > -1 );
 
@@ -914,9 +944,17 @@ bool EarClipTesselator::isTriangleValid( std::list<size_t>::const_iterator u,
     cvf::Vec3d B = ( *m_nodeCoords )[*v];
     cvf::Vec3d C = ( *m_nodeCoords )[*w];
 
-    if ( m_areaTolerance >
+    double mainAxisProjectedArea = ( ( B[m_X] - A[m_X] ) * ( C[m_Y] - A[m_Y] ) ) -
-         ( ( ( B[m_X] - A[m_X] ) * ( C[m_Y] - A[m_Y] ) ) - ( ( B[m_Y] - A[m_Y] ) * ( C[m_X] - A[m_X] ) ) ) )
+                                   ( ( B[m_Y] - A[m_Y] ) * ( C[m_X] - A[m_X] ) );
-        return false;
+    if ( -m_areaTolerance > mainAxisProjectedArea )
+    {
+        // Definite negative triangle
+        return INVALID_TRIANGLE;
+    }
+    else if ( fabs( mainAxisProjectedArea ) < m_areaTolerance )
+    {
+        return NEAR_ZERO_AREA_TRIANGLE;
+    }
 
     std::list<size_t>::const_iterator c;
     std::list<size_t>::const_iterator outside;
@@ -948,15 +986,16 @@ bool EarClipTesselator::isTriangleValid( std::list<size_t>::const_iterator u,
 
         cvf::Vec3d P = ( *m_nodeCoords )[*c];
 
-        if ( isPointInsideTriangle( A, B, C, P ) ) return false;
+        if ( isPointInsideTriangle( A, B, C, P ) ) return INVALID_TRIANGLE;
     }
 
-    return true;
+    return VALID_TRIANGLE;
 }
 
 //--------------------------------------------------------------------------------------------------
 /// Decides if a point P is inside of the triangle defined by A, B, C.
 /// By calculating the "double area" (cross product) of Corner to corner x Corner to point vectors
+/// If cross product is negative, the point is outside
 //--------------------------------------------------------------------------------------------------
 
 bool EarClipTesselator::isPointInsideTriangle( const cvf::Vec3d& A,
@@ -983,7 +1022,8 @@ bool EarClipTesselator::isPointInsideTriangle( const cvf::Vec3d& A,
     double aCROSSbp = ax * bpy - ay * bpx;
     double cCROSSap = cx * apy - cy * apx;
     double bCROSScp = bx * cpy - by * cpx;
-    double tol      = 0;
+    double tol      = -m_areaTolerance;
+
     return ( ( aCROSSbp >= tol ) && ( bCROSScp >= tol ) && ( cCROSSap >= tol ) );
 };
 

ApplicationCode/ReservoirDataModel/cvfGeometryTools.h

@@ -209,9 +209,16 @@ public:
     virtual bool calculateTriangles( std::vector<size_t>* triangles );
 
 protected:
-    bool isTriangleValid( std::list<size_t>::const_iterator u,
+    enum TriangleStatus
-                          std::list<size_t>::const_iterator v,
+    {
-                          std::list<size_t>::const_iterator w ) const;
+        INVALID_TRIANGLE,
+        NEAR_ZERO_AREA_TRIANGLE,
+        VALID_TRIANGLE
+    };
+    TriangleStatus calculateTriangleStatus( std::list<size_t>::const_iterator u,
+                                            std::list<size_t>::const_iterator v,
+                                            std::list<size_t>::const_iterator w ) const;
+
     bool isPointInsideTriangle( const cvf::Vec3d& A, const cvf::Vec3d& B, const cvf::Vec3d& C, const cvf::Vec3d& P ) const;
     double calculateProjectedPolygonArea() const;
 

ApplicationCode/UnitTests/cvfGeometryTools-Test.cpp

@@ -531,3 +531,44 @@ TEST( CellFaceIntersectionTst, PolygonAreaNormal3D )
         EXPECT_DOUBLE_EQ( 0.0, area.z() );
     }
 }
+
+TEST( EarClipTesselator, ErrorTest )
+{
+    std::vector<cvf::Vec3d> remainingPolygon{
+        cvf::Vec3d( 44.66, 20.17, 0 ),
+        cvf::Vec3d( 78.08, 35.26, 0 ),
+        cvf::Vec3d( 93.97, 35.83, 0 ),
+        cvf::Vec3d( 144.95, 44.42, 0 ),
+        cvf::Vec3d( 172.59, 39.73, 0 ),
+        cvf::Vec3d( 227.27, 24.01, 0 ),
+        cvf::Vec3d( 217.46, 45.72, 0 ),
+        cvf::Vec3d( 178.5, 57.61, 0 ),
+        cvf::Vec3d( 141.33, 63.82, 0 ),
+        cvf::Vec3d( 0, 0, 0 ),
+        cvf::Vec3d( 63.77, 0, 0 ),
+    };
+
+    double nativeTriangleArea = 21266;
+
+    cvf::EarClipTesselator tess;
+    tess.setNormal( cvf::Vec3d::Z_AXIS );
+    tess.setMinTriangleArea( 1e-3 * nativeTriangleArea );
+    cvf::Vec3dArray cvfNodes( remainingPolygon );
+    tess.setGlobalNodeArray( cvfNodes );
+
+    std::vector<size_t> polyIndexes;
+    for ( size_t idx = 0; idx < remainingPolygon.size(); ++idx )
+    {
+        polyIndexes.push_back( idx );
+    }
+    tess.setPolygonIndices( polyIndexes );
+
+    std::vector<size_t> triangleIndices;
+    bool                isTesselationOk = tess.calculateTriangles( &triangleIndices );
+
+    // CVF_ASSERT( isTesselationOk );
+    if ( !isTesselationOk )
+    {
+        // continue;
+    }
+}