BoxCornerDomain.h

Go to the documentation of this file.

#ifndef BOXCORNER_DOMAIN_H
#define BOXCORNER_DOMAIN_H
#ifdef HAVE_SAAMG_SOLVER

//#include <vector>
#include <map>
//#include <multimap>
#include <string>
#include <cmath>
#include <iostream>  // Neeeded for stream I/O
#include <fstream>   // Needed for file I/O
//#include <iomanip>   // Needed for I/O manipulators

//#include "Ippl.h"
#include "IrregularDomain.h"


/*

    A_m and B_m are the half apperture of the box


                                     / (A_m,B_m)
                                    /
                                   /
                                  /
    L1_m                         /
------------      --------------+ (-A_m,B_m)
           | L2_m |             |
        C_m|      |             |
           |------|             |      /
         .....                  |     /
(0,0)---.......-----------------+    /
         .....                  |   /
   z                            |  /
   |                            | /
--------------------------------+/ (-A_m,-B_m)

            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;
    else
        b = B_m-C_m;

*/

class BoxCornerDomain : public IrregularDomain {

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,
                    std::string interpl);
    ~BoxCornerDomain();

    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 "BoxCorner";}


    int getNumXY(int z) { return -1;}


    inline double getB(double z) {
      if((z < L1_m) || (z > (L1_m + L2_m)))
            return B_m;
        else
            return B_m - C_m;
    }

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

    //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 -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,
    //int) -> intersection. To simplify things (for now) we use the same
    //structure for X...
    typedef std::multimap< std::pair<int, int>, double > BoxCornerPointList;

    BoxCornerPointList IntersectXDir;

    BoxCornerPointList IntersectYDir;

    std::map<int, int> IdxMap;

    std::map<int, int> CoordMap;

    double A_m;

    double B_m;


    double actBMin_m;

    double actBMax_m;

    double C_m;

    double Length_m;

    double L1_m;

    double L2_m;

    // :FIXME: unused
    //double SemiMajor;
    // :FIXME: unused
    //double SemiMinor;

    int interpolationMethod;
    bool hasGeometryChanged_m;

   std::ofstream os_m;





    inline double getXIntersection(double cx, int z) {
        if(cx < 0)
            return -A_m;
        else
            return  A_m;
    }

    inline double getYIntersection(double cy, int z) {
        if(cy < 0)
            return -B_m;
        else
            return getB(z * hr[2]);
    }

    inline int toCoordIdx(int x, int y, int z) {
        return (z * nr[1] + y) * nr[0] + x;
    }

    /*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;
    }

    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;

    }

    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 BOXCORNER_DOMAIN_H