My Project
Loading...
Searching...
No Matches
GeometryUtil.hpp
1#ifndef GEOMETRYUTIL_H
2#define GEOMETRYUTIL_H
3
4#include <vector>
5#include <opm/common/utility/numeric/VectorUtil.hpp>
6#include <numeric>
7#include <cmath>
8namespace GeometryUtil {
9
10constexpr double epslon = 1e-6;
11// A simple utility function to calculate area of a rectangle
12double calculateRectangleArea(double width, double height);
13
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;
25}
26
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 )
30{
31 constexpr std::array<std::array<int, 3>, 12> faceConfigurations
32 {
33 std::array<int, 3>{0, 1, 5},
34 {1, 5, 4}, // Face 0
35 {0, 4, 6},
36 {4, 6, 2}, // Face 1
37 {2, 3, 7},
38 {3, 7, 6}, // Face 2
39 {1, 3, 7},
40 {3, 7, 5}, // Face 3
41 {0, 1, 3},
42 {1, 3, 2}, // Face 4
43 {4, 5, 7},
44 {5, 7, 6} // Face 5
45 };
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]];
55 }
56 return std::make_tuple(filtered_vectorX,filtered_vectorY,filtered_vectorZ);
57 };
58 // note: some CPG grids may have collapsed faces that are not planar, therefore
59 // the hexadron is subdivided in terahedrons.
60 // calculating the volume of the pyramid with F0 as base and pc as center
61 T totalVolume = 0.0;
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));
65
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));
68 }
69 return totalVolume;
70}
71
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)
76{
77 std::vector<int> in_elements(tpX.size(),0);
78 // check if it is insde or outside boundary box
79 T minX, minY, minZ, maxX, maxY, maxZ;
80 T pcX, pcY, pcZ, element_volume, test_element_volume;
81 bool flag;
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++){
94 // check if center of refined volume is outside the boundary box of a coarse volume.
95 // Only computes volumed base test is this condition is met.
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);
104 }
105 }
106 }
107 }
108 return in_elements;
109}
110
111
112// Example for other utilities...
113
114} // namespace GeometryUtils
115
116#endif // GEOMETRY_UTILS_H