Commit 9dbd7c1c authored by kraus's avatar kraus
Browse files

- remove unused files

- putting own header file first
  o for reason see http://www.cplusplus.com/forum/articles/10627/#msg49679 and
                   http://www.cplusplus.com/forum/articles/10627/#msg49680
parent 814c8c73
......@@ -25,12 +25,13 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include "Elements/OpalOffset/OpalLocalCartesianOffset.h"
#include <string>
#include "AbsBeamline/Offset.h"
#include "Utilities/OpalException.h"
#include "Elements/OpalElement.h"
#include "Elements/OpalOffset/OpalLocalCartesianOffset.h"
#include "Attributes/Attributes.h"
......
......@@ -25,11 +25,12 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include "Elements/OpalOffset/OpalLocalCylindricalOffset.h"
#include <string>
#include "AbsBeamline/Offset.h"
#include "Elements/OpalElement.h"
#include "Elements/OpalOffset/OpalLocalCylindricalOffset.h"
#include "Attributes/Attributes.h"
......
......@@ -25,9 +25,10 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include "Elements/OpalPolynomialTimeDependence.h"
#include <string>
#include "Elements/OpalPolynomialTimeDependence.h"
#include "Algorithms/PolynomialTimeDependence.h"
#include "Attributes/Attributes.h"
......
/*
/*
* Copyright (c) 2012-2014, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#include "Elements/OpalRing.h"
#include <sstream>
#include <limits>
......@@ -37,8 +39,6 @@
#include "Utilities/OpalException.h"
#include "Utilities/OpalRingSection.h"
#include "Elements/OpalRing.h"
// fairly generous tolerance here... maybe too generous? Maybe should be
// user parameter?
const double OpalRing::lengthTolerance_m = 1e-2;
......@@ -136,7 +136,7 @@ bool OpalRing::apply(const Vector_t &R, const Vector_t &centroid,
std::vector<OpalRingSection*> sections = getSectionsAt(R);
bool outOfBounds = true;
// assume field maps don't set B, E to 0...
// assume field maps don't set B, E to 0...
for (size_t i = 0; i < sections.size(); ++i) {
Vector_t B_temp(0.0, 0.0, 0.0);
Vector_t E_temp(0.0, 0.0, 0.0);
......@@ -220,7 +220,7 @@ void OpalRing::rotateToCyclCoordinates(Euclid3D& delta) const {
delta = euclid;
}
Vector_t OpalRing::getNextPosition() const {
if (section_list_m.size() > 0) {
return section_list_m.back()->getEndPosition();
......@@ -314,7 +314,7 @@ void OpalRing::lockRing() {
// Apply symmetry properties; I think it is fastest to just duplicate
// elements rather than try to do some angle manipulations when we do field
// lookups because we do field lookups in Cartesian coordinates in general.
msg << "Applying symmetry to OpalRing" << endl;
msg << "Applying symmetry to OpalRing" << endl;
for (int i = 1; i < symmetry_m; ++i) {
for (size_t j = 0; j < sectionListSize; ++j) {
appendElement(*section_list_m[j]->getComponent());
......@@ -351,7 +351,7 @@ void OpalRing::checkAndClose() {
Vector_t last_pos = section_list_m.back()->getEndPosition();
Vector_t last_norm = section_list_m.back()->getEndNormal();
for (int i = 0; i < 3; ++i) {
if (fabs(first_pos(i) - last_pos(i)) > lengthTolerance_m ||
if (fabs(first_pos(i) - last_pos(i)) > lengthTolerance_m ||
fabs(first_norm(i) - last_norm(i)) > angleTolerance_m)
throw OpalException("OpalRing::lockRing",
"Ring is not closed");
......@@ -385,4 +385,3 @@ bool OpalRing::sectionCompare(OpalRingSection const* const sec1,
OpalRingSection const* const sec2) {
return sec1 > sec2;
}
......@@ -25,9 +25,10 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include "Elements/OpalRingDefinition.h"
#include <limits>
#include "Elements/OpalRingDefinition.h"
#include "Elements/OpalRing.h"
#include "Attributes/Attributes.h"
......
/*
/*
* Copyright (c) 2012, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#include "Elements/OpalVariableRFCavity.h"
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
#include "Attributes/Attributes.h"
#include "Algorithms/AbstractTimeDependence.h"
#include "AbsBeamline/VariableRFCavity.h"
#include "Elements/OpalVariableRFCavity.h"
extern Inform *gmsg;
const std::string OpalVariableRFCavity::doc_string =
const std::string OpalVariableRFCavity::doc_string =
std::string("The \"VARIABLE_RF_CAVITY\" element defines an RF cavity ")+
std::string("with time dependent frequency, phase and amplitude.");
......@@ -117,4 +118,3 @@ void OpalVariableRFCavity::update() {
cavity->setHeight(height);
setElement(cavity->makeAlignWrapper());
}
......@@ -17,8 +17,9 @@
//
// ------------------------------------------------------------------------
#include "AbstractObjects/Attribute.h"
#include "AbstractObjects/Expressions.h"
#include "AbstractObjects/Attribute.h"
#include "AbstractObjects/PlaceRep.h"
#include "AbstractObjects/RangeRep.h"
#include "AbstractObjects/Table.h"
......
......@@ -16,9 +16,10 @@
//
// ------------------------------------------------------------------------
#include "OpalParser/IfStatement.h"
#include "OpalParser/CompoundStatement.h"
#include "AbstractObjects/OpalData.h"
#include "OpalParser/IfStatement.h"
#include "Attributes/Attributes.h"
#include "Parser/Parser.h"
#include "Parser/Token.h"
......
......@@ -4,7 +4,7 @@
// Copyright & License: See Copyright.readme in src directory
// ------------------------------------------------------------------------
// Class ArbitraryDomain
// Interface to iterative solver and boundary geometry
// Interface to iterative solver and boundary geometry
// for space charge calculation
//
// ------------------------------------------------------------------------
......@@ -14,6 +14,8 @@
//#define DEBUG_INTERSECT_RAY_BOUNDARY
#ifdef HAVE_SAAMG_SOLVER
#include "Solvers/ArbitraryDomain.h"
#include <map>
#include <cmath>
#include <iostream>
......@@ -23,12 +25,10 @@
#define MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define MAX2(a,b) (((a) > (b)) ? (a) : (b))
#include "ArbitraryDomain.h"
ArbitraryDomain::ArbitraryDomain(
BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t hr,
BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t hr,
std::string interpl) {
bgeom_m = bgeom;
......@@ -68,7 +68,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
int yGhostOffsetRight = (localId[1].last() == nr[1] - 1) ? 0 : 1;
int xGhostOffsetLeft = (localId[0].first()== 0) ? 0 : 1;
int xGhostOffsetRight = (localId[0].last() == nr[0] - 1) ? 0 : 1;
hasGeometryChanged_m = true;
IntersectLoX.clear();
......@@ -79,7 +79,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
IntersectHiZ.clear();
//calculate intersection
//calculate intersection
Vector_t P, dir, I;
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
P[2] = (idz - (nr[2]-1)/2.0)*hr[2];
......@@ -119,7 +119,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
//number of ghost nodes to the right
//number of ghost nodes to the right
int znumGhostNodesRight = 0;
if(zGhostOffsetRight != 0) {
for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
......@@ -130,7 +130,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
//number of ghost nodes to the left
//number of ghost nodes to the left
int znumGhostNodesLeft = 0;
if(zGhostOffsetLeft != 0) {
for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
......@@ -141,7 +141,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
//number of ghost nodes to the right
//number of ghost nodes to the right
int ynumGhostNodesRight = 0;
if(yGhostOffsetRight != 0) {
for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
......@@ -152,7 +152,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
//number of ghost nodes to the left
//number of ghost nodes to the left
int ynumGhostNodesLeft = 0;
if(yGhostOffsetLeft != 0) {
for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
......@@ -164,7 +164,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
//number of ghost nodes to the right
//number of ghost nodes to the right
int xnumGhostNodesRight = 0;
if(xGhostOffsetRight != 0) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
......@@ -175,7 +175,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
//number of ghost nodes to the left
//number of ghost nodes to the left
int xnumGhostNodesLeft = 0;
if(xGhostOffsetLeft != 0) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
......@@ -186,7 +186,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
//xy points in z plane
int numxy;
int numxy;
int numtotalxy = 0;
numXY.clear();
......@@ -210,7 +210,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
//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++) {
......@@ -229,10 +229,10 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
setHr(hr);
globalMeanR_m = globalMeanR;
globalToLocalQuaternion_m = globalToLocalQuaternion;
globalToLocalQuaternion_m = globalToLocalQuaternion;
localToGlobalQuaternion_m[0] = globalToLocalQuaternion[0];
for (int i=1; i<4; i++)
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion[i];
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion[i];
int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr[2] - 1) ? 0 : 1;
......@@ -240,7 +240,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
int yGhostOffsetRight = (localId[1].last() == nr[1] - 1) ? 0 : 1;
int xGhostOffsetLeft = (localId[0].first()== 0) ? 0 : 1;
int xGhostOffsetRight = (localId[0].last() == nr[0] - 1) ? 0 : 1;
hasGeometryChanged_m = true;
IntersectLoX.clear();
......@@ -250,12 +250,12 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
IntersectLoZ.clear();
IntersectHiZ.clear();
//calculate intersection
//calculate intersection
Vector_t P, saveP, dir, I;
// TODO: Find and set the reference point for any time of geometry
// 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
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);
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
......@@ -264,19 +264,19 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
saveP[0] = (idx - (nr[0]-1)/2.0)*hr[0];
P = saveP;
rotateWithQuaternion(P, localToGlobalQuaternion_m);
rotateWithQuaternion(P, localToGlobalQuaternion_m);
P += globalMeanR_m;
if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
P0 = P;
std::tuple<int, int, int> pos(idx, idy, idz);
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);
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
......@@ -287,19 +287,19 @@ 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);
rotateWithQuaternion(I, globalToLocalQuaternion_m);
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;
#endif
}
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);
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
......@@ -310,7 +310,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);
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
......@@ -322,28 +322,28 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
// setXRangeMax(MIN2(intersectMaxCoords_m(0),I[0]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
#endif
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)){
// setXRangeMin(MAX2(intersectMinCoords_m(0),I[0]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
#endif
}
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
#endif
}
}
}
......@@ -371,7 +371,7 @@ inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
// Conversion from (x,y,z) to index on the 3D grid
int ArbitraryDomain::getIdx(int idx, int idy, int idz) {
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) {
......@@ -395,7 +395,7 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
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);
......@@ -422,12 +422,12 @@ inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
countL = IntersectLoX.count(coordxyz);
if(countH == 1 && countL == 1)
ret = ret && (P[0] <= itrH->second) && (P[0] >= itrL->second);
return ret;
return ret;
}
int ArbitraryDomain::getNumXY(int z) {
return numXY[z];
}
......@@ -471,22 +471,22 @@ void ArbitraryDomain::ConstantInterpolation(int idx, int idy, int idz, double& W
if(!isInside(idx-1,idy,idz))
W = 0.0;
if(!isInside(idx+1,idy,idz))
if(!isInside(idx+1,idy,idz))
E = 0.0;
if(!isInside(idx,idy+1,idz))
N = 0.0;
if(!isInside(idx,idy-1,idz))
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;
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)
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;
......@@ -501,7 +501,7 @@ void ArbitraryDomain::LinearInterpolation(int idx, int idy, int idz, double& W,
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);
......@@ -544,7 +544,7 @@ void ArbitraryDomain::LinearInterpolation(int idx, int idy, int idz, double& W,
// we are a right boundary point
if(dx >= 0 && dx > cx)
{
{
C += 1/((dx-cx)*de);
E = 0.0;
} else {
......@@ -582,7 +582,7 @@ void ArbitraryDomain::LinearInterpolation(int idx, int idy, int idz, double& W,
S = -1/(ds*ds);
}
F = -1/(hr[2]*hr[2]);
F = -1/(hr[2]*hr[2]);
B = -1/(hr[2]*hr[2]);
//XXX: In stand-alone only Dirichlet for validation purposes
......@@ -594,13 +594,13 @@ void ArbitraryDomain::LinearInterpolation(int idx, int idy, int idz, double& W,
//C += 1/hr[2]*1/hr[2];
// case where we are on the Neumann BC in Z-direction
// where we distinguish two cases
// where we distinguish two cases
if(z == 0)
F = 0.0;
else
else
B = 0.0;
//for test no neumann
//for test no neumann
//C += 2/((hr[2]*nr[2]/2.0) * hr[2]);
//
// double d = hr[2]*(nr[2])/2;
......@@ -616,10 +616,10 @@ void ArbitraryDomain::LinearInterpolation(int idx, int idy, int idz, double& W,
scaleFactor *= 0.5;
} else
} else
C += 2*1/hr[2]*1/hr[2];
}
*/
*/
void ArbitraryDomain::getNeighbours(int id, int &W, int &E, int &S, int &N, int &F, int &B) {
......@@ -638,7 +638,7 @@ void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W, int &E, i
F = getIdx(idx, idy, idz - 1);
B = getIdx(idx, idy, idz + 1);
if(!isInside(idx+1,idy,idz))
if(!isInside(idx+1,idy,idz))
E = -1;
if(!isInside(idx-1,idy,idz))
......@@ -647,14 +647,14 @@ void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W, int &E, i
if(!isInside(idx,idy+1,idz))
N = -1;