Skip to content

GitLab

Commit e0cf1208 authored by Andreas Adelmann's avatar Andreas Adelmann
Browse files

Merge branch 'master' of gitorious.psi.ch:opal/src

parents 73416bbc 051beb73
Hide whitespace changes
Inline Side-by-side
Showing with 185 additions and 255 deletions
classic/5.0/src/Algorithms/PartBunch.cpp
View file @ e0cf1208
......@@ -690,11 +690,6 @@ double PartBunch::getMaxdEBins() {
void PartBunch::computeSelfFields(int binNumber) {
IpplTimings::startTimer(selfFieldTimer_m);
// Need to set mean R to shift geometry from file to account for beam movement
// No rotations necessary, i.e. quaternion stays default unit-quatertion
// TODO: Include transversal offsets here!
globalMeanR_m = Vector_t(0.0, 0.0, get_sPos());
/// Set initial charge density to zero. Create image charge
/// potential field.
rho_m = 0.0;
......@@ -947,11 +942,6 @@ void PartBunch::computeSelfFields() {
rho_m = 0.0;
eg_m = Vector_t(0.0);
// Need to set mean R to shift geometry from file to account for beam movement
// No rotations necessary, i.e. quaternion stays default unit-quatertion
// TODO: Include transversal offsets here!
globalMeanR_m = Vector_t(0.0, 0.0, get_sPos());
if(fs_m->hasValidSolver()) {
//use the mesh that is already set
//if (fs_m->getFieldSolverType() == "SAAMG")
......@@ -1216,15 +1206,10 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
* -) Fixed an error where gamma was not taken into account correctly in direction of movement (y in cyclotron)
* -) Comment: There is no account for image charges in the cyclotron tracker (yet?)!
*/
void PartBunch::computeSelfFields_cycl(double gamma,
Vector_t const meanR,
Quaternion_t const quaternion) {
void PartBunch::computeSelfFields_cycl(double gamma) {
IpplTimings::startTimer(selfFieldTimer_m);
globalMeanR_m = meanR;
globalToLocalQuaternion_m = quaternion;
/// set initial charge density to zero.
rho_m = 0.0;
......@@ -1411,12 +1396,9 @@ void PartBunch::computeSelfFields_cycl(double gamma,
* -) Fixed an error where gamma was not taken into account correctly in direction of movement (y in cyclotron)
* -) Comment: There is no account for image charges in the cyclotron tracker (yet?)!
*/
void PartBunch::computeSelfFields_cycl(int bin, Vector_t const meanR, Quaternion_t const quaternion) {
void PartBunch::computeSelfFields_cycl(int bin) {
IpplTimings::startTimer(selfFieldTimer_m);
globalMeanR_m = meanR;
globalToLocalQuaternion_m = quaternion;
/// set initial charge dentsity to zero.
rho_m = 0.0;
......
classic/5.0/src/Algorithms/PartBunch.h
View file @ e0cf1208
......@@ -277,17 +277,8 @@ public:
/** /brief used for self fields with binned distribution */
void computeSelfFields(int b);
//void computeSelfFields_cycl(double gamma);
//void computeSelfFields_cycl(int b);
// Replaced computeSelfFields_cycl() with versions that have meanR and the quaternion of the
// rotation of the particle bunch in order to take into account the rotation
// when finding the boundary conditions for the fieldsolver. -DW
void computeSelfFields_cycl(double gamma, Vector_t const meanR=Vector_t(0.0, 0.0, 0.0),
Quaternion_t const quaternion=Quaternion_t(1.0, 0.0, 0.0, 0.0));
void computeSelfFields_cycl(int b, Vector_t const meanR=Vector_t(0.0, 0.0, 0.0),
Quaternion_t const quaternion=Quaternion_t(1.0, 0.0, 0.0, 0.0));
void computeSelfFields_cycl(double gamma);
void computeSelfFields_cycl(int b);
void resetInterpolationCache(bool clearCache = false);
......@@ -410,8 +401,11 @@ public:
inline void setNumBunch(int n);
inline int getNumBunch() const;
inline void setGlobalMeanR(Vector_t globalMeanR) {globalMeanR_m = globalMeanR;}
inline Vector_t getGlobalMeanR() {return globalMeanR_m;}
inline Vektor<double, 4> getGlobalToLocalQuaternion() {return globalToLocalQuaternion_m;}
inline void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion) {
globalToLocalQuaternion_m = globalToLocalQuaternion;}
inline Quaternion_t getGlobalToLocalQuaternion() {return globalToLocalQuaternion_m;}
void setSteptoLastInj(int n);
int getSteptoLastInj();
......
src/Algorithms/ParallelCyclotronTracker.cpp
View file @ e0cf1208
......@@ -2187,7 +2187,9 @@ void ParallelCyclotronTracker::Tracker_RK4() {
repartition();
itsBunch->computeSelfFields_cycl(temp_meangamma, 0.001 * meanR, quaternionToYAxis);
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(temp_meangamma);
IpplTimings::startTimer(TransformTimer_m);
......@@ -3070,7 +3072,9 @@ void ParallelCyclotronTracker::Tracker_Generic() {
for(int b = 0; b < itsBunch->getLastemittedBin(); b++) {
itsBunch->setBinCharge(b, itsBunch->getChargePerParticle());
itsBunch->computeSelfFields_cycl(b, 0.001 * meanR, quaternionToYAxis);
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(b);
}
itsBunch->Q = itsBunch->getChargePerParticle();
......@@ -3108,7 +3112,9 @@ void ParallelCyclotronTracker::Tracker_Generic() {
repartition();
itsBunch->computeSelfFields_cycl(temp_meangamma, 0.001 * meanR, quaternionToYAxis);
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(temp_meangamma);
IpplTimings::startTimer(TransformTimer_m);
......
src/Algorithms/ParallelTTracker.cpp
View file @ e0cf1208
......@@ -2246,12 +2246,14 @@ void ParallelTTracker::computeSpaceChargeFields() {
binNumber < itsBunch->GetNumberOfEnergyBins(); ++binNumber) {
itsBunch->setBinCharge(binNumber);
itsBunch->setGlobalMeanR(itsBunch->get_centroid());
itsBunch->computeSelfFields(binNumber);
itsBunch->Q = Q_back;
}
} else {
itsBunch->setGlobalMeanR(itsBunch->get_centroid());
itsBunch->computeSelfFields();
}
......@@ -3217,4 +3219,4 @@ Vector_t ParallelTTracker::calcMeanP() const {
}
reduce(meanP, meanP, OpAddAssign());
return meanP / Vector_t(itsBunch->getTotalNum());
}
\ No newline at end of file
}
src/CMakeLists.txt
View file @ e0cf1208
......@@ -107,7 +107,7 @@ add_library( OPALib ${OPAL_SRCS} opal.cpp )
set_target_properties( OPALib PROPERTIES OUTPUT_NAME OPAL )
add_executable( opal Main.cpp )
target_link_libraries( opal OPALib CLASSIC ${OPAL_LIBS} ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES} ${CCSE_LIBRARIES} ${OTHER_CMAKE_EXE_LINKER_FLAGS} -lgfortran)
target_link_libraries( opal OPALib CLASSIC ${OPAL_LIBS} ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES} ${CCSE_LIBRARIES} ${OTHER_CMAKE_EXE_LINKER_FLAGS} -lgfortran ${CMAKE_DL_LIBS})
# build unit tests; unit tests are defined in unit_tests directory, we build the
# executable here so that we can use the OPAL src and lib definitions
......
src/Solvers/ArbitraryDomain.cpp
View file @ e0cf1208
......@@ -25,27 +25,25 @@
#define MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define MAX2(a,b) (((a) > (b)) ? (a) : (b))
ArbitraryDomain::ArbitraryDomain(
BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t hr,
std::string interpl) {
bgeom_m = bgeom;
intersectMinCoords_m = bgeom->getmincoords();
intersectMaxCoords_m = bgeom->getmaxcoords();
setNr(nr);
setHr(hr);
startId = 0;
if(interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "LINEAR")
interpolationMethod = LINEAR;
else if(interpl == "QUADRATIC")
interpolationMethod = QUADRATIC;
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t hr,
std::string interpl){
bgeom_m = bgeom;
minCoords_m = bgeom->getmincoords();
maxCoords_m = bgeom->getmaxcoords();
setNr(nr);
setHr(hr);
startId = 0;
if (interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "LINEAR")
interpolationMethod = LINEAR;
else if(interpl == "QUADRATIC")
interpolationMethod = QUADRATIC;
}
ArbitraryDomain::~ArbitraryDomain() {
......@@ -57,7 +55,7 @@ void ArbitraryDomain::Compute(Vector_t hr) {
setHr(hr);
}
/*
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
setHr(hr);
......@@ -224,15 +222,18 @@ 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){
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
setHr(hr);
globalMeanR_m = globalMeanR;
globalToLocalQuaternion_m = globalToLocalQuaternion;
localToGlobalQuaternion_m[0] = globalToLocalQuaternion[0];
globalMeanR_m = getGlobalMeanR();
globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
for (int i=1; i<4; i++)
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion[i];
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr[2] - 1) ? 0 : 1;
......@@ -252,11 +253,8 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
//calculate intersection
Vector_t P, saveP, dir, I;
// TODO: Find and set the reference point for any time of geometry
//Reference Point inside the boundary for IsoDar Geo
Vector_t P0 = Vector_t(0,0,bgeom_m->getmincoords()[2]+hr[2]);
//Reference Point where the boundary geometry is cylinder
//P0 = Vector_t(0,0,0); // Uncomment for cylinder Benchmarking -DW
//Reference Point
Vector_t P0 = globalMeanR_m;
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
......@@ -274,8 +272,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;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
} else {
......@@ -285,8 +282,7 @@ 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;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
} else {
......@@ -297,8 +293,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;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
} else {
......@@ -308,8 +303,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;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
} else {
......@@ -320,8 +314,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;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
} else {
......@@ -331,8 +324,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;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
} else {
......@@ -387,9 +379,9 @@ inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) {
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
Vector_t P;
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];
bool ret = false;
int countH, countL;
......@@ -448,7 +440,7 @@ 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);
......@@ -485,141 +477,96 @@ void ArbitraryDomain::ConstantInterpolation(int idx, int idy, int idz, double& W
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)
{
scaleFactor = 1.0;
double dx=-1, dy=-1, dz=-1;
scaleFactor = 1;
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;
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
itrH = IntersectHiZ.find(coordxyz);
itrL = IntersectLoZ.find(coordxyz);
countH = IntersectHiZ.count(coordxyz);
countL = IntersectLoZ.count(coordxyz);
if(countH == 1 && countL == 1)
ret = (cz <= itrH->second) && (cz >= itrL->second);
std::multimap< std::pair<int, int>, double >::iterator it;
std::pair< std::multimap< std::pair<int, int>, double>::iterator, std::multimap< std::pair<int, int>, double>::iterator > ret;
std::pair<int, int> coordyz(idy, idz);
ret = IntersectXDir.equal_range(coordyz);
for(it = ret.first; it != ret.second; ++it) {
if(fabs(it->second - cx) < hr[0]) {
dx = it->second;
break;
}
}
std::pair<int, int> coordxz(idx, idz);
ret = IntersectYDir.equal_range(coordxz);
for(it = ret.first; it != ret.second; ++it) {
if(fabs(it->second - cy) < hr[1]) {
dy = it->second;
break;
}
}
double dw=hr[0];
double de=hr[0];
double dn=hr[1];
double ds=hr[1];
double dx_w=hr[0];
double dx_e=hr[0];
double dy_n=hr[1];
double dy_s=hr[1];
double dz_f=hr[2];
double dz_b=hr[2];
C = 0.0;
// we are a right boundary point
if(dx >= 0 && dx > cx)
{
C += 1/((dx-cx)*de);
E = 0.0;
} else {
C += 1/(de*de);
E = -1/(de*de);
}
// we are a left boundary point
if(dx <= 0 && dx < cx)
{
C += 1/((cx-dx)*dw);
W = 0.0;
} else {
C += 1/(dw*dw);
W = -1/(dw*dw);
}
// we are a upper boundary point
if(dy >= 0 && dy > cy)
{
C += 1/((dy-cy)*dn);
N = 0.0;
} else {
C += 1/(dn*dn);
N = -1/(dn*dn);
}
// we are a lower boundary point
if(dy <= 0 && dy < cy)
{
C += 1/((cy-dy)*ds);
S = 0.0;
} else {
C += 1/(ds*ds);
S = -1/(ds*ds);
}
F = -1/(hr[2]*hr[2]);
B = -1/(hr[2]*hr[2]);
//XXX: In stand-alone only Dirichlet for validation purposes
if(z == 0 || z == nr[2]-1) {
// Dirichlet
C += 2/hr[2]*1/hr[2];
//C += 1/hr[2]*1/hr[2];
// case where we are on the Neumann BC in Z-direction
// where we distinguish two cases
if(z == 0)
F = 0.0;
else
B = 0.0;
//for test no neumann
//C += 2/((hr[2]*nr[2]/2.0) * hr[2]);
//
// double d = hr[2]*(nr[2])/2;
// C += 2/(d * hr[2]);
////neumann stuff
//W /= 2.0;
//E /= 2.0;
//N /= 2.0;
//S /= 2.0;
//C /= 2.0;
scaleFactor *= 0.5;
std::tuple<int, int, int> coordxyz(idx, idy, idz);
std::multimap < std::tuple<int, int, int>, double >::iterator itrH, itrL;
} else
C += 2*1/hr[2]*1/hr[2];
if (idx == nr[0]-1)
dx_e = fabs(IntersectHiX.find(coordxyz)->second - cx);
if (idx == 0)
dx_w = fabs(IntersectLoX.find(coordxyz)->second - cx);
if (idy == nr[1]-1)
dy_n = fabs(IntersectHiY.find(coordxyz)->second - cy);
if (idy == 0)
dy_s = fabs(IntersectLoY.find(coordxyz)->second - cy);
if (idz == nr[2]-1)
dz_b = fabs(IntersectHiZ.find(coordxyz)->second - cz);
if (idz == 0)
dz_f = fabs(IntersectLoZ.find(coordxyz)->second - cz);
if(dx_w != 0)
W = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
else
W = 0;
if(dx_e != 0)
E = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
else
E = 0;
if(dy_n != 0)
N = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
else
N = 0;
if(dy_s != 0)
S = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
else
S = 0;
if(dz_f != 0)
F = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
else
F = 0;
if(dz_b != 0)
B = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
else
B = 0;
//RHS scaleFactor for current 3D index
//0.5* comes from discretiztaion
//scaleFactor = 0.5*(dw+de)*(dn+ds)*(df+db);
scaleFactor = 0.5;
if(dx_w + dx_e != 0)
scaleFactor *= (dx_w + dx_e);
if(dy_n + dy_s != 0)
scaleFactor *= (dy_n + dy_s);
if(dz_f + dz_b != 0)
scaleFactor *= (dz_f + dz_b);
//catch the case where a point lies on the boundary
double m1 = dx_w * dx_e;
double m2 = dy_n * dy_s;
if(dx_e == 0)
m1 = dx_w;
if(dx_w == 0)
m1 = dx_e;
if(dy_n == 0)
m2 = dy_s;
if(dy_s == 0)
m2 = dy_n;
C = 2 / hr[2];
if(dx_w != 0 || dx_e != 0)
C *= (dx_w + dx_e);
if(dy_n != 0 || dy_s != 0)
C *= (dy_n + dy_s);
if(dx_w != 0 || dx_e != 0)
C += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
if(dy_n != 0 || dy_s != 0)
C += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
}
*/
void ArbitraryDomain::getNeighbours(int id, int &W, int &E, int &S, int &N, int &F, int &B) {
......@@ -665,7 +612,7 @@ inline void ArbitraryDomain::crossProduct(double A[], double B[], double C[]) {
C[2] = A[0] * B[1] - A[1] * B[0];
}
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates a Vector_t (3 elements) using a quaternion.
// Flip direction of rotation by quaternionVectorcomponent *= -1
......@@ -678,7 +625,7 @@ inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Vektor<double, 4
+ 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}
inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates the positive xaxis using a quaternion.
v(0) = quaternion(0) * quaternion(0)
......@@ -690,7 +637,7 @@ inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Vektor<doub
v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}
inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates the positive yaxis using a quaternion.
v(0) = 2.0 * (quaternion(1) * quaternion(2) - quaternion(0) * quaternion(3));
......@@ -703,7 +650,7 @@ inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Vektor<doub
v(2) = 2.0 * (quaternion(2) * quaternion(3) + quaternion(0) * quaternion(1));
}
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates the positive zaxis using a quaternion.
v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
......
src/Solvers/ArbitraryDomain.h
View file @ e0cf1208
......@@ -20,7 +20,7 @@ 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(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, Vector_t globalMeanR, Quaternion_t globalToLocalQuaternion, std::string interpl);
~ArbitraryDomain();
......@@ -41,24 +41,26 @@ public:
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);
void Compute(Vector_t hr, NDIndex<3> localId);
int getStartId() {return startId;}
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); }