Commit 74dde0f2 authored by frey_m's avatar frey_m
Browse files

SAAMG: use structs for stencil indices and values

parent d8a653c9
...@@ -383,59 +383,56 @@ int ArbitraryDomain::getNumXY(int z) { ...@@ -383,59 +383,56 @@ int ArbitraryDomain::getNumXY(int z) {
return numXY[z]; return numXY[z];
} }
void ArbitraryDomain::getBoundaryStencil(int idx, int idy, int idz, double &W, void ArbitraryDomain::getBoundaryStencil(int idx, int idy, int idz,
double &E, double &S, double &N, double &F, StencilValue_t& value, double &scaleFactor)
double &B, double &C, double &scaleFactor)
{ {
scaleFactor = 1.0; scaleFactor = 1.0;
// determine which interpolation method we use for points near the boundary // determine which interpolation method we use for points near the boundary
switch(interpolationMethod){ switch(interpolationMethod){
case CONSTANT: case CONSTANT:
constantInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor); constantInterpolation(idx,idy,idz,value,scaleFactor);
break; break;
case LINEAR: case LINEAR:
linearInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor); linearInterpolation(idx,idy,idz,value,scaleFactor);
break; break;
case QUADRATIC: case QUADRATIC:
// QuadraticInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor); // QuadraticInterpolation(idx,idy,idz,value,scaleFactor);
break; break;
} }
// stencil center value has to be positive! // stencil center value has to be positive!
assert(C > 0); assert(value.center > 0);
} }
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz, double& W, void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
double& E, double& S, double& N, double& F, StencilValue_t& value, double& /*scaleFactor*/)
double& B, double& C, double& /*scaleFactor*/)
{ {
W = -1/(hr[0]*hr[0]); value.west = -1/(hr[0]*hr[0]);
E = -1/(hr[0]*hr[0]); value.east = -1/(hr[0]*hr[0]);
N = -1/(hr[1]*hr[1]); value.north = -1/(hr[1]*hr[1]);
S = -1/(hr[1]*hr[1]); value.south = -1/(hr[1]*hr[1]);
F = -1/(hr[2]*hr[2]); value.front = -1/(hr[2]*hr[2]);
B = -1/(hr[2]*hr[2]); value.back = -1/(hr[2]*hr[2]);
C = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]); value.center = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
if(!isInside(idx-1,idy,idz)) if(!isInside(idx-1,idy,idz))
W = 0.0; value.west = 0.0;
if(!isInside(idx+1,idy,idz)) if(!isInside(idx+1,idy,idz))
E = 0.0; value.east = 0.0;
if(!isInside(idx,idy+1,idz)) if(!isInside(idx,idy+1,idz))
N = 0.0; value.north = 0.0;
if(!isInside(idx,idy-1,idz)) if(!isInside(idx,idy-1,idz))
S = 0.0; value.south = 0.0;
if(!isInside(idx,idy,idz-1)) if(!isInside(idx,idy,idz-1))
F = 0.0; value.front = 0.0;
if(!isInside(idx,idy,idz+1)) if(!isInside(idx,idy,idz+1))
B = 0.0; value.back = 0.0;
} }
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W, void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
double& E, double& S, double& N, double& F, StencilValue_t& value, double &scaleFactor)
double& B, double& C, double &scaleFactor)
{ {
scaleFactor = 1; scaleFactor = 1;
...@@ -449,7 +446,7 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W, ...@@ -449,7 +446,7 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W,
double dy_s=hr[1]; double dy_s=hr[1];
double dz_f=hr[2]; double dz_f=hr[2];
double dz_b=hr[2]; double dz_b=hr[2];
C = 0.0; value.center = 0.0;
std::tuple<int, int, int> coordxyz(idx, idy, idz); std::tuple<int, int, int> coordxyz(idx, idy, idz);
...@@ -467,29 +464,29 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W, ...@@ -467,29 +464,29 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W,
dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz); dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
if(dx_w != 0) if(dx_w != 0)
W = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w; value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
else else
W = 0; value.west = 0;
if(dx_e != 0) if(dx_e != 0)
E = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e; value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
else else
E = 0; value.east = 0;
if(dy_n != 0) if(dy_n != 0)
N = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n; value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
else else
N = 0; value.north = 0;
if(dy_s != 0) if(dy_s != 0)
S = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s; value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
else else
S = 0; value.south = 0;
if(dz_f != 0) if(dz_f != 0)
F = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f; value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
else else
F = 0; value.front = 0;
if(dz_b != 0) if(dz_b != 0)
B = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b; value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
else else
B = 0; value.back = 0;
//RHS scaleFactor for current 3D index //RHS scaleFactor for current 3D index
//0.5* comes from discretiztaion //0.5* comes from discretiztaion
...@@ -514,44 +511,43 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W, ...@@ -514,44 +511,43 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W,
if(dy_s == 0) if(dy_s == 0)
m2 = dy_n; m2 = dy_n;
C = 2 / hr[2]; value.center = 2 / hr[2];
if(dx_w != 0 || dx_e != 0) if(dx_w != 0 || dx_e != 0)
C *= (dx_w + dx_e); value.center *= (dx_w + dx_e);
if(dy_n != 0 || dy_s != 0) if(dy_n != 0 || dy_s != 0)
C *= (dy_n + dy_s); value.center *= (dy_n + dy_s);
if(dx_w != 0 || dx_e != 0) if(dx_w != 0 || dx_e != 0)
C += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1; value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
if(dy_n != 0 || dy_s != 0) if(dy_n != 0 || dy_s != 0)
C += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2; value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
} }
void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W, void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, StencilIndex_t& index)
int &E, int &S, int &N, int &F, int &B)
{ {
W = getIdx(idx - 1, idy, idz); index.west = getIdx(idx - 1, idy, idz);
E = getIdx(idx + 1, idy, idz); index.east = getIdx(idx + 1, idy, idz);
N = getIdx(idx, idy + 1, idz); index.north = getIdx(idx, idy + 1, idz);
S = getIdx(idx, idy - 1, idz); index.south = getIdx(idx, idy - 1, idz);
F = getIdx(idx, idy, idz - 1); index.front = getIdx(idx, idy, idz - 1);
B = getIdx(idx, idy, idz + 1); index.back = getIdx(idx, idy, idz + 1);
if(!isInside(idx+1,idy,idz)) if(!isInside(idx+1,idy,idz))
E = -1; index.east = -1;
if(!isInside(idx-1,idy,idz)) if(!isInside(idx-1,idy,idz))
W = -1; index.west = -1;
if(!isInside(idx,idy+1,idz)) if(!isInside(idx,idy+1,idz))
N = -1; index.north = -1;
if(!isInside(idx,idy-1,idz)) if(!isInside(idx,idy-1,idz))
S = -1; index.south = -1;
if(!isInside(idx,idy,idz-1)) if(!isInside(idx,idy,idz-1))
F = -1; index.front = -1;
if(!isInside(idx,idy,idz+1)) if(!isInside(idx,idy,idz+1))
B = -1; index.back = -1;
} }
......
...@@ -43,6 +43,8 @@ class BoundaryGeometry; ...@@ -43,6 +43,8 @@ class BoundaryGeometry;
class ArbitraryDomain : public IrregularDomain { class ArbitraryDomain : public IrregularDomain {
public: public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
std::string interpl); std::string interpl);
...@@ -54,13 +56,11 @@ public: ...@@ -54,13 +56,11 @@ public:
~ArbitraryDomain(); ~ArbitraryDomain();
/// returns discretization at (x,y,z) /// returns discretization at (x,y,z)
void getBoundaryStencil(int idx, int idy, int idz, double &W, double &E, void getBoundaryStencil(int idx, int idy, int idz, StencilValue_t& value,
double &S, double &N, double &F, double &B, double &C,
double &scaleFactor); double &scaleFactor);
/// returns index of neighbours at (x,y,z) /// returns index of neighbours at (x,y,z)
void getNeighbours(int idx, int idy, int idz, int &W, void getNeighbours(int idx, int idy, int idz, StencilIndex_t& index);
int &E, int &S, int &N, int &F, int &B);
/// returns type of boundary condition /// returns type of boundary condition
std::string getType() {return "Geometric";} std::string getType() {return "Geometric";}
...@@ -158,17 +158,14 @@ private: ...@@ -158,17 +158,14 @@ private:
} }
// Different interpolation methods for boundary points // Different interpolation methods for boundary points
void constantInterpolation(int idx, int idy, int idz, double &W, double &E, void constantInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &S, double &N, double &F, double &B, double &scaleFactor);
double &C, double &scaleFactor);
void linearInterpolation(int idx, int idy, int idz, double &W, double &E, void linearInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &S, double &N, double &F, double &B, double &scaleFactor);
double &C, double &scaleFactor);
void quadraticInterpolation(int idx, int idy, int idz, double &W, double &E, void quadraticInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &S, double &N, double &F, double &B, double &scaleFactor);
double &C, double &scaleFactor);
// Rotate positive axes with quaternion -DW // Rotate positive axes with quaternion -DW
inline void rotateWithQuaternion(Vector_t &v, Quaternion_t const quaternion); inline void rotateWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
......
...@@ -155,83 +155,84 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){ ...@@ -155,83 +155,84 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
} }
void BoxCornerDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) { void BoxCornerDomain::getBoundaryStencil(int x, int y, int z, StencilValue_t& value, double &scaleFactor) {
// determine which interpolation method we use for points near the boundary // determine which interpolation method we use for points near the boundary
switch(interpolationMethod) { switch(interpolationMethod) {
case CONSTANT: case CONSTANT:
constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor); constantInterpolation(x, y, z, value, scaleFactor);
break; break;
case LINEAR: case LINEAR:
linearInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor); linearInterpolation(x, y, z, value, scaleFactor);
break; break;
case QUADRATIC: case QUADRATIC:
quadraticInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor); quadraticInterpolation(x, y, z, value, scaleFactor);
break; break;
} }
// stencil center value has to be positive! // stencil center value has to be positive!
assert(C > 0); assert(value.center > 0);
} }
void BoxCornerDomain::getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) { void BoxCornerDomain::getNeighbours(int x, int y, int z, StencilIndex_t& index) {
if(x > 0) if(x > 0)
W = getIdx(x - 1, y, z); index.west = getIdx(x - 1, y, z);
else else
W = -1; index.west = -1;
if(x < nr[0] - 1) if(x < nr[0] - 1)
E = getIdx(x + 1, y, z); index.east = getIdx(x + 1, y, z);
else else
E = -1; index.east = -1;
if(y < nr[1] - 1) if(y < nr[1] - 1)
N = getIdx(x, y + 1, z); index.north = getIdx(x, y + 1, z);
else else
N = -1; index.north = -1;
if(y > 0) if(y > 0)
S = getIdx(x, y - 1, z); index.south = getIdx(x, y - 1, z);
else else
S = -1; index.south = -1;
if(z > 0) if(z > 0)
F = getIdx(x, y, z - 1); index.front = getIdx(x, y, z - 1);
else else
F = -1; index.front = -1;
if(z < nr[2] - 1) if(z < nr[2] - 1)
B = getIdx(x, y, z + 1); index.back = getIdx(x, y, z + 1);
else else
B = -1; index.back = -1;
} }
void BoxCornerDomain::constantInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) { void BoxCornerDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
{
scaleFactor = 1.0; scaleFactor = 1.0;
W = -1 / hr[0] * 1 / hr[0]; value.west = -1 / hr[0] * 1 / hr[0];
E = -1 / hr[0] * 1 / hr[0]; value.east = -1 / hr[0] * 1 / hr[0];
N = -1 / hr[1] * 1 / hr[1]; value.north = -1 / hr[1] * 1 / hr[1];
S = -1 / hr[1] * 1 / hr[1]; value.south = -1 / hr[1] * 1 / hr[1];
F = -1 / hr[2] * 1 / hr[2]; value.front = -1 / hr[2] * 1 / hr[2];
B = -1 / hr[2] * 1 / hr[2]; value.back = -1 / hr[2] * 1 / hr[2];
C = 2 / hr[0] * 1 / hr[0] + 2 / hr[1] * 1 / hr[1] + 2 / 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];
// we are a right boundary point // we are a right boundary point
if(!isInside(x + 1, y, z)) if(!isInside(x + 1, y, z))
E = 0.0; value.east = 0.0;
// we are a left boundary point // we are a left boundary point
if(!isInside(x - 1, y, z)) if(!isInside(x - 1, y, z))
W = 0.0; value.west = 0.0;
// we are a upper boundary point // we are a upper boundary point
if(!isInside(x, y + 1, z)) if(!isInside(x, y + 1, z))
N = 0.0; value.north = 0.0;
// we are a lower boundary point // we are a lower boundary point
if(!isInside(x, y - 1, z)) if(!isInside(x, y - 1, z))
S = 0.0; value.south = 0.0;
if(z == 1 || z == nr[2] - 2) { if(z == 1 || z == nr[2] - 2) {
...@@ -240,9 +241,9 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, double &W, doub ...@@ -240,9 +241,9 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, double &W, doub
// IFF: this values should not matter because they // IFF: this values should not matter because they
// never make it into the discretization matrix // never make it into the discretization matrix
if(z == 1) if(z == 1)
F = 0.0; value.front = 0.0;
else else
B = 0.0; value.back = 0.0;
// add contribution of Robin discretization to center point // add contribution of Robin discretization to center point
// d the distance between the center of the bunch and the boundary // d the distance between the center of the bunch and the boundary
...@@ -251,22 +252,23 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, double &W, doub ...@@ -251,22 +252,23 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, double &W, doub
//double cz = hr[2]*(nr[2]-1); //double cz = hr[2]*(nr[2]-1);
//double d = sqrt(cx*cx+cy*cy+cz*cz); //double d = sqrt(cx*cx+cy*cy+cz*cz);
double d = hr[2] * (nr[2] - 1) / 2; double d = hr[2] * (nr[2] - 1) / 2;
C += 2 / (d * hr[2]); value.center += 2 / (d * hr[2]);
//C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]); //value.center += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
// scale all stencil-points in z-plane with 0.5 (Robin discretization) // scale all stencil-points in z-plane with 0.5 (Robin discretization)
W /= 2.0; value.west /= 2.0;
E /= 2.0; value.east /= 2.0;
N /= 2.0; value.north /= 2.0;
S /= 2.0; value.south /= 2.0;
C /= 2.0; value.center /= 2.0;
scaleFactor *= 0.5; scaleFactor *= 0.5;
} }
} }
void BoxCornerDomain::linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) { void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
{
scaleFactor = 1.0; scaleFactor = 1.0;
double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0; double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
...@@ -296,51 +298,51 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, double &W, double ...@@ -296,51 +298,51 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, double &W, double
double de = hr[0]; double de = hr[0];
double dn = hr[1]; double dn = hr[1];
double ds = hr[1]; double ds = hr[1];
C = 0.0; value.center = 0.0;
//we are a right boundary point //we are a right boundary point
if(!isInside(x + 1, y, z)) { if(!isInside(x + 1, y, z)) {
double dx = getXIntersection(cx, z); double dx = getXIntersection(cx, z);
C += 1 / ((dx - cx) * de); value.center += 1 / ((dx - cx) * de);
E = 0.0; value.east = 0.0;
} else { } else {
C += 1 / (de * de); value.center += 1 / (de * de);
E = -1 / (de * de); value.east = -1 / (de * de);
} }
//we are a left boundary point //we are a left boundary point
if(!isInside(x - 1, y, z)) { if(!isInside(x - 1, y, z)) {
double dx = getXIntersection(cx, z); double dx = getXIntersection(cx, z);
C += 1 / ((std::abs(dx) - std::abs(cx)) * dw); value.center += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
W = 0.0; value.west = 0.0;
} else { } else {
C += 1 / (dw * dw); value.center += 1 / (dw * dw);
W = -1 / (dw * dw); value.west = -1 / (dw * dw);
} }
//we are a upper boundary point //we are a upper boundary point
if(!isInside(x, y + 1, z)) {