ArbitraryDomain.h 5.75 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Class ArbitraryDomain
//   Interface to iterative solver and boundary geometry
//   for space charge calculation
//
// Copyright (c) 2008-2020
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
//

13 14
#ifndef ARBITRARY_DOMAIN_H
#define ARBITRARY_DOMAIN_H
15

gsell's avatar
gsell committed
16 17 18
#include <mpi.h>
#include <hdf5.h>
#include "H5hut.h"
19

gsell's avatar
gsell committed
20
#include <map>
21
#include <string>
snuverink_j's avatar
snuverink_j committed
22 23
#include <tuple>
#include <vector>
24
#include "IrregularDomain.h"
snuverink_j's avatar
snuverink_j committed
25 26

class BoundaryGeometry;
gsell's avatar
gsell committed
27 28 29 30 31

class ArbitraryDomain : public IrregularDomain {

public:

32
    ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
33
    ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, Vector_t globalMeanR, Quaternion_t globalToLocalQuaternion, std::string interpl);
34 35

    ~ArbitraryDomain();
gsell's avatar
gsell committed
36 37

    /// returns discretization at (x,y,z)
38
    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
39
    /// returns discretization at 3D index
40
    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
41
    /// returns index of neighbours at (x,y,z)
42
    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
43
    /// returns index of neighbours at 3D index
44
    void getNeighbours(int idxyz, int &W, int &E, int &S, int &N, int &F, int &B);
gsell's avatar
gsell committed
45
    /// returns type of boundary condition
46
    std::string getType() {return "Geometric";}
gsell's avatar
gsell committed
47
    /// queries if a given (x,y,z) coordinate lies inside the domain
48
    bool isInside(int idx, int idy, int idz);
49
    /// returns number of nodes in xy plane
50
    int getNumXY(int idz);
51 52 53 54
    // calculates intersection
    void compute(Vector_t hr);
    // calculates intersection with rotated and shifted geometry
    void compute(Vector_t hr, NDIndex<3> localId);
gsell's avatar
gsell committed
55

56
    int getStartId() {return startId;}
gsell's avatar
gsell committed
57

58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
    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; }
73

74 75

    bool hasGeometryChanged() { return hasGeometryChanged_m; }
gsell's avatar
gsell committed
76 77

private:
78 79 80 81
    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
82 83

    /// all intersection points with gridlines in X direction
84 85 86 87 88
    PointList IntersectHiX, IntersectLoX;

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

gsell's avatar
gsell committed
89
    /// all intersection points with gridlines in Z direction
90
    PointList IntersectHiZ, IntersectLoZ;
91 92

    // meanR to shift from global to local frame
93
    Vector_t globalMeanR_m;
94
    //    Quaternion_t globalToLocalQuaternion_m;  because defined in parent class
95
    Quaternion_t localToGlobalQuaternion_m;
gsell's avatar
gsell committed
96

97
    int startId;
gsell's avatar
gsell committed
98

99
    // Here we store the number of nodes in a xy layer for a given z coordinate
gsell's avatar
gsell committed
100
    std::map<int, int> numXY;
101 102
    std::map<int, int> numYZ;
    std::map<int, int> numXZ;
gsell's avatar
gsell committed
103

104
    // Number of nodes in the xy plane (for this case: independent of the z coordinate)
105
    int nxy_m[1000];
106
    // mapping (x,y,z) -> idxyz
gsell's avatar
gsell committed
107
    std::map<int, int> IdxMap;
108
    // Mapping idxyz -> (x,y,z)
gsell's avatar
gsell committed
109
    std::map<int, int> CoordMap;
110 111 112 113

    // Mapping all cells that are inside the geometry
    std::map<int, bool> IsInsideMap;

114
    // Interpolation type
115
    int interpolationMethod;
116

117
    // Flag indicating if geometry has changed for the current time-step
118 119 120 121
    bool hasGeometryChanged_m;

    Vector_t Geo_nr_m;
    Vector_t Geo_hr_m;
122
    Vector_t geomCentroid_m;
123 124
    Vector_t minCoords_m;
    Vector_t maxCoords_m;
125
    Vector_t globalInsideP0_m;
gsell's avatar
gsell committed
126

127 128 129 130 131
    // 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)
132
    inline void getCoord(int idxyz, int &x, int &y, int &z);
gsell's avatar
gsell committed
133 134 135 136

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

137
    // Different interpolation methods for boundary points
138 139 140
    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);
141 142

    // Rotate positive axes with quaternion -DW
143
    inline void rotateWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
144

145 146 147
    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);
gsell's avatar
gsell committed
148 149
};

150 151 152 153 154 155 156 157 158
#endif

// vi: set et ts=4 sw=4 sts=4:
// Local Variables:
// mode:c
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End: