Commit 559191f1 authored by frey_m's avatar frey_m
Browse files

EllipticDomain: clean up + function for Robin BC

parent 7a6247e0
......@@ -80,10 +80,9 @@ EllipticDomain::~EllipticDomain() {
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
//there is nothing to be done if the mesh spacings have not changed
// there is nothing to be done if the mesh spacings in x and y have not changed
if (hr[0] == getHr()[0] &&
hr[1] == getHr()[1] &&
hr[2] == getHr()[2])
hr[1] == getHr()[1])
{
hasGeometryChanged_m = false;
return;
......@@ -131,18 +130,16 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
double yd = 0.0;
double xd = 0.0;
double pos = 0.0;
double mx = (nr[0] - 1) * hr[0] / 2.0;
double my = (nr[1] - 1) * hr[1] / 2.0;
//calculate intersection with the ellipse
// calculate intersection with the ellipse
for(x = localId[0].first(); x <= localId[0].last(); x++) {
pos = x * hr[0] - mx;
pos = - semiMajor_m + hr[0] * (x + 0.5);
if (std::abs(pos) >= semiMajor_m)
{
intersectYDir_m.insert(std::pair<int, double>(x, 0));
intersectYDir_m.insert(std::pair<int, double>(x, 0));
}else{
yd = std::abs(std::sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
}else {
yd = std::abs(std::sqrt(sminsq - sminsq * pos * pos / smajsq));
intersectYDir_m.insert(std::pair<int, double>(x, yd));
intersectYDir_m.insert(std::pair<int, double>(x, -yd));
}
......@@ -150,13 +147,13 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
}
for(y = localId[0].first(); y < localId[1].last(); y++) {
pos = y * hr[1] - my;
pos = - semiMinor_m + hr[1] * (y + 0.5);
if (std::abs(pos) >= semiMinor_m)
{
intersectXDir_m.insert(std::pair<int, double>(y, 0));
intersectXDir_m.insert(std::pair<int, double>(y, 0));
}else{
xd = std::abs(std::sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
}else {
xd = std::abs(std::sqrt(smajsq - smajsq * pos * pos / sminsq));
intersectXDir_m.insert(std::pair<int, double>(y, xd));
intersectXDir_m.insert(std::pair<int, double>(y, -xd));
}
......@@ -169,9 +166,10 @@ void EllipticDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t&
{
Vector_t mymax = Vector_t(0.0, 0.0, 0.0);
// apply bounding box increment, i.e., "BBOXINCR" input argument
double zmin = std::signbit(rmin[2]) ? rmin[2] * (1.0 + dh) : rmin[2] * (1.0 - dh);
double zmax = std::signbit(rmax[2]) ? rmax[2] * (1.0 - dh) : rmax[2] * (1.0 + dh);
// apply bounding box increment dh, i.e., "BBOXINCR" input argument
double zsize = rmax[2] - rmin[2];
double zmin = rmin[2] - zsize * (1.0 + dh);
double zmax = rmax[2] + zsize * (1.0 + dh);
origin = Vector_t(-semiMajor_m, -semiMinor_m, zmin);
mymax = Vector_t( semiMajor_m, semiMinor_m, zmax);
......@@ -187,7 +185,7 @@ void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W,
{
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_m) {
case CONSTANT:
constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
......@@ -291,23 +289,8 @@ void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV,
if (!isInside(x, y - 1, z))
SV = 0.0;
if (z == 0 || z == nr[2] - 1) {
// case where we are on the Robin BC in Z-direction
// where we distinguish two cases
// IFF: this values should not matter because they
// never make it into the discretization matrix
if (z == 0)
FV = 0.0;
else
BV = 0.0;
// add contribution of Robin discretization to center point
// d the distance between the center of the bunch and the boundary
double d = hr[2] * (nr[2] - 1) / 2.0;
CV += 2.0 / (d * hr[2]);
}
// handle boundary condition in z direction
robinBoundaryStencil(z, FV, BV, CV);
}
void EllipticDomain::linearInterpolation(int x, int y, int z, double &W,
......@@ -317,12 +300,13 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, double &W,
{
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 = - semiMajor_m + hr[0] * (x + 0.5);
double cy = - semiMinor_m + hr[1] * (y + 0.5);
double dx = 0.0;
std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
if(cx < 0)
if (cx < 0)
++it;
dx = it->second;
......@@ -339,70 +323,47 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, double &W,
double ds = hr[1];
C = 0.0;
//we are a right boundary point
// we are a right boundary point
if (!isInside(x + 1, y, z)) {
C += 1 / ((dx - cx) * de);
C += 1.0 / ((dx - cx) * de);
E = 0.0;
} else {
C += 1 / (de * de);
E = -1 / (de * de);
C += 1.0 / (de * de);
E = -1.0 / (de * de);
}
//we are a left boundary point
// we are a left boundary point
if (!isInside(x - 1, y, z)) {
C += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
C += 1.0 / ((std::abs(dx) - std::abs(cx)) * dw);
W = 0.0;
} else {
C += 1 / (dw * dw);
W = -1 / (dw * dw);
C += 1.0 / (dw * dw);
W = -1.0 / (dw * dw);
}
//we are a upper boundary point
// we are a upper boundary point
if (!isInside(x, y + 1, z)) {
C += 1 / ((dy - cy) * dn);
C += 1.0 / ((dy - cy) * dn);
N = 0.0;
} else {
C += 1 / (dn * dn);
N = -1 / (dn * dn);
C += 1.0 / (dn * dn);
N = -1.0 / (dn * dn);
}
//we are a lower boundary point
// we are a lower boundary point
if (!isInside(x, y - 1, z)) {
C += 1 / ((std::abs(dy) - std::abs(cy)) * ds);
C += 1.0 / ((std::abs(dy) - std::abs(cy)) * ds);
S = 0.0;
} else {
C += 1 / (ds * ds);
S = -1 / (ds * ds);
C += 1.0 / (ds * ds);
S = -1.0 / (ds * ds);
}
F = -1 / (hr[2] * hr[2]);
B = -1 / (hr[2] * hr[2]);
C += 2 / (hr[2] * hr[2]);
// handle boundary condition in z direction
/*
if(z == 0 || z == nr[2] - 1) {
// case where we are on the NEUMAN BC in Z-direction
// where we distinguish two cases
if(z == 0)
F = 0.0;
else
B = 0.0;
//hr[2]*(nr2[2]-1)/2 = radius
double d = hr[2] * (nr[2] - 1) / 2;
C += 2 / (d * hr[2]);
W /= 2.0;
E /= 2.0;
N /= 2.0;
S /= 2.0;
C /= 2.0;
scaleFactor *= 0.5;
}
*/
F = -1.0 / (hr[2] * hr[2]);
B = -1.0 / (hr[2] * hr[2]);
C += 2.0 / (hr[2] * hr[2]);
robinBoundaryStencil(z, F, B, C);
}
void EllipticDomain::quadraticInterpolation(int x, int y, int z,
......@@ -410,6 +371,8 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z,
double &N, double &F, double &B, double &C,
double &scaleFactor)
{
scaleFactor = 1.0;
double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
......@@ -431,197 +394,80 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z,
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
W = 1.0;
E = 1.0;
N = 1.0;
S = 1.0;
F = 1.0;
B = 1.0;
W = 0.0;
E = 0.0;
S = 0.0;
N = 0.0;
C = 0.0;
//TODO: = cx+hr[0] > dx && cx > 0
//if((x-nr[0]/2.0+1)*hr[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) {
////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) {
////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) {
////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) {
//we are a right boundary point
// we are a right boundary point
if (!isInside(x + 1, y, z)) {
de = dx - cx;
double s = dx - cx;
C -= 2.0 * (s - 1.0) / (s * de);
E = 0.0;
W += (s - 1.0) / ((s + 1.0) * de);
} else {
C += 1.0 / (de * de);
E = -1.0 / (de * de);
}
//if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
//we are a left boundary point
// we are a left boundary point
if (!isInside(x - 1, y, z)) {
dw = std::abs(dx) - std::abs(cx);
double s = std::abs(dx) - std::abs(cx);
C -= 2.0 * (s - 1.0) / (s * de);
W = 0.0;
E += (s - 1.0) / ((s + 1.0) * de);
} else {
C += 1.0 / (dw * dw);
W = -1.0 / (dw * dw);
}
//if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
//we are a upper boundary point
// we are a upper boundary point
if (!isInside(x, y + 1, z)) {
dn = dy - cy;
double s = dy - cy;
C -= 2.0 * (s - 1.0) / (s * dn);
N = 0.0;
S += (s - 1.0) / ((s + 1.0) * dn);
} else {
C += 1.0 / (dn * dn);
N = -1.0 / (dn * dn);
}
//if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
//we are a lower boundary point
// we are a lower boundary point
if (!isInside(x, y - 1, z)) {
ds = std::abs(dy) - std::abs(cy);
double s = std::abs(dy) - std::abs(cy);
C -= 2.0 * (s - 1.0) / (s * dn);
S = 0.0;
N += (s - 1.0) / ((s + 1.0) * dn);
} else {
C += 1.0 / (ds * ds);
S = -1.0 / (ds * ds);
}
//2/dw*(dw_de)
W *= -1.0 / (dw * (dw + de));
E *= -1.0 / (de * (dw + de));
N *= -1.0 / (dn * (dn + ds));
S *= -1.0 / (ds * (dn + ds));
F = -1 / (hr[2] * (hr[2] + hr[2]));
B = -1 / (hr[2] * (hr[2] + hr[2]));
//TODO: problem when de,dw,dn,ds == 0
//is NOT a regular BOUND PT
C += 1 / de * 1 / dw;
C += 1 / dn * 1 / ds;
C += 1 / hr[2] * 1 / hr[2];
scaleFactor = 0.5;
//for regular gridpoints no problem with symmetry, just boundary
//z direction is right
//implement isLastInside(dir)
//we have LOCAL x,y coordinates!
/*
if(dw != 0 && !wIsB)
W = -1/dw * (dn+ds) * 2*hr[2];
else
W = 0;
if(de != 0 && !eIsB)
E = -1/de * (dn+ds) * 2*hr[2];
else
E = 0;
if(dn != 0 && !nIsB)
N = -1/dn * (dw+de) * 2*hr[2];
else
N = 0;
if(ds != 0 && !sIsB)
S = -1/ds * (dw+de) * 2*hr[2];
else
S = 0;
F = -(dw+de)*(dn+ds)/hr[2];
B = -(dw+de)*(dn+ds)/hr[2];
*/
//if(dw != 0)
//W = -2*hr[2]*(dn+ds)/dw;
//else
//W = 0;
//if(de != 0)
//E = -2*hr[2]*(dn+ds)/de;
//else
//E = 0;
//if(dn != 0)
//N = -2*hr[2]*(dw+de)/dn;
//else
//N = 0;
//if(ds != 0)
//S = -2*hr[2]*(dw+de)/ds;
//else
//S = 0;
//F = -(dw+de)*(dn+ds)/hr[2];
//B = -(dw+de)*(dn+ds)/hr[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]);
// catch the case where a point lies on the boundary
//FIXME: do this more elegant!
//double m1 = dw*de;
//double m2 = dn*ds;
//if(de == 0)
//m1 = dw;
//if(dw == 0)
//m1 = de;
//if(dn == 0)
//m2 = ds;
//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];
//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;
//if(dn != 0 || ds != 0)
//C += (2*hr[2])*(dw+de)*(dn+ds)/m2;
//handle Neumann case
//if(z == 0 || z == nr[2]-1) {
//if(z == 0)
//F = 0.0;
//else
//B = 0.0;
////neumann stuff
//W = W/2.0;
//E = E/2.0;
//N = N/2.0;
//S = S/2.0;
//C /= 2.0;
//scaleFactor /= 2.0;
//}
// handle boundary condition in z direction
F = -1.0 / (hr[2] * hr[2]);
B = -1.0 / (hr[2] * hr[2]);
C += 2.0 / (hr[2] * hr[2]);
robinBoundaryStencil(z, F, B, C);
}
void EllipticDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) {
if (z == 0 || z == nr[2] - 1) {
// case where we are on the NEUMAN BC in Z-direction
// case where we are on the Robin BC in Z-direction
// where we distinguish two cases
// IFF: this values should not matter because they
// never make it into the discretization matrix
if (z == 0)
F = 0.0;
else
B = 0.0;
//C += 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;
C += 2 / (d * hr[2]);
W /= 2.0;
E /= 2.0;
N /= 2.0;
S /= 2.0;
C /= 2.0;
scaleFactor /= 2.0;
// add contribution of Robin discretization to center point
// d the distance between the center of the bunch and the boundary
double d = 0.5 * hr[2] * (nr[2] - 1);
C += 2.0 / (d * hr[2]);
}
}
......
......@@ -75,8 +75,8 @@ public:
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z) {
double xx = - semiMajor_m + hr[0] * (x + 0.5); //(x - (nr[0] - 1) / 2.0) * hr[0];
double yy = - semiMinor_m + hr[1] * (y + 0.5); //(y - (nr[1] - 1) / 2.0) * hr[1];
double xx = - semiMajor_m + hr[0] * (x + 0.5);
double yy = - semiMinor_m + hr[1] * (y + 0.5);
bool isInsideEllipse = (xx * xx / (semiMajor_m * semiMajor_m) +
yy * yy / (semiMinor_m * semiMinor_m) < 1);
......@@ -176,6 +176,8 @@ private:
double &E, double &S, double &N,
double &F, double &B, double &C,
double &scaleFactor);
void robinBoundaryStencil(int z, double &F, double &B, double &C);
};
#endif //#ifdef ELLIPTICAL_DOMAIN_H
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment