Commit 97420544 authored by Tuelin Kaman's avatar Tuelin Kaman

BoundaryGeometry:ray-boundary intersection and quaternion rotations are in field solver

parent 3cff89ba
......@@ -1102,8 +1102,8 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
void PartBunch::computeSelfFields_cycl(double gamma, Vector_t const meanR, Vektor<double, 4> const quaternion) {
IpplTimings::startTimer(selfFieldTimer_m);
globalMeanR = meanR;
globalToLocalQuaternion = quaternion;
globalMeanR_m = meanR;
globalToLocalQuaternion_m = quaternion;
/// set initial charge density to zero.
rho_m = 0.0;
......@@ -1179,8 +1179,8 @@ void PartBunch::computeSelfFields_cycl(double gamma, Vector_t const meanR, Vekto
void PartBunch::computeSelfFields_cycl(int bin, Vector_t const meanR, Vektor<double, 4> const quaternion) {
IpplTimings::startTimer(selfFieldTimer_m);
globalMeanR = meanR;
globalToLocalQuaternion = quaternion;
globalMeanR_m = meanR;
globalToLocalQuaternion_m = quaternion;
/// set initial charge dentsity to zero.
rho_m = 0.0;
......
......@@ -402,8 +402,8 @@ public:
inline void setNumBunch(int n);
inline int getNumBunch() const;
inline Vector_t getGlobalMeanR() {return globalMeanR;}
inline Vektor<double, 4> getGlobalToLocalQuaternion() {return globalToLocalQuaternion;}
inline Vector_t getGlobalMeanR() {return globalMeanR_m;}
inline Vektor<double, 4> getGlobalToLocalQuaternion() {return globalToLocalQuaternion_m;}
void setSteptoLastInj(int n);
int getSteptoLastInj();
......@@ -555,8 +555,8 @@ private:
/// energy spread of the beam in keV
double dE_m;
Vektor<double, 4> globalToLocalQuaternion;
Vector_t globalMeanR;
Vector_t globalMeanR_m;
Vektor<double, 4> globalToLocalQuaternion_m;
/// maximal extend of particles
Vector_t rmax_m;
......
This diff is collapsed.
#ifndef ARBITRARY_DOMAIN_H
#define ARBITRARY_DOMAIN_H
#ifdef HAVE_SAAMG_SOLVER
#include <mpi.h>
......@@ -19,31 +20,32 @@ class ArbitraryDomain : public IrregularDomain {
public:
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion, std::string interpl);
~ArbitraryDomain();
/// returns discretization at (x,y,z)
void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void getBoundaryStencil(int idx, int idy, int idz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
/// returns discretization at 3D index
void getBoundaryStencil(int idxyz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
/// returns index of neighbours at (x,y,z)
void getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B);
void getNeighbours(int idx, int idy, int idz, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns index of neighbours at 3D index
void getNeighbours(int idxyz, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns type of boundary condition
std::string getType() {return "Geometric";}
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z);
inline bool isInside(int idx, int idy, int idz);
/// returns number of nodes in xy plane
int getNumXY(int z);
int getNumXY(int idz);
/// calculates intersection
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
// calculates intersection with rotated and shifted geometry
void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
int getStartId() {return startId;}
/// conversion from (x,y,z) to index on the 3D grid
// int getIdx(int x, int y, int z);
double getXRangeMin() { return Geo_mincoords_m(0); }
double getXRangeMax() { return Geo_maxcoords_m(0); }
......@@ -68,24 +70,29 @@ private:
/// all intersection points with gridlines in Z direction
PointList IntersectHiZ, IntersectLoZ;
// meanR to shift from global to local frame
Vector_t globalMeanR_m;
Vektor<double, 4> globalToLocalQuaternion_m;
Vektor<double, 4> localToGlobalQuaternion_m;
int startId;
/// here we store the number of nodes in a xy layer for a given z coordinate
// Here we store the number of nodes in a xy layer for a given z coordinate
std::map<int, int> numXY;
std::map<int, int> numYZ;
std::map<int, int> numXZ;
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
// Number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
/// mapping (x,y,z) -> idxyz
// mapping (x,y,z) -> idxyz
std::map<int, int> IdxMap;
/// mapping idxyz -> (x,y,z)
// Mapping idxyz -> (x,y,z)
std::map<int, int> CoordMap;
/// interpolation type
// Interpolation type
int interpolationMethod;
/// flag indicating if geometry has changed for the current time-step
// Flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
Vector_t Geo_nr_m;
......@@ -93,20 +100,27 @@ private:
Vector_t Geo_mincoords_m;
Vector_t Geo_maxcoords_m;
/// conversion from (x,y,z) to index in xyz plane
inline int toCoordIdx(int x, int y, int z);
/// conversion from (x,y,z) to index on the 3D grid
int getIdx(int x, int y, int z);
/// conversion from a 3D index to (x,y,z)
// Conversion from (x,y,z) to index in xyz plane
inline int toCoordIdx(int idx, int idy, int idz);
// Conversion from (x,y,z) to index on the 3D grid
int getIdx(int idx, int idy, int idz);
// Conversion from a 3D index to (x,y,z)
inline void getCoord(int idxyz, int &x, int &y, int &z);
inline void crossProduct(double A[], double B[], double C[]);
inline double dotProduct(double v1[], double v2[]) { return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]); }
/// different interpolation methods for boundary points
void ConstantInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void LinearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void QuadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
// Different interpolation methods for boundary points
void ConstantInterpolation(int idx, int idy, int idz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void LinearInterpolation(int idx, int idy, int idz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
void QuadraticInterpolation(int idx, int idy, int idz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
// Rotate positive axes with quaternion -DW
inline void rotateWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateXAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateYAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateZAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
};
#endif //#ifdef HAVE_SAAMG_SOLVER
......
......@@ -55,7 +55,7 @@ BoxCornerDomain::~BoxCornerDomain() {
// for the moment we center the box corner geometry around the center of the grid
// hr holds the grid-spacings (boundary ellipse embedded in hr-grid)
void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId){
void BoxCornerDomain::Compute(Vector_t hr){
//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]) {
......@@ -63,8 +63,6 @@ void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId){
// return;
// }
setHr(hr);
hasGeometryChanged_m = true;
......@@ -132,6 +130,13 @@ void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId){
*/
}
void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion){
}
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) {
// determine which interpolation method we use for points near the boundary
......
......@@ -99,7 +99,10 @@ public:
/// set semi-major
//void setSemiMajor(double sm) {SemiMajor = sm;}
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
double getXRangeMin() { return -A_m; }
double getXRangeMax() { return A_m; }
......
......@@ -36,8 +36,7 @@ EllipticDomain::~EllipticDomain() {
// for this geometry we only have to calculate the intersection on
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
// hr holds the grid-spacings (boundary ellipse embedded in hr-grid)
void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId){
void EllipticDomain::Compute(Vector_t hr){
//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;
......@@ -116,6 +115,11 @@ void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
}
void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion){
}
void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
scaleFactor = 1.0;
......
......@@ -39,7 +39,10 @@ public:
/// set semi-major
void setSemiMajor(double sm) {SemiMajor = sm;}
/// calculates intersection
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
double getXRangeMin() { return -SemiMajor; }
double getXRangeMax() { return SemiMajor; }
......
......@@ -22,8 +22,9 @@ public:
/** method to compute the intersection points with the boundary geometry (stored in some appropriate data structure)
* \param hr updated mesh spacings
*/
virtual void Compute(Vector_t hr) = 0;
virtual void Compute(Vector_t hr, NDIndex<3> localId) = 0;
virtual void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion) = 0;
/** method to get the number of gridpoints in a given z plane
* \param z coordinate of the z plane
* \return int number of grid nodes in the given z plane
......
......@@ -60,7 +60,7 @@ MGPoissonSolver::MGPoissonSolver(PartBunch &beam,Mesh_t *mesh, FieldLayout_t *fl
}
else if(currentGeometry->getTopology() == "BOXCORNER") {
bp = new BoxCornerDomain(currentGeometry->getA(), currentGeometry->getB(), currentGeometry->getC(), currentGeometry->getLength(),currentGeometry->getL1(), currentGeometry->getL2(), orig_nr_m, hr_m, interpl);
// bp->Compute(itsBunch_m->get_hr(), localId);
bp->Compute(itsBunch_m->get_hr());
} else {
ERRORMSG("Geometry not known" << endl);
exit(1);
......@@ -182,11 +182,11 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
bp->setMinMaxZ(itsBunch_m->get_origin()[2], itsBunch_m->get_maxExtend()[2]);
bp->setNr(nr_m);
bp->setHr(hr);
localId = layout_m->getLocalNDIndex();
IpplTimings::startTimer(FunctionTimer1_m);
bp->Compute(hr, localId); // Build the index and coord map
bp->Compute(hr, localId, itsBunch_m->getGlobalMeanR(), itsBunch_m->getGlobalToLocalQuaternion());
IpplTimings::stopTimer(FunctionTimer1_m);
IpplTimings::startTimer(FunctionTimer2_m);
......
......@@ -18,7 +18,7 @@ RectangularDomain::RectangularDomain(double a, double b, Vector_t nr, Vector_t h
nxy_m = nr[0] * nr[1];
}
void RectangularDomain::Compute(Vector_t hr, NDIndex<3> localId){
void RectangularDomain::Compute(Vector_t hr){
setHr(hr);
nxy_m = nr[0] * nr[1];
}
......
......@@ -17,7 +17,7 @@ public:
RectangularDomain(double a, double b, Vector_t nr, Vector_t hr);
/// calculates intersection with the beam pipe
void Compute(Vector_t hr, NDIndex<3> localId);
void Compute(Vector_t hr);
/// returns number of nodes in xy plane (here independent of z coordinate)
int getNumXY(int z);
......
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