Commit 8451cf66 authored by frey_m's avatar frey_m
Browse files

SAAMG: const functions, clean up + nr_m is int vector

parent b6b06cb5
......@@ -38,7 +38,7 @@
#include "Utilities/OpalException.h"
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
Vector_t nr,
IntVector_t nr,
Vector_t hr,
std::string interpl)
: IrregularDomain(nr, hr, interpl)
......@@ -242,23 +242,8 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
INFOMSG(level2 << "* Done." << endl);
}
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
return (idz * nr_m[1] + idy) * nr_m[0] + idx;
}
int ArbitraryDomain::indexAccess(int x, int y, int z) {
return idxMap_m[toCoordIdx(x, y, z)];
}
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
return isInsideMap_m[toCoordIdx(idx, idy, idz)];
}
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
StencilValue_t& value, double& /*scaleFactor*/)
StencilValue_t& value, double& /*scaleFactor*/) const
{
value.west = -1/(hr_m[0]*hr_m[0]);
value.east = -1/(hr_m[0]*hr_m[0]);
......@@ -285,7 +270,7 @@ void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
}
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
StencilValue_t& value, double &scaleFactor)
StencilValue_t& value, double &scaleFactor) const
{
scaleFactor = 1;
......
......@@ -43,16 +43,16 @@ class BoundaryGeometry;
class ArbitraryDomain : public IrregularDomain {
public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
ArbitraryDomain(BoundaryGeometry *bgeom, IntVector_t nr, Vector_t hr,
std::string interpl);
~ArbitraryDomain();
/// queries if a given (x,y,z) coordinate lies inside the domain
bool isInside(int idx, int idy, int idz);
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);
......@@ -82,17 +82,25 @@ private:
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 indexAccess(int x, int y, int z);
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y, z));
}
int coordAccess(int idx) const {
return coordMap_m.at(idx);
}
// Different interpolation methods for boundary points
void constantInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
void linearInterpolation(int idx, int idy, int idz, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
};
#endif
......
......@@ -36,7 +36,7 @@
extern Inform *gmsg;
BoxCornerDomain::BoxCornerDomain(double A, double B, double C, double length,
double L1, double L2, Vector_t nr, Vector_t hr,
double L1, double L2, IntVector_t nr, Vector_t hr,
std::string interpl)
: IrregularDomain(nr, hr, interpl)
{
......@@ -129,7 +129,7 @@ void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
}
void BoxCornerDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
scaleFactor = 1.0;
......@@ -190,7 +190,7 @@ void BoxCornerDomain::constantInterpolation(int x, int y, int z, StencilValue_t&
}
void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
scaleFactor = 1.0;
......@@ -294,7 +294,7 @@ void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& v
//FIXME: this probably needs some cleanup/rewriting
void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
double cx = (x - (nr_m[0] - 1) / 2.0) * hr_m[0];
double cy = (y - (nr_m[1] - 1) / 2.0) * hr_m[1];
......
......@@ -74,9 +74,6 @@ L2 = getZRangeMax() - getZRangeMin
class BoxCornerDomain : public IrregularDomain {
public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
/**
* \param A depth of the box
* \param B maximal height of the box
......@@ -86,12 +83,12 @@ public:
* \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,
double L1, double L2, IntVector_t nr, Vector_t hr,
std::string interpl);
~BoxCornerDomain();
/// as a function of z, determine the hight (B) of the geometry
inline double getB(double z) {
inline double getB(double z) const {
if((z < getZRangeMin()) || (z > getZRangeMax()))
return getYRangeMax();
else
......@@ -99,11 +96,11 @@ public:
}
/// 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 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] - 1);
return (xx < getXRangeMax() && yy < b && z >= 0 && z < nr_m[2]);
}
void compute(Vector_t hr, NDIndex<3> localId);
......@@ -136,33 +133,37 @@ private:
inline double getXIntersection(double cx, int /*z*/) {
inline double getXIntersection(double cx, int /*z*/) const {
return (cx < 0) ? getXRangeMin() : getXRangeMax();
}
inline double getYIntersection(double cy, int z) {
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) {
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
int indexAccess(int x, int y, int z) {
return idxMap_m[toCoordIdx(x, y, z)];
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y, z));
}
int coordAccess(int idx) const {
return coordMap_m.at(idx);
}
/// different interpolation methods for boundary points
void constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
void linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
};
......
......@@ -35,7 +35,7 @@
//FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr, Vector_t hr,
std::string interpl)
: IrregularDomain(nr, hr, interpl)
{
......@@ -51,10 +51,16 @@ EllipticDomain::~EllipticDomain() {
}
int EllipticDomain::getNumXY() const {
return nxy_m;
}
// 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
void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
// there is nothing to be done if the mesh spacings in x and y have not changed
if (hr[0] == getHr()[0] &&
hr[1] == getHr()[1])
......@@ -84,8 +90,8 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
* grid points per plane --> otherwise we might
* get not unique global indices in the Tpetra::CrsMatrix
*/
for (x = 0; x < nr_m[0]; ++x) {
for (y = 0; y < nr_m[1]; ++y) {
for (y = 0; y < nr_m[1]; ++y) {
for (x = 0; x < nr_m[0]; ++x) {
if (isInside(x, y, 1)) {
idxMap_m[toCoordIdx(x, y)] = idx;
coordMap_m[idx++] = toCoordIdx(x, y);
......@@ -155,7 +161,7 @@ void EllipticDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t&
}
void EllipticDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
scaleFactor = 1.0;
......@@ -191,7 +197,7 @@ void EllipticDomain::constantInterpolation(int x, int y, int z, StencilValue_t&
}
void EllipticDomain::linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
scaleFactor = 1.0;
......@@ -199,7 +205,7 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, StencilValue_t& va
double cy = getYRangeMin() + hr_m[1] * (y + 0.5);
double dx = 0.0;
std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
std::multimap<int, double>::const_iterator it = intersectXDir_m.find(y);
if (cx < 0)
++it;
......@@ -263,7 +269,7 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, StencilValue_t& va
void EllipticDomain::quadraticInterpolation(int x, int y, int z,
StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
scaleFactor = 1.0;
......@@ -273,7 +279,7 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z,
// since every vector for elliptic domains has ALWAYS size 2 we
// can catch all cases manually
double dx = 0.0;
std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
std::multimap<int, double>::const_iterator it = intersectXDir_m.find(y);
if (cx < 0)
++it;
dx = it->second;
......@@ -346,7 +352,7 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z,
robinBoundaryStencil(z, value.front, value.back, value.center);
}
void EllipticDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) {
void EllipticDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) const {
if (z == 0 || z == nr_m[2] - 1) {
// case where we are on the Robin BC in Z-direction
......
......@@ -39,30 +39,29 @@
class EllipticDomain : public IrregularDomain {
public:
using IrregularDomain::StencilIndex_t;
using IrregularDomain::StencilValue_t;
EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr,
EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr,
Vector_t hr, std::string interpl);
~EllipticDomain();
int getNumXY() const override;
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z) {
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 / (getXRangeMax() * getXRangeMax()) +
yy * yy / (getYRangeMax() * getYRangeMax()) < 1);
return (isInsideEllipse && z > 0 && z < nr_m[2] - 1);
return (isInsideEllipse && z >= 0 && z < nr_m[2]);
}
/// calculates intersection
void compute(Vector_t hr, NDIndex<3> localId);
void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
const Vector_t& rmax, double dh);
const Vector_t& rmax, double dh) override;
private:
......@@ -76,44 +75,34 @@ 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;
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
/// conversion from (x,y) to index in xy plane
inline int toCoordIdx(int x, int y) { return y * nr_m[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
int indexAccess(int x, int y, int z) {
return idxMap_m[toCoordIdx(x, y)] + (z - 1) * nxy_m;
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y)) + z * nxy_m;
}
/// conversion from a 3D index to (x,y,z)
void getCoord(int idx, int &x, int &y, int &z) override {
int coordAccess(int idx) const {
int ixy = idx % nxy_m;
int xy = coordMap_m[ixy];
x = xy % (int)nr_m[0];
y = (xy - x) / nr_m[0];
z = (idx - ixy) / nxy_m + 1;
return coordMap_m.at(ixy);
}
/// different interpolation methods for boundary points
void constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
void linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) override;
double &scaleFactor) const override;
/// function to handle the open boundary condition in longitudinal direction
void robinBoundaryStencil(int z, double &F, double &B, double &C);
void robinBoundaryStencil(int z, double &F, double &B, double &C) const;
};
#endif //#ifdef ELLIPTICAL_DOMAIN_H
......
......@@ -29,7 +29,7 @@
#include <cassert>
IrregularDomain::IrregularDomain(const Vector_t& nr, const Vector_t& hr,
IrregularDomain::IrregularDomain(const IntVector_t& nr, const Vector_t& hr,
const std::string& interpl)
: nr_m(nr)
, hr_m(hr)
......@@ -46,7 +46,7 @@ IrregularDomain::IrregularDomain(const Vector_t& nr, const Vector_t& hr,
}
void IrregularDomain::getNeighbours(int x, int y, int z, StencilIndex_t& index)
void IrregularDomain::getNeighbours(int x, int y, int z, StencilIndex_t& index) const
{
index.west = -1;
index.east = -1;
......@@ -75,32 +75,28 @@ void IrregularDomain::getNeighbours(int x, int y, int z, StencilIndex_t& index)
}
void IrregularDomain::getNeighbours(int id, StencilIndex_t& index)
{
void IrregularDomain::getNeighbours(int id, StencilIndex_t& index) const {
int x = 0, y = 0, z = 0;
getCoord(id, x, y, z);
getNeighbours(x, y, z, index);
}
void IrregularDomain::getCoord(int idx, int& x, int& y, int& z) {
int id = coordMap_m[idx];
x = id % (int)nr_m[0];
id /= nr_m[0];
y = id % (int)nr_m[1];
id /= nr_m[1];
z = id;
void IrregularDomain::getCoord(int idx, int& x, int& y, int& z) const {
int xy = coordAccess(idx);
x = xy % nr_m[0];
y = (xy - x) / nr_m[0];
z = idx / getNumXY();
}
int IrregularDomain::getIdx(int x, int y, int z) {
int IrregularDomain::getIdx(int x, int y, int z) const {
if (x < 0 || y < 0 || z < 0 || !isInside(x, y, z))
return -1;
return indexAccess(x, y, z);
}
void IrregularDomain::getBoundaryStencil(int x, int y, int z, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
scaleFactor = 1.0;
......@@ -123,7 +119,7 @@ void IrregularDomain::getBoundaryStencil(int x, int y, int z, StencilValue_t& va
void IrregularDomain::getBoundaryStencil(int id, StencilValue_t& value,
double &scaleFactor)
double &scaleFactor) const
{
int idx = 0, idy = 0, idz = 0;
getCoord(id, idx, idy, idz);
......@@ -142,7 +138,7 @@ void IrregularDomain::resizeMesh(Vector_t& origin, Vector_t& hr,
void IrregularDomain::constantInterpolation(int /*idx*/, int /*idy*/, int /*idz*/,
StencilValue_t& /*value*/, double &/*scaleFactor*/)
StencilValue_t& /*value*/, double &/*scaleFactor*/) const
{
throw OpalException("IrregularDomain::constantInterpolation()",
"No constant interpolation method Implemented!");
......@@ -150,7 +146,7 @@ void IrregularDomain::constantInterpolation(int /*idx*/, int /*idy*/, int /*idz*
void IrregularDomain::linearInterpolation(int /*idx*/, int /*idy*/, int /*idz*/,
StencilValue_t& /*value*/, double &/*scaleFactor*/)
StencilValue_t& /*value*/, double &/*scaleFactor*/) const
{
throw OpalException("IrregularDomain::linearInterpolation()",
"No linear interpolation method Implemented!");
......@@ -158,7 +154,7 @@ void IrregularDomain::linearInterpolation(int /*idx*/, int /*idy*/, int /*idz*/,
void IrregularDomain::quadraticInterpolation(int /*idx*/, int /*idy*/, int /*idz*/,
StencilValue_t& /*value*/, double &/*scaleFactor*/)
StencilValue_t& /*value*/, double &/*scaleFactor*/) const
{
throw OpalException("IrregularDomain::quadraticInterpolation()",
"No quadratic interpolation method Implemented!");
......
......@@ -55,8 +55,9 @@ public:
typedef Stencil<int> StencilIndex_t;
typedef Stencil<double> StencilValue_t;
typedef Vektor<int, 3> IntVector_t;
IrregularDomain(const Vector_t& nr,
IrregularDomain(const IntVector_t& nr,
const Vector_t& hr,
const std::string& interpl);
......@@ -75,57 +76,64 @@ public:
/// \param scaleFactor of stencil values
void getBoundaryStencil(int x, int y, int z,
StencilValue_t& value,
double &scaleFactor);
double &scaleFactor) const;
/// method to calculate the stencil at a boundary points
/// \param id index of the current element in the matrix
// \param values of stencil element
/// \param scaleFactor of stencil values
void getBoundaryStencil(int id, StencilValue_t& value,
double &scaleFactor);
double &scaleFactor) const;
/// method to calculate the neighbours in the matrix of the current index (x,y,z)
/// \param x index of the current element in the matrix
/// \param y index of the current element in the matrix
/// \param z index of the current element in the matrix
/// \param index stencil indices of an element
void getNeighbours(int x, int y, int z, StencilIndex_t& index);
void getNeighbours(int x, int y, int z, StencilIndex_t& index) const;
void getNeighbours(int idx, StencilIndex_t& index);
void getNeighbours(int idx, StencilIndex_t& index) const;
/// Conversion from a 3D index to (x,y,z)
virtual void getCoord(int idx, int &x, int &y, int &z);
virtual void getCoord(int idx, int &x, int &y, int &z) const;
/// Conversion from (x,y,z) to index on the 3D grid
int getIdx(int x, int y, int z);
int getIdx(int x, int y, int z) const;
virtual int getNumXY() const { return nr_m[0] * nr_m[1]; }
/// method that checks if a given point lies inside the boundary
/// \param x index of the current element in the matrix
/// \param y index of the current element in the matrix
/// \param z index of the current element in the matrix
/// \return boolean indicating if the point lies inside the boundary
virtual bool isInside(int x, int y, int z) = 0;
virtual bool isInside(int x, int y, int z) const = 0;
IntVector_t getNr() const { return nr_m; }
Vector_t getHr() const { return hr_m; }
Vector_t getNr() { return nr_m; }
Vector_t getHr() { return hr_m; }
void setNr(Vector_t nr) { nr_m = nr; }
void setHr(Vector_t hr) { hr_m = hr; }
void setNr(IntVector_t nr) { nr_m = nr; }
void setHr(Vector_t hr) { hr_m = hr; }
void setMinMaxZ(double minz, double maxz) { zMin_m=minz; zMax_m=maxz; }
double getMinZ() { return zMin_m; }
double getMaxZ() { return zMax_m; }
void setMinMaxZ(double minz, double maxz) {
zMin_m = minz;
zMax_m = maxz;
}
double getXRangeMin();
double getXRangeMax();
double getYRangeMin();
double getYRangeMax();
double getZRangeMin();
double getZRangeMax();
double getMinZ() const { return zMin_m; }
double getMaxZ() const { return zMax_m; }
void setRangeMin(const Vector_t&);
void setRangeMax(const Vector_t&);
double getXRangeMin() const { return min_m(0); }
double getXRangeMax() const { return max_m(0); }
double getYRangeMin() const { return min_m(1); }
double getYRangeMax() const { return max_m(1); }
double getZRangeMin() const { return min_m(2); }
double getZRangeMax() const { return max_m(2); }
bool hasGeometryChanged();
void setRangeMin(const Vector_t& min) { min_m = min; }
void setRangeMax(const Vector_t& max) { max_m = max; }
bool hasGeometryChanged() const { return hasGeometryChanged_m; }
virtual ~IrregularDomain() {};
......@@ -135,22 +143,24 @@ public:
protected:
virtual int indexAccess(int x, int y, int z) = 0;
virtual int indexAccess(int x, int y, int z) const = 0;
virtual int coordAccess(int idx) const = 0;
/// different interpolation methods for boundary points
virtual void constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor);
double &scaleFactor) const;
virtual void linearInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor);
double &scaleFactor) const;