5#include <opm/common/utility/numeric/VectorUtil.hpp>
8namespace GeometryUtil {
10constexpr double epslon = 1e-6;
12double calculateRectangleArea(
double width,
double height);
14template <
typename T =
double>
15T calcTetraVol(
const std::array<T,4>& x,
const std::array<T,4>& y,
const std::array<T,4>& z ){
16 T det = x[0]*y[2]*z[1] - x[0]*y[1]*z[2] + x[1]*y[0]*z[2]
17 - x[1]*y[2]*z[0] - x[2]*y[0]*z[1] + x[2]*y[1]*z[0]
18 + x[0]*y[1]*z[3] - x[0]*y[3]*z[1] - x[1]*y[0]*z[3]
19 + x[1]*y[3]*z[0] + x[3]*y[0]*z[1] - x[3]*y[1]*z[0]
20 - x[0]*y[2]*z[3] + x[0]*y[3]*z[2] + x[2]*y[0]*z[3]
21 - x[2]*y[3]*z[0] - x[3]*y[0]*z[2] + x[3]*y[2]*z[0]
22 + x[1]*y[2]*z[3] - x[1]*y[3]*z[2] - x[2]*y[1]*z[3]
23 + x[2]*y[3]*z[1] + x[3]*y[1]*z[2] - x[3]*y[2]*z[1];
24 return std::abs(det)/6;
27template <
typename T =
double>
28T calcHexaVol(
const std::array<T,8>& x,
const std::array<T,8>& y,
const std::array<T,8>& z,
29 const T& cx,
const T& cy,
const T& cz )
31 constexpr std::array<std::array<int, 3>, 12> faceConfigurations
33 std::array<int, 3>{0, 1, 5},
46 auto getNodes = [](
const std::array<T, 8>& X,
const std::array<T, 8>& Y,
const std::array<T, 8>& Z,
47 const std::array<int,3>& ind){
48 std::array<T, 3> filtered_vectorX;
49 std::array<T, 3> filtered_vectorY;
50 std::array<T, 3> filtered_vectorZ;
51 for (std::size_t index = 0; index < ind.size(); index++) {
52 filtered_vectorX[index] = X[ind[index]];
53 filtered_vectorY[index] = Y[ind[index]];
54 filtered_vectorZ[index] = Z[ind[index]];
56 return std::make_tuple(filtered_vectorX,filtered_vectorY,filtered_vectorZ);
62 for (
size_t i = 0; i < faceConfigurations.size(); i += 2) {
63 auto [fX0, fY0, fZ0] = getNodes(x, y, z, faceConfigurations[i]);
64 totalVolume += std::apply(calcTetraVol<T>, VectorUtil::appendNode<double>(fX0, fY0, fZ0, cx, cy, cz));
66 auto [fX1, fY1, fZ1] = getNodes(x, y, z, faceConfigurations[i + 1]);
67 totalVolume += std::apply(calcTetraVol<T>, VectorUtil::appendNode<double>(fX1, fY1, fZ1, cx, cy, cz));
72template <
typename T =
double>
73std::vector<int> isInsideElement(
const std::vector<T>& tpX,
const std::vector<T>& tpY,
const std::vector<T>& tpZ,
74 const std::vector<std::array<T, 8>>& X,
const std::vector<std::array<T, 8>>& Y,
75 const std::vector<std::array<T, 8>>& Z)
77 std::vector<int> in_elements(tpX.size(),0);
79 T minX, minY, minZ, maxX, maxY, maxZ;
80 T pcX, pcY, pcZ, element_volume, test_element_volume;
82 for (std::size_t outerIndex = 0; outerIndex < X.size(); outerIndex++) {
83 minX = *std::min_element(X[outerIndex].begin(), X[outerIndex].end());
84 minY = *std::min_element(Y[outerIndex].begin(), Y[outerIndex].end());
85 minZ = *std::min_element(Z[outerIndex].begin(), Z[outerIndex].end());
86 maxX = *std::max_element(X[outerIndex].begin(), X[outerIndex].end());
87 maxY = *std::max_element(Y[outerIndex].begin(), Y[outerIndex].end());
88 maxZ = *std::max_element(Z[outerIndex].begin(), Z[outerIndex].end());
89 pcX = std::accumulate(X[outerIndex].begin(), X[outerIndex].end(), 0.0)/8;
90 pcY = std::accumulate(Y[outerIndex].begin(), Y[outerIndex].end(), 0.0)/8;
91 pcZ = std::accumulate(Z[outerIndex].begin(), Z[outerIndex].end(), 0.0)/8;
92 element_volume = calcHexaVol(X[outerIndex],Y[outerIndex],Z[outerIndex], pcX, pcY,pcZ);
93 for (
size_t innerIndex = 0; innerIndex < tpX.size(); innerIndex++){
96 flag = (minX < tpX[innerIndex]) && (maxX > tpX[innerIndex]) &&
97 (minY < tpY[innerIndex]) && (maxY > tpY[innerIndex]) &&
98 (minZ < tpZ[innerIndex]) && (maxZ > tpZ[innerIndex]);
99 if (flag && (in_elements[innerIndex] == 0)) {
100 test_element_volume = calcHexaVol(X[outerIndex],Y[outerIndex],Z[outerIndex],
101 tpX[innerIndex], tpY[innerIndex],tpZ[innerIndex]);
102 if (std::abs(test_element_volume - element_volume) < epslon){
103 in_elements[innerIndex] =
static_cast<int>(outerIndex);