ArbitraryDomain.h 6.93 KB
Newer Older
1 2 3 4 5
//
// Class ArbitraryDomain
//   Interface to iterative solver and boundary geometry
//   for space charge calculation
//
6
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
7 8
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//                      2016, Daniel Winklehner, Massachusetts Institute of Technology
9
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
10
// All rights reserved
11
//
12 13 14 15 16
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
frey_m's avatar
frey_m committed
17 18 19 20 21 22 23 24 25 26
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
27
//
28 29
#ifndef ARBITRARY_DOMAIN_H
#define ARBITRARY_DOMAIN_H
30

gsell's avatar
gsell committed
31 32 33
#include <mpi.h>
#include <hdf5.h>
#include "H5hut.h"
34

gsell's avatar
gsell committed
35
#include <map>
36
#include <string>
snuverink_j's avatar
snuverink_j committed
37 38
#include <tuple>
#include <vector>
39
#include "IrregularDomain.h"
snuverink_j's avatar
snuverink_j committed
40 41

class BoundaryGeometry;
gsell's avatar
gsell committed
42 43 44 45 46

class ArbitraryDomain : public IrregularDomain {

public:

frey_m's avatar
frey_m committed
47 48 49 50 51 52
    ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
                    std::string interpl);

    ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
                    Vector_t globalMeanR, Quaternion_t globalToLocalQuaternion,
                    std::string interpl);
53 54

    ~ArbitraryDomain();
gsell's avatar
gsell committed
55 56

    /// returns discretization at (x,y,z)
frey_m's avatar
frey_m committed
57 58 59 60
    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
61
    /// returns discretization at 3D index
frey_m's avatar
frey_m committed
62 63 64
    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
65
    /// returns index of neighbours at (x,y,z)
frey_m's avatar
frey_m committed
66 67 68
    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
69
    /// returns index of neighbours at 3D index
frey_m's avatar
frey_m committed
70 71 72
    void getNeighbours(int idxyz, int &W, int &E,
                       int &S, int &N, int &F, int &B);

gsell's avatar
gsell committed
73
    /// returns type of boundary condition
74
    std::string getType() {return "Geometric";}
gsell's avatar
gsell committed
75
    /// queries if a given (x,y,z) coordinate lies inside the domain
76
    bool isInside(int idx, int idy, int idz);
77
    /// returns number of nodes in xy plane
78
    int getNumXY(int idz);
79 80
    // calculates intersection with rotated and shifted geometry
    void compute(Vector_t hr, NDIndex<3> localId);
gsell's avatar
gsell committed
81

82
    int getStartId() {return startId;}
gsell's avatar
gsell committed
83

84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
    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; }
99

100 101

    bool hasGeometryChanged() { return hasGeometryChanged_m; }
gsell's avatar
gsell committed
102 103

private:
104 105
    BoundaryGeometry *bgeom_m;

frey_m's avatar
frey_m committed
106 107 108
    /** PointList maps from an (x,z) resp. (y,z) pair to double values
     * (=intersections with boundary)
     */
109
    typedef std::multimap< std::tuple<int, int, int>, double > PointList;
gsell's avatar
gsell committed
110 111

    /// all intersection points with gridlines in X direction
112 113 114 115 116
    PointList IntersectHiX, IntersectLoX;

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

gsell's avatar
gsell committed
117
    /// all intersection points with gridlines in Z direction
118
    PointList IntersectHiZ, IntersectLoZ;
119 120

    // meanR to shift from global to local frame
121
    Vector_t globalMeanR_m;
122
    //    Quaternion_t globalToLocalQuaternion_m;  because defined in parent class
123
    Quaternion_t localToGlobalQuaternion_m;
gsell's avatar
gsell committed
124

125
    int startId;
gsell's avatar
gsell committed
126

127
    // Here we store the number of nodes in a xy layer for a given z coordinate
gsell's avatar
gsell committed
128
    std::map<int, int> numXY;
129 130
    std::map<int, int> numYZ;
    std::map<int, int> numXZ;
gsell's avatar
gsell committed
131

132
    // Number of nodes in the xy plane (for this case: independent of the z coordinate)
133
    int nxy_m[1000];
134
    // mapping (x,y,z) -> idxyz
gsell's avatar
gsell committed
135
    std::map<int, int> IdxMap;
136
    // Mapping idxyz -> (x,y,z)
gsell's avatar
gsell committed
137
    std::map<int, int> CoordMap;
138 139 140 141

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

142
    // Interpolation type
143
    int interpolationMethod;
144

145
    // Flag indicating if geometry has changed for the current time-step
146 147 148 149
    bool hasGeometryChanged_m;

    Vector_t Geo_nr_m;
    Vector_t Geo_hr_m;
150
    Vector_t geomCentroid_m;
151 152
    Vector_t minCoords_m;
    Vector_t maxCoords_m;
153
    Vector_t globalInsideP0_m;
gsell's avatar
gsell committed
154

155 156 157 158 159
    // 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)
160
    inline void getCoord(int idxyz, int &x, int &y, int &z);
gsell's avatar
gsell committed
161 162

    inline void crossProduct(double A[], double B[], double C[]);
frey_m's avatar
frey_m committed
163 164 165 166

    inline double dotProduct(double v1[], double v2[]) {
        return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
    }
gsell's avatar
gsell committed
167

168
    // Different interpolation methods for boundary points
frey_m's avatar
frey_m committed
169 170 171 172 173 174 175 176 177 178 179
    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);
180 181

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

184 185 186
    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
187 188
};

189 190 191 192 193 194 195 196 197
#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: