IrregularDomain.h 6.42 KB
Newer Older
1
//
frey_m's avatar
frey_m committed
2 3
// Class IrregularDomain
//   Defines a common abstract interface for different types of boundaries.
4
//
5
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
6
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
7
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
8
// All rights reserved
9
//
10 11 12 13 14
// 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
15 16 17 18 19 20 21 22 23 24
//
// 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/>.
25
//
gsell's avatar
gsell committed
26 27 28
#ifndef IRREGULAR_DOMAIN_H
#define IRREGULAR_DOMAIN_H

29
#include <map>
gsell's avatar
gsell committed
30
#include <string>
31 32
#include "Algorithms/Vektor.h"
#include "Algorithms/Quaternion.h"
33

gsell's avatar
gsell committed
34 35 36 37 38 39 40
/// enumeration corresponding to different interpolation methods at the boundary
enum {
    CONSTANT,
    LINEAR,
    QUADRATIC
};

41

gsell's avatar
gsell committed
42 43 44 45
class IrregularDomain {

public:

46 47 48 49 50 51 52 53 54 55 56 57 58
    template<typename T>
    struct Stencil {
        T center;   // x,   y,   z
        T west;     // x-1, y,   z
        T east;     // x+1, y,   z
        T north;    // x,   y+1, z
        T south;    // x,   y-1, z
        T front;    // x,   y,   z-1
        T back;     // x,   y,   z+1
    };

    typedef Stencil<int>    StencilIndex_t;
    typedef Stencil<double> StencilValue_t;
59
    typedef Vektor<int, 3>  IntVector_t;
60

61
    IrregularDomain(const IntVector_t& nr,
62 63
                    const Vector_t& hr,
                    const std::string& interpl);
frey_m's avatar
frey_m committed
64 65


frey_m's avatar
frey_m committed
66 67
    /** method to compute the intersection points with the boundary geometry
     * (stored in some appropriate data structure)
gsell's avatar
gsell committed
68 69
     * \param hr updated mesh spacings
     */
70
    virtual void compute(Vector_t hr, NDIndex<3> localId) = 0;
gsell's avatar
gsell committed
71 72 73 74 75

    /// method to calculate the stencil at a boundary points
    /// \param x index of the current element in the matrix
    /// \param y index of the current element in the matrix
    /// \param z index of the current element in the matrix
76 77
    /// \param values of stencil element
    /// \param scaleFactor of stencil values
78
    void getBoundaryStencil(int x, int y, int z,
frey_m's avatar
frey_m committed
79 80
                            StencilValue_t& value,
                            double &scaleFactor) const;
gsell's avatar
gsell committed
81 82

    /// method to calculate the stencil at a boundary points
83
    /// \param id index of the current element in the matrix
84 85 86
    // \param values of stencil element
    /// \param scaleFactor of stencil values
    void getBoundaryStencil(int id, StencilValue_t& value,
87
                            double &scaleFactor) const;
gsell's avatar
gsell committed
88 89 90 91 92

    /// method to calculate the neighbours in the matrix of the current index (x,y,z)
    /// \param x index of the current element in the matrix
    /// \param y index of the current element in the matrix
    /// \param z index of the current element in the matrix
93
    /// \param index stencil indices of an element
94
    void getNeighbours(int x, int y, int z, StencilIndex_t& index) const;
95

96
    void getNeighbours(int idx, StencilIndex_t& index) const;
97

frey_m's avatar
frey_m committed
98
    /// Conversion from a 3D index to (x,y,z)
99
    virtual void getCoord(int idx, int &x, int &y, int &z) const;
100

frey_m's avatar
frey_m committed
101
    /// Conversion from (x,y,z) to index on the 3D grid
102 103 104
    int getIdx(int x, int y, int z) const;

    virtual int getNumXY() const { return nr_m[0] * nr_m[1]; }
gsell's avatar
gsell committed
105 106 107 108 109 110

    /// method that checks if a given point lies inside the boundary
    /// \param x index of the current element in the matrix
    /// \param y index of the current element in the matrix
    /// \param z index of the current element in the matrix
    /// \return boolean indicating if the point lies inside the boundary
111 112 113 114
    virtual bool isInside(int x, int y, int z)  const = 0;

    IntVector_t getNr() const { return nr_m; }
    Vector_t    getHr() const { return hr_m; }
gsell's avatar
gsell committed
115

116 117
    void setNr(IntVector_t nr) { nr_m = nr; }
    void setHr(Vector_t hr)    { hr_m = hr; }
gsell's avatar
gsell committed
118

119 120 121 122
    void setMinMaxZ(double minz, double maxz) {
        zMin_m = minz;
        zMax_m = maxz;
    }
gsell's avatar
gsell committed
123

124 125
    double getMinZ() const { return zMin_m; }
    double getMaxZ() const { return zMax_m; }
126

127 128 129 130 131 132
    double getXRangeMin() const { return min_m(0); }
    double getXRangeMax() const { return max_m(0); }
    double getYRangeMin() const { return min_m(1); }
    double getYRangeMax() const { return max_m(1); }
    double getZRangeMin() const { return min_m(2); }
    double getZRangeMax() const { return max_m(2); }
gsell's avatar
gsell committed
133

134 135 136 137
    void setRangeMin(const Vector_t& min) { min_m = min; }
    void setRangeMax(const Vector_t& max) { max_m = max; }

    bool hasGeometryChanged() const { return hasGeometryChanged_m; }
138

139
    virtual ~IrregularDomain() {};
gsell's avatar
gsell committed
140

141
    virtual void resizeMesh(Vector_t& origin, Vector_t& hr,
frey_m's avatar
frey_m committed
142
                            const Vector_t& /*rmin*/, const Vector_t& /*rmax*/,
143
                            double /*dh*/);
144

gsell's avatar
gsell committed
145 146
protected:

147 148 149
    virtual int indexAccess(int x, int y, int z) const = 0;

    virtual int coordAccess(int idx) const = 0;
frey_m's avatar
frey_m committed
150

151 152
    /// different interpolation methods for boundary points
    virtual void constantInterpolation(int x, int y, int z, StencilValue_t& value,
153
                               double &scaleFactor) const;
154 155

    virtual void linearInterpolation(int x, int y, int z, StencilValue_t& value,
156
                             double &scaleFactor) const;
157 158

    virtual void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
159
                                double &scaleFactor) const;
160 161


gsell's avatar
gsell committed
162 163
    // a irregular domain is always defined on a grid
    /// number of mesh points in each direction
164
    IntVector_t nr_m;
gsell's avatar
gsell committed
165
    /// mesh-spacings in each direction
frey_m's avatar
frey_m committed
166
    Vector_t hr_m;
gsell's avatar
gsell committed
167 168 169 170

    /// min/max of bunch in floor coordinates
    double zMin_m;
    double zMax_m;
171

172 173 174
    Vector_t min_m;
    Vector_t max_m;

175 176 177
    /// flag indicating if geometry has changed for the current time-step
    bool hasGeometryChanged_m;

178 179 180
    /// interpolation type
    int interpolationMethod_m;

181 182 183 184 185 186
    /// mapping (x,y,z) -> idx
    std::map<int, int> idxMap_m;

    /// mapping idx -> (x,y,z)
    std::map<int, int> coordMap_m;

gsell's avatar
gsell committed
187 188
};

189 190
#endif

191

192 193 194 195 196 197 198
// 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: