remove some grid related utilities
moved to opm-grid
This commit is contained in:
Arne Morten Kvarving
parent
de2c3817bc
commit
8240cca69e
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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 <CODE>(c,f)</CODE> 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
|
||||
* <CODE>G->cell_facepos[ G->number_of_cells ]</CODE>.
|
||||
*/
|
||||
template<class Grid>
|
||||
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 <CODE>numFaces(G)</CODE>.
|
||||
*/
|
||||
template<class Grid>
|
||||
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 <CODE>numFaces(G)</CODE>.
|
||||
*/
|
||||
template<class Grid>
|
||||
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 */
|
||||
|
|
@ -1,177 +0,0 @@
|
|||
#include <opm/core/linalg/blas_lapack.h>
|
||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||
#include <opm/core/grid/GridHelpers.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
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<d;++i)
|
||||
out[i]=in[i]*area;
|
||||
return out;
|
||||
}
|
||||
|
||||
inline void maybeFreeFaceNormal(const Dune::CpGrid&, const double* array)
|
||||
{
|
||||
delete[] array;
|
||||
}
|
||||
#endif // HAVE_OPM_GRID
|
||||
|
||||
inline const double* multiplyFaceNormalWithArea(const UnstructuredGrid&, int, const double* in)
|
||||
{
|
||||
return in;
|
||||
}
|
||||
|
||||
inline void maybeFreeFaceNormal(const UnstructuredGrid&, const double*)
|
||||
{}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
template<class Grid>
|
||||
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<Grid>::IteratorType cc = beginCellCentroids(*G);
|
||||
typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
|
||||
typename FaceCellTraits<Grid>::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<Grid>::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<class Grid>
|
||||
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<Grid>::Type c2f = cell2Faces(*G);
|
||||
|
||||
for (int c = 0, i = 0; c < numCells(*G); c++) {
|
||||
typedef typename Cell2FacesTraits<Grid>::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<class Grid>
|
||||
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<Grid>::Type c2f = cell2Faces(*G);
|
||||
|
||||
for (int c = 0, i = 0; c < numCells(*G); c++) {
|
||||
typedef typename Cell2FacesTraits<Grid>::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];
|
||||
}
|
||||
}
|
||||
|
|
@ -1,130 +0,0 @@
|
|||
#include "config.h"
|
||||
#include <assert.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/linalg/blas_lapack.h>
|
||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||
|
||||
|
||||
#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<UnstructuredGrid>(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<UnstructuredGrid>(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<UnstructuredGrid>(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];
|
||||
}
|
||||
}
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_TRANS_TPFA_HEADER_INCLUDED
|
||||
#define OPM_TRANS_TPFA_HEADER_INCLUDED
|
||||
|
||||
/**
|
||||
* \file
|
||||
* Routines to assist in the calculation of two-point transmissibilities.
|
||||
*/
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
#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 <CODE>(c,f)</CODE> 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
|
||||
* <CODE>G->cell_facepos[ G->number_of_cells ]</CODE>.
|
||||
*/
|
||||
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 <CODE>G->number_of_faces</CODE>.
|
||||
*/
|
||||
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 <CODE>G->number_of_faces</CODE>.
|
||||
*/
|
||||
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 */
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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 <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
|
||||
|
||||
#include <cassert>
|
||||
#include <cstddef>
|
||||
#include <memory>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
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 <typename T>
|
||||
struct HasProperty;
|
||||
|
||||
/**
|
||||
* Property value retrieval.
|
||||
*
|
||||
* Supported for types \c int and \c double.
|
||||
*
|
||||
* \tparam T Property element type.
|
||||
*/
|
||||
template <typename T>
|
||||
struct GetProperty;
|
||||
|
||||
/**
|
||||
* Specialization of property existence predicate for
|
||||
* type \c int.
|
||||
*/
|
||||
template <>
|
||||
struct HasProperty<int> {
|
||||
/**
|
||||
* 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 <class PropertyContainer>
|
||||
static bool
|
||||
p(PropertyContainer& ecl,
|
||||
const std::string& kw);
|
||||
};
|
||||
|
||||
template <class PropertyContainer>
|
||||
bool
|
||||
HasProperty<int>::p(PropertyContainer& ecl,
|
||||
const std::string& kw)
|
||||
{
|
||||
return ecl.get3DProperties().hasDeckIntGridProperty(kw);
|
||||
}
|
||||
|
||||
/**
|
||||
* Specialization of property existence predicate for
|
||||
* type \c double.
|
||||
*/
|
||||
template <>
|
||||
struct HasProperty<double> {
|
||||
/**
|
||||
* 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 <class PropertyContainer>
|
||||
static bool
|
||||
p(PropertyContainer& ecl,
|
||||
const std::string& kw);
|
||||
};
|
||||
|
||||
template <class PropertyContainer>
|
||||
bool
|
||||
HasProperty<double>::p(PropertyContainer& ecl,
|
||||
const std::string& kw)
|
||||
{
|
||||
return ecl.get3DProperties().hasDeckDoubleGridProperty(kw);
|
||||
}
|
||||
|
||||
/**
|
||||
* Specialization of property value retrieval for type
|
||||
* \c int.
|
||||
*/
|
||||
template <>
|
||||
struct GetProperty<int> {
|
||||
/**
|
||||
* 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 <class PropertyContainer>
|
||||
static const GridProperty<int>*
|
||||
value(PropertyContainer& ecl,
|
||||
const std::string& kw);
|
||||
};
|
||||
|
||||
template <class PropertyContainer>
|
||||
const GridProperty<int>*
|
||||
GetProperty<int>::value(PropertyContainer& ecl,
|
||||
const std::string& kw)
|
||||
{
|
||||
assert (HasProperty<int>::p(ecl, kw));
|
||||
|
||||
return &ecl.get3DProperties().getIntGridProperty(kw);
|
||||
}
|
||||
|
||||
/**
|
||||
* Specialization of property value retrieval for type
|
||||
* \c double.
|
||||
*/
|
||||
template <>
|
||||
struct GetProperty<double> {
|
||||
/**
|
||||
* 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 <class PropertyContainer>
|
||||
static const GridProperty<double>*
|
||||
value(PropertyContainer& ecl,
|
||||
const std::string& kw);
|
||||
};
|
||||
|
||||
template <class PropertyContainer>
|
||||
const GridProperty<double>*
|
||||
GetProperty<double>::value(PropertyContainer& ecl,
|
||||
const std::string& kw)
|
||||
{
|
||||
assert (HasProperty<double>::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 <typename T>
|
||||
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 <class PropertyContainer>
|
||||
static const GridProperty<T>*
|
||||
value(PropertyContainer& ecl,
|
||||
const std::string& kw);
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
template <class PropertyContainer>
|
||||
const GridProperty<T>*
|
||||
EclipsePropertyArray<T>::value(PropertyContainer& ecl,
|
||||
const std::string& kw)
|
||||
{
|
||||
if (! EclPropImpl::HasProperty<T>::p(ecl, kw)) {
|
||||
return nullptr;
|
||||
}
|
||||
return EclPropImpl::GetProperty<T>::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 <typename T>
|
||||
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 <class PropertyContainer>
|
||||
ExtractFromDeck(PropertyContainer& ecl,
|
||||
const std::string& kw,
|
||||
const T dflt)
|
||||
: x_ (Details::EclipsePropertyArray<T>::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<T>* 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 <typename T>
|
||||
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<void>(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 DataArray, class PropertyTag = Tag::Any>
|
||||
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 */
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED
|
||||
#define OPM_REGIONMAPPING_HEADER_INCLUDED
|
||||
|
||||
#include <boost/range.hpp>
|
||||
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
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<int> >
|
||||
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<CellId>::const_iterator CellIter;
|
||||
|
||||
typedef boost::iterator_range<CellIter> 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<RegionId>&
|
||||
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<CellId>::size_type Pos;
|
||||
|
||||
std::unordered_map<RegionId, Pos> binid;
|
||||
std::vector<RegionId> active;
|
||||
|
||||
std::vector<Pos> p; /**< Region start pointers */
|
||||
std::vector<CellId> 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<Pos>(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
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/utility/VelocityInterpolation.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/linalg/blas_lapack.h>
|
||||
|
||||
#include <iostream>
|
||||
|
||||
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<double> N(dim*dim); // Normals matrix. Fortran ordering!
|
||||
std::vector<double> orig_N(dim*dim); // Normals matrix. Fortran ordering!
|
||||
std::vector<double> f(dim); // Flux vector.
|
||||
std::vector<double> orig_f(dim); // Flux vector.
|
||||
std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve
|
||||
const SparseTable<WachspressCoord::CornerInfo>& all_ci = bcmethod_.cornerInfo();
|
||||
const std::vector<int>& 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<WachspressCoord::CornerInfo>& 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
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_VELOCITYINTERPOLATION_HEADER_INCLUDED
|
||||
#define OPM_VELOCITYINTERPOLATION_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/utility/WachspressCoord.hpp>
|
||||
#include <vector>
|
||||
|
||||
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<double> bary_coord_;
|
||||
std::vector<double> corner_velocity_; // size = dim * #corners
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_VELOCITYINTERPOLATION_HEADER_INCLUDED
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <opm/core/utility/WachspressCoord.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
#include <cmath>
|
||||
#include <map>
|
||||
#include <set>
|
||||
|
||||
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<int> cell_vertices;
|
||||
std::vector<int> cell_faces;
|
||||
std::multimap<int, int> 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<CornerInfo> cell_corner_info;
|
||||
std::set<int>::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<int, int>::const_iterator MMIt;
|
||||
std::pair<MMIt, MMIt> frange = vertex_adj_faces.equal_range(ci.vertex);
|
||||
int fi = 0;
|
||||
std::vector<int> 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<int> 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>& WachspressCoord::cornerInfo() const
|
||||
{
|
||||
return corner_info_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/// The class stores some info for each corner.
|
||||
/// \return The corner info container.
|
||||
const std::vector<int>& 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
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_WACHSPRESSCOORD_HEADER_INCLUDED
|
||||
#define OPM_WACHSPRESSCOORD_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/utility/SparseTable.hpp>
|
||||
#include <vector>
|
||||
|
||||
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>& 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<int>& adjacentFaces() const;
|
||||
|
||||
private:
|
||||
const UnstructuredGrid& grid_;
|
||||
SparseTable<CornerInfo> corner_info_; // Corner info by cell.
|
||||
std::vector<int> adj_faces_; // Set of adjacent faces, by corner id. Contains dim face indices per corner.
|
||||
SparseTable<int> nonadj_faces_; // Set of nonadjacent faces, by corner id.
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_WACHSPRESSCOORD_HEADER_INCLUDED
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
|
||||
#include <opm/core/utility/compressedToCartesian.hpp>
|
||||
#include <numeric>
|
||||
|
||||
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<int> compressedToCartesian(const int num_cells,
|
||||
const int* global_cell)
|
||||
{
|
||||
std::vector<int> 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
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_COMPRESSEDTOCARTESIAN_HEADER_INCLUDED
|
||||
#define OPM_COMPRESSEDTOCARTESIAN_HEADER_INCLUDED
|
||||
|
||||
#include <vector>
|
||||
|
||||
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<int> compressedToCartesian(const int num_cells,
|
||||
const int* global_cell);
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_COMPRESSEDTOCARTESIAN_HEADER_INCLUDED
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#include <config.h>
|
||||
#include <cassert>
|
||||
|
||||
#include "extractPvtTableIndex.hpp"
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
|
||||
const Opm::EclipseState& eclState,
|
||||
size_t numCompressed,
|
||||
const int *compressedToCartesianCellIdx)
|
||||
{
|
||||
//Get the PVTNUM data
|
||||
const std::vector<int>& 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;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef OPM_EXTRACT_PVT_TABLE_INDEX_HPP
|
||||
#define OPM_EXTRACT_PVT_TABLE_INDEX_HPP
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
|
||||
const Opm::EclipseState& eclState,
|
||||
size_t numCompressed,
|
||||
const int *compressedToCartesianCellIdx);
|
||||
|
||||
}
|
||||
|
||||
#endif // OPM_EXTRACT_PVT_TABLE_INDEX_HPP
|
||||
|
|
@ -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
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
|
||||
#if HAVE_DYNAMIC_BOOST_TEST
|
||||
#define BOOST_TEST_DYN_LINK
|
||||
#endif
|
||||
|
||||
#define NVERBOSE
|
||||
|
||||
#define BOOST_TEST_MODULE CompressedPropertyTest
|
||||
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <boost/test/floating_point_comparison.hpp>
|
||||
|
||||
#include <opm/core/utility/CompressedPropertyAccess.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Parser/Parser.hpp>
|
||||
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
struct SetupSimple {
|
||||
SetupSimple() :
|
||||
deck( Opm::Parser{}.parseFile( "compressed_gridproperty.data", Opm::ParseContext{} ) ),
|
||||
ecl( deck )
|
||||
{}
|
||||
|
||||
Opm::Deck deck;
|
||||
Opm::EclipseState ecl;
|
||||
};
|
||||
|
||||
|
||||
template <class Setup>
|
||||
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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<double> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<double> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<int> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<int> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::Constant<double> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::Constant<int> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<double> ECLGlobalDoubleArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ECLGlobalDoubleArray> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<double> ECLGlobalDoubleArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::Tag::NTG NTG;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ECLGlobalDoubleArray, NTG> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<int> ECLGlobalIntArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ECLGlobalIntArray> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::ExtractFromDeck<int> ECLGlobalIntArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ECLGlobalIntArray, RegionID> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::Constant<double> ConstantDoubleArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ConstantDoubleArray> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::Constant<double> ConstantDoubleArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ConstantDoubleArray, MyTag> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::Constant<int> ConstantIntArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ConstantIntArray> 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<SetupSimple>)
|
||||
{
|
||||
typedef Opm::GridPropertyAccess::ArrayPolicy
|
||||
::Constant<int> ConstantIntArray;
|
||||
|
||||
typedef Opm::GridPropertyAccess::
|
||||
Compressed<ConstantIntArray, MyTag> 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()
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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 <opm/common/utility/platform_dependent/disable_warnings.h>
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
||||
|
||||
/* --- our own headers --- */
|
||||
|
||||
#include <opm/core/utility/RegionMapping.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <map>
|
||||
|
||||
BOOST_AUTO_TEST_SUITE ()
|
||||
|
||||
|
||||
BOOST_AUTO_TEST_CASE (Forward)
|
||||
{
|
||||
std::vector<int> 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<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
|
||||
|
||||
Opm::RegionMapping<> rm(regions);
|
||||
|
||||
std::vector<int> 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<int, std::vector<int>>;
|
||||
|
||||
// 0 1 2 3 4 5 6 7 8
|
||||
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
|
||||
|
||||
Opm::RegionMapping<> rm(regions);
|
||||
|
||||
std::vector<int> 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<int, std::vector<int>>;
|
||||
|
||||
// 0 1 2 3 4 5 6 7 8
|
||||
std::vector<int> regions = { 2, 4, 2, 4, 2, 7, 6, 3, 6 };
|
||||
|
||||
Opm::RegionMapping<> rm(regions);
|
||||
|
||||
std::vector<int> 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()
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
|
||||
#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 <boost/test/unit_test.hpp>
|
||||
|
||||
#include <opm/core/utility/VelocityInterpolation.hpp>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
#include <cmath>
|
||||
|
||||
using namespace Opm;
|
||||
|
||||
namespace
|
||||
{
|
||||
|
||||
// Compute flux corresponding to a constant velocity vector v.
|
||||
void computeFlux(const UnstructuredGrid& grid, const std::vector<double>& v, std::vector<double>& 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<double>& v0,
|
||||
const std::vector<double>& v1,
|
||||
const std::vector<double>& x,
|
||||
std::vector<double>& 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<double>& v0,
|
||||
const std::vector<double>& v1,
|
||||
std::vector<double>& flux)
|
||||
{
|
||||
const int dim = v0.size();
|
||||
assert(dim == grid.dimensions);
|
||||
flux.resize(grid.number_of_faces);
|
||||
std::vector<double> x(dim);
|
||||
std::vector<double> 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<double>& v1, const std::vector<double>& 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 <class VelInterp>
|
||||
void testConstantVelRepro2d()
|
||||
{
|
||||
// Set up 2d 1-cell cartesian case.
|
||||
GridManager g(1, 1);
|
||||
const UnstructuredGrid& grid = *g.c_grid();
|
||||
std::vector<double> v(2);
|
||||
v[0] = 0.12345;
|
||||
v[1] = -0.6789;
|
||||
std::vector<double> flux;
|
||||
computeFlux(grid, v, flux);
|
||||
VelInterp vic(grid);
|
||||
vic.setupFluxes(&flux[0]);
|
||||
|
||||
// Test a few points
|
||||
std::vector<double> x(2);
|
||||
x[0] = 0.23456;
|
||||
x[1] = 0.87654;
|
||||
std::vector<double> 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 <class VelInterp>
|
||||
void testConstantVelReproPyramid()
|
||||
{
|
||||
// Set up a 3d 1-cell non-cartesian case (a pyramid).
|
||||
UnstructuredGrid grid = makePyramid();
|
||||
std::vector<double> v(3);
|
||||
v[0] = 0.12345;
|
||||
v[1] = -0.6789;
|
||||
v[2] = 0.3456;
|
||||
std::vector<double> flux;
|
||||
computeFlux(grid, v, flux);
|
||||
VelInterp vic(grid);
|
||||
vic.setupFluxes(&flux[0]);
|
||||
|
||||
// Test a few points
|
||||
std::vector<double> x(3);
|
||||
x[0] = 0.123;
|
||||
x[1] = 0.0123;
|
||||
x[2] = 0.213;
|
||||
std::vector<double> 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 <class VelInterp>
|
||||
void testConstantVelReproIrreg2d()
|
||||
{
|
||||
// Set up a 2d 1-cell non-cartesian case (a pyramid).
|
||||
UnstructuredGrid grid = makeIrreg2d();
|
||||
std::vector<double> v(2);
|
||||
v[0] = 0.12345;
|
||||
v[1] = -0.6789;
|
||||
std::vector<double> flux;
|
||||
computeFlux(grid, v, flux);
|
||||
VelInterp vic(grid);
|
||||
vic.setupFluxes(&flux[0]);
|
||||
|
||||
// Test a few points
|
||||
std::vector<double> x(2);
|
||||
x[0] = 1.2345;
|
||||
x[1] = 2.0123;
|
||||
std::vector<double> 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 <class VelInterp>
|
||||
void testConstantVelReproIrregPrism()
|
||||
{
|
||||
// Set up a 3d 1-cell non-cartesian case (a pyramid).
|
||||
UnstructuredGrid grid = makeIrregPrism();
|
||||
std::vector<double> v(3);
|
||||
v[0] = 0.12345;
|
||||
v[1] = -0.6789;
|
||||
v[2] = 0.3456;
|
||||
std::vector<double> flux;
|
||||
computeFlux(grid, v, flux);
|
||||
VelInterp vic(grid);
|
||||
vic.setupFluxes(&flux[0]);
|
||||
|
||||
// Test a few points
|
||||
std::vector<double> x(3);
|
||||
x[0] = 0.123;
|
||||
x[1] = 0.0123;
|
||||
x[2] = 0.213;
|
||||
std::vector<double> 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 <class VelInterp>
|
||||
void testLinearVelReproIrregPrism()
|
||||
{
|
||||
// Set up a 3d 1-cell non-cartesian case (a pyramid).
|
||||
UnstructuredGrid grid = makeIrregPrism();
|
||||
std::vector<double> 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<double> 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<double> flux;
|
||||
computeFluxLinear(grid, v0, v1, flux);
|
||||
VelInterp vic(grid);
|
||||
vic.setupFluxes(&flux[0]);
|
||||
|
||||
// Test a few points
|
||||
std::vector<double> v(3);
|
||||
std::vector<double> x(3);
|
||||
x[0] = 0.123;
|
||||
x[1] = 0.0123;
|
||||
x[2] = 0.213;
|
||||
computeLinearVec(v0, v1, x, v);
|
||||
std::vector<double> 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<VelocityInterpolationConstant>();
|
||||
testConstantVelReproPyramid<VelocityInterpolationConstant>();
|
||||
testConstantVelReproIrreg2d<VelocityInterpolationConstant>();
|
||||
testConstantVelReproIrregPrism<VelocityInterpolationConstant>();
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE(test_VelocityInterpolationECVI)
|
||||
{
|
||||
testConstantVelRepro2d<VelocityInterpolationECVI>();
|
||||
BOOST_CHECK_THROW(testConstantVelReproPyramid<VelocityInterpolationECVI>(), std::exception);
|
||||
testConstantVelReproIrreg2d<VelocityInterpolationECVI>();
|
||||
testConstantVelReproIrregPrism<VelocityInterpolationECVI>();
|
||||
// Though the interpolation has linear precision, the corner velocity
|
||||
// construction does not, so the below test cannot be expected to succeed.
|
||||
// testLinearVelReproIrregPrism<VelocityInterpolationECVI>();
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
|
||||
#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 <boost/test/unit_test.hpp>
|
||||
|
||||
#include <opm/core/utility/WachspressCoord.hpp>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
#include <cmath>
|
||||
|
||||
using namespace Opm;
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
class Interpolator
|
||||
{
|
||||
public:
|
||||
explicit Interpolator(const UnstructuredGrid& grid)
|
||||
: bcmethod_(grid), grid_(grid)
|
||||
{
|
||||
}
|
||||
template <class Func>
|
||||
double interpolate(const Func& f,
|
||||
const int cell,
|
||||
const std::vector<double>& 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<double> 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<double> 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<double> 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<double> 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<double> 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();
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user