ArbitraryDomain.h 5.61 KB
Newer Older
1 2
#ifndef ARBITRARY_DOMAIN_H
#define ARBITRARY_DOMAIN_H
3

4
#ifdef HAVE_SAAMG_SOLVER
gsell's avatar
gsell committed
5 6 7 8

#include <mpi.h>
#include <hdf5.h>
#include "H5hut.h"
9

gsell's avatar
gsell committed
10 11
#include <vector>
#include <map>
12
#include <string>
gsell's avatar
gsell committed
13
#include <math.h>
14 15 16
#include <cmath>
#include "IrregularDomain.h"
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
17 18 19 20 21

class ArbitraryDomain : public IrregularDomain {

public:

22
    ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
23
    ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion, std::string interpl);
24 25

    ~ArbitraryDomain();
gsell's avatar
gsell committed
26 27

    /// returns discretization at (x,y,z)
28
    void getBoundaryStencil(int idx, int idy, int idz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
gsell's avatar
gsell committed
29
    /// returns discretization at 3D index
30
    void getBoundaryStencil(int idxyz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
gsell's avatar
gsell committed
31
    /// returns index of neighbours at (x,y,z)
32
    void getNeighbours(int idx, int idy, int idz, int &W, int &E, int &S, int &N, int &F, int &B);
gsell's avatar
gsell committed
33
    /// returns index of neighbours at 3D index
34
    void getNeighbours(int idxyz, int &W, int &E, int &S, int &N, int &F, int &B);
gsell's avatar
gsell committed
35
    /// returns type of boundary condition
36
    std::string getType() {return "Geometric";}
gsell's avatar
gsell committed
37
    /// queries if a given (x,y,z) coordinate lies inside the domain
38
    inline bool isInside(int idx, int idy, int idz); 
gsell's avatar
gsell committed
39

40
    /// returns number of nodes in xy plane
41
    int getNumXY(int idz);
42 43 44
    /// calculates intersection 
    void Compute(Vector_t hr);
    void Compute(Vector_t hr, NDIndex<3> localId);
45 46
    // calculates intersection with rotated and shifted geometry 
    void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
gsell's avatar
gsell committed
47

48
    int getStartId() {return startId;}
gsell's avatar
gsell committed
49

50 51 52 53 54 55 56 57 58 59 60 61 62
    double getXRangeMin(){ return intersectMinCoords_m(0); }
    double getXRangeMax(){ return intersectMaxCoords_m(0); }
    double getYRangeMin(){ return intersectMinCoords_m(1); }
    double getYRangeMax(){ return intersectMaxCoords_m(1); }
    double getZRangeMin(){ return intersectMinCoords_m(2); }
    double getZRangeMax(){ return intersectMaxCoords_m(2); }
    void setXRangeMin(double xmin){ intersectMinCoords_m(0) = xmin; }
    void setXRangeMax(double xmax){ intersectMaxCoords_m(0) = xmax; }
    void setYRangeMin(double ymin){ intersectMinCoords_m(1) = ymin; }
    void setYRangeMax(double ymax){ intersectMaxCoords_m(1) = ymax; }
    void setZRangeMin(double zmin){ intersectMinCoords_m(2) = zmin; }
    void setZRangeMax(double zmax){ intersectMaxCoords_m(2) = zmax; }

63 64

    bool hasGeometryChanged() { return hasGeometryChanged_m; }
gsell's avatar
gsell committed
65 66

private:
67 68 69 70
    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;
gsell's avatar
gsell committed
71 72

    /// all intersection points with gridlines in X direction
73 74 75 76 77
    PointList IntersectHiX, IntersectLoX;

    /// all intersection points with gridlines in Y direction
    PointList IntersectHiY, IntersectLoY;

gsell's avatar
gsell committed
78
    /// all intersection points with gridlines in Z direction
79
    PointList IntersectHiZ, IntersectLoZ;
80 81 82 83 84
    
    // meanR to shift from global to local frame 
    Vector_t globalMeanR_m;
    Vektor<double, 4> globalToLocalQuaternion_m;
    Vektor<double, 4> localToGlobalQuaternion_m;
gsell's avatar
gsell committed
85

86
    int startId;
gsell's avatar
gsell committed
87

88
    // Here we store the number of nodes in a xy layer for a given z coordinate
gsell's avatar
gsell committed
89
    std::map<int, int> numXY;
90 91
    std::map<int, int> numYZ;
    std::map<int, int> numXZ;
gsell's avatar
gsell committed
92

93
    // Number of nodes in the xy plane (for this case: independent of the z coordinate)
94
    int nxy_m;
95
    // mapping (x,y,z) -> idxyz
gsell's avatar
gsell committed
96
    std::map<int, int> IdxMap;
97
    // Mapping idxyz -> (x,y,z)
gsell's avatar
gsell committed
98
    std::map<int, int> CoordMap;
99
    // Interpolation type
100 101
    int interpolationMethod;
 
102
    // Flag indicating if geometry has changed for the current time-step
103 104 105 106 107 108
    bool hasGeometryChanged_m;

    Vector_t Geo_nr_m;
    Vector_t Geo_hr_m;
    Vector_t Geo_mincoords_m;
    Vector_t Geo_maxcoords_m;
109 110
    Vector_t intersectMinCoords_m;
    Vector_t intersectMaxCoords_m;
gsell's avatar
gsell committed
111

112 113 114 115 116
    // Conversion from (x,y,z) to index in xyz plane
    inline int toCoordIdx(int idx, int idy, int idz);
    // 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)
117
    inline void getCoord(int idxyz, int &x, int &y, int &z);
gsell's avatar
gsell committed
118 119 120 121

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

122 123 124 125 126 127 128 129 130 131 132
    // 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, Vektor<double, 4> const quaternion);
 	
    inline void rotateXAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
    inline void rotateYAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
    inline void rotateZAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
gsell's avatar
gsell committed
133 134
};

135
#endif //#ifdef HAVE_SAAMG_SOLVER
gsell's avatar
gsell committed
136
#endif //#ifdef ARBITRARY_DOMAIN