OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-01-01 03:37:15 -06:00
Code Issues Packages Projects Releases Wiki Activity
db1fe84c33
ResInsight/ApplicationLibCode/UnitTests/cvfGeometryTools-Test.cpp
Magne Sjaastad cc292b197a
Result Divided by Area: Establish concept used to compute flow velocity and normalized trans (#7349) 
* Geometry Tools : Add convenience functions for polygon area

* #7232 Result Divided by Area: Add cell face result and show in GUI

Native support for flow rate is given by mass rate (mass per time) over a cell face. Add a derived result that takes flow rate divided by cell face area to get velocity (distance per time).

Add support for this concept on relevant native results, and indicate this result type in UI using a "/A" postfix

* Speed up divided-by-area calculations by using openmp

* Some refactoring of result data access.

* Make sure NNC data is scaled correctly in vector flow viz.

Co-authored-by: jonjenssen <jon@soundsoft.no>
2021-02-11 03:01:17 +01:00

601 lines
22 KiB
C++
Raw Blame History

/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) Statoil ASA
// Copyright (C) Ceetron Solutions AS
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "gtest/gtest.h"
#include "cvfLibCore.h"
#include "cvfLibGeometry.h"
#include "cvfLibRender.h"
#include "cvfLibViewing.h"
#include "cvfArrayWrapperConst.h"
#include "cvfArrayWrapperToEdit.h"
#include "cvfBoundingBoxTree.h"
#include "cvfGeometryTools.h"
#include <array>
using namespace cvf;
#if 0
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void ControlVolume::calculateCubeFaceStatus(const cvf::Vec3dArray& nodeCoords, double areaTolerance)
{
int cubeFace;
cvf::uint cubeFaceIndices[4];
for (cubeFace = 0; cubeFace < 6; ++cubeFace)
{
surfaceNodeIndices(static_cast<Defines::CubeFace>(cubeFace), cubeFaceIndices);
std::vector<const brv::Connection*> conns;
connections(static_cast<Defines::CubeFace>(cubeFace), &conns);
if (!conns.size())
{
m_cubeFaceStatus[cubeFace] = FREE_FACE;
}
else
{
double area = 0.5 * (nodeCoords[cubeFaceIndices[1]]-nodeCoords[cubeFaceIndices[0]] ^ nodeCoords[cubeFaceIndices[3]]-nodeCoords[cubeFaceIndices[0]]).length();
area += 0.5 * (nodeCoords[cubeFaceIndices[3]]-nodeCoords[cubeFaceIndices[2]] ^ nodeCoords[cubeFaceIndices[1]]-nodeCoords[cubeFaceIndices[2]]).length();
double totConnectionArea = 0;
size_t i;
for (i = 0; i < conns.size(); ++i)
{
totConnectionArea += conns[i]->brfArea();
}
if ( totConnectionArea < area - areaTolerance )
{
m_cubeFaceStatus[cubeFace] = PARTIALLY_COVERED;
}
else
{
m_cubeFaceStatus[cubeFace] = COMPLETELY_COVERED;
}
}
// Create a polygon to store the complete polygon of the faces
// not completely covered by connections
// This polygon will be filled with nodes later
if (m_cubeFaceStatus[cubeFace] != COMPLETELY_COVERED )
{
m_freeFacePolygons[cubeFace] = new std::list<std::pair<cvf::uint, bool> >;
}
}
}
#endif
template <typename NodeArrayType, typename NodeType, typename IndexType>
NodeType quadNormal( ArrayWrapperConst<NodeArrayType, NodeType> nodeCoords, const IndexType cubeFaceIndices[4] )
{
return ( nodeCoords[cubeFaceIndices[2]] - nodeCoords[cubeFaceIndices[0]] ) ^
( nodeCoords[cubeFaceIndices[3]] - nodeCoords[cubeFaceIndices[1]] );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> createVertices()
{
std::vector<cvf::Vec3d> vxs;
vxs.resize( 14, cvf::Vec3d::ZERO );
// clang-format off
vxs[ 0]= cvf::Vec3d( 0 , 0 , 0 );
vxs[ 1]= cvf::Vec3d( 1 , 0 , 0 );
vxs[ 2]= cvf::Vec3d( 1 , 1 , 0 );
vxs[ 3]= cvf::Vec3d( 0 , 1 , 0 );
vxs[ 4]= cvf::Vec3d(-0.4 ,-0.2 , 0.0 );
vxs[ 5]= cvf::Vec3d( 0.4 , 0.6 , 0.0 );
vxs[ 6]= cvf::Vec3d( 0.8 , 0.2 , 0.0 );
vxs[ 7]= cvf::Vec3d( 0.0 ,-0.6 , 0.0 );
vxs[ 8]= cvf::Vec3d( 1.0 , 1.2 , 0.0 );
vxs[ 9]= cvf::Vec3d( 1.4 , 0.8 , 0.0 );
vxs[10]= cvf::Vec3d( 0.4 ,-0.2 , 0.0 );
vxs[11]= cvf::Vec3d( 1.2 , 0.6 , 0.0 );
vxs[12]= cvf::Vec3d( 1.6 , 0.2 , 0.0 );
vxs[13]= cvf::Vec3d( 0.8 ,-0.6 , 0.0 );
// clang-format on
return vxs;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::array<cvf::uint, 4>> getCubeFaces()
{
std::vector<std::array<cvf::uint, 4>> cubeFaces;
cubeFaces.resize( 4 );
cubeFaces[0] = { 0, 1, 2, 3 };
cubeFaces[1] = { 4, 5, 6, 7 };
cubeFaces[2] = { 5, 8, 9, 6 };
cubeFaces[3] = { 10, 11, 12, 13 };
return cubeFaces;
}
std::ostream& operator<<( std::ostream& stream, std::vector<cvf::uint> v )
{
for ( size_t i = 0; i < v.size(); ++i )
{
stream << v[i] << " ";
}
return stream;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( CellFaceIntersectionTst, Intersection1 )
{
std::vector<cvf::Vec3d> nodes = createVertices();
std::vector<cvf::Vec3d> additionalVertices;
std::vector<std::vector<cvf::uint>> overlapPolygons;
auto faces = getCubeFaces();
EdgeIntersectStorage<cvf::uint> edgeIntersectionStorage;
edgeIntersectionStorage.setVertexCount( nodes.size() );
{
std::vector<cvf::uint> polygon;
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
faces[0].data(),
faces[1].data(),
1e-6 );
EXPECT_EQ( (size_t)5, polygon.size() );
EXPECT_EQ( (size_t)2, additionalVertices.size() );
EXPECT_TRUE( isOk );
overlapPolygons.push_back( polygon );
std::cout << polygon << std::endl;
}
{
std::vector<cvf::uint> polygon;
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
faces[0].data(),
faces[2].data(),
1e-6 );
EXPECT_EQ( (size_t)5, polygon.size() );
EXPECT_EQ( (size_t)4, additionalVertices.size() );
EXPECT_TRUE( isOk );
overlapPolygons.push_back( polygon );
std::cout << polygon << std::endl;
}
{
std::vector<cvf::uint> polygon;
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
faces[0].data(),
faces[3].data(),
1e-6 );
EXPECT_EQ( (size_t)3, polygon.size() );
EXPECT_EQ( (size_t)6, additionalVertices.size() );
EXPECT_TRUE( isOk );
overlapPolygons.push_back( polygon );
std::cout << polygon << std::endl;
}
nodes.insert( nodes.end(), additionalVertices.begin(), additionalVertices.end() );
std::vector<cvf::uint> basePolygon;
basePolygon.insert( basePolygon.begin(), faces[0].data(), &( faces[0].data()[4] ) );
for ( cvf::uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx )
{
GeometryTools::insertVertexInPolygon( &basePolygon, wrapArrayConst( &nodes ), vxIdx, 1e-6 );
}
EXPECT_EQ( (size_t)8, basePolygon.size() );
std::cout << "Bp: " << basePolygon << std::endl;
for ( size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx )
{
for ( cvf::uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx )
{
GeometryTools::insertVertexInPolygon( &overlapPolygons[pIdx], wrapArrayConst( &nodes ), vxIdx, 1e-6 );
}
if ( pIdx == 0 )
{
EXPECT_EQ( (size_t)5, overlapPolygons[pIdx].size() );
}
if ( pIdx == 1 )
{
EXPECT_EQ( (size_t)5, overlapPolygons[pIdx].size() );
}
if ( pIdx == 2 )
{
EXPECT_EQ( (size_t)4, overlapPolygons[pIdx].size() );
}
std::cout << "Op" << pIdx << ":" << overlapPolygons[pIdx] << std::endl;
}
Vec3d normal = quadNormal( wrapArrayConst( &nodes ), faces[0].data() );
std::vector<bool> faceOverlapPolygonWindingSameAsCubeFaceFlags;
faceOverlapPolygonWindingSameAsCubeFaceFlags.resize( overlapPolygons.size(), true );
{
std::vector<cvf::uint> freeFacePolygon;
bool hasHoles = false;
std::vector<std::vector<cvf::uint>*> overlapPolygonPtrs;
for ( size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx )
{
overlapPolygonPtrs.push_back( &( overlapPolygons[pIdx] ) );
}
GeometryTools::calculatePartiallyFreeCubeFacePolygon( wrapArrayConst( &nodes ),
wrapArrayConst( &basePolygon ),
normal,
overlapPolygonPtrs,
faceOverlapPolygonWindingSameAsCubeFaceFlags,
&freeFacePolygon,
&hasHoles );
EXPECT_EQ( (size_t)4, freeFacePolygon.size() );
EXPECT_FALSE( hasHoles );
std::cout << "FF1: " << freeFacePolygon << std::endl;
}
{
std::vector<cvf::uint> freeFacePolygon;
bool hasHoles = false;
std::vector<std::vector<cvf::uint>*> overlapPolygonPtrs;
for ( size_t pIdx = 0; pIdx < 1; ++pIdx )
{
overlapPolygonPtrs.push_back( &( overlapPolygons[pIdx] ) );
}
GeometryTools::calculatePartiallyFreeCubeFacePolygon( wrapArrayConst( &nodes ),
wrapArrayConst( &basePolygon ),
normal,
overlapPolygonPtrs,
faceOverlapPolygonWindingSameAsCubeFaceFlags,
&freeFacePolygon,
&hasHoles );
EXPECT_EQ( (size_t)9, freeFacePolygon.size() );
EXPECT_FALSE( hasHoles );
std::cout << "FF2: " << freeFacePolygon << std::endl;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( CellFaceIntersectionTst, Intersection )
{
std::vector<cvf::Vec3d> additionalVertices;
cvf::Vec3dArray nodes;
std::vector<size_t> polygon;
cvf::Array<size_t> ids;
size_t cv1CubeFaceIndices[4] = { 0, 1, 2, 3 };
size_t cv2CubeFaceIndices[4] = { 4, 5, 6, 7 };
nodes.resize( 8 );
nodes.setAll( cvf::Vec3d( 0, 0, 0 ) );
EdgeIntersectStorage<size_t> edgeIntersectionStorage;
edgeIntersectionStorage.setVertexCount( nodes.size() );
// Face 1
nodes[0] = cvf::Vec3d( 0, 0, 0 );
nodes[1] = cvf::Vec3d( 1, 0, 0 );
nodes[2] = cvf::Vec3d( 1, 1, 0 );
nodes[3] = cvf::Vec3d( 0, 1, 0 );
// Face 2
nodes[4] = cvf::Vec3d( 0, 0, 0 );
nodes[5] = cvf::Vec3d( 1, 0, 0 );
nodes[6] = cvf::Vec3d( 1, 1, 0 );
nodes[7] = cvf::Vec3d( 0, 1, 0 );
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
cv1CubeFaceIndices,
cv2CubeFaceIndices,
1e-6 );
EXPECT_EQ( (size_t)4, polygon.size() );
EXPECT_EQ( (size_t)0, additionalVertices.size() );
EXPECT_TRUE( isOk );
// Face 1
nodes[0] = cvf::Vec3d( 0, 0, 0 );
nodes[1] = cvf::Vec3d( 1, 0, 0 );
nodes[2] = cvf::Vec3d( 1, 1, 0 );
nodes[3] = cvf::Vec3d( 0, 1, 0 );
// Face 2
nodes[4] = cvf::Vec3d( 0.5, -0.25, 0 );
nodes[5] = cvf::Vec3d( 1.25, 0.5, 0 );
nodes[6] = cvf::Vec3d( 0.5, 1.25, 0 );
nodes[7] = cvf::Vec3d( -0.25, 0.5, 0 );
polygon.clear();
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
cv1CubeFaceIndices,
cv2CubeFaceIndices,
1e-6 );
EXPECT_EQ( (size_t)8, polygon.size() );
EXPECT_EQ( (size_t)8, additionalVertices.size() );
EXPECT_TRUE( isOk );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( CellFaceIntersectionTst, FreeFacePolygon )
{
std::vector<cvf::Vec3d> additionalVertices;
cvf::Vec3dArray nodes;
std::vector<size_t> polygon;
cvf::Array<size_t> ids;
size_t cv1CubeFaceIndices[4] = { 0, 1, 2, 3 };
size_t cv2CubeFaceIndices[4] = { 4, 5, 6, 7 };
nodes.resize( 8 );
nodes.setAll( cvf::Vec3d( 0, 0, 0 ) );
EdgeIntersectStorage<size_t> edgeIntersectionStorage;
edgeIntersectionStorage.setVertexCount( nodes.size() );
// Face 1
nodes[0] = cvf::Vec3d( 0, 0, 0 );
nodes[1] = cvf::Vec3d( 1, 0, 0 );
nodes[2] = cvf::Vec3d( 1, 1, 0 );
nodes[3] = cvf::Vec3d( 0, 1, 0 );
// Face 2
nodes[4] = cvf::Vec3d( 0, 0, 0 );
nodes[5] = cvf::Vec3d( 1, 0, 0 );
nodes[6] = cvf::Vec3d( 1, 1, 0 );
nodes[7] = cvf::Vec3d( 0, 1, 0 );
bool isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
cv1CubeFaceIndices,
cv2CubeFaceIndices,
1e-6 );
EXPECT_EQ( (size_t)4, polygon.size() );
EXPECT_EQ( (size_t)0, additionalVertices.size() );
EXPECT_TRUE( isOk );
std::vector<bool> faceOverlapPolygonWinding;
std::vector<std::vector<size_t>*> faceOverlapPolygons;
faceOverlapPolygons.push_back( &polygon );
faceOverlapPolygonWinding.push_back( true );
std::vector<size_t> partialFacePolygon;
bool hasHoles = false;
GeometryTools::calculatePartiallyFreeCubeFacePolygon( wrapArrayConst( &nodes ),
wrapArrayConst( cv1CubeFaceIndices, 4 ),
Vec3d( 0, 0, 1 ),
faceOverlapPolygons,
faceOverlapPolygonWinding,
&partialFacePolygon,
&hasHoles );
// Face 1
nodes[0] = cvf::Vec3d( 0, 0, 0 );
nodes[1] = cvf::Vec3d( 1, 0, 0 );
nodes[2] = cvf::Vec3d( 1, 1, 0 );
nodes[3] = cvf::Vec3d( 0, 1, 0 );
// Face 2
nodes[4] = cvf::Vec3d( 0.5, -0.25, 0 );
nodes[5] = cvf::Vec3d( 1.25, 0.5, 0 );
nodes[6] = cvf::Vec3d( 0.5, 1.25, 0 );
nodes[7] = cvf::Vec3d( -0.25, 0.5, 0 );
polygon.clear();
isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( &polygon,
&additionalVertices,
&edgeIntersectionStorage,
wrapArrayConst( &nodes ),
cv1CubeFaceIndices,
cv2CubeFaceIndices,
1e-6 );
EXPECT_EQ( (size_t)8, polygon.size() );
EXPECT_EQ( (size_t)8, additionalVertices.size() );
EXPECT_TRUE( isOk );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST( CellFaceIntersectionTst, PolygonAreaNormal3D )
{
// Test special cases with zero area
{
std::vector<cvf::Vec3d> vxs;
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_TRUE( area == cvf::Vec3d::ZERO );
}
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_TRUE( area == cvf::Vec3d::ZERO );
}
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 1 } );
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_TRUE( area == cvf::Vec3d::ZERO );
}
// Three points
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 1 } );
vxs.push_back( { 0, 1, 1 } );
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_DOUBLE_EQ( -0.5, area.x() );
EXPECT_DOUBLE_EQ( 0.0, area.y() );
EXPECT_DOUBLE_EQ( 0.0, area.z() );
}
// four identical points
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 0 } );
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_TRUE( area == cvf::Vec3d::ZERO );
}
// Square of four points
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 1 } );
vxs.push_back( { 0, 1, 1 } );
vxs.push_back( { 0, 1, 0 } );
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_DOUBLE_EQ( -1.0, area.x() );
EXPECT_DOUBLE_EQ( 0.0, area.y() );
EXPECT_DOUBLE_EQ( 0.0, area.z() );
}
// Square of four points + one point in center of square
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 1 } );
vxs.push_back( { 0, 1, 1 } );
vxs.push_back( { 0, 1, 0 } );
vxs.push_back( { 0, 0.5, 0.5 } ); // center of square
cvf::Vec3d area = GeometryTools::polygonAreaNormal3D( vxs );
EXPECT_DOUBLE_EQ( -0.75, area.x() );
EXPECT_DOUBLE_EQ( 0.0, area.y() );
EXPECT_DOUBLE_EQ( 0.0, area.z() );
}
// Area (float)
{
std::vector<cvf::Vec3f> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 2 } );
vxs.push_back( { 0, 2, 2 } );
vxs.push_back( { 0, 2, 0 } );
auto area = GeometryTools::polygonArea( vxs );
EXPECT_FLOAT_EQ( 4.0, area );
}
// Area (double)
{
std::vector<cvf::Vec3d> vxs;
vxs.push_back( { 0, 0, 0 } );
vxs.push_back( { 0, 0, 2 } );
vxs.push_back( { 0, 2, 2 } );
vxs.push_back( { 0, 2, 0 } );
auto area = GeometryTools::polygonArea( vxs );
EXPECT_DOUBLE_EQ( 4.0, area );
}
}
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;
}
}
Reference in New Issue View Git Blame Copy Permalink