Commit 382323d7 authored by frey_m's avatar frey_m
Browse files

SAAMG: nr --> nr_m and hr --> hr_m

parent 8e0d416c
......@@ -39,8 +39,10 @@
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t /*hr*/,
std::string interpl){
Vector_t hr,
std::string interpl)
: IrregularDomain(nr, hr)
{
bgeom_m = bgeom;
setRangeMin(bgeom->getmincoords());
......@@ -54,10 +56,6 @@ ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
"ArbitraryDomain::ArbitraryDomain()",
"No point inside geometry found/set!");
}
setNr(nr);
for(int i=0; i<3; i++)
Geo_hr_m[i] = (max_m[i] - min_m[i])/nr[i];
setHr(Geo_hr_m);
startId = 0;
......@@ -87,11 +85,11 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr[2] - 1) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr_m[2] - 1) ? 0 : 1;
int yGhostOffsetLeft = (localId[1].first()== 0) ? 0 : 1;
int yGhostOffsetRight = (localId[1].last() == nr[1] - 1) ? 0 : 1;
int yGhostOffsetRight = (localId[1].last() == nr_m[1] - 1) ? 0 : 1;
int xGhostOffsetLeft = (localId[0].first()== 0) ? 0 : 1;
int xGhostOffsetRight = (localId[0].last() == nr[0] - 1) ? 0 : 1;
int xGhostOffsetRight = (localId[0].last() == nr_m[0] - 1) ? 0 : 1;
hasGeometryChanged_m = true;
......@@ -112,17 +110,17 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
// example (-0.13 to +0.025). -DW
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
//saveP_old[2] = (idz - (nr[2]-1)/2.0)*hr[2];
//saveP_old[2] = (idz - (nr_m[2]-1)/2.0)*hr[2];
P[2] = min_m[2] + (idz + 0.5) * hr[2];
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
//saveP_old[1] = (idy - (nr[1]-1)/2.0)*hr[1];
//saveP_old[1] = (idy - (nr_m[1]-1)/2.0)*hr[1];
P[1] = min_m[1] + (idy + 0.5) * hr[1];
for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
//saveP_old[0] = (idx - (nr[0]-1)/2.0)*hr[0];
//saveP_old[0] = (idx - (nr_m[0]-1)/2.0)*hr[0];
P[0] = min_m[0] + (idx + 0.5) * hr[0];
// *gmsg << "Now working on point " << saveP << " (original was " << saveP_old << ")" << endl;
......@@ -244,8 +242,8 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
//number of ghost nodes to the left
int numGhostNodesLeft = 0;
if(localId[2].first() != 0) {
for(int idx = 0; idx < nr[0]; idx++) {
for(int idy = 0; idy < nr[1]; idy++) {
for(int idx = 0; idx < nr_m[0]; idx++) {
for(int idy = 0; idy < nr_m[1]; idy++) {
if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
numGhostNodesLeft++;
}
......@@ -259,8 +257,8 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
numXY.clear();
for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
int numxy = 0;
for(int idx = 0; idx < nr[0]; idx++) {
for(int idy = 0; idy < nr[1]; idy++) {
for(int idx = 0; idx < nr_m[0]; idx++) {
for(int idy = 0; idy < nr_m[1]; idy++) {
if(isInside(idx, idy, idz))
numxy++;
}
......@@ -281,8 +279,8 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
INFOMSG(level2 << "* Building up index and coordinate map..." << endl);
for(int idz = localId[2].first() - zGhostOffsetLeft; idz <= localId[2].last() + zGhostOffsetRight; idz++) {
for(int idy = 0; idy < nr[1]; idy++) {
for(int idx = 0; idx < nr[0]; idx++) {
for(int idy = 0; idy < nr_m[1]; idy++) {
for(int idx = 0; idx < nr_m[0]; idx++) {
if(isInside(idx, idy, idz)) {
IdxMap[toCoordIdx(idx, idy, idz)] = index;
CoordMap[index] = toCoordIdx(idx, idy, idz);
......@@ -297,7 +295,7 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
return (idz * nr[1] + idy) * nr[0] + idx;
return (idz * nr_m[1] + idy) * nr_m[0] + idx;
}
// Conversion from (x,y,z) to index on the 3D grid
......@@ -312,10 +310,10 @@ int ArbitraryDomain::getIdx(int idx, int idy, int idz) {
// Conversion from a 3D index to (x,y,z)
inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) {
int id = CoordMap[idxyz];
idx = id % (int)nr[0];
id /= nr[0];
idy = id % (int)nr[1];
id /= nr[1];
idx = id % (int)nr_m[0];
id /= nr_m[0];
idy = id % (int)nr_m[1];
id /= nr_m[1];
idz = id;
}
......@@ -328,9 +326,9 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
Vector_t P;
P[0] = min_m[0] + (idx + 0.5) * hr[0];
P[1] = min_m[1] + (idy + 0.5) * hr[1];
P[2] = min_m[2] + (idz + 0.5) * hr[2];
P[0] = min_m[0] + (idx + 0.5) * hr_m[0];
P[1] = min_m[1] + (idy + 0.5) * hr_m[1];
P[2] = min_m[2] + (idz + 0.5) * hr_m[2];
return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
}
......@@ -340,9 +338,9 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
Vector_t P;
P[0] = (idx - (nr[0]-1)/2.0) * hr[0];
P[1] = (idy - (nr[1]-1)/2.0) * hr[1];
P[2] = (idz - (nr[2]-1)/2.0) * hr[2];
P[0] = (idx - (nr_m[0]-1)/2.0) * hr_m[0];
P[1] = (idy - (nr_m[1]-1)/2.0) * hr_m[1];
P[2] = (idz - (nr_m[2]-1)/2.0) * hr_m[2];
bool ret = false;
int countH, countL;
......@@ -409,13 +407,13 @@ void ArbitraryDomain::getBoundaryStencil(int idx, int idy, int idz,
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
StencilValue_t& value, double& /*scaleFactor*/)
{
value.west = -1/(hr[0]*hr[0]);
value.east = -1/(hr[0]*hr[0]);
value.north = -1/(hr[1]*hr[1]);
value.south = -1/(hr[1]*hr[1]);
value.front = -1/(hr[2]*hr[2]);
value.back = -1/(hr[2]*hr[2]);
value.center = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
value.west = -1/(hr_m[0]*hr_m[0]);
value.east = -1/(hr_m[0]*hr_m[0]);
value.north = -1/(hr_m[1]*hr_m[1]);
value.south = -1/(hr_m[1]*hr_m[1]);
value.front = -1/(hr_m[2]*hr_m[2]);
value.back = -1/(hr_m[2]*hr_m[2]);
value.center = 2/(hr_m[0]*hr_m[0]) + 2/(hr_m[1]*hr_m[1]) + 2/(hr_m[2]*hr_m[2]);
if(!isInside(idx-1,idy,idz))
value.west = 0.0;
......@@ -438,29 +436,29 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
{
scaleFactor = 1;
double cx = (idx - (nr[0]-1)/2.0)*hr[0];
double cy = (idy - (nr[1]-1)/2.0)*hr[1];
double cz = (idz - (nr[2]-1)/2.0)*hr[2];
double cx = (idx - (nr_m[0]-1)/2.0)*hr_m[0];
double cy = (idy - (nr_m[1]-1)/2.0)*hr_m[1];
double cz = (idz - (nr_m[2]-1)/2.0)*hr_m[2];
double dx_w=hr[0];
double dx_e=hr[0];
double dy_n=hr[1];
double dy_s=hr[1];
double dz_f=hr[2];
double dz_b=hr[2];
double dx_w=hr_m[0];
double dx_e=hr_m[0];
double dy_n=hr_m[1];
double dy_s=hr_m[1];
double dz_f=hr_m[2];
double dz_b=hr_m[2];
value.center = 0.0;
std::tuple<int, int, int> coordxyz(idx, idy, idz);
if (idx == nr[0]-1)
if (idx == nr_m[0]-1)
dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
if (idx == 0)
dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
if (idy == nr[1]-1)
if (idy == nr_m[1]-1)
dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
if (idy == 0)
dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
if (idz == nr[2]-1)
if (idz == nr_m[2]-1)
dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
if (idz == 0)
dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
......@@ -513,7 +511,7 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
if(dy_s == 0)
m2 = dy_n;
value.center = 2 / hr[2];
value.center = 2 / hr_m[2];
if(dx_w != 0 || dx_e != 0)
value.center *= (dx_w + dx_e);
if(dy_n != 0 || dy_s != 0)
......
......@@ -110,8 +110,6 @@ private:
// Interpolation type
int interpolationMethod;
Vector_t Geo_nr_m;
Vector_t Geo_hr_m;
Vector_t geomCentroid_m;
Vector_t minCoords_m;
Vector_t maxCoords_m;
......
......@@ -38,15 +38,13 @@ extern Inform *gmsg;
BoxCornerDomain::BoxCornerDomain(double A, double B, double C, double length,
double L1, double L2, Vector_t nr, Vector_t hr,
std::string interpl)
: IrregularDomain(nr, hr)
{
setRangeMin(Vector_t(-A, -B, L1));
setRangeMax(Vector_t( A, B, L1 + L2));
C_m = C;
length_m = length;
setNr(nr);
setHr(hr);
if(interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "LINEAR")
......@@ -78,7 +76,7 @@ BoxCornerDomain::~BoxCornerDomain() {
void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
//there is nothing to be done if the mesh spacings have not changed
// if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
// if(hr_m[0] == getHr()[0] && hr_m[1] == getHr()[1] && hr_m[2] == getHr()[2]) {
// hasGeometryChanged_m = false;
// return;
// }
......@@ -104,9 +102,9 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
// build a index and coordinate map
int idx = 0;
int x, y, z;
for(x = 0; x < nr[0]; x++) {
for(y = 0; y < nr[1]; y++) {
for(z = 0; z < nr[2]; z++) {
for(x = 0; x < nr_m[0]; x++) {
for(y = 0; y < nr_m[1]; y++) {
for(z = 0; z < nr_m[2]; z++) {
if(isInside(x, y, z)) {
IdxMap[toCoordIdx(x, y, z)] = idx;
CoordMap[idx++] = toCoordIdx(x, y, z);
......@@ -126,18 +124,18 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
// calculate intersection
for(int z = 0; z < nr[2]; z++) {
for(int z = 0; z < nr_m[2]; z++) {
for(int x = 0; x < nr[0]; x++) {
for(int x = 0; x < nr_m[0]; x++) {
// the x coordinate does not change in the CornerBox geometry
std::pair<int, int> pos(x, z);
IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, 0.5*getXRangeMax()));
IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, 0.5*getXRangeMin()));
}
for(int y = 0; y < nr[1]; y++) {
for(int y = 0; y < nr_m[1]; y++) {
std::pair<int, int> pos(y, z);
double yt = getB(z*hr[2]);
double yt = getB(z*hr_m[2]);
double yb = 0.5*getYRangeMin();
IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, yt));
IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, yb));
......@@ -172,13 +170,13 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, StencilValue_t&
{
scaleFactor = 1.0;
value.west = -1 / hr[0] * 1 / hr[0];
value.east = -1 / hr[0] * 1 / hr[0];
value.north = -1 / hr[1] * 1 / hr[1];
value.south = -1 / hr[1] * 1 / hr[1];
value.front = -1 / hr[2] * 1 / hr[2];
value.back = -1 / hr[2] * 1 / hr[2];
value.center = 2 / hr[0] * 1 / hr[0] + 2 / hr[1] * 1 / hr[1] + 2 / hr[2] * 1 / hr[2];
value.west = -1 / hr_m[0] * 1 / hr_m[0];
value.east = -1 / hr_m[0] * 1 / hr_m[0];
value.north = -1 / hr_m[1] * 1 / hr_m[1];
value.south = -1 / hr_m[1] * 1 / hr_m[1];
value.front = -1 / hr_m[2] * 1 / hr_m[2];
value.back = -1 / hr_m[2] * 1 / hr_m[2];
value.center = 2 / hr_m[0] * 1 / hr_m[0] + 2 / hr_m[1] * 1 / hr_m[1] + 2 / hr_m[2] * 1 / hr_m[2];
// we are a right boundary point
if(!isInside(x + 1, y, z))
......@@ -196,7 +194,7 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, StencilValue_t&
if(!isInside(x, y - 1, z))
value.south = 0.0;
if(z == 1 || z == nr[2] - 2) {
if(z == 1 || z == nr_m[2] - 2) {
// case where we are on the Robin BC in Z-direction
// where we distinguish two cases
......@@ -209,13 +207,13 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, StencilValue_t&
// add contribution of Robin discretization to center point
// d the distance between the center of the bunch and the boundary
//double cx = (x-(nr[0]-1)/2)*hr[0];
//double cy = (y-(nr[1]-1)/2)*hr[1];
//double cz = hr[2]*(nr[2]-1);
//double cx = (x-(nr_m[0]-1)/2)*hr_m[0];
//double cy = (y-(nr_m[1]-1)/2)*hr_m[1];
//double cz = hr_m[2]*(nr_m[2]-1);
//double d = sqrt(cx*cx+cy*cy+cz*cz);
double d = hr[2] * (nr[2] - 1) / 2;
value.center += 2 / (d * hr[2]);
//value.center += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
double d = hr_m[2] * (nr_m[2] - 1) / 2;
value.center += 2 / (d * hr_m[2]);
//value.center += 2/((hr_m[2]*(nr_m[2]-1)/2.0) * hr_m[2]);
// scale all stencil-points in z-plane with 0.5 (Robin discretization)
value.west /= 2.0;
......@@ -233,8 +231,8 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& v
{
scaleFactor = 1.0;
double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
double cy = y * hr[1] - (nr[1] - 1) * hr[1] / 2.0;
double cx = x * hr_m[0] - (nr_m[0] - 1) * hr_m[0] / 2.0;
double cy = y * hr_m[1] - (nr_m[1] - 1) * hr_m[1] / 2.0;
//XXX: calculate intersection on the fly
/*
......@@ -256,10 +254,10 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& v
dy = it->second;
*/
double dw = hr[0];
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
double dw = hr_m[0];
double de = hr_m[0];
double dn = hr_m[1];
double ds = hr_m[1];
value.center = 0.0;
//we are a right boundary point
......@@ -302,12 +300,12 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& v
value.south = -1 / (ds * ds);
}
value.front = -1 / (hr[2] * hr[2]);
value.back = -1 / (hr[2] * hr[2]);
value.center += 2 / (hr[2] * hr[2]);
value.front = -1 / (hr_m[2] * hr_m[2]);
value.back = -1 / (hr_m[2] * hr_m[2]);
value.center += 2 / (hr_m[2] * hr_m[2]);
// handle boundary condition in z direction
if(z == 0 || z == nr[2] - 1) {
if(z == 0 || z == nr_m[2] - 1) {
// case where we are on the NEUMAN BC in Z-direction
// where we distinguish two cases
......@@ -316,9 +314,9 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& v
else
value.back = 0.0;
//hr[2]*(nr2[2]-1)/2 = radius
double d = hr[2] * (nr[2] - 1) / 2;
value.center += 2 / (d * hr[2]);
//hr_m[2]*(nr_m2[2]-1)/2 = radius
double d = hr_m[2] * (nr_m[2] - 1) / 2;
value.center += 2 / (d * hr_m[2]);
value.west /= 2.0;
value.east /= 2.0;
......@@ -335,8 +333,8 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& v
void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
{
double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
double cx = (x - (nr_m[0] - 1) / 2.0) * hr_m[0];
double cy = (y - (nr_m[1] - 1) / 2.0) * hr_m[1];
double dx = getXIntersection(cx, z);
double dy = getYIntersection(cy, z);
......@@ -361,10 +359,10 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
dy = it->second;
*/
double dw = hr[0];
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
double dw = hr_m[0];
double de = hr_m[0];
double dn = hr_m[1];
double ds = hr_m[1];
value.west = 1.0;
value.east = 1.0;
value.north = 1.0;
......@@ -373,54 +371,54 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
value.back = 1.0;
value.center = 0.0;
//TODO: = cx+hr[0] > dx && cx > 0
//if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
//TODO: = cx+hr_m[0] > dx && cx > 0
//if((x-nr_m[0]/2.0+1)*hr_m[0] > dx && cx > 0) {
////we are a right boundary point
////if(!isInside(x+1,y,z)) {
//de = dx-cx;
//}
//if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
//if((x-nr_m[0]/2.0-1)*hr_m[0] < dx && cx < 0) {
////we are a left boundary point
////if(!isInside(x-1,y,z)) {
//dw = std::abs(dx)-std::abs(cx);
//}
//if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
//if((y-nr_m[1]/2.0+1)*hr_m[1] > dy && cy > 0) {
////we are a upper boundary point
////if(!isInside(x,y+1,z)) {
//dn = dy-cy;
//}
//if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
//if((y-nr_m[1]/2.0-1)*hr_m[1] < dy && cy < 0) {
////we are a lower boundary point
////if(!isInside(x,y-1,z)) {
//ds = std::abs(dy)-std::abs(cy);
//}
//TODO: = cx+hr[0] > dx && cx > 0
//if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
//TODO: = cx+hr_m[0] > dx && cx > 0
//if((x-nr_m[0]/2.0+1)*hr_m[0] > dx && cx > 0) {
//we are a right boundary point
if(!isInside(x + 1, y, z)) {
de = dx - cx;
value.east = 0.0;
}
//if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
//if((x-nr_m[0]/2.0-1)*hr_m[0] < dx && cx < 0) {
//we are a left boundary point
if(!isInside(x - 1, y, z)) {
dw = std::abs(dx) - std::abs(cx);
value.west = 0.0;
}
//if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
//if((y-nr_m[1]/2.0+1)*hr_m[1] > dy && cy > 0) {
//we are a upper boundary point
if(!isInside(x, y + 1, z)) {
dn = dy - cy;
value.north = 0.0;
}
//if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
//if((y-nr_m[1]/2.0-1)*hr_m[1] < dy && cy < 0) {
//we are a lower boundary point
if(!isInside(x, y - 1, z)) {
ds = std::abs(dy) - std::abs(cy);
......@@ -432,14 +430,14 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
value.east *= -1.0 / (de * (dw + de));
value.north *= -1.0 / (dn * (dn + ds));
value.south *= -1.0 / (ds * (dn + ds));
value.front = -1 / (hr[2] * (hr[2] + hr[2]));
value.back = -1 / (hr[2] * (hr[2] + hr[2]));
value.front = -1 / (hr_m[2] * (hr_m[2] + hr_m[2]));
value.back = -1 / (hr_m[2] * (hr_m[2] + hr_m[2]));
//TODO: problem when de,dw,dn,ds == 0
//is NOT a regular BOUND PT
value.center += 1 / de * 1 / dw;
value.center += 1 / dn * 1 / ds;
value.center += 1 / hr[2] * 1 / hr[2];
value.center += 1 / hr_m[2] * 1 / hr_m[2];
scaleFactor = 0.5;
......@@ -450,47 +448,47 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
/*
if(dw != 0 && !wIsB)
W = -1/dw * (dn+ds) * 2*hr[2];
W = -1/dw * (dn+ds) * 2*hr_m[2];
else
W = 0;
if(de != 0 && !eIsB)
E = -1/de * (dn+ds) * 2*hr[2];
E = -1/de * (dn+ds) * 2*hr_m[2];
else
E = 0;
if(dn != 0 && !nIsB)
N = -1/dn * (dw+de) * 2*hr[2];
N = -1/dn * (dw+de) * 2*hr_m[2];
else
N = 0;
if(ds != 0 && !sIsB)
S = -1/ds * (dw+de) * 2*hr[2];
S = -1/ds * (dw+de) * 2*hr_m[2];
else
S = 0;
F = -(dw+de)*(dn+ds)/hr[2];
B = -(dw+de)*(dn+ds)/hr[2];
F = -(dw+de)*(dn+ds)/hr_m[2];
B = -(dw+de)*(dn+ds)/hr_m[2];
*/
//if(dw != 0)
//W = -2*hr[2]*(dn+ds)/dw;
//W = -2*hr_m[2]*(dn+ds)/dw;
//else
//W = 0;
//if(de != 0)
//E = -2*hr[2]*(dn+ds)/de;
//E = -2*hr_m[2]*(dn+ds)/de;
//else
//E = 0;
//if(dn != 0)
//N = -2*hr[2]*(dw+de)/dn;
//N = -2*hr_m[2]*(dw+de)/dn;
//else
//N = 0;
//if(ds != 0)
//S = -2*hr[2]*(dw+de)/ds;
//S = -2*hr_m[2]*(dw+de)/ds;
//else
//S = 0;
//F = -(dw+de)*(dn+ds)/hr[2];
//B = -(dw+de)*(dn+ds)/hr[2];
//F = -(dw+de)*(dn+ds)/hr_m[2];
//B = -(dw+de)*(dn+ds)/hr_m[2];
//// RHS scaleFactor for current 3D index
//// Factor 0.5 results from the SW/quadratic extrapolation
//scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr[2]);
//scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr_m[2]);
// catch the case where a point lies on the boundary
//FIXME: do this more elegant!
......@@ -505,19 +503,19 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
//if(ds == 0)
//m2 = dn;
////XXX: dn+ds || dw+de can be 0
////C = 2*(dn+ds)*(dw+de)/hr[2];
//C = 2/hr[2];
////C = 2*(dn+ds)*(dw+de)/hr_m[2];
//C = 2/hr_m[2];
//if(dw != 0 || de != 0)
//C *= (dw+de);
//if(dn != 0 || ds != 0)
//C *= (dn+ds);
//if(dw != 0 || de != 0)
//C += (2*hr[2])*(dn+ds)*(dw+de)/m1;
//C += (2*hr_m[2])*(dn+ds)*(dw+de)/m1;
//if(dn != 0 || ds != 0)
//C += (2*hr[2])*(dw+de)*(dn+ds)/m2;
//C += (2*hr_m[2])*(dw+de)*(dn+ds)/m2;
//handle Neumann case
//if(z == 0 || z == nr[2]-1) {
//if(z == 0 || z == nr_m[2]-1) {
//if(z == 0)
//F = 0.0;
......@@ -535,7 +533,7 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
//}
// handle boundary condition in z direction
if(z == 0 || z == nr[2] - 1) {
if(z == 0 || z == nr_m[2] - 1) {
// case where we are on the NEUMAN BC in Z-direction
// where we distinguish two cases
......@@ -544,10 +542,10 @@ void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t
else
value.back = 0.0;
//value.center += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
//hr[2]*(nr2[2]-1)/2 = radius
double d = hr[2] * (nr[2] - 1) / 2;
value.center += 2 / (d * hr[2]);
//value.center += 2/((hr_m[2]*(nr_m[2]-1)/2.0) * hr_m[2]);
//hr_m[2]*(nr_m2[2]-1)/2 = radius
double d = hr_m[2] * (nr_m[2] - 1) / 2;
value.center += 2 / (d * hr_m[2]);
value.west /= 2.0;
value.east /= 2.0;
......
......@@ -107,10 +107,10 @@ public: