1 #ifdef HAVE_SAAMG_SOLVER
19 BoxCornerDomain::BoxCornerDomain(
double A,
double B,
double C,
double Length,
double L1,
double L2,
Vector_t nr,
Vector_t hr, std::string interpl) {
30 if(interpl ==
"CONSTANT")
31 interpolationMethod = CONSTANT;
32 else if(interpl ==
"LINEAR")
33 interpolationMethod = LINEAR;
34 else if(interpl ==
"QUADRATIC")
35 interpolationMethod = QUADRATIC;
38 *gmsg <<
" Write BoxCorner data to file boxcorner.dat" <<
endl;
39 std::string file(
"boxcorner.dat");
40 os_m.open(file.c_str());
42 *gmsg <<
"Unable to open output file " << file <<
endl;
48 BoxCornerDomain::~BoxCornerDomain() {
58 void BoxCornerDomain::compute(
Vector_t hr){
67 hasGeometryChanged_m =
true;
69 double bL= getB(getMinZ());
70 double bH= getB(getMaxZ());
75 INFOMSG(
" BoxCorner L= " << Length_m <<
" L1= " << L1_m <<
" L2= " << L2_m <<
" A= " << A_m <<
" B= " << B_m <<
" C= " << C_m
76 <<
" bL= " << bL <<
" bH= " << bH <<
" actBMin= " << actBMin_m <<
" actBMax=max(bL,bH)= " << actBMax_m << endl);
84 IntersectYDir.clear();
85 IntersectXDir.clear();
90 for(x = 0; x < nr[0]; x++) {
91 for(y = 0; y < nr[1]; y++) {
92 for(z = 0; z < nr[2]; z++) {
93 if(isInside(x, y, z)) {
94 IdxMap[toCoordIdx(x, y, z)] = idx;
95 CoordMap[idx++] = toCoordIdx(x, y, z);
136 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) {
139 switch(interpolationMethod) {
141 constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
144 linearInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
147 quadraticInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
155 void BoxCornerDomain::getBoundaryStencil(
int idx,
double &W,
double &E,
double &S,
double &N,
double &F,
double &B,
double &C,
double &scaleFactor) {
158 int x = 0, y = 0, z = 0;
159 getCoord(idx, x, y, z);
160 getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
164 void BoxCornerDomain::getNeighbours(
int idx,
int &W,
int &E,
int &S,
int &N,
int &F,
int &B) {
166 int x = 0, y = 0, z = 0;
167 getCoord(idx, x, y, z);
168 getNeighbours(x, y, z, W, E, S, N, F, B);
172 void BoxCornerDomain::getNeighbours(
int x,
int y,
int z,
int &W,
int &E,
int &S,
int &N,
int &F,
int &B) {
175 W = getIdx(x - 1, y, z);
179 E = getIdx(x + 1, y, z);
184 N = getIdx(x, y + 1, z);
188 S = getIdx(x, y - 1, z);
193 F = getIdx(x, y, z - 1);
197 B = getIdx(x, y, z + 1);
203 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) {
207 W = -1 / hr[0] * 1 / hr[0];
208 E = -1 / hr[0] * 1 / hr[0];
209 N = -1 / hr[1] * 1 / hr[1];
210 S = -1 / hr[1] * 1 / hr[1];
211 F = -1 / hr[2] * 1 / hr[2];
212 B = -1 / hr[2] * 1 / hr[2];
213 C = 2 / hr[0] * 1 / hr[0] + 2 / hr[1] * 1 / hr[1] + 2 / hr[2] * 1 / hr[2];
216 if(!isInside(x + 1, y, z))
220 if(!isInside(x - 1, y, z))
224 if(!isInside(x, y + 1, z))
228 if(!isInside(x, y - 1, z))
231 if(z == 1 || z == nr[2] - 2) {
248 double d = hr[2] * (nr[2] - 1) / 2;
249 C += 2 / (d * hr[2]);
263 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) {
267 double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
268 double cy = y * hr[1] - (nr[1] - 1) * hr[1] / 2.0;
297 if(!isInside(x + 1, y, z)) {
298 double dx = getXIntersection(cx, z);
299 C += 1 / ((dx - cx) * de);
307 if(!isInside(x - 1, y, z)) {
308 double dx = getXIntersection(cx, z);
317 if(!isInside(x, y + 1, z)) {
318 double dy = getYIntersection(cy, z);
319 C += 1 / ((dy - cy) * dn);
327 if(!isInside(x, y - 1, z)) {
328 double dy = getYIntersection(cy, z);
336 F = -1 / (hr[2] * hr[2]);
337 B = -1 / (hr[2] * hr[2]);
338 C += 2 / (hr[2] * hr[2]);
341 if(z == 0 || z == nr[2] - 1) {
351 double d = hr[2] * (nr[2] - 1) / 2;
352 C += 2 / (d * hr[2]);
366 void BoxCornerDomain::quadraticInterpolation(
int x,
int y,
int z,
double &W,
double &E,
double &S,
double &N,
double &F,
double &B,
double &C,
double &scaleFactor) {
368 double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
369 double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
371 double dx = getXIntersection(cx, z);
372 double dy = getYIntersection(cy, z);
434 if(!isInside(x + 1, y, z)) {
441 if(!isInside(x - 1, y, z)) {
448 if(!isInside(x, y + 1, z)) {
455 if(!isInside(x, y - 1, z)) {
461 W *= -1.0 / (dw * (dw + de));
462 E *= -1.0 / (de * (dw + de));
463 N *= -1.0 / (dn * (dn + ds));
464 S *= -1.0 / (ds * (dn + ds));
465 F = -1 / (hr[2] * (hr[2] + hr[2]));
466 B = -1 / (hr[2] * (hr[2] + hr[2]));
470 C += 1 / de * 1 / dw;
471 C += 1 / dn * 1 / ds;
472 C += 1 / hr[2] * 1 / hr[2];
568 if(z == 0 || z == nr[2] - 1) {
579 double d = hr[2] * (nr[2] - 1) / 2;
580 C += 2 / (d * hr[2]);
593 #endif //#ifdef HAVE_SAAMG_SOLVER
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Inform & endl(Inform &inf)