OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

riTRANXYZ: Paving the way: Renaming, Removing obsolete code, Store fault index for each cell face, Added prototype code for the calculation

Browse Source
This commit is contained in:
Jacob Støren 2014-08-18 18:30:52 +02:00
parent ab483fee81
commit 147fbffe55
12 changed files with 241 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
ApplicationCode/ReservoirDataModel/RigCell.cpp
View File

@ -243,16 +243,20 @@ cvf::Vec3d RigCell::faceCenter(cvf::StructGridInterface::FaceType face) const
}
//--------------------------------------------------------------------------------------------------
///
/// Returns an area vector for the cell face. The direction is the face normal, and the length is
/// equal to the face area (projected to the plane represented by the diagonal in case of warp)
/// The components of this area vector are equal to the area of the face projection onto
/// the corresponding plane.
/// See http://geomalgorithms.com/a01-_area.html
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigCell::faceNormal(cvf::StructGridInterface::FaceType face) const
cvf::Vec3d RigCell::faceNormalWithAreaLenght(cvf::StructGridInterface::FaceType face) const
{
cvf::ubyte faceVertexIndices[4];
cvf::StructGridInterface::cellFaceVertexIndices(face, faceVertexIndices);
const std::vector<cvf::Vec3d>& nodeCoords = m_hostGrid->mainGrid()->nodes();
return ( nodeCoords[m_cornerIndices[faceVertexIndices[2]]] - nodeCoords[m_cornerIndices[faceVertexIndices[0]]]) ^
( nodeCoords[m_cornerIndices[faceVertexIndices[3]]] - nodeCoords[m_cornerIndices[faceVertexIndices[1]]]);
return 0.5*( nodeCoords[m_cornerIndices[faceVertexIndices[2]]] - nodeCoords[m_cornerIndices[faceVertexIndices[0]]]) ^
( nodeCoords[m_cornerIndices[faceVertexIndices[3]]] - nodeCoords[m_cornerIndices[faceVertexIndices[1]]]);
}
//--------------------------------------------------------------------------------------------------

2
ApplicationCode/ReservoirDataModel/RigCell.h
View File

@ -63,7 +63,7 @@ public:
cvf::Vec3d center() const;
cvf::Vec3d faceCenter(cvf::StructGridInterface::FaceType face) const;
cvf::Vec3d faceNormal(cvf::StructGridInterface::FaceType face) const;
cvf::Vec3d faceNormalWithAreaLenght(cvf::StructGridInterface::FaceType face) const;
int firstIntersectionPoint(const cvf::Ray& ray, cvf::Vec3d* intersectionPoint) const;
bool isLongPyramidCell(double maxHeightFactor = 5, double nodeNearTolerance = 1e-3 ) const;

144
ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.cpp Normal file
View File

@ -0,0 +1,144 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) Statoil ASA, 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 "RigCombRiTransResultAccessor.h"
#include "RigGridBase.h"
#include <cmath>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigCombRiTransResultAccessor::RigCombRiTransResultAccessor(const RigGridBase* grid)
: m_grid(grid)
{
m_cdarchy = 0.008527; // (ECLIPSE 100) (METRIC)
}
//--------------------------------------------------------------------------------------------------
/// Only sensible to provide the positive values, as the negative ones will never be used.
/// The negative faces gets their value from the neighbor cell in that direction
//--------------------------------------------------------------------------------------------------
void RigCombRiTransResultAccessor::setPermResultAccessors( RigResultAccessor* xPermAccessor,
RigResultAccessor* yPermAccessor,
RigResultAccessor* zPermAccessor)
{
m_xPermAccessor = xPermAccessor;
m_yPermAccessor = yPermAccessor;
m_zPermAccessor = zPermAccessor;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCombRiTransResultAccessor::setNTGResultAccessor(RigResultAccessor* ntgAccessor)
{
m_ntgAccessor = ntgAccessor;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigCombRiTransResultAccessor::cellScalar(size_t gridLocalCellIndex) const
{
CVF_TIGHT_ASSERT(false);
return HUGE_VAL;
}
double RigCombRiTransResultAccessor::getPermValue(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId)
{
switch (faceId)
{
case cvf::StructGridInterface::POS_I:
case cvf::StructGridInterface::NEG_I:
{
return m_xPermAccessor->cellScalar(gridLocalCellIndex);
}
break;
case cvf::StructGridInterface::POS_J:
case cvf::StructGridInterface::NEG_J:
{
return m_yPermAccessor->cellScalar(gridLocalCellIndex);
}
break;
case cvf::StructGridInterface::POS_K:
case cvf::StructGridInterface::NEG_K:
{
return zPermAccessor->cellScalar(gridLocalCellIndex);
}
break;
}
return HUGE_VAL;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigCombRiTransResultAccessor::calculateHalfCellTrans( size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId)
{
double perm = getPermValue(gridLocalCellIndex, faceId);
double ntg = 1.0;
if (faceId != cvf::StructGridInterface::POS_K && faceId != cvf::StructGridInterface::NEG_K)
{
m_ntgAccessor->cellScalar(gridLocalCellIndex);
}
RigCell& cell = m_grid->cell(gridLocalCellIndex);
cvf::Vec3d centerToFace = cell.faceCenter(faceId) - cell.center();
cvf::Vec3d faceAreaVec = cell.faceNormalWithAreaLenght(faceId);
double halfCellTrans = perm*ntg*(faceAreaVec*centerToFace) / (centerToFace*centerToFace);
return halfCellTrans;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigCombRiTransResultAccessor::cellFaceScalar(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const
{
size_t reservoirCellIndex = m_grid->reservoirCellIndex(gridLocalCellIndex);
RigFault* fault = m_grid->mainGrid()->findFaultFromCellIndexAndCellFace(reservoirCellIndex, faceId);
if (fault)
{
}
else
{
size_t i, j, k, neighborGridCellIdx;
m_grid->ijkFromCellIndex(gridLocalCellIndex, &i, &j, &k);
if (m_grid->cellIJKNeighbor(i, j, k, faceId, &neighborGridCellIdx))
{
double halfCellTrans = calculateHalfCellTrans(gridLocalCellIndex, faceId);
double neighborHalfCellTrans = calculateHalfCellTrans(neighborGridCellIdx, cvf::StructGridInterface::oppositeFace(faceId));
return m_cdarchy/ ( ( 1 / halfCellTrans) + (1 / neighborHalfCellTrans) );
}
}
return HUGE_VAL;
}

62
ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.h Normal file
View File

@ -0,0 +1,62 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) Statoil ASA, 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.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include "RigResultAccessor.h"
#include "cvfCollection.h"
class RigGridBase;
//==================================================================================================
///
//==================================================================================================
class RigCombRiTransResultAccessor : public RigResultAccessor
{
public:
RigCombRiTransResultAccessor(const RigGridBase* grid);
void setPermResultAccessors(RigResultAccessor* xPermAccessor,
RigResultAccessor* yPermAccessor,
RigResultAccessor* zPermAccessor);
void setNTGResultAccessor( RigResultAccessor* ntgAccessor);
/// CDARCY Darcy’s constant
/// = 0.00852702 (E300); 0.008527 (ECLIPSE 100) (METRIC)
/// = 0.00112712 (E300); 0.001127 (ECLIPSE 100) (FIELD)
/// = 3.6 (LAB)
/// = 0.00864 (PVT-M)
void setCDARCHY(double cDarchy) { m_cdarchy = cDarchy;}
virtual double cellScalar(size_t gridLocalCellIndex) const;
virtual double cellFaceScalar(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const;
private:
cvf::ref<RigResultAccessor> m_xPermAccessor;
cvf::ref<RigResultAccessor> m_yPermAccessor;
cvf::ref<RigResultAccessor> m_zPermAccessor;
cvf::ref<RigResultAccessor> m_ntgAccessor;
double m_cdarchy;
const RigGridBase* m_grid;
};

5
ApplicationCode/ReservoirDataModel/RigFault.h
View File

@ -46,7 +46,7 @@ public:
m_faultIdxForCellFace.resize(reservoirCellCount, initVal);
}
inline int faultIdx(size_t reservoirCellIndex, cvf::StructGridInterface::FaceType face)
inline int faultIdx(size_t reservoirCellIndex, cvf::StructGridInterface::FaceType face) const
{
return m_faultIdxForCellFace[reservoirCellIndex][face];
}
@ -94,9 +94,6 @@ public:
std::vector<size_t>& connectionIndices() { return m_connectionIndices; }
const std::vector<size_t>& connectionIndices() const { return m_connectionIndices; }
static RigFaultsPrCellAccumulator* faultsPrCellAccumulator() { CVF_ASSERT(m_faultsPrCellAcc.notNull()); return m_faultsPrCellAcc.p();}
static void initFaultsPrCellAccumulator(size_t reservoirCellCount) { m_faultsPrCellAcc = new RigFaultsPrCellAccumulator(reservoirCellCount); }
private:
QString m_name;

15
ApplicationCode/ReservoirDataModel/RigGridBase.cpp
View File

@ -490,18 +490,3 @@ bool RigGridCellFaceVisibilityFilter::isFaceVisible(size_t i, size_t j, size_t k
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFaultFaceVisibilityFilter::isFaceVisible(size_t i, size_t j, size_t k, cvf::StructGridInterface::FaceType face, const cvf::UByteArray* cellVisibility) const
{
size_t cellIndex = m_grid->cellIndexFromIJK(i, j, k);
if (m_grid->cell(cellIndex).isCellFaceFault(face))
{
return true;
}
return false;
}

14
ApplicationCode/ReservoirDataModel/RigGridBase.h
View File

@ -133,17 +133,3 @@ public:
private:
const RigGridBase* m_grid;
};
class RigFaultFaceVisibilityFilter : public cvf::CellFaceVisibilityFilter
{
public:
RigFaultFaceVisibilityFilter(const RigGridBase* grid)
: m_grid(grid)
{
}
virtual bool isFaceVisible( size_t i, size_t j, size_t k, cvf::StructGridInterface::FaceType face, const cvf::UByteArray* cellVisibility ) const;
private:
const RigGridBase* m_grid;
};

31
ApplicationCode/ReservoirDataModel/RigMainGrid.cpp
View File

@ -20,6 +20,8 @@
#include "cvfAssert.h"
#include "RimDefines.h"
#include "RigFault.h"
RigMainGrid::RigMainGrid(void)
: RigGridBase(this)
@ -213,13 +215,12 @@ void RigMainGrid::setFaults(const cvf::Collection<RigFault>& faults)
//--------------------------------------------------------------------------------------------------
void RigMainGrid::calculateFaults()
{
//RigFault::initFaultsPrCellAccumulator(m_cells.size());
cvf::ref<RigFaultsPrCellAccumulator> faultsPrCellAcc = new RigFaultsPrCellAccumulator(m_cells.size());
m_faultsPrCellAcc = new RigFaultsPrCellAccumulator(m_cells.size());
// Spread fault idx'es on the cells from the faults
for (size_t fIdx = 0 ; fIdx < m_faults.size(); ++fIdx)
{
m_faults[fIdx]->accumulateFaultsPrCell(faultsPrCellAcc.p(), static_cast<int>(fIdx));
m_faults[fIdx]->accumulateFaultsPrCell(m_faultsPrCellAcc.p(), static_cast<int>(fIdx));
}
// Find the geometrical faults that is in addition
@ -244,7 +245,7 @@ void RigMainGrid::calculateFaults()
{
cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::FaceType(faceIdx);
if (faultsPrCellAcc->faultIdx(gcIdx, face) == RigFaultsPrCellAccumulator::NO_FAULT)
if (m_faultsPrCellAcc->faultIdx(gcIdx, face) == RigFaultsPrCellAccumulator::NO_FAULT)
{
// Find neighbor cell
if (firstNO_FAULTFaceForCell) // To avoid doing this for every face, and only when detecting a NO_FAULT
@ -288,11 +289,8 @@ void RigMainGrid::calculateFaults()
// To avoid doing this calculation for the opposite face
faultsPrCellAcc->setFaultIdx(gcIdx, face, unNamedFaultIdx);
faultsPrCellAcc->setFaultIdx(neighborReservoirCellIdx, StructGridInterface::oppositeFace(face), unNamedFaultIdx);
//m_cells[gcIdx].setCellFaceFault(face);
//m_cells[neighborReservoirCellIdx].setCellFaceFault(StructGridInterface::oppositeFace(face));
m_faultsPrCellAcc->setFaultIdx(gcIdx, face, unNamedFaultIdx);
m_faultsPrCellAcc->setFaultIdx(neighborReservoirCellIdx, StructGridInterface::oppositeFace(face), unNamedFaultIdx);
// Add as fault face only if the grid index is less than the neighbors
@ -330,8 +328,8 @@ void RigMainGrid::calculateFaults()
if (conn.m_c1Face != StructGridInterface::NO_FACE)
{
fIdx1 = faultsPrCellAcc->faultIdx(conn.m_c1GlobIdx, conn.m_c1Face);
fIdx2 = faultsPrCellAcc->faultIdx(conn.m_c2GlobIdx, StructGridInterface::oppositeFace(conn.m_c1Face));
fIdx1 = m_faultsPrCellAcc->faultIdx(conn.m_c1GlobIdx, conn.m_c1Face);
fIdx2 = m_faultsPrCellAcc->faultIdx(conn.m_c2GlobIdx, StructGridInterface::oppositeFace(conn.m_c1Face));
}
if (fIdx1 < 0 && fIdx2 < 0)
@ -373,8 +371,15 @@ bool RigMainGrid::faceNormalsIsOutwards() const
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const RigFault* RigMainGrid::findFaultFromCellIndexAndCellFace(size_t cellIndex, cvf::StructGridInterface::FaceType face) const
const RigFault* RigMainGrid::findFaultFromCellIndexAndCellFace(size_t reservoirCellIndex, cvf::StructGridInterface::FaceType face) const
{
int faultIdx = m_faultsPrCellAcc->faultIdx(reservoirCellIndex, face);
if (faultIdx != RigFaultsPrCellAccumulator::NO_FAULT )
{
return m_faults.at(faultIdx);
}
#if 0
for (size_t i = 0; i < m_faults.size(); i++)
{
const RigFault* rigFault = m_faults.at(i);
@ -399,6 +404,6 @@ const RigFault* RigMainGrid::findFaultFromCellIndexAndCellFace(size_t cellIndex,
}
}
}
#endif
return NULL;
}

3
ApplicationCode/ReservoirDataModel/RigMainGrid.h
View File

@ -51,7 +51,7 @@ public:
void setFaults(const cvf::Collection<RigFault>& faults);
const cvf::Collection<RigFault>& faults() { return m_faults; }
void calculateFaults();
const RigFault* findFaultFromCellIndexAndCellFace(size_t cellIndex, cvf::StructGridInterface::FaceType face) const;
const RigFault* findFaultFromCellIndexAndCellFace(size_t reservoirCellIndex, cvf::StructGridInterface::FaceType face) const;
bool faceNormalsIsOutwards() const;
void computeCachedData();
@ -76,6 +76,7 @@ private:
cvf::Collection<RigFault> m_faults;
cvf::ref<RigNNCData> m_nncData;
cvf::ref<RigFaultsPrCellAccumulator> m_faultsPrCellAcc;
cvf::Vec3d m_displayModelOffset;

2
ApplicationCode/ReservoirDataModel/RigNNCData.cpp
View File

@ -95,7 +95,7 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid)
cvf::Vec3d fc1 = c1.faceCenter((cvf::StructGridInterface::FaceType)(fIdx));
cvf::Vec3d fc2 = c2.faceCenter(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)));
cvf::Vec3d fc1ToFc2 = fc2 - fc1;
normal = c1.faceNormal((cvf::StructGridInterface::FaceType)(fIdx));
normal = c1.faceNormalWithAreaLenght((cvf::StructGridInterface::FaceType)(fIdx));
normal.normalize();
// Check that face centers are approx in the face plane
if (normal.dot(fc1ToFc2) < 0.01*fc1ToFc2.length())

6
ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp
View File

@ -599,7 +599,7 @@ bool EarClipTesselator::calculateTriangles( std::vector<size_t>* triangleIndices
// We want m_polygonIndices to be a counter-clockwise polygon to make the validation test work
if (calculatePolygonArea() < 0 )
if (calculateProjectedPolygonArea() < 0 )
{
m_polygonIndices.reverse();
}
@ -732,7 +732,7 @@ bool EarClipTesselator::isPointInsideTriangle(const cvf::Vec3d& A, const cvf::Ve
/// Computes area of the currently stored 2D polygon/contour
//--------------------------------------------------------------------------------------------------
double EarClipTesselator::calculatePolygonArea() const
double EarClipTesselator::calculateProjectedPolygonArea() const
{
CVF_ASSERT(m_X > -1 && m_Y > -1);
@ -834,7 +834,7 @@ bool FanEarClipTesselator::calculateTriangles(std::vector<size_t>* triangles)
// We want m_polygonIndices to be a counter-clockwise polygon to make the validation test work
if (calculatePolygonArea() < 0 )
if (calculateProjectedPolygonArea() < 0 )
{
m_polygonIndices.reverse();
}

2
ApplicationCode/ReservoirDataModel/cvfGeometryTools.h
View File

@ -138,7 +138,7 @@ public:
protected:
bool isTriangleValid( 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 calculatePolygonArea() const;
double calculateProjectedPolygonArea() const;
protected:
std::list<size_t> m_polygonIndices;