mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-10 15:36:09 -06:00
209 lines
8.1 KiB
C++
209 lines
8.1 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) 2018- Equinor ASA
|
|
//
|
|
// 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 "RiaCellDividingTools.h"
|
|
|
|
#include "cvfAssert.h"
|
|
|
|
#include <cmath>
|
|
#include <limits>
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Internal functions
|
|
//--------------------------------------------------------------------------------------------------
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Splits a line in a number of equal parts
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<cvf::Vec3d> splitLine(cvf::Vec3d ptStart, cvf::Vec3d ptEnd, size_t partCount);
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Calculates all points on a face described by edge points from all four edges.
|
|
/// The result is a grid of points including the given edge points
|
|
///
|
|
/// edgeXPtsHigh
|
|
/// |-------------|
|
|
/// | |
|
|
/// edgeYPtsLow | | edgeYPtsHigh
|
|
/// | |
|
|
/// |-------------|
|
|
/// edgeXPtsLow
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<std::vector<cvf::Vec3d>> calcFacePoints(const std::vector<cvf::Vec3d> edgeXPtsLow,
|
|
const std::vector<cvf::Vec3d> edgeXPtsHigh,
|
|
const std::vector<cvf::Vec3d> edgeYPtsLow,
|
|
const std::vector<cvf::Vec3d> edgeYPtsHigh);
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<cvf::Vec3d>
|
|
RiaCellDividingTools::createHexCornerCoords(std::array<cvf::Vec3d, 8> mainCellCorners, size_t nx, size_t ny, size_t nz)
|
|
{
|
|
std::array<std::pair<size_t, size_t>, 12> edgeCorners = {
|
|
std::make_pair(0, 1),
|
|
std::make_pair(3, 2),
|
|
std::make_pair(4, 5),
|
|
std::make_pair(7, 6), // X
|
|
std::make_pair(0, 3),
|
|
std::make_pair(4, 7),
|
|
std::make_pair(1, 2),
|
|
std::make_pair(5, 6), // Y
|
|
std::make_pair(0, 4),
|
|
std::make_pair(1, 5),
|
|
std::make_pair(3, 7),
|
|
std::make_pair(2, 6), // Z
|
|
};
|
|
|
|
std::array<size_t, 3> nxyz = {nx, ny, nz};
|
|
std::array<std::vector<cvf::Vec3d>, 12> edgePoints;
|
|
|
|
for (int i = 0; i < 12; i++)
|
|
{
|
|
int partCountsIndex = i / 4;
|
|
edgePoints[i] =
|
|
splitLine(mainCellCorners[edgeCorners[i].first], mainCellCorners[edgeCorners[i].second], nxyz[partCountsIndex]);
|
|
}
|
|
|
|
// lowIJ, highIJ, lowJK, highKJ,
|
|
std::vector<std::vector<std::vector<cvf::Vec3d>>> nodes;
|
|
nodes.reserve((nx + 1) * (ny + 1) * (nz + 1));
|
|
|
|
auto xyFacePtsLow = calcFacePoints(edgePoints[0], edgePoints[1], edgePoints[4], edgePoints[6]);
|
|
auto xyFacePtsHigh = calcFacePoints(edgePoints[2], edgePoints[3], edgePoints[5], edgePoints[7]);
|
|
auto yzFacePtsLow = calcFacePoints(edgePoints[4], edgePoints[5], edgePoints[8], edgePoints[10]);
|
|
auto yzFacePtsHigh = calcFacePoints(edgePoints[6], edgePoints[7], edgePoints[9], edgePoints[11]);
|
|
auto xzFacePtsLow = calcFacePoints(edgePoints[0], edgePoints[2], edgePoints[8], edgePoints[9]);
|
|
auto xzFacePtsHigh = calcFacePoints(edgePoints[1], edgePoints[3], edgePoints[10], edgePoints[11]);
|
|
|
|
nodes.push_back(xyFacePtsLow);
|
|
|
|
for (size_t z = 1; z < nz; z++)
|
|
{
|
|
auto xyFacePoints = calcFacePoints(xzFacePtsLow[z], xzFacePtsHigh[z], yzFacePtsLow[z], yzFacePtsHigh[z]);
|
|
nodes.push_back(xyFacePoints);
|
|
}
|
|
|
|
nodes.push_back(xyFacePtsHigh);
|
|
|
|
std::vector<cvf::Vec3d> coords;
|
|
coords.reserve(nx * ny * nz * 8);
|
|
|
|
for (size_t z = 1; z < nz + 1; z++)
|
|
{
|
|
for (size_t y = 1; y < ny + 1; y++)
|
|
{
|
|
for (size_t x = 1; x < nx + 1; x++)
|
|
{
|
|
std::array<cvf::Vec3d, 8> cs;
|
|
|
|
cs[0] = nodes[z - 1][y - 1][x - 1];
|
|
cs[1] = nodes[z - 1][y - 1][x];
|
|
cs[2] = nodes[z - 1][y][x];
|
|
cs[3] = nodes[z - 1][y][x - 1];
|
|
|
|
cs[4] = nodes[z][y - 1][x - 1];
|
|
cs[5] = nodes[z][y - 1][x];
|
|
cs[6] = nodes[z][y][x];
|
|
cs[7] = nodes[z][y][x - 1];
|
|
|
|
coords.insert(coords.end(), cs.begin(), cs.end());
|
|
}
|
|
}
|
|
}
|
|
return coords;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RiaCellDividingTools::computeFlowDistance(const std::array<cvf::Vec3d, 8>& cellVertices,
|
|
const cvf::Vec3d& areaCenter)
|
|
{
|
|
auto subCellCorners = createHexCornerCoords(cellVertices, 2, 2, 2);
|
|
|
|
double weightedDistanceTotal = 0.0;
|
|
double weightTotal = 0.0;
|
|
|
|
for (size_t c = 0; c < 8; c++)
|
|
{
|
|
double weight = 1.0;
|
|
weightTotal += weight;
|
|
|
|
cvf::Vec3d centerOfSubCell = cvf::Vec3d::ZERO;
|
|
{
|
|
cvf::Vec3d vertexSum = cvf::Vec3d::ZERO;
|
|
for (size_t v = 0; v < 8; v++) vertexSum += subCellCorners[c * 8 + v];
|
|
centerOfSubCell = vertexSum / 8;
|
|
}
|
|
|
|
auto dist = (centerOfSubCell - areaCenter).length();
|
|
weightedDistanceTotal += (dist * weight);
|
|
}
|
|
|
|
return weightedDistanceTotal / weightTotal;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<cvf::Vec3d> splitLine(cvf::Vec3d ptStart, cvf::Vec3d ptEnd, size_t partCount)
|
|
{
|
|
std::vector<cvf::Vec3d> pts = {ptStart};
|
|
|
|
for (size_t i = 1; i < partCount; i++)
|
|
{
|
|
pts.push_back(cvf::Vec3d(ptStart.x() + (ptEnd.x() - ptStart.x()) * i / partCount,
|
|
ptStart.y() + (ptEnd.y() - ptStart.y()) * i / partCount,
|
|
ptStart.z() + (ptEnd.z() - ptStart.z()) * i / partCount));
|
|
}
|
|
pts.push_back(ptEnd);
|
|
return pts;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<std::vector<cvf::Vec3d>> calcFacePoints(const std::vector<cvf::Vec3d> edgeXPtsLow,
|
|
const std::vector<cvf::Vec3d> edgeXPtsHigh,
|
|
const std::vector<cvf::Vec3d> edgeYPtsLow,
|
|
const std::vector<cvf::Vec3d> edgeYPtsHigh)
|
|
{
|
|
CVF_ASSERT(edgeXPtsLow.size() == edgeXPtsHigh.size() && edgeYPtsLow.size() == edgeYPtsHigh.size());
|
|
|
|
size_t xSize = edgeXPtsLow.size();
|
|
size_t ySize = edgeYPtsLow.size();
|
|
|
|
std::vector<std::vector<cvf::Vec3d>> pts;
|
|
|
|
// Add low edge points
|
|
pts.push_back(edgeXPtsLow);
|
|
|
|
// Interior points
|
|
for (size_t y = 1; y < ySize - 1; y++)
|
|
{
|
|
auto interiorPts = splitLine(edgeYPtsLow[y], edgeYPtsHigh[y], xSize - 1);
|
|
pts.push_back(interiorPts);
|
|
}
|
|
|
|
// Add low edge points
|
|
pts.push_back(edgeXPtsHigh);
|
|
return pts;
|
|
}
|