IrregularDomain.h 7.16 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 29 30
#ifndef IRREGULAR_DOMAIN_H
#define IRREGULAR_DOMAIN_H

#include <vector>
#include <string>
31
#include "Algorithms/PBunchDefs.h"
32 33
#include "Algorithms/Quaternion.h"

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:

frey_m's avatar
frey_m committed
46 47
    /** method to compute the intersection points with the boundary geometry
     * (stored in some appropriate data structure)
gsell's avatar
gsell committed
48 49
     * \param hr updated mesh spacings
     */
50
    virtual void compute(Vector_t hr, NDIndex<3> localId) = 0;
gsell's avatar
gsell committed
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    /** method to get the number of gridpoints in a given z plane
     * \param z coordinate of the z plane
     * \return int number of grid nodes in the given z plane
     */
    virtual int getNumXY(int z) = 0;

    /// 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
    /// \param W stencil value of the element in the west of idx: (x-1)
    /// \param E stencil value of the element in the east of idx: (x+1)
    /// \param S stencil value of the element in the south of idx: (y-1)
    /// \param N stencil value of the element in the north of idx: (y+1)
    /// \param F stencil value of the element in front of idx: (z-1)
    /// \param B stencil value of the element in the back of idx: (z+1)
    /// \param C stencil value of the element in the center
frey_m's avatar
frey_m committed
68 69 70 71
    virtual void getBoundaryStencil(int x, int y, int z,
                                    double &W, double &E, double &S,
                                    double &N, double &F, double &B,
                                    double &C, double &scaleFactor) = 0;
gsell's avatar
gsell committed
72 73 74 75 76 77 78 79 80 81

    /// method to calculate the stencil at a boundary points
    /// \param idx index of the current element in the matrix
    /// \param W stencil value of the element in the west of idx: (x-1)
    /// \param E stencil value of the element in the east of idx: (x+1)
    /// \param S stencil value of the element in the south of idx: (y-1)
    /// \param N stencil value of the element in the north of idx: (y+1)
    /// \param F stencil value of the element in front of idx: (z-1)
    /// \param B stencil value of the element in the back of idx: (z+1)
    /// \param C stencil value of the element in the center
frey_m's avatar
frey_m committed
82 83 84
    virtual void getBoundaryStencil(int idx, double &W, double &E, double &S,
                                    double &N, double &F, double &B, double &C,
                                    double &scaleFactor) = 0;
gsell's avatar
gsell committed
85 86 87 88 89 90 91 92 93 94 95

    /// 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
    /// \param W stencil index of the element in the west of idx: (x-1)
    /// \param E stencil index of the element in the east of idx: (x+1)
    /// \param S stencil index of the element in the south of idx: (y-1)
    /// \param N stencil index of the element in the north of idx: (y+1)
    /// \param F stencil index of the element in front of idx: (z-1)
    /// \param B stencil index of the element in the back of idx: (z+1)
frey_m's avatar
frey_m committed
96 97 98 99 100
    virtual void getNeighbours(int x, int y, int z, int &W, int &E, int &S,
                               int &N, int &F, int &B) = 0;

    virtual void getNeighbours(int idx, int &W, int &E, int &S, int &N,
                               int &F, int &B) = 0;
gsell's avatar
gsell committed
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

    /// method that identifies a specialized boundary geometry
    /// \return std::string containing a description of the boundary geometry used
    virtual std::string getType() = 0;

    /// 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
    virtual bool isInside(int x, int y, int z) = 0;

    Vector_t getNr() { return nr; }
    Vector_t getHr() { return hr; }
    void setNr(Vector_t nri) { nr = nri; }
    void setHr(Vector_t hri) { hr = hri; }

    void setMinMaxZ(double minz, double maxz) { zMin_m=minz; zMax_m=maxz; }
    double getMinZ() { return zMin_m; }
    double getMaxZ() { return zMax_m; }

122 123
    void setGlobalMeanR(Vector_t rmean) { rMean_m = rmean;}
    Vector_t getGlobalMeanR() { return rMean_m; }
124

125 126
    void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion){
        globalToLocalQuaternion_m = globalToLocalQuaternion;}
frey_m's avatar
frey_m committed
127

128
    Quaternion_t getGlobalToLocalQuaternion() { return globalToLocalQuaternion_m;}
129

gsell's avatar
gsell committed
130 131 132 133
    virtual double getXRangeMin() = 0;
    virtual double getXRangeMax() = 0;
    virtual double getYRangeMin() = 0;
    virtual double getYRangeMax() = 0;
134 135
    virtual double getZRangeMin() = 0;
    virtual double getZRangeMax() = 0;
gsell's avatar
gsell committed
136 137 138

    virtual int getIdx(int x, int y, int z) = 0;
    virtual bool hasGeometryChanged() = 0;
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 143 144
                            const Vector_t& /*rmin*/, const Vector_t& /*rmax*/,
                            double /*dh*/)
    {
145 146 147 148 149 150 151 152 153 154 155
        double xmin = getXRangeMin();
        double xmax = getXRangeMax();
        double ymin = getYRangeMin();
        double ymax = getYRangeMax();
        double zmin = getZRangeMin();
        double zmax = getZRangeMax();

        origin = Vector_t(xmin, ymin , zmin);
        Vector_t mymax = Vector_t(xmax, ymax , zmax);

        for (int i = 0; i < 3; i++)
frey_m's avatar
frey_m committed
156
            hr[i] = (mymax[i] - origin[i]) / nr[i];
157 158
    };

gsell's avatar
gsell committed
159 160 161 162 163 164 165 166 167 168 169
protected:

    // a irregular domain is always defined on a grid
    /// number of mesh points in each direction
    Vector_t nr;
    /// mesh-spacings in each direction
    Vector_t hr;

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

    /// mean position of bunch (m)
    Vector_t rMean_m;
173
    Quaternion_t globalToLocalQuaternion_m;
gsell's avatar
gsell committed
174 175
};

176 177 178 179 180 181 182 183 184
#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: