Merge pull request #4393 from arturcastiel/lgr-refactor

LGR-EGrid: Rebased Version of the PR#4354

Supersedes and closes #4354

CMakeLists_files.cmake

@@ -44,6 +44,8 @@ list (APPEND MAIN_SOURCE_FILES
       opm/common/utility/parameters/ParameterRequirement.cpp
       opm/common/utility/parameters/ParameterTools.cpp
       opm/common/utility/numeric/calculateCellVol.cpp
+      opm/common/utility/numeric/GeometryUtil.cpp
+      opm/common/utility/numeric/VectorUtil.cpp
       opm/common/utility/numeric/MonotCubicInterpolator.cpp
       opm/common/utility/numeric/RootFinders.cpp
       opm/material/common/Spline.cpp
@@ -736,7 +738,11 @@ if(ENABLE_ECL_OUTPUT)
           tests/msim/action_count_no_run_function.py
           tests/msim/open_well_past.py
           tests/msim/open_well_too_late.py
-          tests/VFP_CASE.DATA)
+          tests/VFP_CASE.DATA
+          tests/CARFIN-COLUMN.EGRID
+          tests/CARFIN-DOUBLE.EGRID
+          tests/CARFIN-NESTED.EGRID
+          tests/CARFIN5.EGRID)
 endif()
 
 list (APPEND EXAMPLE_SOURCE_FILES
@@ -825,6 +831,8 @@ list( APPEND PUBLIC_HEADER_FILES
       opm/common/utility/numeric/cmp.hpp
       opm/common/utility/numeric/blas_lapack.h
       opm/common/utility/numeric/calculateCellVol.hpp
+      opm/common/utility/numeric/GeometryUtil.hpp
+      opm/common/utility/numeric/VectorUtil.hpp
       opm/common/utility/numeric/buildUniformMonotoneTable.hpp
       opm/common/utility/numeric/linearInterpolation.hpp
       opm/common/utility/numeric/MonotCubicInterpolator.hpp

opm/common/utility/numeric/GeometryUtil.cpp

@@ -0,0 +1,5 @@
+#include <opm/common/utility/numeric/GeometryUtil.hpp>
+
+namespace GeometryUtil {
+
+} // namespace GeometryUtils

opm/common/utility/numeric/GeometryUtil.hpp

@@ -0,0 +1,116 @@
+#ifndef GEOMETRYUTIL_H
+#define GEOMETRYUTIL_H
+
+#include <vector>
+#include <opm/common/utility/numeric/VectorUtil.hpp>
+#include <numeric>
+#include <cmath> 
+namespace GeometryUtil {
+
+constexpr double epslon = 1e-6;         
+// A simple utility function to calculate area of a rectangle
+double calculateRectangleArea(double width, double height);
+
+template <typename T = double>
+T  calcTetraVol(const std::array<T,4>& x, const std::array<T,4>& y, const std::array<T,4>& z ){
+    T det =	x[0]*y[2]*z[1] - x[0]*y[1]*z[2] + x[1]*y[0]*z[2]
+            - x[1]*y[2]*z[0] - x[2]*y[0]*z[1] + x[2]*y[1]*z[0]
+            + x[0]*y[1]*z[3] - x[0]*y[3]*z[1] - x[1]*y[0]*z[3]
+            + x[1]*y[3]*z[0] + x[3]*y[0]*z[1] - x[3]*y[1]*z[0]
+            - x[0]*y[2]*z[3] + x[0]*y[3]*z[2] + x[2]*y[0]*z[3]
+            - x[2]*y[3]*z[0] - x[3]*y[0]*z[2] + x[3]*y[2]*z[0]
+            + x[1]*y[2]*z[3] - x[1]*y[3]*z[2] - x[2]*y[1]*z[3]
+            + x[2]*y[3]*z[1] + x[3]*y[1]*z[2] - x[3]*y[2]*z[1];
+            return std::abs(det)/6; 
+}
+
+template <typename T = double>
+T calcHexaVol(const std::array<T,8>& x, const std::array<T,8>& y, const std::array<T,8>& z, 
+              const T& cx,  const T& cy, const T& cz )
+{
+    constexpr std::array<std::array<int, 3>, 12> faceConfigurations
+    {
+        std::array<int, 3>{0, 1, 5}, 
+                          {1, 5, 4},     // Face 0
+                          {0, 4, 6},
+                          {4, 6, 2},     // Face 1
+                          {2, 3, 7}, 
+                          {3, 7, 6},     // Face 2
+                          {1, 3, 7}, 
+                          {3, 7, 5},     // Face 3
+                          {0, 1, 3}, 
+                          {1, 3, 2},     // Face 4
+                          {4, 5, 7}, 
+                          {5, 7, 6}     // Face 5
+    };
+    auto getNodes = [](const std::array<T, 8>& X, const std::array<T, 8>& Y, const std::array<T, 8>& Z,
+                                    const std::array<int,3>&  ind){
+        std::array<T, 3> filtered_vectorX;
+        std::array<T, 3> filtered_vectorY;
+        std::array<T, 3> filtered_vectorZ;
+        for (std::size_t index = 0; index < ind.size(); index++) {
+            filtered_vectorX[index] = X[ind[index]];
+            filtered_vectorY[index] = Y[ind[index]];
+            filtered_vectorZ[index] = Z[ind[index]];
+        }
+        return std::make_tuple(filtered_vectorX,filtered_vectorY,filtered_vectorZ);
+    };
+    // note: some CPG grids may have collapsed faces that are not planar, therefore
+    // the hexadron is subdivided in terahedrons.
+    // calculating the volume of the pyramid with F0 as base and pc as center                         
+    T totalVolume = 0.0;
+    for (size_t i = 0; i < faceConfigurations.size(); i += 2) {
+        auto [fX0, fY0, fZ0] = getNodes(x, y, z, faceConfigurations[i]);
+        totalVolume += std::apply(calcTetraVol<T>, VectorUtil::appendNode<double>(fX0, fY0, fZ0, cx, cy, cz));
+
+        auto [fX1, fY1, fZ1] = getNodes(x, y, z, faceConfigurations[i + 1]);
+        totalVolume += std::apply(calcTetraVol<T>, VectorUtil::appendNode<double>(fX1, fY1, fZ1, cx, cy, cz));
+    }
+    return totalVolume;
+};
+
+template <typename T = double>
+std::vector<int> isInsideElement(const std::vector<T>& tpX, const std::vector<T>& tpY, const std::vector<T>& tpZ,  
+                                                const std::vector<std::array<T, 8>>& X, const std::vector<std::array<T, 8>>& Y,
+                                                const std::vector<std::array<T, 8>>& Z)
+{
+    std::vector<int> in_elements(tpX.size(),0);
+    // check if it is insde or outside boundary box
+    T minX, minY, minZ, maxX, maxY, maxZ;
+    T pcX, pcY, pcZ, element_volume, test_element_volume;
+    bool flag;
+    for (std::size_t outerIndex = 0; outerIndex < X.size(); outerIndex++) {
+        minX = *std::min_element(X[outerIndex].begin(), X[outerIndex].end());
+        minY = *std::min_element(Y[outerIndex].begin(), Y[outerIndex].end());
+        minZ = *std::min_element(Z[outerIndex].begin(), Z[outerIndex].end());
+        maxX = *std::max_element(X[outerIndex].begin(), X[outerIndex].end());
+        maxY = *std::max_element(Y[outerIndex].begin(), Y[outerIndex].end());
+        maxZ = *std::max_element(Z[outerIndex].begin(), Z[outerIndex].end());
+        pcX = std::accumulate(X[outerIndex].begin(), X[outerIndex].end(), 0.0)/8;
+        pcY = std::accumulate(Y[outerIndex].begin(), Y[outerIndex].end(), 0.0)/8;            
+        pcZ = std::accumulate(Z[outerIndex].begin(), Z[outerIndex].end(), 0.0)/8;
+        element_volume = calcHexaVol(X[outerIndex],Y[outerIndex],Z[outerIndex], pcX, pcY,pcZ);                         
+        for (size_t innerIndex  = 0; innerIndex < tpX.size(); innerIndex++){
+            // check if center of refined volume is outside the boundary box of a coarse volume.
+            // Only computes volumed base test is this condition is met.
+            flag  = (minX < tpX[innerIndex]) && (maxX > tpX[innerIndex]) &&
+                    (minY < tpY[innerIndex]) && (maxY > tpY[innerIndex]) &&
+                    (minZ < tpZ[innerIndex]) && (maxZ > tpZ[innerIndex]);                            
+            if (flag && (in_elements[innerIndex] == 0)) {
+                test_element_volume = calcHexaVol(X[outerIndex],Y[outerIndex],Z[outerIndex], 
+                                                tpX[innerIndex], tpY[innerIndex],tpZ[innerIndex]);                         
+                if (std::abs(test_element_volume - element_volume) < epslon){
+                        in_elements[innerIndex] = static_cast<int>(outerIndex);
+                }                             
+            }
+        }
+    }
+    return in_elements;
+};
+
+
+// Example for other utilities...
+
+} // namespace GeometryUtils
+
+#endif // GEOMETRY_UTILS_H

opm/common/utility/numeric/VectorUtil.cpp

@@ -0,0 +1,6 @@
+#include <opm/common/utility/numeric/VectorUtil.hpp>
+
+namespace VectorUtil {
+
+
+} // namespace VectorUtil

opm/common/utility/numeric/VectorUtil.hpp

@@ -0,0 +1,114 @@
+#ifndef VECTORUTIL_H
+#define VECTORUTIL_H
+
+#include <tuple>
+#include <vector>
+#include <functional>
+#include <stdexcept> 
+#include <array>
+#include <tuple>
+
+
+namespace VectorUtil {
+
+
+template <typename T = double>
+std::tuple<std::array<T,4>, std::array<T,4>, std::array<T,4>> 
+appendNode(const std::array<T,3>& X, const std::array<T,3>& Y, const std::array<T,3>& Z, 
+           const T& xc, const T& yc,  const T& zc) 
+{
+    std::array<T,4> tX;
+    std::array<T,4> tY;
+    std::array<T,4> tZ;
+    std::copy(X.begin(), X.end(), tX.begin());
+    tX[3]= xc;
+    std::copy(Y.begin(), Y.end(), tY.begin());
+    tY[3]= yc;            
+    std::copy(Z.begin(), Z.end(), tZ.begin());
+    tZ[3]= zc; 
+    return std::make_tuple(tX,tY,tZ);
+}
+
+// Implementation of generation General Operation between two vectors of the same type
+template <typename T, typename Operation>
+std::vector<T> vectorOperation(const std::vector<T>& vecA, const std::vector<T>& vecB, Operation op) {
+    if (vecA.size() != vecB.size()) {
+        throw std::invalid_argument("Error: Vectors must have the same size!"); // Throwing exception
+    }
+    std::vector<T> result;
+    result.reserve(vecA.size());
+    // Use std::transform with the passed operation
+    std::transform(vecA.begin(), vecA.end(), vecB.begin(), std::back_inserter(result), op);
+    return result;
+}
+
+
+// A simple utility function to calculate area of a rectangle
+std::tuple<std::array<double,4>, std::array<double,4>, std::array<double,4>> 
+appendNode(const std::array<double,3>&, const std::array<double,3>&, 
+           const std::array<double,3>&, const double&, const double&, 
+           const double&);
+template <typename T>
+std::vector<T> filterArray(const std::vector<std::size_t>& X, const std::vector<int>& ind){
+    std::vector<T> filtered_vectorX(ind.size(),0);
+    for (std::size_t index = 0; index < ind.size(); index++) {
+        filtered_vectorX[index] = X[ind[index]];
+    }
+    return filtered_vectorX;
+};       
+
+// T type of the object
+// Rout type of the object output 
+// Method type of method
+template <typename T,  typename Rout, typename Rin = std::size_t, typename Method, typename... Args>
+std::vector<Rout> callMethodForEachInputOnObject(const T& obj, Method mtd, const std::vector<Rin>& input_vector, Args&&... args) {
+    std::vector<Rout> result;
+    // Reserve space for each vector in the tuple
+    result.reserve(input_vector.size());
+    // Iterate over the input_vector and fill the tuple's vectors
+    for (const auto& element : input_vector) {
+        Rout value = (obj.*mtd)(element, std::forward<Args>(args)...);
+        result.push_back(std::move(value));
+    }
+    return result;
+}
+
+
+
+template <typename T>
+std::tuple<std::vector<T>,std::vector<T>, std::vector<T>> splitXYZ(std::vector<std::array<T,3>>& input_vector ){
+std::vector<T> X, Y, Z;
+    X.reserve(input_vector.size());
+    Y.reserve(input_vector.size());
+    Z.reserve(input_vector.size());
+    for (auto& element : input_vector) {
+            X.push_back(std::move(element[0]));  // Add to first vector
+            Y.push_back(std::move(element[1]));    // Add to second vector
+            Z.push_back(std::move(element[2]));   // Add to third vector
+    }
+    return std::make_tuple(X,Y,Z);
+}
+
+
+template <typename T,  typename Rout, typename Rin = std::size_t, typename Method, typename... Args>
+auto callMethodForEachInputOnObjectXYZ(const T& obj, Method mtd, const std::vector<Rin>& input_vector, Args&&... args) {
+    using X = typename Rout::value_type;
+    auto result = callMethodForEachInputOnObject<T, Rout, Rin, Method, Args...>(obj, mtd, input_vector, std::forward<Args>(args)...);
+    return splitXYZ<X>(result);
+}
+
+template <typename T>
+    void test(std::vector<T> vec){
+    std::size_t index = 1;
+    for (T el : vec){
+        index++;
+    }
+}
+
+
+
+// Example for other utilities...
+
+} // namespace VectorUtil
+
+#endif // VECTORUTIL_H