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

Fix support for velocity interpolation in 2d.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2012-10-17 21:30:53 +02:00
parent f9efd72ecc
commit dcee95ec96
1 changed files with 21 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
opm/core/utility/VelocityInterpolation.cpp
View File

@@ -86,6 +86,15 @@ namespace Opm
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,
@@ -97,6 +106,16 @@ namespace Opm
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
@@ -139,8 +158,7 @@ namespace Opm
CornerInfo ci;
ci.corner_id = corner_id_count++;;
ci.vertex = *it;
typedef double* DblPtr;
DblPtr fnorm[3] = { 0, 0, 0 };
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;
@@ -155,7 +173,7 @@ namespace Opm
}
ASSERT(fi == dim);
adj_faces_.insert(adj_faces_.end(), vert_adj_faces.begin(), vert_adj_faces.end());
const double corner_vol = std::fabs(determinantOf(fnorm[0], fnorm[1], fnorm[2]));
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());