OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added class WachspressCoord.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2012-10-29 14:55:44 +01:00
parent 3288c7d505
commit 1b31b207fe
5 changed files with 678 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@@ -122,6 +122,7 @@ opm/core/transport/transport_source.c \
opm/core/utility/MonotCubicInterpolator.cpp \
opm/core/utility/StopWatch.cpp \
opm/core/utility/VelocityInterpolation.cpp \
opm/core/utility/WachspressCoord.cpp \
opm/core/utility/miscUtilities.cpp \
opm/core/utility/miscUtilitiesBlackoil.cpp \
opm/core/utility/parameters/Parameter.cpp \
@@ -258,6 +259,7 @@ opm/core/utility/StopWatch.hpp \
opm/core/utility/UniformTableLinear.hpp \
opm/core/utility/Units.hpp \
opm/core/utility/VelocityInterpolation.hpp \
opm/core/utility/WachspressCoord.hpp \
opm/core/utility/buildUniformMonotoneTable.hpp \
opm/core/utility/initState.hpp \
opm/core/utility/initState_impl.hpp \

228
opm/core/utility/WachspressCoord.cpp Normal file
View File

@@ -0,0 +1,228 @@
/*
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 <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) {
THROW("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) {
THROW("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_;
}
/// 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

81
opm/core/utility/WachspressCoord.hpp Normal file
View File

@@ -0,0 +1,81 @@
/*
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.
/// \return The corner info container.
const SparseTable<CornerInfo>& cornerInfo() 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

4
tests/Makefile.am
View File

@@ -25,6 +25,7 @@ test_read_vag \
test_readpolymer \
test_sf2p \
test_velocityinterpolation \
test_wachspresscoord \
test_wells \
test_writeVtkData \
unit_test
@@ -63,6 +64,9 @@ test_sf2p_SOURCES = test_sf2p.cpp
test_velocityinterpolation_SOURCES = test_velocityinterpolation.cpp
test_velocityinterpolation_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)
test_wachspresscoord_SOURCES = test_wachspresscoord.cpp
test_wachspresscoord_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)
test_writeVtkData_SOURCES = test_writeVtkData.cpp
test_wells_SOURCES = test_wells.cpp

363
tests/test_wachspresscoord.cpp Normal file
View File

@@ -0,0 +1,363 @@
/*
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/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();
}