EllipticDomain.h

Go to the documentation of this file.

#ifndef ELLIPTICAL_DOMAIN_H
#define ELLIPTICAL_DOMAIN_H
#ifdef HAVE_SAAMG_SOLVER

#include <vector>
#include <map>
#include <string>
#include <cmath>
#include "IrregularDomain.h"
#include "Structure/BoundaryGeometry.h"

class EllipticDomain : public IrregularDomain {

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, Vector_t hr, std::string interpl);
    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);

    ~EllipticDomain();

    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, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
    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 &W, int &E, int &S, int &N, int &F, int &B);
    std::string getType() {return "Elliptic";}
    inline bool isInside(int x, int y, int z) {
        double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
        double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
        return ((xx * xx / (SemiMajor * SemiMajor) + yy * yy / (SemiMinor * SemiMinor) < 1) && z != 0 && z != nr[2] - 1);
    }

    int getNumXY(int z) { return nxy_m; }
    void setSemiMinor(double sm) {SemiMinor = sm;}
    void setSemiMajor(double sm) {SemiMajor = sm;}

    void compute(Vector_t hr);
    void compute(Vector_t hr, NDIndex<3> localId);

    double getXRangeMin() { return -SemiMajor; }
    double getXRangeMax() { return SemiMajor;  }
    double getYRangeMin() { return -SemiMinor; }
    double getYRangeMax() { return SemiMinor;  }
    double getZRangeMin() { return zMin_m; }
    double getZRangeMax() { return zMax_m; }


    //TODO: ?
    int getStartIdx() {return 0;}

    bool hasGeometryChanged() { return hasGeometryChanged_m; }

private:

    typedef std::multimap<int, double> EllipticPointList;

    EllipticPointList IntersectXDir;

    EllipticPointList IntersectYDir;

    std::map<int, int> IdxMap;

    std::map<int, int> CoordMap;

    double SemiMajor;
    double SemiMinor;
    int nxy_m;
    int interpolationMethod;
    bool hasGeometryChanged_m;

    inline int toCoordIdx(int x, int y) { return y * nr[0] + x; }
    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 - 1) * nxy_m;
        else
            return -1;
    }
    inline void getCoord(int idx, int &x, int &y, int &z) {
        int ixy = idx % nxy_m;
        int xy = CoordMap[ixy];
        int inr = (int)nr[0];
        x = xy % inr;
        y = (xy - x) / nr[0];
        z = (idx - ixy) / nxy_m + 1;
    }

    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);

};

#endif //#ifdef HAVE_SAAMG_SOLVER
#endif //#ifdef ELLIPTICAL_DOMAIN_H