Commit 8e0d416c authored by frey_m's avatar frey_m
Browse files

SAAMG: cleanup get*RangeMin() + get*RangeMax() functions

parent 35fd36af
......@@ -42,9 +42,11 @@ ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
Vector_t /*hr*/,
std::string interpl){
bgeom_m = bgeom;
minCoords_m = bgeom->getmincoords();
maxCoords_m = bgeom->getmaxcoords();
geomCentroid_m = (minCoords_m + maxCoords_m)/2.0;
setRangeMin(bgeom->getmincoords());
setRangeMax(bgeom->getmaxcoords());
geomCentroid_m = (min_m + max_m)/2.0;
bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m);
if (have_inside_pt == false) {
......@@ -54,7 +56,7 @@ ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
}
setNr(nr);
for(int i=0; i<3; i++)
Geo_hr_m[i] = (maxCoords_m[i] - minCoords_m[i])/nr[i];
Geo_hr_m[i] = (max_m[i] - min_m[i])/nr[i];
setHr(Geo_hr_m);
startId = 0;
......@@ -111,17 +113,17 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
//saveP_old[2] = (idz - (nr[2]-1)/2.0)*hr[2];
P[2] = minCoords_m[2] + (idz + 0.5) * 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];
P[1] = minCoords_m[1] + (idy + 0.5) * 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];
P[0] = minCoords_m[0] + (idx + 0.5) * hr[0];
P[0] = min_m[0] + (idx + 0.5) * hr[0];
// *gmsg << "Now working on point " << saveP << " (original was " << saveP_old << ")" << endl;
......@@ -326,9 +328,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] = minCoords_m[0] + (idx + 0.5) * hr[0];
P[1] = minCoords_m[1] + (idy + 0.5) * hr[1];
P[2] = minCoords_m[2] + (idz + 0.5) * hr[2];
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];
return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
}
......
......@@ -68,22 +68,6 @@ public:
int getStartId() {return startId;}
double getXRangeMin(){ return minCoords_m(0); }
double getYRangeMin(){ return minCoords_m(1); }
double getZRangeMin(){ return minCoords_m(2); }
double getXRangeMax(){ return maxCoords_m(0); }
double getYRangeMax(){ return maxCoords_m(1); }
double getZRangeMax(){ return maxCoords_m(2); }
void setXRangeMin(double xmin){ minCoords_m(0) = xmin; }
void setYRangeMin(double ymin){ minCoords_m(1) = ymin; }
void setZRangeMin(double zmin){ minCoords_m(2) = zmin; }
void setXRangeMax(double xmax){ maxCoords_m(0) = xmax; }
void setYRangeMax(double ymax){ maxCoords_m(1) = ymax; }
void setZRangeMax(double zmax){ maxCoords_m(2) = zmax; }
private:
BoundaryGeometry *bgeom_m;
......
......@@ -35,18 +35,14 @@
extern Inform *gmsg;
BoxCornerDomain::BoxCornerDomain(Vector_t nr, Vector_t hr) {
setNr(nr);
setHr(hr);
}
BoxCornerDomain::BoxCornerDomain(double A, double B, double C, double Length, double L1, double L2, Vector_t nr, Vector_t hr, std::string interpl) {
A_m = A;
B_m = B;
BoxCornerDomain::BoxCornerDomain(double A, double B, double C, double length,
double L1, double L2, Vector_t nr, Vector_t hr,
std::string interpl)
{
setRangeMin(Vector_t(-A, -B, L1));
setRangeMax(Vector_t( A, B, L1 + L2));
C_m = C;
Length_m = Length;
L1_m = L1;
L2_m = L2;
length_m = length;
setNr(nr);
setHr(hr);
......@@ -93,12 +89,9 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
double bL= getB(getMinZ());
double bH= getB(getMaxZ());
actBMin_m = -B_m;
actBMin_m = getYRangeMin();
actBMax_m = std::max(bL,bH);
INFOMSG(" BoxCorner L= " << Length_m << " L1= " << L1_m << " L2= " << L2_m << " A= " << A_m << " B= " << B_m << " C= " << C_m
<< " bL= " << bL << " bH= " << bH << " actBMin= " << actBMin_m << " actBMax=max(bL,bH)= " << actBMax_m << endl);
//reset number of points inside domain
// clear previous coordinate maps
......@@ -138,14 +131,14 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
for(int x = 0; x < nr[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*A_m));
IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, -0.5*A_m));
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++) {
std::pair<int, int> pos(y, z);
double yt = getB(z*hr[2]);
double yb = -0.5*B_m;
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));
}
......
......@@ -36,34 +36,39 @@
/*
A_m and B_m are the half apperture of the box
A and B are the half apperture of the box
/ (A_m,B_m)
/ (A,B)
/
/
/
L1_m /
------------ --------------+ (-A_m,B_m)
| L2_m | |
C_m| | |
L1 /
------------ --------------+ (-A,B)
| L2 | |
C| | |
|------| | /
..... | /
(0,0)---.......-----------------+ /
..... | /
z | /
| | /
--------------------------------+/ (-A_m,-B_m)
--------------------------------+/ (-A,-B)
Length_m
Test in which of the 3 parts of the geometry we are in.
if((z < L1_m) || (z > (L1_m + L2_m)))
b = B_m;
if((z < L1) || (z > (L1 + L2)))
b = B;
else
b = B_m-C_m;
b = B-C;
A = max_m(0)
B = max_m(1)
L1 = min_m(2)
L2 = max_m(2) - min_m(2)
*/
class BoxCornerDomain : public IrregularDomain {
......@@ -72,8 +77,16 @@ public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
BoxCornerDomain(Vector_t nr, Vector_t hr);
BoxCornerDomain(double A, double B, double C, double Length, double L1, double L2, Vector_t nr, Vector_t hr,
/**
* \param A depth of the box
* \param B maximal height of the box
* \param C height of the corner
* \param length of the structure
* \param L1 length of the first part of the structure
* \param L2 length of the corner
*/
BoxCornerDomain(double A, double B, double C, double length,
double L1, double L2, Vector_t nr, Vector_t hr,
std::string interpl);
~BoxCornerDomain();
......@@ -86,10 +99,10 @@ public:
/// as a function of z, determine the hight (B) of the geometry
inline double getB(double z) {
if((z < L1_m) || (z > (L1_m + L2_m)))
return B_m;
if((z < min_m(2)) || (z > max_m(2)))
return max_m(1);
else
return B_m - C_m;
return max_m(1) - C_m;
}
/// queries if a given (x,y,z) coordinate lies inside the domain
......@@ -97,24 +110,11 @@ public:
const double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
const double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
const double b = getB(z * hr[2]);
return (xx < A_m && yy < b && z != 0 && z != nr[2] - 1);
return (xx < getXRangeMax() && yy < b && z != 0 && z != nr[2] - 1);
}
/// set semi-minor
//void setSemiMinor(double sm) {SemiMinor = sm;}
/// set semi-major
//void setSemiMajor(double sm) {SemiMajor = sm;}
void compute(Vector_t hr, NDIndex<3> localId);
double getXRangeMin() { return -A_m; }
double getXRangeMax() { return A_m; }
double getYRangeMin() { return -B_m;} // actBMin_m; }
double getYRangeMax() { return B_m; } // actBMax_m; }
double getZRangeMin() { return L1_m;}
double getZRangeMax() { return L1_m+L2_m; }
//TODO: ?
int getStartIdx() {return 0;}
......@@ -139,59 +139,31 @@ private:
/// mapping idx -> (x,y,z)
std::map<int, int> CoordMap;
/// depth of the box
double A_m;
/// the maximal hight of the box
double B_m;
/// because the geometry can change in the y direction
double actBMin_m;
double actBMax_m;
/// hight of the corner
double C_m;
/// lenght of the structure
double Length_m;
/// lenght of the first part of the structure
double L1_m;
/// lenght of the corner
double L2_m;
/// semi-major of the ellipse
// :FIXME: unused
//double SemiMajor;
/// semi-minor of the ellipse
// :FIXME: unused
//double SemiMinor;
/// length of the structure
double length_m;
/// interpolation type
int interpolationMethod;
/// for debug reasons
std::ofstream os_m;
std::ofstream os_m;
/// height of the corner
double C_m;
inline double getXIntersection(double cx, int /*z*/) {
if(cx < 0)
return -A_m;
else
return A_m;
return (cx < 0) ? min_m(0) : max_m(0);
}
inline double getYIntersection(double cy, int z) {
if(cy < 0)
return -B_m;
else
return getB(z * hr[2]);
return (cy < 0) ? min_m(1) : getB(z * hr[2]);
}
/// conversion from (x,y,z) to index in xyz plane
......
......@@ -35,16 +35,14 @@
//FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)
EllipticDomain::EllipticDomain(Vector_t nr, Vector_t hr) {
setNr(nr);
setHr(hr);
}
EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr,
Vector_t hr, std::string interpl)
: semiMajor_m(semimajor)
, semiMinor_m(semiminor)
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
std::string interpl)
{
Vector_t min(-bgeom->getA(), -bgeom->getB(), bgeom->getS());
Vector_t max( bgeom->getA(), bgeom->getB(), bgeom->getS() + bgeom->getLength());
setRangeMin(min);
setRangeMax(max);
setMinMaxZ(min[2], max[2]);
setNr(nr);
setHr(hr);
......@@ -56,21 +54,6 @@ EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr,
interpolationMethod_m = QUADRATIC;
}
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
std::string interpl)
: EllipticDomain(bgeom->getA(), bgeom->getB(), nr, hr, interpl)
{
setMinMaxZ(bgeom->getS(), bgeom->getLength());
}
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl)
: EllipticDomain(bgeom, nr,
Vector_t(2.0 * bgeom->getA() / nr[0],
2.0 * bgeom->getB() / nr[1],
(bgeom->getLength() - bgeom->getS()) / nr[2]),
interpl)
{ }
EllipticDomain::~EllipticDomain() {
//nothing so far
}
......@@ -125,16 +108,16 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
case LINEAR:
case QUADRATIC:
double smajsq = semiMajor_m * semiMajor_m;
double sminsq = semiMinor_m * semiMinor_m;
double smajsq = getXRangeMax() * getXRangeMax();
double sminsq = getYRangeMax() * getYRangeMax();
double yd = 0.0;
double xd = 0.0;
double pos = 0.0;
// calculate intersection with the ellipse
for (x = localId[0].first(); x <= localId[0].last(); x++) {
pos = - semiMajor_m + hr[0] * (x + 0.5);
if (std::abs(pos) >= semiMajor_m)
pos = getXRangeMin() + hr[0] * (x + 0.5);
if (std::abs(pos) >= getXRangeMax())
{
intersectYDir_m.insert(std::pair<int, double>(x, 0));
intersectYDir_m.insert(std::pair<int, double>(x, 0));
......@@ -147,8 +130,8 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
}
for (y = localId[0].first(); y < localId[1].last(); y++) {
pos = - semiMinor_m + hr[1] * (y + 0.5);
if (std::abs(pos) >= semiMinor_m)
pos = getYRangeMin() + hr[1] * (y + 0.5);
if (std::abs(pos) >= getYRangeMax())
{
intersectXDir_m.insert(std::pair<int, double>(y, 0));
intersectXDir_m.insert(std::pair<int, double>(y, 0));
......@@ -172,8 +155,8 @@ void EllipticDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t&
setMinMaxZ(rmin[2] - zsize * (1.0 + dh),
rmax[2] + zsize * (1.0 + dh));
origin = Vector_t(-semiMajor_m, -semiMinor_m, getMinZ());
mymax = Vector_t( semiMajor_m, semiMinor_m, getMaxZ());
origin = Vector_t(getXRangeMin(), getYRangeMin(), getMinZ());
mymax = Vector_t(getXRangeMax(), getYRangeMax(), getMaxZ());
for (int i = 0; i < 3; ++i)
hr[i] = (mymax[i] - origin[i]) / nr[i];
......@@ -242,8 +225,8 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, StencilValue_t& va
{
scaleFactor = 1.0;
double cx = - semiMajor_m + hr[0] * (x + 0.5);
double cy = - semiMinor_m + hr[1] * (y + 0.5);
double cx = getXRangeMin() + hr[0] * (x + 0.5);
double cy = getYRangeMin() + hr[1] * (y + 0.5);
double dx = 0.0;
std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
......
......@@ -42,16 +42,9 @@ public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
EllipticDomain(Vector_t nr, Vector_t hr);
EllipticDomain(double semimajor, double semiminor, Vector_t nr,
Vector_t hr, std::string interpl);
EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr,
Vector_t hr, std::string interpl);
EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
~EllipticDomain();
/// returns discretization at (x,y,z)
......@@ -60,11 +53,11 @@ 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);
double yy = - semiMinor_m + hr[1] * (y + 0.5);
double xx = getXRangeMin() + hr[0] * (x + 0.5);
double yy = getYRangeMin() + hr[1] * (y + 0.5);
bool isInsideEllipse = (xx * xx / (semiMajor_m * semiMajor_m) +
yy * yy / (semiMinor_m * semiMinor_m) < 1);
bool isInsideEllipse = (xx * xx / (getXRangeMax() * getXRangeMax()) +
yy * yy / (getYRangeMax() * getYRangeMax()) < 1);
return (isInsideEllipse && z > 0 && z < nr[2] - 1);
}
......@@ -75,13 +68,6 @@ public:
/// calculates intersection
void compute(Vector_t hr, NDIndex<3> localId);
double getXRangeMin() { return -semiMajor_m; }
double getXRangeMax() { return semiMajor_m; }
double getYRangeMin() { return -semiMinor_m; }
double getYRangeMax() { return semiMinor_m; }
double getZRangeMin() { return zMin_m; }
double getZRangeMax() { return zMax_m; }
void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
const Vector_t& rmax, double dh);
......@@ -103,12 +89,6 @@ private:
/// mapping idx -> (x,y,z)
std::map<int, int> coordMap_m;
/// semi-major of the ellipse
double semiMajor_m;
/// semi-minor of the ellipse
double semiMinor_m;
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
......
......@@ -69,4 +69,14 @@ void IrregularDomain::getBoundaryStencil(int id, StencilValue_t& value,
int idx = 0, idy = 0, idz = 0;
getCoord(id, idx, idy, idz);
getBoundaryStencil(idx, idy, idz, value, scaleFactor);
}
\ No newline at end of file
}
void IrregularDomain::resizeMesh(Vector_t& origin, Vector_t& hr,
const Vector_t& /*rmin*/, const Vector_t& /*rmax*/,
double /*dh*/)
{
origin = min_m;
for (int i = 0; i < 3; i++)
hr[i] = (max_m[i] - min_m[i]) / nr[i];
};
\ No newline at end of file
......@@ -121,12 +121,15 @@ public:
Quaternion_t getGlobalToLocalQuaternion() { return globalToLocalQuaternion_m;}
virtual double getXRangeMin() = 0;
virtual double getXRangeMax() = 0;
virtual double getYRangeMin() = 0;
virtual double getYRangeMax() = 0;
virtual double getZRangeMin() = 0;
virtual double getZRangeMax() = 0;
double getXRangeMin();
double getXRangeMax();
double getYRangeMin();
double getYRangeMax();
double getZRangeMin();
double getZRangeMax();
void setRangeMin(const Vector_t&);
void setRangeMax(const Vector_t&);
virtual int getIdx(int x, int y, int z) = 0;
......@@ -136,21 +139,7 @@ public:
virtual void resizeMesh(Vector_t& origin, Vector_t& hr,
const Vector_t& /*rmin*/, const Vector_t& /*rmax*/,
double /*dh*/)
{
double xmin = getXRangeMin();
double xmax = getXRangeMax();
double ymin = getYRangeMin();
double ymax = getYRangeMax();
double zmin = getZRangeMin();
double zmax = getZRangeMax();
origin = Vector_t(xmin, ymin , zmin);
Vector_t mymax = Vector_t(xmax, ymax , zmax);
for (int i = 0; i < 3; i++)
hr[i] = (mymax[i] - origin[i]) / nr[i];
};
double /*dh*/);
protected:
......@@ -164,6 +153,9 @@ protected:
double zMin_m;
double zMax_m;
Vector_t min_m;
Vector_t max_m;
/// mean position of bunch (m)
Vector_t rMean_m;
Quaternion_t globalToLocalQuaternion_m;
......@@ -178,6 +170,38 @@ inline bool IrregularDomain::hasGeometryChanged() {
return hasGeometryChanged_m;
}
inline double IrregularDomain::getXRangeMin() {
return min_m(0);
}
inline double IrregularDomain::getXRangeMax() {
return max_m(0);
}
inline double IrregularDomain::getYRangeMin() {
return min_m(1);
}
inline double IrregularDomain::getYRangeMax() {
return max_m(1);
}
inline double IrregularDomain::getZRangeMin() {
return min_m(2);
}
inline double IrregularDomain::getZRangeMax() {
return max_m(2);
}
inline void IrregularDomain::setRangeMin(const Vector_t& min) {
min_m = min;
}
inline void IrregularDomain::setRangeMax(const Vector_t& max) {
max_m = max;
}
#endif
......
......@@ -28,17 +28,9 @@
#include "Solvers/RectangularDomain.h"
#include "Utilities/OpalException.h"
RectangularDomain::RectangularDomain(Vector_t nr, Vector_t hr) {
setNr(nr);
setHr(hr);
a_m = 0.1;
b_m = 0.1;
nxy_m = nr[0] * nr[1];
}
RectangularDomain::RectangularDomain(double a, double b, Vector_t nr, Vector_t hr) {
a_m = a;
b_m = b;
setRangeMin(Vector_t(-a, -b, getMinZ()));
setRangeMax(Vector_t( a, b, getMax/()));
setNr(nr);
setHr(hr);
nxy_m = nr[0] * nr[1];
......
......@@ -36,9 +36,11 @@ public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
/// constructor
RectangularDomain(Vector_t nr, Vector_t hr);
/// constructor