Commit 03d4f489 authored by Tuelin Kaman's avatar Tuelin Kaman
Browse files

SAAMG test cases for elliptical and arbitrary geometries; parallel decomposition and mapping

parent e49b8232
......@@ -824,8 +824,8 @@ void PartBunch::resizeMesh() {
double xmax = fs_m->solver_m->getXRangeMax();
double ymin = fs_m->solver_m->getYRangeMin();
double ymax = fs_m->solver_m->getYRangeMax();
double zmin = fs_m->solver_m->getZRangeMin();
double zmax = fs_m->solver_m->getZRangeMax();
double zmin = rmin_m[2]; //fs_m->solver_m->getZRangeMin();
double zmax = rmax_m[2]; //fs_m->solver_m->getZRangeMax();
if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
ymin > rmin_m[1] || ymax < rmax_m[1]) {
......@@ -846,9 +846,9 @@ void PartBunch::resizeMesh() {
boundp();
get_bounds(rmin_m, rmax_m);
}
Vector_t mymin = Vector_t(xmin, ymin , zmin);
Vector_t mymax = Vector_t(xmax, ymax , zmax);
// extend domain with extra "ghost" point
Vector_t mymin = Vector_t(xmin, ymin , zmin-hr_m[2]);
Vector_t mymax = Vector_t(xmax, ymax , zmax+hr_m[2]);
for(int i = 0; i < 3; i++)
hr_m[i] = (mymax[i] - mymin[i]) / nr_m[i];
......@@ -877,8 +877,8 @@ void PartBunch::computeSelfFields() {
if(fs_m->hasValidSolver()) {
// if(fs_m->getFieldSolverType() == "SAAMG" && !isGridFixed())
// resizeMesh();
if (fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
INFOMSG("after resizeMesh" << hr_m << endl);
//scatter charges onto grid
......@@ -926,7 +926,9 @@ void PartBunch::computeSelfFields() {
INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
#endif
// charge density is in rho_m
IpplTimings::startTimer(compPotenTimer_m);
fs_m->solver_m->computePotential(rho_m, hr_scaled);
IpplTimings::stopTimer(compPotenTimer_m);
//do the multiplication of the grid-cube volume coming
//from the discretization of the convolution integral.
......@@ -1126,7 +1128,9 @@ void PartBunch::computeSelfFields_cycl(double gamma, Vector_t const meanR, Vekto
/// now charge density is in rho_m
/// calculate Possion equation (without coefficient: -1/(eps))
IpplTimings::startTimer(compPotenTimer_m);
fs_m->solver_m->computePotential(rho_m, hr_scaled);
IpplTimings::stopTimer(compPotenTimer_m);
/// additional work of FFT solver
/// now the scalar potential is given back in rho_m
......
......@@ -1718,7 +1718,7 @@ void ParallelCyclotronTracker::Tracker_RK4() {
// main integration loop
*gmsg << "* ---------------------------- Start tracking ----------------------------" << endl;
for(; step_m < maxSteps_m; step_m++) {
// *gmsg << "step_m=" << step_m << endl;
*gmsg << "step_m= " << step_m << endl;
bool dumpEachTurn = false;
if(initialTotalNum_m > 2) {
// single particle dumping
......@@ -4060,9 +4060,9 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
globalToLocal(itsBunch->R, phi, meanR);
itsBunch->R *= Vector_t(0.001); // mm --> m
itsBunch->boundp();
checkNumPart(std::string("* Before repartation: "));
checkNumPart(std::string("* Before repartition: "));
repartition();
checkNumPart(std::string("* After repartation: "));
checkNumPart(std::string("* After repartition: "));
itsBunch->R *= Vector_t(1000.0); // m --> mm
localToGlobal(itsBunch->R, phi, meanR);
}
......
......@@ -11,7 +11,7 @@
// $Author: kaman $
// $Date: 2014 $
// ------------------------------------------------------------------------
#define DEBUG_INTERSECT_RAY_BOUNDARY
//#define DEBUG_INTERSECT_RAY_BOUNDARY
#ifdef HAVE_SAAMG_SOLVER
#include <map>
......@@ -19,6 +19,10 @@
#include <iostream>
#include <assert.h>
#include <math.h>
#define MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define MAX2(a,b) (((a) > (b)) ? (a) : (b))
#include "ArbitraryDomain.h"
ArbitraryDomain::ArbitraryDomain(
......@@ -28,8 +32,8 @@ ArbitraryDomain::ArbitraryDomain(
std::string interpl) {
bgeom_m = bgeom;
Geo_mincoords_m = bgeom->getmincoords();
Geo_maxcoords_m = bgeom->getmaxcoords();
intersectMinCoords_m = bgeom->getmincoords();
intersectMaxCoords_m = bgeom->getmaxcoords();
setNr(nr);
setHr(hr);
......@@ -57,7 +61,7 @@ void ArbitraryDomain::Compute(Vector_t hr) {
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
setHr(hr);
setNr(nr);
int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr[2] - 1) ? 0 : 1;
int yGhostOffsetLeft = (localId[1].first()== 0) ? 0 : 1;
......@@ -224,7 +228,6 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion){
setHr(hr);
globalMeanR_m = globalMeanR;
globalToLocalQuaternion_m = globalToLocalQuaternion;
localToGlobalQuaternion_m[0] = globalToLocalQuaternion[0];
......@@ -249,7 +252,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
//calculate intersection
Vector_t P, saveP, dir, I;
Vector_t P0 = Vector_t(0,0,Geo_mincoords_m[2]+hr[2]); //Reference Point inside the boundary
Vector_t P0 = Vector_t(0,0,bgeom_m->getmincoords()[2]+hr[2]); //Reference Point inside the boundary
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
......@@ -268,6 +271,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
setYRangeMax(MIN2(intersectMaxCoords_m(2),I[2]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
......@@ -278,9 +282,10 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
setYRangeMin(MAX2(intersectMinCoords_m(2),I[2]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
......@@ -289,6 +294,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
setZRangeMax(MIN2(intersectMaxCoords_m(1),I[1]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
......@@ -299,6 +305,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
setZRangeMin(MAX2(intersectMinCoords_m(1),I[1]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
......@@ -310,6 +317,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
setXRangeMax(MIN2(intersectMaxCoords_m(0),I[0]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
......@@ -320,6 +328,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)){
setXRangeMin(MAX2(intersectMinCoords_m(0),I[0]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
......@@ -328,134 +337,43 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
*gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
}
} else
continue;
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
}
}
}
}
//number of ghost nodes to the right
int znumGhostNodesRight = 0;
if (zGhostOffsetRight != 0) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
if (isInside(idx, idy, localId[2].last() + zGhostOffsetRight))
znumGhostNodesRight++;
}
}
}
//number of ghost nodes to the left
int znumGhostNodesLeft = 0;
if (zGhostOffsetLeft != 0) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
if (isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
znumGhostNodesLeft++;
}
}
}
//number of ghost nodes to the right
int ynumGhostNodesRight = 0;
if (yGhostOffsetRight != 0) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
if (isInside(idx, localId[1].last() + yGhostOffsetRight, idz))
ynumGhostNodesRight++;
}
}
}
//number of ghost nodes to the left
int ynumGhostNodesLeft = 0;
if (yGhostOffsetLeft != 0) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
if (isInside(idx, localId[1].first() - yGhostOffsetLeft, idz))
ynumGhostNodesLeft++;
}
}
}
//number of ghost nodes to the right
int xnumGhostNodesRight = 0;
if (xGhostOffsetRight != 0) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
if (isInside(localId[0].last() + xGhostOffsetRight, idy, idz))
xnumGhostNodesRight++;
}
}
}
//number of ghost nodes to the left
int xnumGhostNodesLeft = 0;
if (xGhostOffsetLeft != 0) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
if (isInside(localId[0].first() - xGhostOffsetLeft, idy, idz))
xnumGhostNodesLeft++;
}
}
}
//xy points in z plane
int numxy;
int numtotalxy = 0;
numXY.clear();
for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
numxy =0;
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx =localId[0].first(); idx <= localId[0].last(); idx++) {
if (isInside(idx, idy, idz))
numxy++;
}
}
numtotalxy += numxy;
}
startId = 0;
MPI_Scan(&numtotalxy, &startId, 1, MPI_INTEGER, MPI_SUM, Ippl::getComm());
startId -= numtotalxy;
//build up index and coord map
IdxMap.clear();
CoordMap.clear();
register int id = startId - xnumGhostNodesLeft - ynumGhostNodesLeft - znumGhostNodesLeft;
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
if (isInside(idx, idy, idz)) {
IdxMap[toCoordIdx(idx, idy, idz)] = id;
CoordMap[id++] = toCoordIdx(idx, idy, idz);
}
}
register int id=0;
for (int idz = 0; idz < nr[2]; idz++) {
for (int idy = 0; idy < nr[1]; idy++) {
for (int idx = 0; idx < nr[0]; idx++) {
IdxMap[toCoordIdx(idx, idy, idz)] = id;
CoordMap[id++] = toCoordIdx(idx, idy, idz);
}
}
}
}
// 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[1] + idy) * nr[0] + idx;
}
// Conversion from (x,y,z) to index on the 3D grid
int ArbitraryDomain::getIdx(int idx, int idy, int idz) {
if (isInside(idx, idy, idz) && idx>=0 && idy >=0 && idz >=0 )
return IdxMap[toCoordIdx(idx, idy, idz)];
else
return -1;
return IdxMap[toCoordIdx(idx, idy, 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];
......@@ -464,27 +382,15 @@ inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) {
}
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
/* Expensive computation to check */
Vector_t P0, P;
Vector_t P;
P0 = Vector_t(0, 0, Geo_mincoords_m[2]+hr[2]); //Reference Point inside the boundary
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[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];
rotateWithQuaternion(P, localToGlobalQuaternion_m);
P += globalMeanR_m;
return (bgeom_m->fastIsInside(P0, P) % 2 == 0);
/*
bool ret = false;
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];
int countH, countL;
int countH, countL;
std::multimap < std::tuple<int, int, int>, double >::iterator itrH, itrL;
std::tuple<int, int, int> coordxyz(idx, idy, idz);
//check if z is inside with x,y coords
......@@ -493,9 +399,8 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
countH = IntersectHiZ.count(coordxyz);
countL = IntersectLoZ.count(coordxyz);
if(countH == 1 && countL == 1)
ret = (cz <= itrH->second) && (cz >= itrL->second);
ret = (P[2] <= itrH->second) && (P[2] >= itrL->second);
//check if y is inside with x,z coords
itrH = IntersectHiY.find(coordxyz);
......@@ -503,9 +408,8 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
countH = IntersectHiY.count(coordxyz);
countL = IntersectLoY.count(coordxyz);
if(countH == 1 && countL == 1)
ret = ret && (cy <= itrH->second) && (cy >= itrL->second);
ret = ret && (P[1] <= itrH->second) && (P[1] >= itrL->second);
//check if x is inside with y,z coords
itrH = IntersectHiX.find(coordxyz);
......@@ -513,12 +417,10 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
countH = IntersectHiX.count(coordxyz);
countL = IntersectLoX.count(coordxyz);
if(countH == 1 && countL == 1)
ret = ret && (cx <= itrH->second) && (cx >= itrL->second);
ret = ret && (P[0] <= itrH->second) && (P[0] >= itrL->second);
return ret;
*/
}
int ArbitraryDomain::getNumXY(int z) {
......@@ -543,10 +445,10 @@ void ArbitraryDomain::getBoundaryStencil(int idx, int idy, int idz, double &W, d
ConstantInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor);
break;
case LINEAR:
// LinearInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor);
// LinearInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor);
break;
case QUADRATIC:
// QuadraticInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor);
// QuadraticInterpolation(idx,idy,idz,W,E,S,N,F,B,C,scaleFactor);
break;
}
......@@ -564,25 +466,22 @@ void ArbitraryDomain::ConstantInterpolation(int idx, int idy, int idz, double& W
B = -1/(hr[2]*hr[2]);
C = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
if(!isInside(idx+1,idy,idz))
E = 0.0;
if(!isInside(idx-1,idy,idz))
W = 0.0;
if(!isInside(idx+1,idy,idz))
E = 0.0;
if(!isInside(idx,idy+1,idz))
N = 0.0;
if(!isInside(idx,idy-1,idz))
S = 0.0;
if(!isInside(idx,idy,idz-1))
F = 0.0;
if(!isInside(idx,idy,idz+1))
B = 0.0;
}
/*
void ArbitraryDomain::LinearInterpolation(int idx, int idy, int idz, double& W, double& E, double& S, double& N, double& F, double& B, double& C, double &scaleFactor)
{
......@@ -729,35 +628,30 @@ void ArbitraryDomain::getNeighbours(int id, int &W, int &E, int &S, int &N, int
void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W, int &E, int &S, int &N, int &F, int &B) {
if(idx > 0)
W = getIdx(idx - 1, idy, idz);
else
W = -1;
W = getIdx(idx - 1, idy, idz);
E = getIdx(idx + 1, idy, idz);
N = getIdx(idx, idy + 1, idz);
S = getIdx(idx, idy - 1, idz);
F = getIdx(idx, idy, idz - 1);
B = getIdx(idx, idy, idz + 1);
if(idx < nr[0] - 1)
E = getIdx(idx + 1, idy, idz);
else
if(!isInside(idx+1,idy,idz))
E = -1;
if(idy < nr[1] - 1)
N = getIdx(idx, idy + 1, idz);
else
if(!isInside(idx-1,idy,idz))
W = -1;
if(!isInside(idx,idy+1,idz))
N = -1;
if(idy > 0)
S = getIdx(idx, idy - 1, idz);
else
if(!isInside(idx,idy-1,idz))
S = -1;
if(!isInside(idx,idy,idz-1))
F = -1;
if(idz > 0)
F = getIdx(idx, idy, idz - 1);
else
F = -1;
if(idz < nr[2] - 1)
B = getIdx(idx, idy, idz + 1);
else
B = -1;
if(!isInside(idx,idy,idz+1))
B = -1;
}
......
......@@ -47,12 +47,19 @@ public:
int getStartId() {return startId;}
double getXRangeMin() { return Geo_mincoords_m(0); }
double getXRangeMax() { return Geo_maxcoords_m(0); }
double getYRangeMin() { return Geo_mincoords_m(1); }
double getYRangeMax() { return Geo_maxcoords_m(1); }
double getZRangeMin() { return Geo_mincoords_m(2); }
double getZRangeMax() { return Geo_maxcoords_m(2); }
double getXRangeMin(){ return intersectMinCoords_m(0); }
double getXRangeMax(){ return intersectMaxCoords_m(0); }
double getYRangeMin(){ return intersectMinCoords_m(1); }
double getYRangeMax(){ return intersectMaxCoords_m(1); }
double getZRangeMin(){ return intersectMinCoords_m(2); }
double getZRangeMax(){ return intersectMaxCoords_m(2); }
void setXRangeMin(double xmin){ intersectMinCoords_m(0) = xmin; }
void setXRangeMax(double xmax){ intersectMaxCoords_m(0) = xmax; }
void setYRangeMin(double ymin){ intersectMinCoords_m(1) = ymin; }
void setYRangeMax(double ymax){ intersectMaxCoords_m(1) = ymax; }
void setZRangeMin(double zmin){ intersectMinCoords_m(2) = zmin; }
void setZRangeMax(double zmax){ intersectMaxCoords_m(2) = zmax; }
bool hasGeometryChanged() { return hasGeometryChanged_m; }
......@@ -99,6 +106,8 @@ private:
Vector_t Geo_hr_m;
Vector_t Geo_mincoords_m;
Vector_t Geo_maxcoords_m;
Vector_t intersectMinCoords_m;
Vector_t intersectMaxCoords_m;
// Conversion from (x,y,z) to index in xyz plane
inline int toCoordIdx(int idx, int idy, int idz);
......
......@@ -60,9 +60,7 @@ void EllipticDomain::Compute(Vector_t hr){
for(x = 0; x < nr[0]; x++) {
for(y = 0; y < nr[1]; y++) {
if(isInside(x, y, 1)) {
//IdxMap[toCoordIdx(x, y)] = idx++;
IdxMap[toCoordIdx(x, y)] = idx;
CoordMap[idx++] = toCoordIdx(x, y);
nxy_m++;
......@@ -116,10 +114,86 @@ void EllipticDomain::Compute(Vector_t hr){
}
void EllipticDomain::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]) {
hasGeometryChanged_m = false;
return;
}
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
IdxMap.clear();
CoordMap.clear();
//clear previous intersection points
IntersectYDir.clear();
IntersectXDir.clear();
// build a index and coordinate map
register int idx = 0;
register int x, y;
for(x = localId[0].first(); x<= localId[0].last(); x++) {
for(y = localId[1].first(); y <= localId[1].last(); y++) {
if(isInside(x, y, 1)) {
IdxMap[toCoordIdx(x, y)] = idx;
CoordMap[idx++] = toCoordIdx(x, y);
nxy_m++;
}
}
}
switch(interpolationMethod) {
case CONSTANT:
break;
case LINEAR:
case QUADRATIC:
double smajsq = SemiMajor * SemiMajor;
double sminsq = SemiMinor * SemiMinor;
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
for(x = localId[0].first(); x <= localId[0].last(); x++) {
pos = x * hr[0] - mx;
if (pos <= -SemiMajor || pos >= SemiMajor)
{
IntersectYDir.insert(std::pair<int, double>(x, 0));
IntersectYDir.insert(std::pair<int, double>(x, 0));
}else{
yd = std::abs(sqrt(sminsq</