EllipticDomain.h 6.36 KB
Newer Older
frey_m's avatar
frey_m committed
1 2
//
// Class EllipticDomain
3 4 5
//   This class provides an elliptic beam pipe. The mesh adapts to the bunch size
//   in the longitudinal direction. At the intersection of the mesh with the
//   beam pipe, three stencil interpolation methods are available.
frey_m's avatar
frey_m committed
6
//
7
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
8
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
9
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
10 11
// All rights reserved
//
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 27
//
// 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/>.
//
gsell's avatar
gsell committed
28 29 30 31 32 33 34 35
#ifndef ELLIPTICAL_DOMAIN_H
#define ELLIPTICAL_DOMAIN_H

#include <vector>
#include <map>
#include <string>
#include <cmath>
#include "IrregularDomain.h"
36
#include "Structure/BoundaryGeometry.h"
37
#include "Utilities/OpalException.h"
gsell's avatar
gsell committed
38 39 40 41 42 43

class EllipticDomain : public IrregularDomain {

public:

    EllipticDomain(Vector_t nr, Vector_t hr);
44 45 46 47 48 49 50

    EllipticDomain(double semimajor, double semiminor, Vector_t nr,
                   Vector_t hr, std::string interpl);

    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr,
                   Vector_t hr, std::string interpl);

51 52
    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);

gsell's avatar
gsell committed
53 54 55
    ~EllipticDomain();

    /// returns discretization at (x,y,z)
56 57 58 59 60
    void getBoundaryStencil(int x, int y, int z,
                            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
62 63 64 65
    void getBoundaryStencil(int idx, double &W, double &E,
                            double &S, double &N, double &F,
                            double &B, double &C, double &scaleFactor);

gsell's avatar
gsell committed
66
    /// returns index of neighbours at (x,y,z)
67 68 69
    void getNeighbours(int x, int y, int z, int &W, int &E,
                       int &S, int &N, int &F, int &B);

gsell's avatar
gsell committed
70 71
    /// returns index of neighbours at 3D index
    void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
72

gsell's avatar
gsell committed
73 74
    /// returns type of boundary condition
    std::string getType() {return "Elliptic";}
75

gsell's avatar
gsell committed
76 77
    /// queries if a given (x,y,z) coordinate lies inside the domain
    inline bool isInside(int x, int y, int z) {
78 79
        double xx = - semiMajor_m + hr[0] * (x + 0.5);
        double yy = - semiMinor_m + hr[1] * (y + 0.5);
80 81 82 83 84

        bool isInsideEllipse = (xx * xx / (semiMajor_m * semiMajor_m) +
                                yy * yy / (semiMinor_m * semiMinor_m) < 1);

        return (isInsideEllipse && z > 0 && z < nr[2] - 1);
gsell's avatar
gsell committed
85 86
    }

87
    int getNumXY(int /*z*/) { return nxy_m; }
gsell's avatar
gsell committed
88

89 90

    /// calculates intersection
91
    void compute(Vector_t hr, NDIndex<3> localId);
gsell's avatar
gsell committed
92

93 94 95 96
    double getXRangeMin() { return -semiMajor_m; }
    double getXRangeMax() { return semiMajor_m;  }
    double getYRangeMin() { return -semiMinor_m; }
    double getYRangeMax() { return semiMinor_m;  }
97 98 99
    double getZRangeMin() { return zMin_m; }
    double getZRangeMax() { return zMax_m; }

gsell's avatar
gsell committed
100 101
    bool hasGeometryChanged() { return hasGeometryChanged_m; }

102

103 104
    void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
                    const Vector_t& rmax, double dh);
105

gsell's avatar
gsell committed
106 107 108 109
private:

    /// Map from a single coordinate (x or y) to a list of intersection values with
    /// boundary.
110
    typedef std::multimap<int, double> EllipticPointList_t;
gsell's avatar
gsell committed
111 112

    /// all intersection points with grid lines in X direction
113
    EllipticPointList_t intersectXDir_m;
gsell's avatar
gsell committed
114 115

    /// all intersection points with grid lines in Y direction
116
    EllipticPointList_t intersectYDir_m;
gsell's avatar
gsell committed
117 118

    /// mapping (x,y,z) -> idx
119
    std::map<int, int> idxMap_m;
gsell's avatar
gsell committed
120 121

    /// mapping idx -> (x,y,z)
122
    std::map<int, int> coordMap_m;
gsell's avatar
gsell committed
123 124

    /// semi-major of the ellipse
125 126
    double semiMajor_m;

gsell's avatar
gsell committed
127
    /// semi-minor of the ellipse
128 129
    double semiMinor_m;

gsell's avatar
gsell committed
130 131
    /// number of nodes in the xy plane (for this case: independent of the z coordinate)
    int nxy_m;
132

gsell's avatar
gsell committed
133
    /// interpolation type
134 135
    int interpolationMethod_m;

gsell's avatar
gsell committed
136 137 138 139 140
    /// flag indicating if geometry has changed for the current time-step
    bool hasGeometryChanged_m;

    /// conversion from (x,y) to index in xy plane
    inline int toCoordIdx(int x, int y) { return y * nr[0] + x; }
141

gsell's avatar
gsell committed
142 143
    /// conversion from (x,y,z) to index on the 3D grid
    inline int getIdx(int x, int y, int z) {
144 145
        if (isInside(x, y, z))
            return idxMap_m[toCoordIdx(x, y)] + (z - 1) * nxy_m;
gsell's avatar
gsell committed
146 147 148
        else
            return -1;
    }
149

gsell's avatar
gsell committed
150 151 152
    /// conversion from a 3D index to (x,y,z)
    inline void getCoord(int idx, int &x, int &y, int &z) {
        int ixy = idx % nxy_m;
153
        int xy = coordMap_m[ixy];
gsell's avatar
gsell committed
154 155 156 157 158 159 160
        int inr = (int)nr[0];
        x = xy % inr;
        y = (xy - x) / nr[0];
        z = (idx - ixy) / nxy_m + 1;
    }

    /// different interpolation methods for boundary points
161 162 163 164 165 166 167 168 169 170 171 172 173 174
    void constantInterpolation(int x, int y, int z,
                               double &W, double &E, double &S,
                               double &N, double &F, double &B,
                               double &C, double &scaleFactor);

    void linearInterpolation(int x, int y, int z, double &W,
                             double &E, double &S, double &N,
                             double &F, double &B, double &C,
                             double &scaleFactor);

    void quadraticInterpolation(int x, int y, int z, double &W,
                                double &E, double &S, double &N,
                                double &F, double &B, double &C,
                                double &scaleFactor);
175

frey_m's avatar
frey_m committed
176
    /// function to handle the open boundary condition in longitudinal direction
177
    void robinBoundaryStencil(int z, double &F, double &B, double &C);
gsell's avatar
gsell committed
178 179 180
};

#endif //#ifdef ELLIPTICAL_DOMAIN_H
181 182 183 184 185 186 187 188

// 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: