Implemented 2d quadrature rules for degree 1 and 2 (order 2 and 3).
Also added cartesian 2d cell testcase.
This commit is contained in:
Atgeirr Flø Rasmussen
parent
91f9b3c2de
commit
ff1bfc810a
@ -54,6 +54,18 @@ namespace Opm
|
||||
double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
|
||||
return std::fabs(determinantOf(a, b, c) / 6.0);
|
||||
}
|
||||
|
||||
/// Calculates the area of a triangle consisting of 3 vertices
|
||||
/// with 2-dimensional coordinates
|
||||
inline double triangleArea2d(const double* p0,
|
||||
const double* p1,
|
||||
const double* p2)
|
||||
{
|
||||
double a[2] = { p1[0] - p0[0], p1[1] - p0[1] };
|
||||
double b[2] = { p2[0] - p0[0], p2[1] - p0[1] };
|
||||
double a_cross_b = a[0]*b[1] - a[1]*b[0];
|
||||
return 0.5*std::fabs(a_cross_b);
|
||||
}
|
||||
} // anonymous namespace
|
||||
|
||||
|
||||
@ -61,10 +73,35 @@ namespace Opm
|
||||
/// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
|
||||
/// Note that this class does multiply weights by cell volume,
|
||||
/// so weights always sum to cell volume.
|
||||
///
|
||||
/// Degree 1 method:
|
||||
/// Midpoint (centroid) method.
|
||||
/// n = 1, w_0 = cell volume, x_0 = cell centroid
|
||||
/// Degree 2 method:
|
||||
///
|
||||
/// Degree 2 method for 2d (but see the note):
|
||||
/// Based on subdivision of the cell into triangles,
|
||||
/// with the centroid as a common vertex, and the triangle
|
||||
/// edge midpoint rule.
|
||||
/// Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
|
||||
/// Its area is A_i.
|
||||
/// n = 2 * nn (nn = num nodes in face)
|
||||
/// For i = 0..(nn-1):
|
||||
/// w_i = 1/3 A_i.
|
||||
/// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
|
||||
/// x_i = (N_i + N_{i+1})/2
|
||||
/// x_{nn+i} = (C + N_i)/2
|
||||
/// All N and A indices are interpreted cyclic, modulus nn.
|
||||
/// Note: for simplicity of implementation, we currently use
|
||||
/// n = 3 * nn
|
||||
/// For i = 0..(nn-1):
|
||||
/// w_{3*i + {0,1,2}} = 1/3 A_i
|
||||
/// x_{3*i} = (N_i + N_{i+1})/2
|
||||
/// x_{3*i + {1,2}} = (C + N_{i,i+1})/2
|
||||
/// This is simpler, because we can implement it easily
|
||||
/// based on iteration over faces without requiring any
|
||||
/// particular (cyclic) ordering.
|
||||
///
|
||||
/// Degree 2 method for 3d:
|
||||
/// Based on subdivision of each cell face into triangles
|
||||
/// with the face centroid as a common vertex, and then
|
||||
/// subdividing the cell into tetrahedra with the cell
|
||||
@ -87,8 +124,8 @@ namespace Opm
|
||||
const int degree)
|
||||
: grid_(grid), cell_(cell), degree_(degree)
|
||||
{
|
||||
if (grid.dimensions != 3) {
|
||||
THROW("CellQuadrature only implemented for 3D case.");
|
||||
if (grid.dimensions > 3) {
|
||||
THROW("CellQuadrature only implemented for up to 3 dimensions.");
|
||||
}
|
||||
if (degree > 2) {
|
||||
THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
|
||||
@ -97,10 +134,14 @@ namespace Opm
|
||||
|
||||
int numQuadPts() const
|
||||
{
|
||||
if (degree_ < 2) {
|
||||
if (degree_ < 2 || grid_.dimensions == 1) {
|
||||
return 1;
|
||||
}
|
||||
// Degree 2 case.
|
||||
if (grid_.dimensions == 2) {
|
||||
return 3*(grid_.cell_facepos[cell_ + 1] - grid_.cell_facepos[cell_]);
|
||||
}
|
||||
ASSERT(grid_.dimensions == 3);
|
||||
int sumnodes = 0;
|
||||
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
|
||||
const int face = grid_.cell_faces[hf];
|
||||
@ -118,6 +159,29 @@ namespace Opm
|
||||
return;
|
||||
}
|
||||
// Degree 2 case.
|
||||
if (dim == 2) {
|
||||
if (index % 3 == 0) {
|
||||
// Boundary midpoint. This is the face centroid.
|
||||
const int hface = grid_.cell_facepos[cell_] + index/3;
|
||||
const int face = grid_.cell_faces[hface];
|
||||
const double* fc = grid_.face_centroids + dim*face;
|
||||
std::copy(fc, fc + dim, coord);
|
||||
} else {
|
||||
// Interiour midpoint. This is the average of the
|
||||
// cell centroid and a face node (they should
|
||||
// always have two nodes in 2d).
|
||||
const int hface = grid_.cell_facepos[cell_] + index/3;
|
||||
const int face = grid_.cell_faces[hface];
|
||||
const int nodeoff = (index % 3) - 1; // == 0 or 1
|
||||
const int node = grid_.face_nodes[grid_.face_nodepos[face] + nodeoff];
|
||||
const double* nc = grid_.node_coordinates + dim*node;
|
||||
for (int dd = 0; dd < dim; ++dd) {
|
||||
coord[dd] = 0.5*(nc[dd] + cc[dd]);
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
ASSERT(dim == 3);
|
||||
int tetindex = index / 4;
|
||||
const int subindex = index % 4;
|
||||
const double* nc = grid_.node_coordinates;
|
||||
@ -155,6 +219,15 @@ namespace Opm
|
||||
// Degree 2 case.
|
||||
const int dim = grid_.dimensions;
|
||||
const double* cc = grid_.cell_centroids + dim*cell_;
|
||||
if (dim == 2) {
|
||||
const int hface = grid_.cell_facepos[cell_] + index/3;
|
||||
const int face = grid_.cell_faces[hface];
|
||||
const int* nptr = grid_.face_nodes + grid_.face_nodepos[face];
|
||||
const double* nc0 = grid_.node_coordinates + dim*nptr[0];
|
||||
const double* nc1 = grid_.node_coordinates + dim*nptr[1];
|
||||
return triangleArea2d(nc0, nc1, cc)/3.0;
|
||||
}
|
||||
ASSERT(dim == 3);
|
||||
int tetindex = index / 4;
|
||||
const double* nc = grid_.node_coordinates;
|
||||
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
|
||||
|
||||
@ -57,10 +57,15 @@ namespace Opm
|
||||
/// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
|
||||
/// Note that this class does multiply weights by face area,
|
||||
/// so weights always sum to face area.
|
||||
///
|
||||
/// Degree 1 method:
|
||||
/// Midpoint (centroid) method.
|
||||
/// n = 1, w_0 = face area, x_0 = face centroid
|
||||
/// Degree 2 method:
|
||||
///
|
||||
/// Degree 2 method for 2d:
|
||||
/// Simpson's method (actually this is degree 3).
|
||||
///
|
||||
/// Degree 2 method for 3d:
|
||||
/// Based on subdivision of the face into triangles,
|
||||
/// with the centroid as a common vertex, and the triangle
|
||||
/// edge midpoint rule.
|
||||
@ -81,8 +86,8 @@ namespace Opm
|
||||
const int degree)
|
||||
: grid_(grid), face_(face), degree_(degree)
|
||||
{
|
||||
if (grid_.dimensions != 3) {
|
||||
THROW("FaceQuadrature only implemented for 3D case.");
|
||||
if (grid_.dimensions > 3) {
|
||||
THROW("FaceQuadrature only implemented for up to 3 dimensions.");
|
||||
}
|
||||
if (degree_ > 2) {
|
||||
THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
|
||||
@ -91,18 +96,22 @@ namespace Opm
|
||||
|
||||
int numQuadPts() const
|
||||
{
|
||||
if (degree_ < 2) {
|
||||
if (degree_ < 2 || grid_.dimensions < 2) {
|
||||
return 1;
|
||||
}
|
||||
// Degree 2 case.
|
||||
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
|
||||
if (grid_.dimensions == 2) {
|
||||
return 3;
|
||||
} else {
|
||||
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
|
||||
}
|
||||
}
|
||||
|
||||
void quadPtCoord(const int index, double* coord) const
|
||||
{
|
||||
const int dim = grid_.dimensions;
|
||||
const double* fc = grid_.face_centroids + dim*face_;
|
||||
if (degree_ < 2) {
|
||||
if (degree_ < 2 || dim < 2) {
|
||||
std::copy(fc, fc + dim, coord);
|
||||
return;
|
||||
}
|
||||
@ -110,6 +119,13 @@ namespace Opm
|
||||
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
|
||||
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
|
||||
const double* nc = grid_.node_coordinates;
|
||||
if (dim == 2) {
|
||||
ASSERT(nn == 2);
|
||||
const double* pa[3] = { nc + dim*fnodes[0], fc, nc + dim*fnodes[1] };
|
||||
std::copy(pa[index], pa[index] + dim, coord);
|
||||
return;
|
||||
}
|
||||
ASSERT(dim == 3);
|
||||
if (index < nn) {
|
||||
// Boundary edge midpoint.
|
||||
const int node0 = fnodes[index];
|
||||
@ -134,6 +150,11 @@ namespace Opm
|
||||
}
|
||||
// Degree 2 case.
|
||||
const int dim = grid_.dimensions;
|
||||
if (dim == 2) {
|
||||
const double simpsonw[3] = { 1.0/6.0, 4.0/6.0, 1.0/6.0 };
|
||||
return simpsonw[index];
|
||||
}
|
||||
ASSERT(dim == 3);
|
||||
const double* fc = grid_.face_centroids + dim*face_;
|
||||
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
|
||||
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
|
||||
|
||||
@ -66,22 +66,6 @@ namespace
|
||||
mutable std::vector<double> pt;
|
||||
};
|
||||
|
||||
struct LinearFunc
|
||||
{
|
||||
double operator()(const double* x) const
|
||||
{
|
||||
return 1.0*x[0] + 2.0*x[1] + 1.0*x[2] + 3.0;
|
||||
}
|
||||
};
|
||||
|
||||
struct QuadraticFunc
|
||||
{
|
||||
double operator()(const double* x) const
|
||||
{
|
||||
return 1.0*x[0]*x[1] + 2.0*x[1] + 4.0*x[2] + 3.0;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template <class Quadrature, class Func>
|
||||
void testSingleCase(const UnstructuredGrid& grid,
|
||||
@ -98,45 +82,124 @@ namespace
|
||||
} // anonymous namespace
|
||||
|
||||
|
||||
|
||||
static void test3dCart()
|
||||
namespace cart2d
|
||||
{
|
||||
// Set up 3d 1-cell cartesian case.
|
||||
GridManager g(1, 1, 1);
|
||||
const UnstructuredGrid& grid = *g.c_grid();
|
||||
struct ConstantFunc
|
||||
{
|
||||
double operator()(const double*) const
|
||||
{
|
||||
return 1.234;
|
||||
}
|
||||
};
|
||||
|
||||
// CellQuadrature tests.
|
||||
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 1, 5.0);
|
||||
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 2, 5.0);
|
||||
testSingleCase<CellQuadrature, QuadraticFunc>(grid, 0, 2, 6.25);
|
||||
struct LinearFunc
|
||||
{
|
||||
double operator()(const double* x) const
|
||||
{
|
||||
return 1.0*x[0] + 2.0*x[1] + 3.0;
|
||||
}
|
||||
};
|
||||
|
||||
// FaceQuadrature tests, degree 1 precision.
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 1, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 1, 5.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 1, 4.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 1, 6.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 4, 1, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 5, 1, 5.5);
|
||||
struct QuadraticFunc
|
||||
{
|
||||
double operator()(const double* x) const
|
||||
{
|
||||
return 3.0*x[0]*x[0] + 1.0*x[0]*x[1] + 2.0*x[1] + 3.0;
|
||||
}
|
||||
};
|
||||
|
||||
// FaceQuadrature tests, degree 2 precision.
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 2, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 2, 5.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 2, 4.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 2, 6.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 4, 2, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 5, 2, 5.5);
|
||||
static void test()
|
||||
{
|
||||
// Set up 2d 1-cell cartesian case.
|
||||
GridManager g(1, 1);
|
||||
const UnstructuredGrid& grid = *g.c_grid();
|
||||
|
||||
// FaceQuadrature tests, quadratic function, degree 2 precision.
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 0, 2, 6.0);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 1, 2, 6.5);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 2, 2, 5.0);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 3, 2, 7.5);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 4, 2, 4.25);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 5, 2, 8.25);
|
||||
}
|
||||
// CellQuadrature tests.
|
||||
testSingleCase<CellQuadrature, ConstantFunc>(grid, 0, 1, 1.234);
|
||||
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 1, 4.5);
|
||||
testSingleCase<CellQuadrature, ConstantFunc>(grid, 0, 2, 1.234);
|
||||
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 2, 4.5);
|
||||
testSingleCase<CellQuadrature, QuadraticFunc>(grid, 0, 2, 5.25);
|
||||
|
||||
// FaceQuadrature tests, degree 1 precision.
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 1, 4);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 1, 5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 1, 3.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 1, 5.5);
|
||||
|
||||
// FaceQuadrature tests, degree 2 precision.
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 2, 4);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 2, 5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 2, 3.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 2, 5.5);
|
||||
|
||||
// FaceQuadrature tests, quadratic function, degree 2 precision.
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 0, 2, 4.0);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 1, 2, 7.5);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 2, 2, 4.0);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 3, 2, 6.5);
|
||||
}
|
||||
} // namespace cart2d
|
||||
|
||||
|
||||
namespace cart3d
|
||||
{
|
||||
struct LinearFunc
|
||||
{
|
||||
double operator()(const double* x) const
|
||||
{
|
||||
return 1.0*x[0] + 2.0*x[1] + 1.0*x[2] + 3.0;
|
||||
}
|
||||
};
|
||||
|
||||
struct QuadraticFunc
|
||||
{
|
||||
double operator()(const double* x) const
|
||||
{
|
||||
return 1.0*x[0]*x[1] + 2.0*x[1] + 4.0*x[2] + 3.0;
|
||||
}
|
||||
};
|
||||
|
||||
static void test()
|
||||
{
|
||||
// Set up 3d 1-cell cartesian case.
|
||||
GridManager g(1, 1, 1);
|
||||
const UnstructuredGrid& grid = *g.c_grid();
|
||||
|
||||
// CellQuadrature tests.
|
||||
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 1, 5.0);
|
||||
testSingleCase<CellQuadrature, LinearFunc>(grid, 0, 2, 5.0);
|
||||
testSingleCase<CellQuadrature, QuadraticFunc>(grid, 0, 2, 6.25);
|
||||
|
||||
// FaceQuadrature tests, degree 1 precision.
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 1, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 1, 5.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 1, 4.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 1, 6.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 4, 1, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 5, 1, 5.5);
|
||||
|
||||
// FaceQuadrature tests, degree 2 precision.
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 0, 2, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 1, 2, 5.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 2, 2, 4.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 3, 2, 6.0);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 4, 2, 4.5);
|
||||
testSingleCase<FaceQuadrature, LinearFunc>(grid, 5, 2, 5.5);
|
||||
|
||||
// FaceQuadrature tests, quadratic function, degree 2 precision.
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 0, 2, 6.0);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 1, 2, 6.5);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 2, 2, 5.0);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 3, 2, 7.5);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 4, 2, 4.25);
|
||||
testSingleCase<FaceQuadrature, QuadraticFunc>(grid, 5, 2, 8.25);
|
||||
}
|
||||
|
||||
} // namespace cart3d
|
||||
|
||||
BOOST_AUTO_TEST_CASE(test_quadratures)
|
||||
{
|
||||
test3dCart();
|
||||
cart2d::test();
|
||||
cart3d::test();
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user