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

White space changes

Browse Source
This commit is contained in:
Joakim Hove
2019-12-20 12:23:14 +01:00
parent 067fcd636c
commit 3fd9b20cd7
1 changed files with 57 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

116
src/opm/common/utility/numeric/calculateCellVol.cpp
View File

@@ -22,14 +22,14 @@
#include <opm/common/utility/numeric/calculateCellVol.hpp>
#include <opm/common/ErrorMacros.hpp>
/*
/*
Cell volume calculation based on following publication:
D. K Pointing, Corner Point Geometry in Reservoir Simulation ,
ECMOR I - 1st European Conference on the Mathematics of Oil Recovery,
1989
D. K Pointing, Corner Point Geometry in Reservoir Simulation ,
ECMOR I - 1st European Conference on the Mathematics of Oil Recovery,
1989
*/
/*
The expressions {C(0,0,0),..C(1,1,1)} have a nice interpretation in
@@ -44,107 +44,105 @@ double C(const double* r, int i1, int i2, int i3){
double temp;
temp = 0.0;
if ((i1 == 0) && (i2 == 0) && (i3 == 0)){
temp = r[0];
temp = r[0];
}
if ((i1 == 1) && (i2 == 0) && (i3 == 0)){
temp = r[1] - r[0];
temp = r[1] - r[0];
}
if ((i1 == 0) && (i2 == 1) && (i3 == 0)){
temp = r[2] - r[0];
temp = r[2] - r[0];
}
if ((i1 == 0) && (i2 == 0) && (i3 == 1)){
temp = r[4]-r[0];
temp = r[4]-r[0];
}
if ((i1 == 1) && (i2 == 1) && (i3 == 0)){
temp = r[3] + r[0] - r[2] - r[1];
temp = r[3] + r[0] - r[2] - r[1];
}
if ((i1 == 0) && (i2 == 1) && (i3 == 1)){
temp = r[6] + r[0] - r[4] - r[2];
temp = r[6] + r[0] - r[4] - r[2];
}
if ((i1 == 1) && (i2 == 0) && (i3 == 1)){
temp = r[5] + r[0] - r[4] - r[1];
temp = r[5] + r[0] - r[4] - r[1];
}
if ((i1 == 1) && (i2 == 1) && (i3 == 1)){
temp = r[7] + r[4] + r[2] + r[1] - r[6] - r[5] - r[3] - r[0];
temp = r[7] + r[4] + r[2] + r[1] - r[6] - r[5] - r[3] - r[0];
}
return temp;
}
double perm123sign(int i1, int i2, int i3){
double temp;
if ((i1 ==1 ) && (i2==2) && (i3 == 3)){
temp = 1.0;
temp = 1.0;
} else if ((i1 == 1) && (i2==3) && (i3 == 2)){
temp = -1.0;
temp = -1.0;
} else if ((i1 == 2) && (i2 == 1) && (i3 == 3)){
temp = -1.0;
temp = -1.0;
} else if ((i1 == 2) && (i2 == 3) && (i3 == 1)){
temp = 1.0;
temp = 1.0;
} else if ((i1 == 3) && (i2 == 1) && (i3 == 2)){
temp = 1.0;
temp = 1.0;
} else if ((i1 == 3) && (i2 == 2) && (i3 == 1)){
temp = -1.0;
temp = -1.0;
} else {
OPM_THROW(std::logic_error, "Wrong indices in perm123sign");
OPM_THROW(std::logic_error, "Wrong indices in perm123sign");
}
return temp;
}
double calculateCellVol(const std::array<double,8>& X, const std::array<double,8>& Y, const std::array<double,8>& Z){
double volume = 0.0;
const double* vect[3];
const double* vect[3];
int permutation[] = {1, 2, 3};
do {
for (int j = 0; j < 3; ++j){
if (permutation[j] == 1){
vect[j] = X.data();
} else if (permutation[j] == 2){
vect[j] = Y.data();
} else if (permutation[j] == 3){
vect[j] = Z.data();
} else {
// this condition can never happen, since all values in 'permutation'
// is covered, but compiler analysis may not be deep enough to not give
// warnings about 'vect' being uninitialized further down.
assert(false);
vect[j] = 0;
}
}
for (int j = 0; j < 3; ++j){
if (permutation[j] == 1){
vect[j] = X.data();
} else if (permutation[j] == 2){
vect[j] = Y.data();
} else if (permutation[j] == 3){
vect[j] = Z.data();
} else {
// this condition can never happen, since all values in 'permutation'
// is covered, but compiler analysis may not be deep enough to not give
// warnings about 'vect' being uninitialized further down.
assert(false);
vect[j] = 0;
}
}
for (int pb = 0; pb < 2; ++pb){
for (int pg = 0; pg < 2; ++pg){
for (int qa = 0; qa < 2; ++qa){
for (int qg = 0; qg < 2; ++qg){
for (int ra = 0; ra < 2; ++ra){
for (int rb = 0; rb < 2; ++rb){
const double cprod = C(vect[0], 1, pb, pg)*C(vect[1], qa, 1, qg)*C(vect[2], ra, rb, 1);
const double denom = (qa+ra+1) * (pb+rb+1) * (pg+qg+1);
volume += perm123sign(permutation[0], permutation[1], permutation[2]) * cprod / denom;
}
}
}
}
}
}
for (int pb = 0; pb < 2; ++pb){
for (int pg = 0; pg < 2; ++pg){
for (int qa = 0; qa < 2; ++qa){
for (int qg = 0; qg < 2; ++qg){
for (int ra = 0; ra < 2; ++ra){
for (int rb = 0; rb < 2; ++rb){
const double cprod = C(vect[0], 1, pb, pg)*C(vect[1], qa, 1, qg)*C(vect[2], ra, rb, 1);
const double denom = (qa+ra+1) * (pb+rb+1) * (pg+qg+1);
volume += perm123sign(permutation[0], permutation[1], permutation[2]) * cprod / denom;
}
}
}
}
}
}
} while (std::next_permutation(std::begin(permutation), std::end(permutation)));
return std::fabs(volume);