Commit 123cd63c authored by frey_m's avatar frey_m

Merge branch '565-follow-up-from-resolve-saamg-enable-rectangulardomain' into 'master'

Resolve "Code duplication in Domains"

Closes #565

See merge request !396
parents 99361fd6 f48154af
This diff is collapsed.
......@@ -43,123 +43,64 @@ class BoundaryGeometry;
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, Quaternion_t globalToLocalQuaternion, std::string interpl);
ArbitraryDomain(BoundaryGeometry *bgeom, IntVector_t nr, Vector_t hr,
std::string interpl);
~ArbitraryDomain();
/// returns discretization at (x,y,z)
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 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
bool isInside(int idx, int idy, int idz);
/// returns number of nodes in xy plane
int getNumXY(int idz);
// calculates intersection
void compute(Vector_t hr);
bool isInside(int idx, int idy, int idz) const {
return isInsideMap_m.at(toCoordIdx(idx, idy, idz));
}
// calculates intersection with rotated and shifted geometry
void compute(Vector_t hr, NDIndex<3> localId);
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; }
bool hasGeometryChanged() { return hasGeometryChanged_m; }
private:
BoundaryGeometry *bgeom_m;
/// PointList maps from an (x,z) resp. (y,z) pair to double values (=intersections with boundary)
typedef std::multimap< std::tuple<int, int, int>, double > PointList;
/** PointList_t maps from an (x,z) resp. (y,z) pair to double values
* (=intersections with boundary)
*/
typedef std::multimap< std::tuple<int, int, int>, double > PointList_t;
/// all intersection points with gridlines in X direction
PointList IntersectHiX, IntersectLoX;
PointList_t intersectHiX_m, intersectLoX_m;
/// all intersection points with gridlines in Y direction
PointList IntersectHiY, IntersectLoY;
PointList_t intersectHiY_m, intersectLoY_m;
/// all intersection points with gridlines in Z direction
PointList IntersectHiZ, IntersectLoZ;
// meanR to shift from global to local frame
Vector_t globalMeanR_m;
// Quaternion_t globalToLocalQuaternion_m; because defined in parent class
Quaternion_t localToGlobalQuaternion_m;
int startId;
PointList_t intersectHiZ_m, intersectLoZ_m;
// 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)
int nxy_m[1000];
// mapping (x,y,z) -> idxyz
std::map<int, int> IdxMap;
// Mapping idxyz -> (x,y,z)
std::map<int, int> CoordMap;
std::map<int, int> numXY_m;
// Mapping all cells that are inside the geometry
std::map<int, bool> IsInsideMap;
std::map<int, bool> isInsideMap_m;
// Interpolation type
int interpolationMethod;
// Flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
Vector_t Geo_nr_m;
Vector_t Geo_hr_m;
Vector_t geomCentroid_m;
Vector_t minCoords_m;
Vector_t maxCoords_m;
Vector_t globalInsideP0_m;
// Conversion from (x,y,z) to index in xyz plane
inline int toCoordIdx(int idx, int idy, int idz);
int toCoordIdx(int idx, int idy, int idz) const {
return (idz * nr_m[1] + idy) * nr_m[0] + idx;
}
// 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);
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y, 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]); }
int coordAccess(int idx) const {
return coordMap_m.at(idx);
}
// 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, Quaternion_t const quaternion);
void constantInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &scaleFactor) const override;
inline void rotateXAxisWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
inline void rotateYAxisWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
inline void rotateZAxisWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
void linearInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &scaleFactor) const override;
};
#endif
......
This diff is collapsed.
......@@ -29,111 +29,80 @@
#include <map>
#include <string>
#include <cmath>
#include <iostream> // Neeeded for stream I/O
#include <fstream> // Needed for file I/O
#include "IrregularDomain.h"
#include "Solvers/RegularDomain.h"
/*
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 = getXRangeMax()
B = getYRangeMax()
L1 = getZRangeMin()
L2 = getZRangeMax() - getZRangeMin
*/
class BoxCornerDomain : public IrregularDomain {
class BoxCornerDomain : public RegularDomain {
public:
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, IntVector_t nr, Vector_t hr,
std::string interpl);
~BoxCornerDomain();
/// 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);
/// returns discretization at 3D index
void getBoundaryStencil(int idx, 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);
/// returns index of neighbours at 3D index
void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns type of boundary condition
std::string getType() {return "BoxCorner";}
/// we do not need this
int getNumXY(int /*z*/) { return -1;}
/// 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;
inline double getB(double z) const {
if((z < getZRangeMin()) || (z > getZRangeMax()))
return getYRangeMax();
else
return B_m - C_m;
return getYRangeMax() - C_m;
}
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z) {
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);
inline bool isInside(int x, int y, int z) const {
const double xx = (x - (nr_m[0] - 1) / 2.0) * hr_m[0];
const double yy = (y - (nr_m[1] - 1) / 2.0) * hr_m[1];
const double b = getB(z * hr_m[2]);
return (xx < getXRangeMax() && yy < b && z >= 0 && z < nr_m[2]);
}
/// set semi-minor
//void setSemiMinor(double sm) {SemiMinor = sm;}
/// set semi-major
//void setSemiMajor(double sm) {SemiMajor = sm;}
void compute(Vector_t hr);
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;}
bool hasGeometryChanged() { return hasGeometryChanged_m; }
private:
//XXX: since the Y coorindate is dependent on the Z value we need (int,
......@@ -149,100 +118,47 @@ private:
/// all intersection points with grid lines in Y direction
BoxCornerPointList IntersectYDir;
/// mapping (x,y,z) -> idx
std::map<int, int> IdxMap;
/// 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;
/// interpolation type
int interpolationMethod;
/// flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
/// for debug reasons
std::ofstream os_m;
/// length of the structure
double length_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;
inline double getXIntersection(double cx, int /*z*/) const {
return (cx < 0) ? getXRangeMin() : getXRangeMax();
}
inline double getYIntersection(double cy, int z) {
if(cy < 0)
return -B_m;
else
return getB(z * hr[2]);
inline double getYIntersection(double cy, int z) const {
return (cy < 0) ? getYRangeMin() : getB(z * hr_m[2]);
}
/// conversion from (x,y,z) to index in xyz plane
inline int toCoordIdx(int x, int y, int z) {
return (z * nr[1] + y) * nr[0] + x;
inline int toCoordIdx(int x, int y, int z) const {
return (z * nr_m[1] + y) * nr_m[0] + x;
}
/// conversion from (x,y,z) to index on the 3D grid
/*inline*/
inline int getIdx(int x, int y, int z) {
if(isInside(x, y, z) && x >= 0 && y >= 0 && z >= 0)
return IdxMap[toCoordIdx(x, y, z)];
else
return -1;
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y, z));
}
/// conversion from a 3D index to (x,y,z)
inline void getCoord(int idx, int &x, int &y, int &z) {
int idxx = CoordMap[idx];
x = idxx % (int)nr[0];
idxx /= nr[0];
y = idxx % (int)nr[1];
idxx /= nr[1];
z = idxx;
int coordAccess(int idx) const {
return coordMap_m.at(idx);
}
/// 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);
void linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) const override;
void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) const override;
};
......
......@@ -7,7 +7,6 @@ set (_SRCS
set (HDRS
FFTBoxPoissonSolver.h
FFTPoissonSolver.h
IrregularDomain.h
P3MPoissonSolver.h
PoissonSolver.h
)
......@@ -17,16 +16,20 @@ if (ENABLE_SAAMG_SOLVER)
ArbitraryDomain.cpp
BoxCornerDomain.cpp
EllipticDomain.cpp
IrregularDomain.cpp
MGPoissonSolver.cpp
RectangularDomain.cpp
RegularDomain.cpp
)
list (APPEND HDRS
ArbitraryDomain.h
BoxCornerDomain.h
EllipticDomain.h
IrregularDomain.h
MGPoissonSolver.h
RectangularDomain.h
RegularDomain.h
)
endif ()
......
This diff is collapsed.
......@@ -32,81 +32,32 @@
#include <map>
#include <string>
#include <cmath>
#include "IrregularDomain.h"
#include "Solvers/RegularDomain.h"
#include "Structure/BoundaryGeometry.h"
#include "Utilities/OpalException.h"
class EllipticDomain : public IrregularDomain {
class EllipticDomain : public RegularDomain {
public:
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,
EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr,
Vector_t hr, std::string interpl);
EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
~EllipticDomain();
/// 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);
/// returns discretization at 3D index
void getBoundaryStencil(int idx, 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);
/// returns index of neighbours at 3D index
void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns type of boundary condition
std::string getType() {return "Elliptic";}
/// 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);
bool isInside(int x, int y, int z) const {
double xx = getXRangeMin() + hr_m[0] * (x + 0.5);
double yy = getYRangeMin() + hr_m[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);
return (isInsideEllipse && z >= 0 && z < nr_m[2]);
}
int getNumXY(int /*z*/) { return nxy_m; }
/// calculates intersection
void compute(Vector_t /*hr*/) {
throw OpalException("EllipticDomain::compute()", "This function is not available.");
}
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; }
bool hasGeometryChanged() { return hasGeometryChanged_m; }
void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
const Vector_t& rmax, double dh);
private:
/// Map from a single coordinate (x or y) to a list of intersection values with
......@@ -119,66 +70,25 @@ private:
/// all intersection points with grid lines in Y direction
EllipticPointList_t intersectYDir_m;
/// mapping (x,y,z) -> idx
std::map<int, int> idxMap_m;
/// 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;
/// interpolation type
int interpolationMethod_m;
/// flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
/// conversion from (x,y) to index in xy plane
inline int toCoordIdx(int x, int y) { return y * nr[0] + x; }
int toCoordIdx(int x, int y) const { return y * nr_m[0] + x; }
/// conversion from (x,y,z) to index on the 3D grid
inline int getIdx(int x, int y, int z) {
if (isInside(x, y, z))
return idxMap_m[toCoordIdx(x, y)] + (z - 1) * nxy_m;
else
return -1;
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y)) + z * getNumXY();
}
/// conversion from a 3D index to (x,y,z)
inline void getCoord(int idx, int &x, int &y, int &z) {
int ixy = idx % nxy_m;
int xy = coordMap_m[ixy];
int inr = (int)nr[0];
x = xy % inr;
y = (xy - x) / nr[0];
z = (idx - ixy) / nxy_m + 1;
int coordAccess(int idx) const {
int ixy = idx % getNumXY();
return coordMap_m.at(ixy);
}
/// 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);
/// function to handle the open boundary condition in longitudinal direction
void robinBoundaryStencil(int z, double &F, double &B, double &C);
void linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) const override;
void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) const override;
};
#endif //#ifdef ELLIPTICAL_DOMAIN_H
......
//
// Class IrregularDomain
// Defines a common abstract interface for different types of boundaries.
//
// Copyright (c) 2008, Yves Ineichen, ETH Zürich,
// 2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
// 2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Solvers/IrregularDomain.h"
#include "Utilities/OpalException.h"
#include <cassert>
IrregularDomain::IrregularDomain(const IntVector_t& nr, const Vector_t& hr,
const std::string& interpl)