From 8240cca69e85f2b46bfb57ea2012960665f1e2dd Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 16 Jan 2018 16:30:57 +0100 Subject: [PATCH] remove some grid related utilities moved to opm-grid --- CMakeLists_files.cmake | 19 - opm/core/pressure/tpfa/TransTpfa.hpp | 105 ---- opm/core/pressure/tpfa/TransTpfa_impl.hpp | 177 ------ opm/core/pressure/tpfa/trans_tpfa.c | 130 ----- opm/core/pressure/tpfa/trans_tpfa.h | 110 ---- opm/core/utility/CompressedPropertyAccess.hpp | 504 ------------------ opm/core/utility/RegionMapping.hpp | 181 ------- opm/core/utility/VelocityInterpolation.cpp | 185 ------- opm/core/utility/VelocityInterpolation.hpp | 122 ----- opm/core/utility/WachspressCoord.cpp | 238 --------- opm/core/utility/WachspressCoord.hpp | 85 --- opm/core/utility/compressedToCartesian.cpp | 48 -- opm/core/utility/compressedToCartesian.hpp | 40 -- opm/core/utility/extractPvtTableIndex.cpp | 47 -- opm/core/utility/extractPvtTableIndex.hpp | 33 -- tests/compressed_gridproperty.data | 48 -- tests/test_compressedpropertyaccess.cpp | 347 ------------ tests/test_regionmapping.cpp | 155 ------ tests/test_velocityinterpolation.cpp | 487 ----------------- tests/test_wachspresscoord.cpp | 363 ------------- 20 files changed, 3424 deletions(-) delete mode 100644 opm/core/pressure/tpfa/TransTpfa.hpp delete mode 100644 opm/core/pressure/tpfa/TransTpfa_impl.hpp delete mode 100644 opm/core/pressure/tpfa/trans_tpfa.c delete mode 100644 opm/core/pressure/tpfa/trans_tpfa.h delete mode 100644 opm/core/utility/CompressedPropertyAccess.hpp delete mode 100644 opm/core/utility/RegionMapping.hpp delete mode 100644 opm/core/utility/VelocityInterpolation.cpp delete mode 100644 opm/core/utility/VelocityInterpolation.hpp delete mode 100644 opm/core/utility/WachspressCoord.cpp delete mode 100644 opm/core/utility/WachspressCoord.hpp delete mode 100644 opm/core/utility/compressedToCartesian.cpp delete mode 100644 opm/core/utility/compressedToCartesian.hpp delete mode 100644 opm/core/utility/extractPvtTableIndex.cpp delete mode 100644 opm/core/utility/extractPvtTableIndex.hpp delete mode 100644 tests/compressed_gridproperty.data delete mode 100644 tests/test_compressedpropertyaccess.cpp delete mode 100644 tests/test_regionmapping.cpp delete mode 100644 tests/test_velocityinterpolation.cpp delete mode 100644 tests/test_wachspresscoord.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 87e1ace2..23805bc2 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -51,7 +51,6 @@ list (APPEND MAIN_SOURCE_FILES opm/core/pressure/msmfem/partition.c opm/core/pressure/tpfa/cfs_tpfa_residual.c opm/core/pressure/tpfa/ifs_tpfa.c - opm/core/pressure/tpfa/trans_tpfa.c opm/core/props/BlackoilPropertiesBasic.cpp opm/core/props/BlackoilPropertiesFromDeck.cpp opm/core/props/IncompPropertiesBasic.cpp @@ -74,10 +73,6 @@ list (APPEND MAIN_SOURCE_FILES opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp opm/core/transport/reorder/reordersequence.cpp opm/core/transport/reorder/tarjan.c - opm/core/utility/VelocityInterpolation.cpp - opm/core/utility/WachspressCoord.cpp - opm/core/utility/compressedToCartesian.cpp - opm/core/utility/extractPvtTableIndex.cpp opm/core/utility/miscUtilities.cpp opm/core/utility/miscUtilitiesBlackoil.cpp opm/core/wells/InjectionSpecification.cpp @@ -92,19 +87,15 @@ list (APPEND MAIN_SOURCE_FILES # originally generated with the command: # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort list (APPEND TEST_SOURCE_FILES - tests/test_compressedpropertyaccess.cpp tests/test_dgbasis.cpp tests/test_flowdiagnostics.cpp tests/test_parallelistlinformation.cpp - tests/test_velocityinterpolation.cpp tests/test_wells.cpp - tests/test_wachspresscoord.cpp tests/test_linearsolver.cpp tests/test_parallel_linearsolver.cpp tests/test_satfunc.cpp tests/test_shadow.cpp tests/test_equil.cpp - tests/test_regionmapping.cpp tests/test_blackoilstate.cpp tests/test_wellsmanager.cpp tests/test_wellcontrols.cpp @@ -124,7 +115,6 @@ list (APPEND TEST_DATA_FILES tests/capillary.DATA tests/capillary_overlap.DATA tests/capillarySwatinit.DATA - tests/compressed_gridproperty.data tests/deadfluids.DATA tests/equil_livegas.DATA tests/equil_liveoil.DATA @@ -198,13 +188,10 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/pressure/mimetic/mimetic.h opm/core/pressure/msmfem/dfs.h opm/core/pressure/msmfem/partition.h - opm/core/pressure/tpfa/TransTpfa.hpp - opm/core/pressure/tpfa/TransTpfa_impl.hpp opm/core/pressure/tpfa/cfs_tpfa_residual.h opm/core/pressure/tpfa/compr_quant_general.h opm/core/pressure/tpfa/compr_source.h opm/core/pressure/tpfa/ifs_tpfa.h - opm/core/pressure/tpfa/trans_tpfa.h opm/core/props/BlackoilPhases.hpp opm/core/props/BlackoilPropertiesBasic.hpp opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -247,13 +234,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp opm/core/transport/reorder/reordersequence.h opm/core/transport/reorder/tarjan.h - opm/core/utility/CompressedPropertyAccess.hpp opm/core/utility/initHydroCarbonState.hpp - opm/core/utility/RegionMapping.hpp - opm/core/utility/VelocityInterpolation.hpp - opm/core/utility/WachspressCoord.hpp - opm/core/utility/compressedToCartesian.hpp - opm/core/utility/extractPvtTableIndex.hpp opm/core/utility/miscUtilities.hpp opm/core/utility/miscUtilitiesBlackoil.hpp opm/core/utility/miscUtilities_impl.hpp diff --git a/opm/core/pressure/tpfa/TransTpfa.hpp b/opm/core/pressure/tpfa/TransTpfa.hpp deleted file mode 100644 index 9169130f..00000000 --- a/opm/core/pressure/tpfa/TransTpfa.hpp +++ /dev/null @@ -1,105 +0,0 @@ -/* - Copyright 2010 SINTEF ICT, Applied Mathematics. - Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_TRANSTPFA_HEADER_INCLUDED -#define OPM_TRANSTPFA_HEADER_INCLUDED - -/** - * \file - * Routines to assist in the calculation of two-point transmissibilities. - */ - -/** - * Calculate static, one-sided transmissibilities for use in the two-point flux - * approximation method. - * - * The one-sided transmissibilities are defined by the formula - * \f[ - * t_i = \frac{\vec{n}_f \mathsf{K}_c \vec{c}_c}{\lVert \vec{c}_c \rVert^2} - * \f] - * in which @c i is the half-face index corresponding to the cell-face index - * pair (c,f) and \f$\vec{c}_{cf} = \Bar{x}_f - \Bar{x}_c\f$ is the - * centroid difference vector. - * - * @param[in] G Grid. - * @param[in] perm Permeability. One symmetric, positive definite tensor - * per grid cell. - * @param[out] htrans One-sided transmissibilities. Array of size at least - * G->cell_facepos[ G->number_of_cells ]. - */ -template -void -tpfa_htrans_compute(const Grid *G , - const double *perm , - double *htrans); - -/** - * Compute two-point transmissibilities from one-sided transmissibilities. - * - * The two-point transmissibilities are given by the simple, near-harmonic - * average (save a factor of two) - * \f[ - * \mathsf{T}_f = \big(\frac{1}{t_1} + \frac{1}{t_2}\big)^{-1} - * = \frac{t_1t_2}{t_1 + t_2} - * \f] - * in which \f$t_1\f$ and \f$t_2\f$ are the one-sided transmissibilities that - * connect the neighbouring cells of face @c f. - * - * @param[in] G Grid. - * @param[in] htrans One-sided transmissibilities as defined by function - * tpfa_htrans_compute(). - * @param[out] trans Interface, two-point transmissibilities. Array of size - * at least numFaces(G). - */ -template -void -tpfa_trans_compute(const Grid *G , - const double *htrans, - double *trans ); - -/** - * Calculate effective two-point transmissibilities from one-sided, total - * mobility weighted, transmissibilities. - * - * Specifically, compute the following product - * \f[ - * \mathsf{T}_f = \big(\frac{1}{\lambda_1t_1} + \frac{1}{\lambda_2t_2}\big)^{-1} - * = \lambda_1\lambda_2 \frac{t_1t_2}{t_1 + t_2} - * \f] - * in which \f$t_1\f$ and \f$t_2\f$ are the one-sided, static transmissibility - * values connecting the cells of face @c f and \f$\lambda_1\f$ and - * \f$\lambda_2\f$ denote the total mobilities of the respective cells. - * - * @param[in] G Grid. - * @param[in] totmob Total mobilities. One positive scalar value for each cell. - * @param[in] htrans One-sided transmissibilities as defined by function - * tpfa_htrans_compute(). - * @param[out] trans Effective, two-point transmissibilities. Array of size at - * least numFaces(G). - */ -template -void -tpfa_eff_trans_compute(const Grid *G , - const double *totmob, - const double *htrans, - double *trans ); - -#include "TransTpfa_impl.hpp" -#endif /* OPM_TRANS_TPFA_HEADER_INCLUDED */ diff --git a/opm/core/pressure/tpfa/TransTpfa_impl.hpp b/opm/core/pressure/tpfa/TransTpfa_impl.hpp deleted file mode 100644 index 3a51e3ce..00000000 --- a/opm/core/pressure/tpfa/TransTpfa_impl.hpp +++ /dev/null @@ -1,177 +0,0 @@ -#include -#include -#include - -#include - -namespace Dune -{ -class CpGrid; -} - -namespace Opm -{ -namespace UgGridHelpers -{ -int dimensions(const Dune::CpGrid&); - -double faceArea(const Dune::CpGrid&, int); -} -} - -namespace -{ -#ifdef HAVE_OPM_GRID -inline const double* multiplyFaceNormalWithArea(const Dune::CpGrid& grid, int face_index, const double* in) -{ - int d=Opm::UgGridHelpers::dimensions(grid); - double* out=new double[d]; - double area=Opm::UgGridHelpers::faceArea(grid, face_index); - - for(int i=0;i -void -tpfa_htrans_compute(const Grid* G, const double *perm, double *htrans) -/* ---------------------------------------------------------------------- */ -{ - using namespace Opm::UgGridHelpers; - int d, j; - double s, dist, denom; - - double Kn[3]; - typename CellCentroidTraits::IteratorType cc = beginCellCentroids(*G); - typename Cell2FacesTraits::Type c2f = cell2Faces(*G); - typename FaceCellTraits::Type face_cells = faceCells(*G); - - const double *n; - const double *K; - - MAT_SIZE_T nrows, ncols, ldA, incx, incy; - double a1, a2; - - d = dimensions(*G); - - nrows = ncols = ldA = d; - incx = incy = 1 ; - a1 = 1.0; a2 = 0.0 ; - - for (int c =0, i = 0; c < numCells(*G); c++) { - K = perm + (c * d * d); - - typedef typename Cell2FacesTraits::Type::row_type FaceRow; - FaceRow faces = c2f[c]; - - for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end(); - f!=end; ++f, ++i) - { - s = 2.0*(face_cells(*f, 0) == c) - 1.0; - n = faceNormal(*G, *f); - const double* nn=multiplyFaceNormalWithArea(*G, *f, n); - const double* fc = &(faceCentroid(*G, *f)[0]); - dgemv_("No Transpose", &nrows, &ncols, - &a1, K, &ldA, nn, &incx, &a2, &Kn[0], &incy); - maybeFreeFaceNormal(*G, nn); - - htrans[i] = denom = 0.0; - for (j = 0; j < d; j++) { - dist = fc[j] - getCoordinate(cc, j); - - htrans[i] += s * dist * Kn[j]; - denom += dist * dist; - } - - assert (denom > 0); - htrans[i] /= denom; - htrans[i] = std::abs(htrans[i]); - } - // Move to next cell centroid. - cc = increment(cc, 1, d); - } -} - - -/* ---------------------------------------------------------------------- */ -template -void -tpfa_trans_compute(const Grid* G, const double *htrans, double *trans) -/* ---------------------------------------------------------------------- */ -{ - using namespace Opm::UgGridHelpers; - - for (int f = 0; f < numFaces(*G); f++) { - trans[f] = 0.0; - } - - typename Cell2FacesTraits::Type c2f = cell2Faces(*G); - - for (int c = 0, i = 0; c < numCells(*G); c++) { - typedef typename Cell2FacesTraits::Type::row_type FaceRow; - FaceRow faces = c2f[c]; - - for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end(); - f!=end; ++f, ++i) - { - trans[*f] += 1.0 / htrans[i]; - } - } - - for (int f = 0; f < numFaces(*G); f++) { - trans[f] = 1.0 / trans[f]; - } -} - - -/* ---------------------------------------------------------------------- */ - template -void -tpfa_eff_trans_compute(const Grid* G, - const double *totmob, - const double *htrans, - double *trans) -/* ---------------------------------------------------------------------- */ -{ - using namespace Opm::UgGridHelpers; - - for (int f = 0; f < numFaces(*G); f++) { - trans[f] = 0.0; - } - - typename Cell2FacesTraits::Type c2f = cell2Faces(*G); - - for (int c = 0, i = 0; c < numCells(*G); c++) { - typedef typename Cell2FacesTraits::Type::row_type FaceRow; - FaceRow faces = c2f[c]; - - for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end(); - f!=end; ++f, ++i) - { - trans[*f] += 1.0 / (totmob[c] * htrans[i]); - } - } - - for (int f = 0; f < numFaces(*G); f++) { - trans[f] = 1.0 / trans[f]; - } -} diff --git a/opm/core/pressure/tpfa/trans_tpfa.c b/opm/core/pressure/tpfa/trans_tpfa.c deleted file mode 100644 index 7134be54..00000000 --- a/opm/core/pressure/tpfa/trans_tpfa.c +++ /dev/null @@ -1,130 +0,0 @@ -#include "config.h" -#include -#include -#include -#include - -#include -#include - - -#ifdef __cplusplus -#include "TransTpfa.hpp" -#endif - -/* ---------------------------------------------------------------------- */ -/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */ -/* ---------------------------------------------------------------------- */ -void -tpfa_htrans_compute(struct UnstructuredGrid *G, const double *perm, double *htrans) -/* ---------------------------------------------------------------------- */ -{ - #ifdef __cplusplus - return tpfa_htrans_compute(G, totmob, htrans, trans); - #endif - - int c, d, f, i, j; - double s, dist, denom; - - double Kn[3]; - double *cc, *fc, *n; - const double *K; - - MAT_SIZE_T nrows, ncols, ldA, incx, incy; - double a1, a2; - - d = G->dimensions; - - nrows = ncols = ldA = d; - incx = incy = 1 ; - a1 = 1.0; a2 = 0.0 ; - - for (c = i = 0; c < G->number_of_cells; c++) { - K = perm + (c * d * d); - cc = G->cell_centroids + (c * d); - - for (; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; - s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0; - - n = G->face_normals + (f * d); - fc = G->face_centroids + (f * d); - - dgemv_("No Transpose", &nrows, &ncols, - &a1, K, &ldA, n, &incx, &a2, &Kn[0], &incy); - - htrans[i] = denom = 0.0; - for (j = 0; j < d; j++) { - dist = fc[j] - cc[j]; - - htrans[i] += s * dist * Kn[j]; - denom += dist * dist; - } - - assert (denom > 0); - htrans[i] /= denom; - htrans[i] = fabs(htrans[i]); - } - } -} - - -/* ---------------------------------------------------------------------- */ -void -tpfa_trans_compute(struct UnstructuredGrid *G, const double *htrans, double *trans) -/* ---------------------------------------------------------------------- */ -{ - #ifdef __cplusplus - return tpfa_trans_compute(G, totmob, htrans, trans); - #endif - - int c, i, f; - - for (f = 0; f < G->number_of_faces; f++) { - trans[f] = 0.0; - } - - for (c = i = 0; c < G->number_of_cells; c++) { - for (; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; - - trans[f] += 1.0 / htrans[i]; - } - } - - for (f = 0; f < G->number_of_faces; f++) { - trans[f] = 1.0 / trans[f]; - } -} - - -/* ---------------------------------------------------------------------- */ -void -tpfa_eff_trans_compute(struct UnstructuredGrid *G, - const double *totmob, - const double *htrans, - double *trans) -/* ---------------------------------------------------------------------- */ -{ - #ifdef __cplusplus - return tpfa_eff_trans_compute(G, totmob, htrans, trans); - #endif - - int c, i, f; - - for (f = 0; f < G->number_of_faces; f++) { - trans[f] = 0.0; - } - - for (c = i = 0; c < G->number_of_cells; c++) { - for (; i < G->cell_facepos[c + 1]; i++) { - f = G->cell_faces[i]; - - trans[f] += 1.0 / (totmob[c] * htrans[i]); - } - } - - for (f = 0; f < G->number_of_faces; f++) { - trans[f] = 1.0 / trans[f]; - } -} diff --git a/opm/core/pressure/tpfa/trans_tpfa.h b/opm/core/pressure/tpfa/trans_tpfa.h deleted file mode 100644 index 5da1490f..00000000 --- a/opm/core/pressure/tpfa/trans_tpfa.h +++ /dev/null @@ -1,110 +0,0 @@ -/* - Copyright 2010 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_TRANS_TPFA_HEADER_INCLUDED -#define OPM_TRANS_TPFA_HEADER_INCLUDED - -/** - * \file - * Routines to assist in the calculation of two-point transmissibilities. - */ - -#include - -#ifdef __cplusplus -extern "C" { -#endif - -/** - * Calculate static, one-sided transmissibilities for use in the two-point flux - * approximation method. - * - * The one-sided transmissibilities are defined by the formula - * \f[ - * t_i = \frac{\vec{n}_f \mathsf{K}_c \vec{c}_c}{\lVert \vec{c}_c \rVert^2} - * \f] - * in which @c i is the half-face index corresponding to the cell-face index - * pair (c,f) and \f$\vec{c}_{cf} = \Bar{x}_f - \Bar{x}_c\f$ is the - * centroid difference vector. - * - * @param[in] G Grid. - * @param[in] perm Permeability. One symmetric, positive definite tensor - * per grid cell. - * @param[out] htrans One-sided transmissibilities. Array of size at least - * G->cell_facepos[ G->number_of_cells ]. - */ -void -tpfa_htrans_compute(struct UnstructuredGrid *G , - const double *perm , - double *htrans); - -/** - * Compute two-point transmissibilities from one-sided transmissibilities. - * - * The two-point transmissibilities are given by the simple, near-harmonic - * average (save a factor of two) - * \f[ - * \mathsf{T}_f = \big(\frac{1}{t_1} + \frac{1}{t_2}\big)^{-1} - * = \frac{t_1t_2}{t_1 + t_2} - * \f] - * in which \f$t_1\f$ and \f$t_2\f$ are the one-sided transmissibilities that - * connect the neighbouring cells of face @c f. - * - * @param[in] G Grid. - * @param[in] htrans One-sided transmissibilities as defined by function - * tpfa_htrans_compute(). - * @param[out] trans Interface, two-point transmissibilities. Array of size - * at least G->number_of_faces. - */ -void -tpfa_trans_compute(struct UnstructuredGrid *G , - const double *htrans, - double *trans ); - -/** - * Calculate effective two-point transmissibilities from one-sided, total - * mobility weighted, transmissibilities. - * - * Specifically, compute the following product - * \f[ - * \mathsf{T}_f = \big(\frac{1}{\lambda_1t_1} + \frac{1}{\lambda_2t_2}\big)^{-1} - * = \lambda_1\lambda_2 \frac{t_1t_2}{t_1 + t_2} - * \f] - * in which \f$t_1\f$ and \f$t_2\f$ are the one-sided, static transmissibility - * values connecting the cells of face @c f and \f$\lambda_1\f$ and - * \f$\lambda_2\f$ denote the total mobilities of the respective cells. - * - * @param[in] G Grid. - * @param[in] totmob Total mobilities. One positive scalar value for each cell. - * @param[in] htrans One-sided transmissibilities as defined by function - * tpfa_htrans_compute(). - * @param[out] trans Effective, two-point transmissibilities. Array of size at - * least G->number_of_faces. - */ -void -tpfa_eff_trans_compute(struct UnstructuredGrid *G , - const double *totmob, - const double *htrans, - double *trans ); - -#ifdef __cplusplus -} -#endif - -#endif /* OPM_TRANS_TPFA_HEADER_INCLUDED */ diff --git a/opm/core/utility/CompressedPropertyAccess.hpp b/opm/core/utility/CompressedPropertyAccess.hpp deleted file mode 100644 index 562557b8..00000000 --- a/opm/core/utility/CompressedPropertyAccess.hpp +++ /dev/null @@ -1,504 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2014 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_COMPRESSEDPROPERTYACCESS_HPP_HEADER -#define OPM_COMPRESSEDPROPERTYACCESS_HPP_HEADER - -/** - * \file - * - * Facility for accessing active subset of data arrays defined for all - * global cells. The main component is class template \code - * GridPropertyAccess::Compressed<> \endcode which encapsulates and - * provides read-only access to a data array and while translating - * active cell indices to "global" cell indices. The data array is a - * policy parameter for which preexisting implementations "constant" - * and "extract from ECLIPSE input" are defined in this module. Data - * values in the array must be defined for all global cells. - */ - -#include - -#include -#include -#include -#include -#include - -namespace Opm { - /** - * Nested name-space that serves no other purpose than to - * contextualise \c Compressed class name. - */ - namespace GridPropertyAccess { - /** - * Glue code for interacting with ECLIPSE input decks as - * defined by module opm-parser. - */ - namespace Details { - /** - * Implementation of property query and retrieval from an - * ECLIPSE property container. - */ - namespace EclPropImpl { - /** - * Property existence predicate. - * - * Supported for types \c int and \c double. - * - * \tparam T Property element type. - */ - template - struct HasProperty; - - /** - * Property value retrieval. - * - * Supported for types \c int and \c double. - * - * \tparam T Property element type. - */ - template - struct GetProperty; - - /** - * Specialization of property existence predicate for - * type \c int. - */ - template <> - struct HasProperty { - /** - * Existence predicate implementation. - * - * \tparam PropertyContainer Pointer type - * representing collection of (global) grid - * properties. Must implement method \c - * hasDeckIntGridProperty. - * - * \param[in] ecl Property container. - * - * \param[in] kw ECLIPSE property keyword. - * - * \return Whether or not property \c kw exists in - * the container \c ecl. - */ - template - static bool - p(PropertyContainer& ecl, - const std::string& kw); - }; - - template - bool - HasProperty::p(PropertyContainer& ecl, - const std::string& kw) - { - return ecl.get3DProperties().hasDeckIntGridProperty(kw); - } - - /** - * Specialization of property existence predicate for - * type \c double. - */ - template <> - struct HasProperty { - /** - * Existence predicate implementation. - * - * \tparam PropertyContainer Pointer type - * representing collection of (global) grid - * properties. Must implement method \c - * hasDeckDoubleGridProperty. - * - * \param[in] ecl Property container. - * - * \param[in] kw ECLIPSE property keyword. - * - * \return Whether or not property \c kw exists in - * the container \c ecl. - */ - template - static bool - p(PropertyContainer& ecl, - const std::string& kw); - }; - - template - bool - HasProperty::p(PropertyContainer& ecl, - const std::string& kw) - { - return ecl.get3DProperties().hasDeckDoubleGridProperty(kw); - } - - /** - * Specialization of property value retrieval for type - * \c int. - */ - template <> - struct GetProperty { - /** - * Property value retrieval implementation. - * - * \tparam PropertyContainer Pointer type - * representing collection of (global) grid - * properties. Must implement method \c - * getIntGridProperty. - * - * \param[in] ecl Property container. - * - * \param[in] kw ECLIPSE property keyword. - * - * \return Data values for property \c kw. - */ - template - static const GridProperty* - value(PropertyContainer& ecl, - const std::string& kw); - }; - - template - const GridProperty* - GetProperty::value(PropertyContainer& ecl, - const std::string& kw) - { - assert (HasProperty::p(ecl, kw)); - - return &ecl.get3DProperties().getIntGridProperty(kw); - } - - /** - * Specialization of property value retrieval for type - * \c double. - */ - template <> - struct GetProperty { - /** - * Property value retrieval implementation. - * - * \tparam PropertyContainer Pointer type - * representing collection of (global) grid - * properties. Must implement method \c - * getDoubleGridProperty. - * - * \param[in] ecl Property container. - * - * \param[in] kw ECLIPSE property keyword. - * - * \return Data values for property \c kw. - */ - template - static const GridProperty* - value(PropertyContainer& ecl, - const std::string& kw); - }; - - template - const GridProperty* - GetProperty::value(PropertyContainer& ecl, - const std::string& kw) - { - assert (HasProperty::p(ecl, kw)); - - return &ecl.get3DProperties().getDoubleGridProperty(kw); - } - } // namespace EclPropImpl - - /** - * Conditional retrieval of property values from an - * ECLIPSE input deck. - * - * Supported for types \c int and \c double. - * - * \tparam T Property element type. - */ - template - struct EclipsePropertyArray { - /** - * Retrieve property values if present in container. - * - * \tparam PropertyContainer Pointer type representing - * collection of (global) grid properties. - * - * \param[in] ecl Property container. - * - * \param[in] kw ECLIPSE property keyword. - * - * \return Data values for property \c kw if present, - * an empty \code shared_ptr<> \endcode if not. - */ - template - static const GridProperty* - value(PropertyContainer& ecl, - const std::string& kw); - }; - - template - template - const GridProperty* - EclipsePropertyArray::value(PropertyContainer& ecl, - const std::string& kw) - { - if (! EclPropImpl::HasProperty::p(ecl, kw)) { - return nullptr; - } - return EclPropImpl::GetProperty::value(ecl, kw); - } - } // namespace Details - - /** - * Predefined data array policies for use with class template - * \code Compressed<> \endcode. - */ - namespace ArrayPolicy { - /** - * Data array policy that extracts the array values from - * an ECLIPSE input deck or returns a user specified - * default value if the data vector is not present in a - * particular input deck. - * - * Provides read-only access to the underlying data. - * - * \tparam T Array element type. Must be \c int or \c - * double. - */ - template - class ExtractFromDeck { - public: - /** - * Constructor. - * - * \tparam PropertyContainer Pointer type representing - * collection of (global) grid properties. Typically - * \c EclipseStatePtr or \c EclipseStateConstPtr. - * Must implement methods \c hasDeckIntGridProperty and \c - * getIntGridProperty if \c T is \c int, or \c - * hasDeckDoubleGridProperty and \c getDoubleGridProperty - * if \c T is \c double. - * - * \param[in] ecl Property container. - * - * \param[in] kw ECLIPSE keyword from which to extract - * data array. - * - * \param[in] dlft Default/fall-back data array value - * if \c kw is not defined. - */ - template - ExtractFromDeck(PropertyContainer& ecl, - const std::string& kw, - const T dflt) - : x_ (Details::EclipsePropertyArray::value(ecl, kw)) - , dflt_(dflt) - {} - - /** - * Publicly accessible data array element type. - */ - typedef T value_type; - - /** - * Index type for accessing data array. - */ - typedef std::size_t size_type; - - /** - * Read-only data array access. - * - * \param[in] i Array index. Assumed to identify a - * global (uncompressed) cell. - * - * \return Data array element at global index \c i if - * present in input or user specified fall-back value - * if not. - */ - value_type - operator[](const size_type i) const - { - if (x_) { - return x_->iget(i); - } - else { - return dflt_; - } - } - - private: - /** - * Grid property handle. - * - * Null if data not defined. - */ - const GridProperty* x_ = nullptr; - - /** - * Fall-back data element value if data not defined. - */ - T dflt_; - }; - - /** - * Data array policy that returns a single, constant user - * specified value for every global cell. - * - * Provides read-only access to the underlying data. - * - * \tparam T Array element type. - */ - template - class Constant { - public: - /** - * Constructor - * - * \param[in] c Constant property value used for all - * global cells. - */ - Constant(const T c) - : c_(c) - {} - - /** - * Publicly accessible data array element type. - */ - typedef T value_type; - - /** - * Index type for accessing data array. - */ - typedef std::size_t size_type; - - /** - * Read-only data array access. - * - * \param[in] i Array index. Assumed to identify a - * global (uncompressed) cell. Unused. - * - * \return User specified constant value for every - * (global) cell. - */ - value_type - operator[](const size_type i) const - { - static_cast(i); // Suppress "unused parameter" - - return c_; - } - - private: - /** - * Constant, user specified property value. - */ - T c_; - }; - } // namespace ArrayPolicy - - /** - * Collection of tags to help enforce semantic type checks - * when using class \code Compressed<> \endcode. - */ - namespace Tag { - /** - * Default tag that implies no restriction. - */ - struct Any {}; - - /** - * Tag that restricts usage to NTG (net-to-gross) - * contexts. - */ - struct NTG : public Any {}; - } // namespace Tag - - /** - * Provide compressed (active cell) read-only access to - * globally defined data array. - * - * \tparam DataArray Type representing an array of data - * values, one value for each global (uncompressed) cell in a - * model. Must implement value semantics. Typically one of - * the array policies of name space \c ArrayPolicy. Must - * provide public type \c value_type to infer the data element - * type and \code operator[](i) \endcode to access the - * property value of the \c i'th global cell. - * - * \tparam PropertyTag Type tag that can be used to restrict - * applicability of the resulting \c Compressed array, e.g., - * to enforce net-to-gross ratios only. Default: No - * restriction. - */ - template - class Compressed { - public: - /** - * Constructor - * - * \param[in] x Preconfigured global property value array. - * The \c Compressed array creates a private copy of this - * object. - * - * \param[in] gc Compressed-to-global cell map. Typically - * the \c global_cell field of an \c UnstructuredGrid or - * something very similar. If null, interpreted as - * identity mapping, i.e., as if all cells are active. - */ - Compressed(const DataArray& x, - const int* gc) - : x_ (x) - , gc_(gc) - {} - - /** - * Property value type. - */ - typedef typename DataArray::value_type value_type; - - /** - * Read-only data array access. - * - * \param[in] c Active cell index. - * - * \return Property value in active cell \c c. - */ - value_type - operator[](const int c) const - { - return x_[ (gc_ == 0) ? c : gc_[c] ]; - } - - private: - /** - * Global property value array. - * - * Value semantics to support putting \c Compressed arrays - * into standard containers. - */ - DataArray x_; - - /** - * Compressed-to-global cell index map. \c Null if all - * cells active. - */ - const int* gc_; - }; - } // namespace GridPropertyAccess -} // namespace Opm - -#endif /* OPM_COMPRESSEDPROPERTYACCESS_HPP_HEADER */ diff --git a/opm/core/utility/RegionMapping.hpp b/opm/core/utility/RegionMapping.hpp deleted file mode 100644 index 63d412b0..00000000 --- a/opm/core/utility/RegionMapping.hpp +++ /dev/null @@ -1,181 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED -#define OPM_REGIONMAPPING_HEADER_INCLUDED - -#include - -#include -#include - -namespace Opm -{ - - /** - * Forward and reverse mappings between cells and - * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', - * 'PVTNUM', or 'EQUILNUM' arrays). - * - * \tparam Region Type of a forward region mapping. Expected - * to provide indexed access through - * operator[]() as well as inner types - * 'value_type', 'size_type', and - * 'const_iterator'. - */ - template < class Region = std::vector > - class RegionMapping { - public: - /** - * Constructor. - * - * \param[in] reg Forward region mapping, restricted to - * active cells only. - */ - explicit - RegionMapping(const Region& reg) - : reg_(reg) - { - rev_.init(reg_); - } - - /** - * Type of forward (cell-to-region) mapping result. - * Expected to be an integer. - */ - typedef typename Region::value_type RegionId; - - /** - * Type of reverse (region-to-cell) mapping (element) - * result. - */ - typedef typename Region::size_type CellId; - - /** - * Type of reverse region-to-cell range bounds and - * iterators. - */ - typedef typename std::vector::const_iterator CellIter; - - typedef boost::iterator_range Range; - - /** - * Compute region number of given active cell. - * - * \param[in] c Active cell - * \return Region to which @c c belongs. - */ - RegionId - region(const CellId c) const { return reg_[c]; } - - const std::vector& - activeRegions() const - { - return rev_.active; - } - - /** - * Extract active cells in particular region. - * - * \param[in] r Region number - * - * \return Range of active cells in region @c r. Empty if @c r is - * not an active region. - */ - Range - cells(const RegionId r) const { - const auto id = rev_.binid.find(r); - - if (id == rev_.binid.end()) { - // Region 'r' not an active region. Return empty. - return Range(rev_.c.end(), rev_.c.end()); - } - - const auto i = id->second; - - return Range(rev_.c.begin() + rev_.p[i + 0], - rev_.c.begin() + rev_.p[i + 1]); - } - - private: - /** - * Copy of forward region mapping (cell-to-region). - */ - Region reg_; - - /** - * Reverse mapping (region-to-cell). - */ - struct { - typedef typename std::vector::size_type Pos; - - std::unordered_map binid; - std::vector active; - - std::vector p; /**< Region start pointers */ - std::vector c; /**< Region cells */ - - /** - * Compute reverse mapping. Standard linear insertion - * sort algorithm. - */ - void - init(const Region& reg) - { - binid.clear(); - for (const auto& r : reg) { - ++binid[r]; - } - - p .clear(); p.emplace_back(0); - active.clear(); - { - Pos n = 0; - for (auto& id : binid) { - active.push_back(id.first); - p .push_back(id.second); - - id.second = n++; - } - } - - for (decltype(p.size()) i = 1, sz = p.size(); i < sz; ++i) { - p[0] += p[i]; - p[i] = p[0] - p[i]; - } - - assert (p[0] == static_cast(reg.size())); - - c.resize(reg.size()); - { - CellId i = 0; - for (const auto& r : reg) { - auto& pos = p[ binid[r] + 1 ]; - c[ pos++ ] = i++; - } - } - - p[0] = 0; - } - } rev_; /**< Reverse mapping instance */ - }; - -} // namespace Opm - -#endif // OPM_REGIONMAPPING_HEADER_INCLUDED diff --git a/opm/core/utility/VelocityInterpolation.cpp b/opm/core/utility/VelocityInterpolation.cpp deleted file mode 100644 index 48296584..00000000 --- a/opm/core/utility/VelocityInterpolation.cpp +++ /dev/null @@ -1,185 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include - -#include - -namespace Opm -{ - - // -------- Methods of class VelocityInterpolationInterface -------- - - - VelocityInterpolationInterface::~VelocityInterpolationInterface() - { - } - - - - // -------- Methods of class VelocityInterpolationConstant -------- - - /// Constructor. - /// \param[in] grid A grid. - VelocityInterpolationConstant::VelocityInterpolationConstant(const UnstructuredGrid& grid) - : grid_(grid) - { - } - - /// Set up fluxes for interpolation. - /// \param[in] flux One signed flux per face in the grid. - void VelocityInterpolationConstant::setupFluxes(const double* flux) - { - flux_ = flux; - } - - /// Interpolate velocity. - /// \param[in] cell Cell in which to interpolate. - /// \param[in] x Coordinates of point at which to interpolate. - /// Must be array of length grid.dimensions. - /// \param[out] v Interpolated velocity. - /// Must be array of length grid.dimensions. - void VelocityInterpolationConstant::interpolate(const int cell, - const double* /*x*/, - double* v) const - { - const int dim = grid_.dimensions; - std::fill(v, v + dim, 0.0); - const double* cc = grid_.cell_centroids + cell*dim; - for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { - const int face = grid_.cell_faces[hface]; - const double* fc = grid_.face_centroids + face*dim; - double face_flux = 0.0; - if (cell == grid_.face_cells[2*face]) { - face_flux = flux_[face]; - } else { - assert(cell == grid_.face_cells[2*face + 1]); - face_flux = -flux_[face]; - } - for (int dd = 0; dd < dim; ++dd) { - v[dd] += face_flux * (fc[dd] - cc[dd]) / grid_.cell_volumes[cell]; - } - } - } - - - // -------- Methods of class VelocityInterpolationECVI -------- - - - - /// Constructor. - /// \param[in] grid A grid. - VelocityInterpolationECVI::VelocityInterpolationECVI(const UnstructuredGrid& grid) - : bcmethod_(grid), grid_(grid) - { - } - - /// Set up fluxes for interpolation. - /// Computes the corner velocities. - /// \param[in] flux One signed flux per face in the grid. - void VelocityInterpolationECVI::setupFluxes(const double* flux) - { - // We must now update the velocity member of the CornerInfo - // for each corner. - const int dim = grid_.dimensions; - std::vector N(dim*dim); // Normals matrix. Fortran ordering! - std::vector orig_N(dim*dim); // Normals matrix. Fortran ordering! - std::vector f(dim); // Flux vector. - std::vector orig_f(dim); // Flux vector. - std::vector piv(dim); // For LAPACK solve - const SparseTable& all_ci = bcmethod_.cornerInfo(); - const std::vector& adj_faces = bcmethod_.adjacentFaces(); - corner_velocity_.resize(dim*all_ci.dataSize()); - const int num_cells = grid_.number_of_cells; - for (int cell = 0; cell < num_cells; ++cell) { - const int num_cell_corners = bcmethod_.numCorners(cell); - for (int cell_corner = 0; cell_corner < num_cell_corners; ++cell_corner) { - const int cid = all_ci[cell][cell_corner].corner_id; - for (int adj_ix = 0; adj_ix < dim; ++adj_ix) { - const int face = adj_faces[dim*cid + adj_ix]; - const double* fn = grid_.face_normals + dim*face; - for (int dd = 0; dd < dim; ++dd) { - N[adj_ix + dd*dim] = fn[dd]; // Row adj_ix, column dd - } - f[adj_ix] = flux[face]; - } - // Now we have built N and f. Solve Nv = f. - // Note that the face orientations do not matter, - // as changing an orientation would negate both a - // row in N and the corresponding element of f. - // Solving linear equation with LAPACK. - MAT_SIZE_T n = dim; - MAT_SIZE_T nrhs = 1; - MAT_SIZE_T lda = n; - MAT_SIZE_T ldb = n; - MAT_SIZE_T info = 0; - orig_N = N; - orig_f = f; - dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info); - if (info != 0) { - // Print the local matrix and rhs. - std::cerr << "Failed solving single-cell system Nv = f in cell " << cell - << " with N = \n"; - for (int row = 0; row < n; ++row) { - for (int col = 0; col < n; ++col) { - std::cerr << " " << orig_N[row + n*col]; - } - std::cerr << '\n'; - } - std::cerr << "and f = \n"; - for (int row = 0; row < n; ++row) { - std::cerr << " " << orig_f[row] << '\n'; - } - OPM_THROW(std::runtime_error, "Lapack error: " << info << " encountered in cell " << cell); - } - // The solution ends up in f, so we must copy it. - std::copy(f.begin(), f.end(), corner_velocity_.begin() + dim*cid); - } - } - } - - /// Interpolate velocity. - /// \param[in] cell Cell in which to interpolate. - /// \param[in] x Coordinates of point at which to interpolate. - /// Must be array of length grid.dimensions. - /// \param[out] v Interpolated velocity. - /// Must be array of length grid.dimensions. - void VelocityInterpolationECVI::interpolate(const int cell, - const double* x, - double* v) const - { - const int n = bcmethod_.numCorners(cell); - const int dim = grid_.dimensions; - bary_coord_.resize(n); - bcmethod_.cartToBary(cell, x, &bary_coord_[0]); - std::fill(v, v + dim, 0.0); - const SparseTable& all_ci = bcmethod_.cornerInfo(); - for (int i = 0; i < n; ++i) { - const int cid = all_ci[cell][i].corner_id; - for (int dd = 0; dd < dim; ++dd) { - v[dd] += corner_velocity_[dim*cid + dd] * bary_coord_[i]; - } - } - } - - -} // namespace Opm diff --git a/opm/core/utility/VelocityInterpolation.hpp b/opm/core/utility/VelocityInterpolation.hpp deleted file mode 100644 index 3b7a5bc4..00000000 --- a/opm/core/utility/VelocityInterpolation.hpp +++ /dev/null @@ -1,122 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_VELOCITYINTERPOLATION_HEADER_INCLUDED -#define OPM_VELOCITYINTERPOLATION_HEADER_INCLUDED - -#include -#include - -struct UnstructuredGrid; - -namespace Opm -{ - - /// Abstract interface for velocity interpolation method classes. - class VelocityInterpolationInterface - { - public: - virtual ~VelocityInterpolationInterface(); - - /// Set up fluxes for interpolation. - /// \param[in] flux One signed flux per face in the grid. - virtual void setupFluxes(const double* flux) = 0; - - /// Interpolate velocity. - /// \param[in] cell Cell in which to interpolate. - /// \param[in] x Coordinates of point at which to interpolate. - /// Must be array of length grid.dimensions. - /// \param[out] v Interpolated velocity. - /// Must be array of length grid.dimensions. - virtual void interpolate(const int cell, - const double* x, - double* v) const = 0; - }; - - - - - /// A constant velocity approximation. Will not actually interpolate - /// unless the fluxes are consistent with a constant velocity. - class VelocityInterpolationConstant : public VelocityInterpolationInterface - { - public: - /// Constructor. - /// \param[in] grid A grid. - explicit VelocityInterpolationConstant(const UnstructuredGrid& grid); - - /// Set up fluxes for interpolation. - /// \param[in] flux One signed flux per face in the grid. - virtual void setupFluxes(const double* flux); - - /// Interpolate velocity. - /// \param[in] cell Cell in which to interpolate. - /// \param[in] x Coordinates of point at which to interpolate. - /// Must be array of length grid.dimensions. - /// \param[out] v Interpolated velocity. - /// Must be array of length grid.dimensions. - virtual void interpolate(const int cell, - const double* x, - double* v) const; - private: - const UnstructuredGrid& grid_; - const double* flux_; - }; - - - - - - /// Interpolate velocity using the extended CVI scheme: - /// compute a corner velocity for each cell corner that - /// is consistent with fluxes of adjacent faces, then - /// interpolate with generalized barycentric coordinates. - class VelocityInterpolationECVI : public VelocityInterpolationInterface - { - public: - /// Constructor. - /// \param[in] grid A grid. - explicit VelocityInterpolationECVI(const UnstructuredGrid& grid); - - /// Set up fluxes for interpolation. - /// \param[in] flux One signed flux per face in the grid. - virtual void setupFluxes(const double* flux); - - /// Interpolate velocity. - /// \param[in] cell Cell in which to interpolate. - /// \param[in] x Coordinates of point at which to interpolate. - /// Must be array of length grid.dimensions. - /// \param[out] v Interpolated velocity. - /// Must be array of length grid.dimensions. - virtual void interpolate(const int cell, - const double* x, - double* v) const; - private: - WachspressCoord bcmethod_; - const UnstructuredGrid& grid_; - mutable std::vector bary_coord_; - std::vector corner_velocity_; // size = dim * #corners - }; - - - - -} // namespace Opm - -#endif // OPM_VELOCITYINTERPOLATION_HEADER_INCLUDED diff --git a/opm/core/utility/WachspressCoord.cpp b/opm/core/utility/WachspressCoord.cpp deleted file mode 100644 index 399bc9bf..00000000 --- a/opm/core/utility/WachspressCoord.cpp +++ /dev/null @@ -1,238 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include -#include -#include - -namespace Opm -{ - - - // -------- Helper methods for class WachspressCoord -------- - - namespace - { - /// Calculates the determinant of a 2 x 2 matrix, represented as - /// two two-dimensional arrays. - double determinantOf(const double* a0, - const double* a1) - { - return - a0[0] * a1[1] - a0[1] * a1[0]; - } - - /// Calculates the determinant of a 3 x 3 matrix, represented as - /// three three-dimensional arrays. - double determinantOf(const double* a0, - const double* a1, - const double* a2) - { - return - a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - - a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + - a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); - } - - /// Calculates the volume of the parallelepiped given by - /// the vectors n[i] for i = 0..(dim-1), each n[i] is of size dim. - double cornerVolume(double** n, const int dim) - { - assert(dim == 2 || dim == 3); - double det = (dim == 2) ? determinantOf(n[0], n[1]) : determinantOf(n[0], n[1], n[2]); - return std::fabs(det); - } - - } // anonymous namespace - - - - // -------- Methods of class WachspressCoord -------- - - // The formula used is a modification of the formula given in: - // M. Meyer, A. Barr, H. Lee, and M. Desbrun. - // Generalized barycentric coordinates on irregular poly- - // gons. Journal of Graphics Tools, 7(1):13–22, 2002. - // - // The formula given there is, for a corner i, - // b_i = w_i / sum_{k} w_k - // w_i = V_i / (prod_{j \in adjacent faces} n_j * (x_i - x) ) - // ^^^ ^^^ ^^^ - // corner "volume" normal corner coordinates - // V_i = |Det({n_j}_{j \in adjacent faces})| - // The corner coordinate x_i above can be replaced with any point on face j - // without changing the value of w_i, and we replace it with c_j, the face - // centroid. - // However, this formula has the problem that the denominator of w_i becomes zero - // close to the boundary. Our solution is to multiply all w_i by - /// prod_{all j} n_j * (c_j - x), resulting in the formula: - // w_i = V_i * (prod_{j \in nonadjacent faces} n_j * (c_j - x) ). - // Another implementation note is that the above formulas assumes that - // the normals have length 1 (are unit normals). It is easy to see that this - // can be relaxed, since each normal occurs once in the formula for w_i, and - // all w_i will be scaled by the same number. In our implementation we therefore - // use the area-scaled normals directly as provided by the UnstructuredGrid. - - /// Constructor. - /// \param[in] grid A grid. - WachspressCoord::WachspressCoord(const UnstructuredGrid& grid) - : grid_(grid) - { - enum { Maxdim = 3 }; - const int dim = grid.dimensions; - if (dim > Maxdim) { - OPM_THROW(std::runtime_error, "Grid has more than " << Maxdim << " dimensions."); - } - // Compute static data for each corner. - const int num_cells = grid.number_of_cells; - int corner_id_count = 0; - for (int cell = 0; cell < num_cells; ++cell) { - std::set cell_vertices; - std::vector cell_faces; - std::multimap vertex_adj_faces; - for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell + 1]; ++hface) { - const int face = grid.cell_faces[hface]; - cell_faces.push_back(face); - const int fn0 = grid.face_nodepos[face]; - const int fn1 = grid.face_nodepos[face + 1]; - cell_vertices.insert(grid.face_nodes + fn0, grid.face_nodes + fn1); - for (int fn = fn0; fn < fn1; ++fn) { - const int vertex = grid.face_nodes[fn]; - vertex_adj_faces.insert(std::make_pair(vertex, face)); - } - } - std::sort(cell_faces.begin(), cell_faces.end()); // set_difference requires sorted ranges - std::vector cell_corner_info; - std::set::const_iterator it = cell_vertices.begin(); - for (; it != cell_vertices.end(); ++it) { - CornerInfo ci; - ci.corner_id = corner_id_count++;; - ci.vertex = *it; - double* fnorm[Maxdim] = { 0 }; - typedef std::multimap::const_iterator MMIt; - std::pair frange = vertex_adj_faces.equal_range(ci.vertex); - int fi = 0; - std::vector vert_adj_faces(dim); - for (MMIt face_it = frange.first; face_it != frange.second; ++face_it, ++fi) { - if (fi >= dim) { - OPM_THROW(std::runtime_error, "In cell " << cell << ", vertex " << ci.vertex << " has " - << " more than " << dim << " adjacent faces."); - } - fnorm[fi] = grid_.face_normals + dim*(face_it->second); - vert_adj_faces[fi] = face_it->second; - } - assert(fi == dim); - adj_faces_.insert(adj_faces_.end(), vert_adj_faces.begin(), vert_adj_faces.end()); - const double corner_vol = cornerVolume(fnorm, dim); - ci.volume = corner_vol; - cell_corner_info.push_back(ci); - std::sort(vert_adj_faces.begin(), vert_adj_faces.end()); - std::vector vert_nonadj_faces(cell_faces.size() - vert_adj_faces.size()); - std::set_difference(cell_faces.begin(), cell_faces.end(), - vert_adj_faces.begin(), vert_adj_faces.end(), - vert_nonadj_faces.begin()); - nonadj_faces_.appendRow(vert_nonadj_faces.begin(), vert_nonadj_faces.end()); - } - corner_info_.appendRow(cell_corner_info.begin(), cell_corner_info.end()); - } - assert(corner_id_count == corner_info_.dataSize()); - } - - - - - /// Count of vertices adjacent to a call. - /// \param[in] cell A cell index. - /// \return Number of corners of cell. - int WachspressCoord::numCorners(const int cell) const - { - return corner_info_[cell].size(); - } - - - /// The class stores some info for each corner. - /// \return The corner info container. - const SparseTable& WachspressCoord::cornerInfo() const - { - return corner_info_; - } - - - - /// The class stores some info for each corner. - /// \return The corner info container. - const std::vector& WachspressCoord::adjacentFaces() const - { - return adj_faces_; - } - - - - /// Compute generalized barycentric coordinates for some point x - /// with respect to the vertices of a grid cell. - /// \param[in] cell Cell in which to compute coordinates. - /// \param[in] x Coordinates of point in cartesian coordinates. - /// Must be array of length grid.dimensions. - /// \param[out] xb Coordinates of point in barycentric coordinates. - /// Must be array of length numCorners(cell). - void WachspressCoord::cartToBary(const int cell, - const double* x, - double* xb) const - { - // Note: - // A possible optimization is: compute all n_j * (c_j - x) factors - // once, instead of repeating computation for all corners (for - // which j is a nonadjacent face). - const int n = numCorners(cell); - const int dim = grid_.dimensions; - double totw = 0.0; - for (int i = 0; i < n; ++i) { - const CornerInfo& ci = corner_info_[cell][i]; - // Weight (unnormalized) is equal to: - // V_i * (prod_{j \in nonadjacent faces} n_j * (c_j - x) ) - // ^^^ ^^^ ^^^ - // corner "volume" normal centroid - xb[i] = ci.volume; - const int num_nonadj_faces = nonadj_faces_[ci.corner_id].size(); - for (int j = 0; j < num_nonadj_faces; ++j) { - const int face = nonadj_faces_[ci.corner_id][j]; - double factor = 0.0; - for (int dd = 0; dd < dim; ++dd) { - factor += grid_.face_normals[dim*face + dd]*(grid_.face_centroids[dim*face + dd] - x[dd]); - } - // Assumes outward-pointing normals, so negate factor if necessary. - if (grid_.face_cells[2*face] != cell) { - assert(grid_.face_cells[2*face + 1] == cell); - factor = -factor; - } - xb[i] *= factor; - } - totw += xb[i]; - } - for (int i = 0; i < n; ++i) { - xb[i] /= totw; - } - } - - - -} // namespace Opm diff --git a/opm/core/utility/WachspressCoord.hpp b/opm/core/utility/WachspressCoord.hpp deleted file mode 100644 index 73312ef6..00000000 --- a/opm/core/utility/WachspressCoord.hpp +++ /dev/null @@ -1,85 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_WACHSPRESSCOORD_HEADER_INCLUDED -#define OPM_WACHSPRESSCOORD_HEADER_INCLUDED - -#include -#include - -struct UnstructuredGrid; - -namespace Opm -{ - - /// Class capable of computing Wachspress coordinates in 2d and 3d. - /// The formula used is a modification of the formula given in: - /// M. Meyer, A. Barr, H. Lee, and M. Desbrun. - /// Generalized barycentric coordinates on irregular poly- - /// gons. Journal of Graphics Tools, 7(1):13–22, 2002. - class WachspressCoord - { - public: - /// Constructor. - /// \param[in] grid A grid. - explicit WachspressCoord(const UnstructuredGrid& grid); - - /// Count of vertices adjacent to a call. - /// \param[in] cell A cell index. - /// \return Number of corners of cell. - int numCorners(const int cell) const; - - /// Compute generalized barycentric coordinates for some point x - /// with respect to the vertices of a grid cell. - /// \param[in] cell Cell in which to compute coordinates. - /// \param[in] x Coordinates of point in cartesian coordinates. - /// Must be array of length grid.dimensions. - /// \param[out] xb Coordinates of point in barycentric coordinates. - /// Must be array of length numCorners(cell). - void cartToBary(const int cell, - const double* x, - double* xb) const; - - // A corner is here defined as a {cell, vertex} pair where the - // vertex is adjacent to the cell. - struct CornerInfo - { - int corner_id; // Unique for each corner. - int vertex; // Shared between corners belonging to different cells. - double volume; // Defined as det(N) where N is the matrix of adjacent face normals. - }; - - /// The class stores some info for each corner, made accessible for user convenience. - /// \return The corner info container. - const SparseTable& cornerInfo() const; - - /// The class stores adjacent faces for each corner, made accessible for user convenience. - /// \return The vector of adjacent faces. Size = dim * #corners. - const std::vector& adjacentFaces() const; - - private: - const UnstructuredGrid& grid_; - SparseTable corner_info_; // Corner info by cell. - std::vector adj_faces_; // Set of adjacent faces, by corner id. Contains dim face indices per corner. - SparseTable nonadj_faces_; // Set of nonadjacent faces, by corner id. - }; - -} // namespace Opm - -#endif // OPM_WACHSPRESSCOORD_HEADER_INCLUDED diff --git a/opm/core/utility/compressedToCartesian.cpp b/opm/core/utility/compressedToCartesian.cpp deleted file mode 100644 index 641ee817..00000000 --- a/opm/core/utility/compressedToCartesian.cpp +++ /dev/null @@ -1,48 +0,0 @@ -/* - Copyright 2015 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - - -#include -#include - -namespace Opm -{ - - // Construct explicit mapping from active/compressed to logical cartesian - // indices, either as given in global_cell or as { 0, 1, 2, ....} if null. - // \param[in] num_cells The number of active cells. - // \param[in] global_cell Either null, or an array of size num_cells. - // \return A vector containing the same data as global_cell, - // or the sequence { 0, 1, ... , num_cells - 1 } if - // global_cell was null. - std::vector compressedToCartesian(const int num_cells, - const int* global_cell) - { - std::vector retval; - if (global_cell) { - retval.assign(global_cell, global_cell + num_cells); - } else { - retval.resize(num_cells); - std::iota(retval.begin(), retval.end(), 0); - } - return retval; - } - - -} // namespace Opm diff --git a/opm/core/utility/compressedToCartesian.hpp b/opm/core/utility/compressedToCartesian.hpp deleted file mode 100644 index e0db9076..00000000 --- a/opm/core/utility/compressedToCartesian.hpp +++ /dev/null @@ -1,40 +0,0 @@ -/* - Copyright 2015 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_COMPRESSEDTOCARTESIAN_HEADER_INCLUDED -#define OPM_COMPRESSEDTOCARTESIAN_HEADER_INCLUDED - -#include - -namespace Opm -{ - - // Construct explicit mapping from active/compressed to logical cartesian - // indices, either as given in global_cell or as { 0, 1, 2, ....} if null. - // \param[in] num_cells The number of active cells. - // \param[in] global_cell Either null, or an array of size num_cells. - // \return A vector containing the same data as global_cell, - // or the sequence { 0, 1, ... , num_cells - 1 } if - // global_cell was null. - std::vector compressedToCartesian(const int num_cells, - const int* global_cell); - -} // namespace Opm - -#endif // OPM_COMPRESSEDTOCARTESIAN_HEADER_INCLUDED diff --git a/opm/core/utility/extractPvtTableIndex.cpp b/opm/core/utility/extractPvtTableIndex.cpp deleted file mode 100644 index 738039dd..00000000 --- a/opm/core/utility/extractPvtTableIndex.cpp +++ /dev/null @@ -1,47 +0,0 @@ -/* - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#include -#include - -#include "extractPvtTableIndex.hpp" - -#include -#include - -#include - -namespace Opm { - -void extractPvtTableIndex(std::vector &pvtTableIdx, - const Opm::EclipseState& eclState, - size_t numCompressed, - const int *compressedToCartesianCellIdx) -{ - //Get the PVTNUM data - const std::vector& pvtnumData = eclState.get3DProperties().getIntGridProperty("PVTNUM").getData(); - // Convert this into an array of compressed cells - // Eclipse uses Fortran-style indices which start at 1 - // instead of 0, we subtract 1. - pvtTableIdx.resize(numCompressed); - for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) { - size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx; - assert(cartesianCellIdx < pvtnumData.size()); - pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1; - } -} - -} diff --git a/opm/core/utility/extractPvtTableIndex.hpp b/opm/core/utility/extractPvtTableIndex.hpp deleted file mode 100644 index 72447046..00000000 --- a/opm/core/utility/extractPvtTableIndex.hpp +++ /dev/null @@ -1,33 +0,0 @@ -/* - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#ifndef OPM_EXTRACT_PVT_TABLE_INDEX_HPP -#define OPM_EXTRACT_PVT_TABLE_INDEX_HPP - -#include - -#include - -namespace Opm { - -void extractPvtTableIndex(std::vector &pvtTableIdx, - const Opm::EclipseState& eclState, - size_t numCompressed, - const int *compressedToCartesianCellIdx); - -} - -#endif // OPM_EXTRACT_PVT_TABLE_INDEX_HPP diff --git a/tests/compressed_gridproperty.data b/tests/compressed_gridproperty.data deleted file mode 100644 index 575d108d..00000000 --- a/tests/compressed_gridproperty.data +++ /dev/null @@ -1,48 +0,0 @@ -RUNSPEC - -DIMENS - 1 1 4 -/ - -GRID - -COORD - 0 0 0 0 0 0 - 1 0 0 1 0 0 - 0 1 0 0 1 0 - 1 1 0 1 1 0 -/ - -ZCORN - 0 0 0 0 - 1 1 1 1 - - 1 1 1 1 - 2 2 2 2 - - 2 2 2 2 - 3 3 3 3 - - 3 3 3 3 - 4 4 4 4 -/ - -ACTNUM - 0 1 0 1 -/ - -NTG - 0.1 0.2 0.3 0.4 -/ - -REGIONS - -SATNUM - 4 3 2 1 -/ - --- Note: Dummy 'PROPS' and 'SCHEDULE' sections only -PROPS -SCHEDULE - -END diff --git a/tests/test_compressedpropertyaccess.cpp b/tests/test_compressedpropertyaccess.cpp deleted file mode 100644 index b45715f2..00000000 --- a/tests/test_compressedpropertyaccess.cpp +++ /dev/null @@ -1,347 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2014 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include - -#if HAVE_DYNAMIC_BOOST_TEST -#define BOOST_TEST_DYN_LINK -#endif - -#define NVERBOSE - -#define BOOST_TEST_MODULE CompressedPropertyTest - -#include -#include - -#include - -#include -#include -#include -#include - -#include -#include - -struct SetupSimple { - SetupSimple() : - deck( Opm::Parser{}.parseFile( "compressed_gridproperty.data", Opm::ParseContext{} ) ), - ecl( deck ) - {} - - Opm::Deck deck; - Opm::EclipseState ecl; -}; - - -template -struct TestFixture : public Setup -{ - TestFixture() - : Setup () - , grid (ecl.getInputGrid()) - , reltol(1.0e-10) - { - } - - using Setup::deck; - using Setup::ecl; - - Opm::GridManager grid; - const double reltol; -}; - - -BOOST_AUTO_TEST_SUITE(CompressedPropertyHandling) - -// Construct global array by extracting an "undefined" (unspecified) -// double vector from the input deck. -BOOST_FIXTURE_TEST_CASE(APExtractDoubleUndefined, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalDoubleArray; - - ECLGlobalDoubleArray multpv(ecl, "MULTPV", 1.0); - - BOOST_CHECK_CLOSE(multpv[0], 1.0, reltol); - BOOST_CHECK_CLOSE(multpv[1], 1.0, reltol); - BOOST_CHECK_CLOSE(multpv[2], 1.0, reltol); - BOOST_CHECK_CLOSE(multpv[3], 1.0, reltol); -} - - -// Construct global array by extracting a fully defined (specified) -// double vector from the input deck. -BOOST_FIXTURE_TEST_CASE(APExtractDoubleDefined, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalDoubleArray; - - ECLGlobalDoubleArray ntg(ecl, "NTG", 1.0); - - BOOST_CHECK_CLOSE(ntg[0], 0.1, reltol); - BOOST_CHECK_CLOSE(ntg[1], 0.2, reltol); - BOOST_CHECK_CLOSE(ntg[2], 0.3, reltol); - BOOST_CHECK_CLOSE(ntg[3], 0.4, reltol); -} - - -// Construct global array by extracting an undefined (unspecified) -// integer (int) vector from the input deck. -BOOST_FIXTURE_TEST_CASE(APExtractIntUndefined, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalIntArray; - - ECLGlobalIntArray imbnum(ecl, "IMBNUM", 1); - - BOOST_CHECK_EQUAL(imbnum[0], 1); - BOOST_CHECK_EQUAL(imbnum[1], 1); - BOOST_CHECK_EQUAL(imbnum[2], 1); - BOOST_CHECK_EQUAL(imbnum[3], 1); -} - - -// Construct global array by extracting a fully defined (specified) -// integer (int) vector from the input deck. -BOOST_FIXTURE_TEST_CASE(APExtractIntDefined, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalIntArray; - - ECLGlobalIntArray satnum(ecl, "SATNUM", 1); - - BOOST_CHECK_EQUAL(satnum[0], 4); - BOOST_CHECK_EQUAL(satnum[1], 3); - BOOST_CHECK_EQUAL(satnum[2], 2); - BOOST_CHECK_EQUAL(satnum[3], 1); -} - - -// Construct global, infinitely sized, double array by specifying a -// single constant for all cells. -BOOST_FIXTURE_TEST_CASE(APConstantDouble, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::Constant ConstantDoubleArray; - - const double c = 1.234e5; - - ConstantDoubleArray x(c); - - BOOST_CHECK_CLOSE(x[ 0], c, reltol); - BOOST_CHECK_CLOSE(x[ 1], c, reltol); - BOOST_CHECK_CLOSE(x[ 2], c, reltol); - BOOST_CHECK_CLOSE(x[10000], c, reltol); // Infinite size array -} - - -// Construct global, infinitely sized, integer (int) array by -// specifying a single constant for all cells. -BOOST_FIXTURE_TEST_CASE(APConstantInt, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::Constant ConstantIntArray; - - const int i = 12345; - - ConstantIntArray a(i); - - BOOST_CHECK_EQUAL(a[ 0], i); - BOOST_CHECK_EQUAL(a[ 1], i); - BOOST_CHECK_EQUAL(a[ 2], i); - BOOST_CHECK_EQUAL(a[10001], i); // Inifinite size array -} - - -// Construct compressed double array based on global, undefined array -// extracted from input deck. Default ("any") type-check tag. -BOOST_FIXTURE_TEST_CASE(CAExtractDoubleUndefinedAny, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalDoubleArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - ECLGlobalDoubleArray multpv_glob(ecl, "MULTPV", 1.0); - CompressedArray multpv(multpv_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_CLOSE(multpv[0], 1.0, reltol); - BOOST_CHECK_CLOSE(multpv[1], 1.0, reltol); -} - - -// Construct compressed double array based on global, fully specified -// array extracted from input deck. Type-check tag: NTG. -BOOST_FIXTURE_TEST_CASE(CAExtractDoubleDefinedNTG, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalDoubleArray; - - typedef Opm::GridPropertyAccess::Tag::NTG NTG; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - ECLGlobalDoubleArray ntg_glob(ecl, "NTG", 1.0); - CompressedArray ntg(ntg_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_CLOSE(ntg[0], 0.2, reltol); - BOOST_CHECK_CLOSE(ntg[1], 0.4, reltol); -} - - -// Construct compressed integer (int) array based on global, undefined -// (unspecified) array extracted from input deck. Default ("any") -// type-check tag. -BOOST_FIXTURE_TEST_CASE(CAExtractIntUndefinedAny, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalIntArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - ECLGlobalIntArray imbnum_glob(ecl, "IMBNUM", 1); - CompressedArray imbnum(imbnum_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_EQUAL(imbnum[0], 1); - BOOST_CHECK_EQUAL(imbnum[1], 1); -} - - -// Construct compressed integer (int) array based on global, fully -// specified array extracted from input deck. Custom type-check tag. -struct RegionID : public Opm::GridPropertyAccess::Tag::Any {}; - -BOOST_FIXTURE_TEST_CASE(CAExtractIntDefinedCustom, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::ExtractFromDeck ECLGlobalIntArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - ECLGlobalIntArray satnum_glob(ecl, "SATNUM", 1); - CompressedArray satnum(satnum_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_EQUAL(satnum[0], 3); - BOOST_CHECK_EQUAL(satnum[1], 1); -} - - -// Construct compressed double array based on global constant value -// for all cells. Default ("any") type-check tag. -BOOST_FIXTURE_TEST_CASE(CAConstantDoubleAny, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::Constant ConstantDoubleArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - const double c = 1.234e5; - - ConstantDoubleArray x_glob(c); - CompressedArray x(x_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_CLOSE(x[0], c, reltol); - BOOST_CHECK_CLOSE(x[1], c, reltol); -} - - -// Construct compressed double array based on global constant value -// for all cells. Custom type-check tag. -struct MyTag : public Opm::GridPropertyAccess::Tag::Any {}; - -BOOST_FIXTURE_TEST_CASE(CAConstantDoubleCustom, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::Constant ConstantDoubleArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - const double c = 1.234e5; - - ConstantDoubleArray x_glob(c); - CompressedArray x(x_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_CLOSE(x[0], c, reltol); - BOOST_CHECK_CLOSE(x[1], c, reltol); -} - - -// Construct compressed integer (int) array based on global constant -// value for all cells. Custom type-check tag. -BOOST_FIXTURE_TEST_CASE(CAConstantIntAny, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::Constant ConstantIntArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - const int i = 12345; - - ConstantIntArray x_glob(i); - CompressedArray x(x_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_EQUAL(x[0], i); - BOOST_CHECK_EQUAL(x[1], i); -} - - -// Construct compressed integer (int) array based on global constant -// value for all cells. Custom type-check tag. -BOOST_FIXTURE_TEST_CASE(CAConstantIntCustom, - TestFixture) -{ - typedef Opm::GridPropertyAccess::ArrayPolicy - ::Constant ConstantIntArray; - - typedef Opm::GridPropertyAccess:: - Compressed CompressedArray; - - const int i = 12345; - - ConstantIntArray x_glob(i); - CompressedArray x(x_glob, grid.c_grid()->global_cell); - - BOOST_CHECK_EQUAL(x[0], i); - BOOST_CHECK_EQUAL(x[1], i); -} - -BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/test_regionmapping.cpp b/tests/test_regionmapping.cpp deleted file mode 100644 index 760a719f..00000000 --- a/tests/test_regionmapping.cpp +++ /dev/null @@ -1,155 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" - -/* --- Boost.Test boilerplate --- */ -#if HAVE_DYNAMIC_BOOST_TEST -#define BOOST_TEST_DYN_LINK -#endif - -#define NVERBOSE // Suppress own messages when throw()ing - -#define BOOST_TEST_MODULE RegionMapping - -#include -#include -#include - -/* --- our own headers --- */ - -#include - -#include -#include - -BOOST_AUTO_TEST_SUITE () - - -BOOST_AUTO_TEST_CASE (Forward) -{ - std::vector regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 }; - - Opm::RegionMapping<> rm(regions); - - for (decltype(regions.size()) i = 0, n = regions.size(); i < n; ++i) { - BOOST_CHECK_EQUAL(rm.region(i), regions[i]); - } -} - - -BOOST_AUTO_TEST_CASE (ActiveRegions) -{ - // 0 1 2 3 4 5 6 7 8 - std::vector regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 }; - - Opm::RegionMapping<> rm(regions); - - std::vector region_ids = { 2, 3, 4, 5, 6, 7 }; - - auto active = [®ion_ids](const int reg) - { - auto b = region_ids.begin(); - auto e = region_ids.end(); - - return std::find(b, e, reg) != e; - }; - - BOOST_CHECK_EQUAL(rm.activeRegions().size(), region_ids.size()); - - for (const auto& reg : rm.activeRegions()) { - BOOST_CHECK(active(reg)); - } -} - - -BOOST_AUTO_TEST_CASE (Consecutive) -{ - using RegionCells = std::map>; - - // 0 1 2 3 4 5 6 7 8 - std::vector regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 }; - - Opm::RegionMapping<> rm(regions); - - std::vector region_ids = { 2, 3, 4, 5, 6, 7 }; - RegionCells region_cells; - { - using VT = RegionCells::value_type; - - region_cells.insert(VT(2, { 0, 2, 4 })); - region_cells.insert(VT(3, { 7 })); - region_cells.insert(VT(4, { 3 })); - region_cells.insert(VT(5, { 1 })); - region_cells.insert(VT(6, { 6, 8 })); - region_cells.insert(VT(7, { 5 })); - } - - for (const auto& reg : region_ids) { - const auto& cells = rm.cells(reg); - const auto& expect = region_cells[reg]; - - BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(), - expect.begin(), expect.end()); - } - - // Verify that there are no cells in unused regions 0 and 1. - for (const auto& r : { 0, 1 }) { - BOOST_CHECK(rm.cells(r).empty()); - } -} - - -BOOST_AUTO_TEST_CASE (NonConsecutive) -{ - using RegionCells = std::map>; - - // 0 1 2 3 4 5 6 7 8 - std::vector regions = { 2, 4, 2, 4, 2, 7, 6, 3, 6 }; - - Opm::RegionMapping<> rm(regions); - - std::vector region_ids = { 2, 3, 4, 6, 7 }; - RegionCells region_cells; - { - using VT = RegionCells::value_type; - - region_cells.insert(VT(2, { 0, 2, 4 })); - region_cells.insert(VT(3, { 7 })); - region_cells.insert(VT(4, { 1, 3 })); - region_cells.insert(VT(6, { 6, 8 })); - region_cells.insert(VT(7, { 5 })); - } - - for (const auto& reg : region_ids) { - const auto& cells = rm.cells(reg); - const auto& expect = region_cells[reg]; - - BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(), - expect.begin(), expect.end()); - } - - // Verify that there are no cells in unused regions 0, 1, and 5. - for (const auto& r : { 0, 1, 5 }) { - BOOST_CHECK(rm.cells(r).empty()); - } -} - - -BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/test_velocityinterpolation.cpp b/tests/test_velocityinterpolation.cpp deleted file mode 100644 index f4b4c3f9..00000000 --- a/tests/test_velocityinterpolation.cpp +++ /dev/null @@ -1,487 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include - -#if HAVE_DYNAMIC_BOOST_TEST -#define BOOST_TEST_DYN_LINK -#endif -#define NVERBOSE // to suppress our messages when throwing - -#define BOOST_TEST_MODULE VelocityInterpolationTest -#include - -#include -#include -#include -#include - -using namespace Opm; - -namespace -{ - - // Compute flux corresponding to a constant velocity vector v. - void computeFlux(const UnstructuredGrid& grid, const std::vector& v, std::vector& flux) - { - const int dim = v.size(); - assert(dim == grid.dimensions); - flux.resize(grid.number_of_faces); - for (int face = 0; face < grid.number_of_faces; ++face) { - flux[face] = std::inner_product(v.begin(), v.end(), grid.face_normals + face*dim, 0.0); - } - } - - // Compute a linear vector function v = v0 + x*v1. - void computeLinearVec(const std::vector& v0, - const std::vector& v1, - const std::vector& x, - std::vector& v) - { - assert(v0.size() == v1.size() && v0.size() == x.size()); - const int dim = v0.size(); - v.resize(dim); - for (int dd = 0; dd < dim; ++dd) { - v[dd] = v0[dd] + x[dd]*v1[dd]; - } - } - -#ifdef __clang__ -#pragma clang diagnostic push -#pragma clang diagnostic ignored "-Wunneeded-internal-declaration" -#endif /* __clang__ */ - // Compute flux corresponding to a velocity vector v = v0 + x*v1. - void computeFluxLinear(const UnstructuredGrid& grid, - const std::vector& v0, - const std::vector& v1, - std::vector& flux) - { - const int dim = v0.size(); - assert(dim == grid.dimensions); - flux.resize(grid.number_of_faces); - std::vector x(dim); - std::vector v(dim); - for (int face = 0; face < grid.number_of_faces; ++face) { - const double* fc = grid.face_centroids + face*dim; - std::copy(fc, fc + dim, x.begin()); - computeLinearVec(v0, v1, x, v); - flux[face] = std::inner_product(v.begin(), v.end(), grid.face_normals + face*dim, 0.0); - } - } -#ifdef __clang__ -#pragma clang diagnostic pop -#endif /* __clang__ */ - - double vectorDiff2(const std::vector& v1, const std::vector& v2) - { - assert(v1.size() == v2.size()); - const int sz = v1.size(); - double vdiff = 0.0; - for (int i = 0; i < sz; ++i) { - vdiff += (v1[i] - v2[i])*(v1[i] - v2[i]); - } - return vdiff; - } - -} // anonymous namespace - - -template -void testConstantVelRepro2d() -{ - // Set up 2d 1-cell cartesian case. - GridManager g(1, 1); - const UnstructuredGrid& grid = *g.c_grid(); - std::vector v(2); - v[0] = 0.12345; - v[1] = -0.6789; - std::vector flux; - computeFlux(grid, v, flux); - VelInterp vic(grid); - vic.setupFluxes(&flux[0]); - - // Test a few points - std::vector x(2); - x[0] = 0.23456; - x[1] = 0.87654; - std::vector v_interp(2); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 1.0; - x[1] = 0.5; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 1.0; - x[1] = 1.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); -} - - -namespace -{ - - // Data for a pyramid. Node numbering goes - // lexicographic on bottom, then top. - // Face numbering goes xmin, xmax, ymin, ymax, bottom. - namespace Pyramid - { - static int face_nodes[] = { 0, 4, 2, 3, 4, 1, 0, 1, 4, 4, 3, 2, 0, 2, 3, 1, }; - static int face_nodepos[] = { 0, 3, 6, 9, 12, 16 }; - static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; - static int cell_faces[] = { 0, 1, 2, 3, 4 }; - static int cell_facepos[] = { 0, 5 }; - static double node_coordinates[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; - static double face_centroids[] = { 0, 1.0/3.0, 1.0/3.0, - 2.0/3.0, 1.0/3.0, 1.0/3.0, - 1.0/3.0, 0, 1.0/3.0, - 1.0/3.0, 2.0/3.0, 1.0/3.0, - 0.5, 0.5, 0 }; - static double face_areas[] = { 0.5, std::sqrt(2.0), 0.5, std::sqrt(2.0), 1.0 }; - static double face_normals[] = { -0.5000, 0, 0, - 0.5000, 0, 0.5000, - 0, -0.5000, 0, - 0, 0.5000, 0.5000, - 0, 0, -1.0000 }; - static double cell_centroids[] = { 0.375, 0.375, 0.25 }; - static double cell_volumes[] = { 1.0/3.0 }; - - } // namespace Pyramid - - UnstructuredGrid makePyramid() - { - // Make a 3d 1-cell grid, where the single cell is a pyramid. - UnstructuredGrid grid; - grid.dimensions = 3; - grid.number_of_cells = 1; - grid.number_of_faces = 5; - grid.number_of_nodes = 5; - grid.face_nodes = Pyramid::face_nodes; - grid.face_nodepos = Pyramid::face_nodepos; - grid.face_cells = Pyramid::face_cells; - grid.cell_faces = Pyramid::cell_faces; - grid.cell_facepos = Pyramid::cell_facepos; - grid.node_coordinates = Pyramid::node_coordinates; - grid.face_centroids = Pyramid::face_centroids; - grid.face_areas = Pyramid::face_areas; - grid.face_normals = Pyramid::face_normals; - grid.cell_centroids = Pyramid::cell_centroids; - grid.cell_volumes = Pyramid::cell_volumes; - return grid; - } - -} // anonymous namespace - - -template -void testConstantVelReproPyramid() -{ - // Set up a 3d 1-cell non-cartesian case (a pyramid). - UnstructuredGrid grid = makePyramid(); - std::vector v(3); - v[0] = 0.12345; - v[1] = -0.6789; - v[2] = 0.3456; - std::vector flux; - computeFlux(grid, v, flux); - VelInterp vic(grid); - vic.setupFluxes(&flux[0]); - - // Test a few points - std::vector x(3); - x[0] = 0.123; - x[1] = 0.0123; - x[2] = 0.213; - std::vector v_interp(3); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - x[2] = 1.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.1; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); -} - - -namespace -{ - // Data for an irregular 2d polygon. - namespace Irreg2d - { - static int face_nodes[] = { 0, 1, 1, 2, 2, 3, 3, 4, 4, 0 }; - static int face_nodepos[] = { 0, 2, 4, 6, 8, 10 }; - static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; - static int cell_faces[] = { 0, 1, 2, 3, 4 }; - static int cell_facepos[] = { 0, 5 }; - static double node_coordinates[] = { 0, 0, 3, 0, 3, 2, 1, 3, 0, 2 }; - static double face_centroids[] = { 1.5, 0, 3, 1, 2, 2.5, 0.5, 2.5, 0, 1 }; - static double face_areas[] = { 3, 2, std::sqrt(5.0), std::sqrt(2.0), 2 }; - static double face_normals[] = { 0, -3, 2, 0, 1, 2, -1, 1, -2, 0 }; - static double cell_centroids[] = { 22.0/15.0, 19.0/15.0 }; - static double cell_volumes[] = { 7.5 }; - } // namespace Irreg2d - - UnstructuredGrid makeIrreg2d() - { - // Make a 2d 1-cell grid, where the single cell is a polygon. - UnstructuredGrid grid; - grid.dimensions = 2; - grid.number_of_cells = 1; - grid.number_of_faces = 5; - grid.number_of_nodes = 5; - grid.face_nodes = Irreg2d::face_nodes; - grid.face_nodepos = Irreg2d::face_nodepos; - grid.face_cells = Irreg2d::face_cells; - grid.cell_faces = Irreg2d::cell_faces; - grid.cell_facepos = Irreg2d::cell_facepos; - grid.node_coordinates = Irreg2d::node_coordinates; - grid.face_centroids = Irreg2d::face_centroids; - grid.face_areas = Irreg2d::face_areas; - grid.face_normals = Irreg2d::face_normals; - grid.cell_centroids = Irreg2d::cell_centroids; - grid.cell_volumes = Irreg2d::cell_volumes; - return grid; - } - -} // anonymous namespace - - -template -void testConstantVelReproIrreg2d() -{ - // Set up a 2d 1-cell non-cartesian case (a pyramid). - UnstructuredGrid grid = makeIrreg2d(); - std::vector v(2); - v[0] = 0.12345; - v[1] = -0.6789; - std::vector flux; - computeFlux(grid, v, flux); - VelInterp vic(grid); - vic.setupFluxes(&flux[0]); - - // Test a few points - std::vector x(2); - x[0] = 1.2345; - x[1] = 2.0123; - std::vector v_interp(2); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 1.0; - x[1] = 3.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 3.0; - x[1] = 1.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); -} - - - -namespace -{ - // Data for an irregular 3d prism. - namespace IrregPrism - { - static int face_nodes[] = { 0, 4, 2, 1, 3, 5, 0, 1, 5, 4, 2, 4, 5, 3, 2, 3, 0, 1}; - static int face_nodepos[] = { 0, 3, 6, 10, 14, 18 }; - static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; - static int cell_faces[] = { 0, 1, 2, 3, 4 }; - static int cell_facepos[] = { 0, 5 }; - static double node_coordinates[] = { 0, 0, 0, - 2, 0, 0, - 0, 1, 0, - 2, 1, 0, - 0, 0, 1, - 1, 0, 1 }; - static double face_centroids[] = { 0, 1.0/3.0, 1.0/3.0, - 5.0/3.0, 1.0/3.0, 1.0/3.0, - 7.0/9.0, 0, 4.0/9.0, - 7.0/9.0, 5.0/9.0, 4.0/9.0, - 1, 0.5, 0 }; - static double face_areas[] = { 0.500000000000000, - 0.707106781186548, - 1.500000000000000, - 2.121320343559642, - 2.000000000000000 }; - static double face_normals[] = { -0.500000000000000, 0, 0, - 0.500000000000000, 0.000000000000000, 0.500000000000000, - 0, -1.500000000000000, 0, - 0, 1.500000000000000, 1.500000000000000, - 0, 0, -2.000000000000000 }; - static double cell_centroids[] = { 0.85, 0.35, 0.3 }; - static double cell_volumes[] = { 5.0/6.0 }; - } // namespace IrregPrism - - UnstructuredGrid makeIrregPrism() - { - // Make a 3d 1-cell grid, where the single cell is a prism. - UnstructuredGrid grid; - grid.dimensions = 3; - grid.number_of_cells = 1; - grid.number_of_faces = 5; - grid.number_of_nodes = 6; - grid.face_nodes = IrregPrism::face_nodes; - grid.face_nodepos = IrregPrism::face_nodepos; - grid.face_cells = IrregPrism::face_cells; - grid.cell_faces = IrregPrism::cell_faces; - grid.cell_facepos = IrregPrism::cell_facepos; - grid.node_coordinates = IrregPrism::node_coordinates; - grid.face_centroids = IrregPrism::face_centroids; - grid.face_areas = IrregPrism::face_areas; - grid.face_normals = IrregPrism::face_normals; - grid.cell_centroids = IrregPrism::cell_centroids; - grid.cell_volumes = IrregPrism::cell_volumes; - return grid; - } - -} // anonymous namespace - -template -void testConstantVelReproIrregPrism() -{ - // Set up a 3d 1-cell non-cartesian case (a pyramid). - UnstructuredGrid grid = makeIrregPrism(); - std::vector v(3); - v[0] = 0.12345; - v[1] = -0.6789; - v[2] = 0.3456; - std::vector flux; - computeFlux(grid, v, flux); - VelInterp vic(grid); - vic.setupFluxes(&flux[0]); - - // Test a few points - std::vector x(3); - x[0] = 0.123; - x[1] = 0.0123; - x[2] = 0.213; - std::vector v_interp(3); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - x[2] = 1.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 1.0; - x[1] = 0.0; - x[2] = 1.0; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.5; - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); -} - -template -void testLinearVelReproIrregPrism() -{ - // Set up a 3d 1-cell non-cartesian case (a pyramid). - UnstructuredGrid grid = makeIrregPrism(); - std::vector v0(3); - // v0[0] = 0.12345; - // v0[1] = -0.6789; - // v0[2] = 0.423; - v0[0] = 0.0; - v0[1] = 0.0; - v0[2] = 0.0; - std::vector v1(3); - // v1[0] = -0.1; - // v1[1] = 0.454; - // v1[2] = 0.21; - v1[0] = 0.0; - v1[1] = 0.0; - v1[2] = 1.0; - std::vector flux; - computeFluxLinear(grid, v0, v1, flux); - VelInterp vic(grid); - vic.setupFluxes(&flux[0]); - - // Test a few points - std::vector v(3); - std::vector x(3); - x[0] = 0.123; - x[1] = 0.0123; - x[2] = 0.213; - computeLinearVec(v0, v1, x, v); - std::vector v_interp(3); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - x[2] = 1.0; - computeLinearVec(v0, v1, x, v); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 1.0; - x[1] = 0.0; - x[2] = 1.0; - computeLinearVec(v0, v1, x, v); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.5; - computeLinearVec(v0, v1, x, v); - vic.interpolate(0, &x[0], &v_interp[0]); - BOOST_CHECK(vectorDiff2(v, v_interp) < 1e-12); -} - - -BOOST_AUTO_TEST_CASE(test_VelocityInterpolationConstant) -{ - testConstantVelRepro2d(); - testConstantVelReproPyramid(); - testConstantVelReproIrreg2d(); - testConstantVelReproIrregPrism(); -} - -BOOST_AUTO_TEST_CASE(test_VelocityInterpolationECVI) -{ - testConstantVelRepro2d(); - BOOST_CHECK_THROW(testConstantVelReproPyramid(), std::exception); - testConstantVelReproIrreg2d(); - testConstantVelReproIrregPrism(); - // Though the interpolation has linear precision, the corner velocity - // construction does not, so the below test cannot be expected to succeed. - // testLinearVelReproIrregPrism(); -} - - diff --git a/tests/test_wachspresscoord.cpp b/tests/test_wachspresscoord.cpp deleted file mode 100644 index 6b8b144f..00000000 --- a/tests/test_wachspresscoord.cpp +++ /dev/null @@ -1,363 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include - -#if HAVE_DYNAMIC_BOOST_TEST -#define BOOST_TEST_DYN_LINK -#endif -#define NVERBOSE // to suppress our messages when throwing - -#define BOOST_TEST_MODULE WachspressCoordTest -#include - -#include -#include -#include -#include - -using namespace Opm; - - -namespace -{ - class Interpolator - { - public: - explicit Interpolator(const UnstructuredGrid& grid) - : bcmethod_(grid), grid_(grid) - { - } - template - double interpolate(const Func& f, - const int cell, - const std::vector& x) const - { - const int ncor = bcmethod_.numCorners(cell); - bary_coord_.resize(ncor); - bcmethod_.cartToBary(cell, &x[0], &bary_coord_[0]); - double val = 0.0; - for (int cor = 0; cor < ncor; ++cor) { - const int vertex = bcmethod_.cornerInfo()[cell][cor].vertex; - const double vval = f(grid_.node_coordinates + grid_.dimensions*vertex); - val += vval*bary_coord_[cor]; - } - return val; - } - private: - WachspressCoord bcmethod_; - const UnstructuredGrid grid_; - mutable std::vector bary_coord_; - }; - -} // anonymous namespace - -struct LinearFunc -{ - double operator()(const double* x) const - { - return 1.0*x[0] + 2.0*x[1] + 3.0; - } -}; - -static void test2dCart() -{ - // Set up 2d 1-cell cartesian case. - GridManager g(1, 1); - const UnstructuredGrid& grid = *g.c_grid(); - Interpolator interp(grid); - LinearFunc f; - - // Test a few points - std::vector x(2); - x[0] = 0.23456; - x[1] = 0.87654; - double val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 1.0; - x[1] = 0.5; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 1.0; - x[1] = 1.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); -} - - -namespace -{ - - // Data for a pyramid. Node numbering goes - // lexicographic on bottom, then top. - // Face numbering goes xmin, xmax, ymin, ymax, bottom. - namespace Pyramid - { - static int face_nodes[] = { 0, 4, 2, 3, 4, 1, 0, 1, 4, 4, 3, 2, 0, 2, 3, 1, }; - static int face_nodepos[] = { 0, 3, 6, 9, 12, 16 }; - static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; - static int cell_faces[] = { 0, 1, 2, 3, 4 }; - static int cell_facepos[] = { 0, 5 }; - static double node_coordinates[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; - static double face_centroids[] = { 0, 1.0/3.0, 1.0/3.0, - 2.0/3.0, 1.0/3.0, 1.0/3.0, - 1.0/3.0, 0, 1.0/3.0, - 1.0/3.0, 2.0/3.0, 1.0/3.0, - 0.5, 0.5, 0 }; - static double face_areas[] = { 0.5, std::sqrt(2.0), 0.5, std::sqrt(2.0), 1.0 }; - static double face_normals[] = { -0.5000, 0, 0, - 0.5000, 0, 0.5000, - 0, -0.5000, 0, - 0, 0.5000, 0.5000, - 0, 0, -1.0000 }; - static double cell_centroids[] = { 0.375, 0.375, 0.25 }; - static double cell_volumes[] = { 1.0/3.0 }; - - } // namespace Pyramid - - UnstructuredGrid makePyramid() - { - // Make a 3d 1-cell grid, where the single cell is a pyramid. - UnstructuredGrid grid; - grid.dimensions = 3; - grid.number_of_cells = 1; - grid.number_of_faces = 5; - grid.number_of_nodes = 5; - grid.face_nodes = Pyramid::face_nodes; - grid.face_nodepos = Pyramid::face_nodepos; - grid.face_cells = Pyramid::face_cells; - grid.cell_faces = Pyramid::cell_faces; - grid.cell_facepos = Pyramid::cell_facepos; - grid.node_coordinates = Pyramid::node_coordinates; - grid.face_centroids = Pyramid::face_centroids; - grid.face_areas = Pyramid::face_areas; - grid.face_normals = Pyramid::face_normals; - grid.cell_centroids = Pyramid::cell_centroids; - grid.cell_volumes = Pyramid::cell_volumes; - return grid; - } - -} // anonymous namespace - - -static void testPyramid() -{ - // Set up a 3d 1-cell non-cartesian case (a pyramid). - UnstructuredGrid grid = makePyramid(); - Interpolator interp(grid); - LinearFunc f; - - // Test a few points - std::vector x(3); - x[0] = 0.123; - x[1] = 0.0123; - x[2] = 0.213; - double val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - x[2] = 1.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.1; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); -} - - -namespace -{ - // Data for an irregular 2d polygon. - namespace Irreg2d - { - static int face_nodes[] = { 0, 1, 1, 2, 2, 3, 3, 4, 4, 0 }; - static int face_nodepos[] = { 0, 2, 4, 6, 8, 10 }; - static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; - static int cell_faces[] = { 0, 1, 2, 3, 4 }; - static int cell_facepos[] = { 0, 5 }; - static double node_coordinates[] = { 0, 0, 3, 0, 3, 2, 1, 3, 0, 2 }; - static double face_centroids[] = { 1.5, 0, 3, 1, 2, 2.5, 0.5, 2.5, 0, 1 }; - static double face_areas[] = { 3, 2, std::sqrt(5.0), std::sqrt(2.0), 2 }; - static double face_normals[] = { 0, -3, 2, 0, 1, 2, -1, 1, -2, 0 }; - static double cell_centroids[] = { 22.0/15.0, 19.0/15.0 }; - static double cell_volumes[] = { 7.5 }; - } // namespace Irreg2d - - UnstructuredGrid makeIrreg2d() - { - // Make a 2d 1-cell grid, where the single cell is a polygon. - UnstructuredGrid grid; - grid.dimensions = 2; - grid.number_of_cells = 1; - grid.number_of_faces = 5; - grid.number_of_nodes = 5; - grid.face_nodes = Irreg2d::face_nodes; - grid.face_nodepos = Irreg2d::face_nodepos; - grid.face_cells = Irreg2d::face_cells; - grid.cell_faces = Irreg2d::cell_faces; - grid.cell_facepos = Irreg2d::cell_facepos; - grid.node_coordinates = Irreg2d::node_coordinates; - grid.face_centroids = Irreg2d::face_centroids; - grid.face_areas = Irreg2d::face_areas; - grid.face_normals = Irreg2d::face_normals; - grid.cell_centroids = Irreg2d::cell_centroids; - grid.cell_volumes = Irreg2d::cell_volumes; - return grid; - } - -} // anonymous namespace - - -static void testIrreg2d() -{ - // Set up a 2d 1-cell where the single cell is a polygon. - UnstructuredGrid grid = makeIrreg2d(); - Interpolator interp(grid); - LinearFunc f; - - // Test a few points - std::vector x(2); - x[0] = 1.2345; - x[1] = 2.0123; - double val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 1.0; - x[1] = 3.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 3.0; - x[1] = 1.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); -} - - -namespace -{ - // Data for an irregular 3d prism. - namespace IrregPrism - { - static int face_nodes[] = { 0, 4, 2, 1, 3, 5, 0, 1, 5, 4, 2, 4, 5, 3, 2, 3, 0, 1}; - static int face_nodepos[] = { 0, 3, 6, 10, 14, 18 }; - static int face_cells[] = { 0, -1, 0, -1, 0, -1, 0, -1, 0, -1 }; - static int cell_faces[] = { 0, 1, 2, 3, 4 }; - static int cell_facepos[] = { 0, 5 }; - static double node_coordinates[] = { 0, 0, 0, - 2, 0, 0, - 0, 1, 0, - 2, 1, 0, - 0, 0, 1, - 1, 0, 1 }; - static double face_centroids[] = { 0, 1.0/3.0, 1.0/3.0, - 5.0/3.0, 1.0/3.0, 1.0/3.0, - 7.0/9.0, 0, 4.0/9.0, - 7.0/9.0, 5.0/9.0, 4.0/9.0, - 1, 0.5, 0 }; - static double face_areas[] = { 0.500000000000000, - 0.707106781186548, - 1.500000000000000, - 2.121320343559642, - 2.000000000000000 }; - static double face_normals[] = { -0.500000000000000, 0, 0, - 0.500000000000000, 0.000000000000000, 0.500000000000000, - 0, -1.500000000000000, 0, - 0, 1.500000000000000, 1.500000000000000, - 0, 0, -2.000000000000000 }; - static double cell_centroids[] = { 0.85, 0.35, 0.3 }; - static double cell_volumes[] = { 5.0/6.0 }; - } // namespace IrregPrism - - UnstructuredGrid makeIrregPrism() - { - // Make a 3d 1-cell grid, where the single cell is a prism. - UnstructuredGrid grid; - grid.dimensions = 3; - grid.number_of_cells = 1; - grid.number_of_faces = 5; - grid.number_of_nodes = 6; - grid.face_nodes = IrregPrism::face_nodes; - grid.face_nodepos = IrregPrism::face_nodepos; - grid.face_cells = IrregPrism::face_cells; - grid.cell_faces = IrregPrism::cell_faces; - grid.cell_facepos = IrregPrism::cell_facepos; - grid.node_coordinates = IrregPrism::node_coordinates; - grid.face_centroids = IrregPrism::face_centroids; - grid.face_areas = IrregPrism::face_areas; - grid.face_normals = IrregPrism::face_normals; - grid.cell_centroids = IrregPrism::cell_centroids; - grid.cell_volumes = IrregPrism::cell_volumes; - return grid; - } - -} // anonymous namespace - -static void testIrregPrism() -{ - // Set up a 3d 1-cell non-cartesian case (a prism). - UnstructuredGrid grid = makeIrregPrism(); - Interpolator interp(grid); - LinearFunc f; - - // Test a few points - std::vector x(3); - x[0] = 0.123; - x[1] = 0.0123; - x[2] = 0.213; - double val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.0; - x[1] = 0.0; - x[2] = 1.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 1.0; - x[1] = 0.0; - x[2] = 1.0; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); - x[0] = 0.5; - x[1] = 0.5; - x[2] = 0.5; - val = interp.interpolate(f, 0, x); - BOOST_CHECK(std::fabs(val - f(&x[0])) < 1e-12); -} - - -BOOST_AUTO_TEST_CASE(test_WachspressCoord) -{ - test2dCart(); - BOOST_CHECK_THROW(testPyramid(), std::exception); - testIrreg2d(); - testIrregPrism(); -}