BoxCornerDomain.h 4.82 KB
Newer Older
1 2 3 4
//
// Class BoxCornerDomain
//   :FIXME: add brief description
//
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
#ifndef BOXCORNER_DOMAIN_H
27
#define BOXCORNER_DOMAIN_H
gsell's avatar
gsell committed
28 29 30

#include <map>
#include <string>
31
#include <utility>
gsell's avatar
gsell committed
32

33
#include "Solvers/RegularDomain.h"
gsell's avatar
gsell committed
34 35 36

/*

37
    A and B are the half aperture of the box
gsell's avatar
gsell committed
38

39
                                     / (A,B)
gsell's avatar
gsell committed
40 41 42
                                    /
                                   /
                                  /
43 44 45 46
    L1                         /
------------      --------------+ (-A,B)
           | L2 |             |
        C|      |             |
gsell's avatar
gsell committed
47
           |------|             |      /
kraus's avatar
kraus committed
48
         .....                  |     /
gsell's avatar
gsell committed
49 50 51
(0,0)---.......-----------------+    /
         .....                  |   /
   z                            |  /
kraus's avatar
kraus committed
52
   |                            | /
53
--------------------------------+/ (-A,-B)
gsell's avatar
gsell committed
54

55
            Length_m
gsell's avatar
gsell committed
56 57 58

Test in which of the 3 parts of the geometry we are in.

59 60
    if((z < L1) || (z > (L1 + L2)))
        b = B;
gsell's avatar
gsell committed
61
    else
62
        b = B-C;
gsell's avatar
gsell committed
63

64

frey_m's avatar
frey_m committed
65 66 67 68
A  = getXRangeMax()
B  = getYRangeMax()
L1 = getZRangeMin()
L2 = getZRangeMax() - getZRangeMin
gsell's avatar
gsell committed
69 70
*/

71
class BoxCornerDomain : public RegularDomain {
gsell's avatar
gsell committed
72 73

public:
74 75 76 77 78 79 80 81
    /**
     * \param A depth of the box
     * \param B maximal height of the box
     * \param C height of the corner
     * \param length of the structure
     * \param L1 length of the first part of the structure
     * \param L2 length of the corner
     */
82
    BoxCornerDomain(double A, double B, double C,
83
                    double L1, double L2, IntVector_t nr, Vector_t hr,
gsell's avatar
gsell committed
84 85 86
                    std::string interpl);
    ~BoxCornerDomain();

87
    /// as a function of z, determine the height (B) of the geometry
88
    inline double getB(double z) const {
frey_m's avatar
frey_m committed
89 90
      if((z < getZRangeMin()) || (z > getZRangeMax()))
            return getYRangeMax();
gsell's avatar
gsell committed
91
        else
frey_m's avatar
frey_m committed
92
            return getYRangeMax() - C_m;
gsell's avatar
gsell committed
93 94 95
    }

    /// queries if a given (x,y,z) coordinate lies inside the domain
96
    inline bool isInside(int x, int y, int z) const {
frey_m's avatar
frey_m committed
97 98 99
        const double xx = (x - (nr_m[0] - 1) / 2.0) * hr_m[0];
        const double yy = (y - (nr_m[1] - 1) / 2.0) * hr_m[1];
        const double b = getB(z * hr_m[2]);
100
        return (xx < getXRangeMax() && yy < b && z >= 0 && z < nr_m[2]);
gsell's avatar
gsell committed
101 102
    }

103
    void compute(Vector_t hr, NDIndex<3> localId);
gsell's avatar
gsell committed
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

private:

    //XXX: since the Y coorindate is dependent on the Z value we need (int,
    //int) -> intersection. To simplify things (for now) we use the same
    //structure for X...
    /// Map from a ([(x or y], z) to a list of intersection values with
    /// boundary.
    typedef std::multimap< std::pair<int, int>, double > BoxCornerPointList;

    /// all intersection points with grid lines in X direction
    BoxCornerPointList IntersectXDir;

    /// all intersection points with grid lines in Y direction
    BoxCornerPointList IntersectYDir;

    /// because the geometry can change in the y direction
    double actBMin_m;

    double actBMax_m;

125 126
    /// height of the corner
    double C_m;
gsell's avatar
gsell committed
127

128
    inline double getXIntersection(double cx, int /*z*/) const {
frey_m's avatar
frey_m committed
129
        return (cx < 0) ? getXRangeMin() : getXRangeMax();
gsell's avatar
gsell committed
130 131
    }

132
    inline double getYIntersection(double cy, int z) const {
frey_m's avatar
frey_m committed
133
        return (cy < 0) ? getYRangeMin() : getB(z * hr_m[2]);
gsell's avatar
gsell committed
134 135 136
    }

    /// conversion from (x,y,z) to index in xyz plane
137
    inline int toCoordIdx(int x, int y, int z) const {
frey_m's avatar
frey_m committed
138
        return (z * nr_m[1] + y) * nr_m[0] + x;
gsell's avatar
gsell committed
139 140 141
    }

    /// conversion from (x,y,z) to index on the 3D grid
142 143 144 145 146 147
    int indexAccess(int x, int y, int z) const {
        return idxMap_m.at(toCoordIdx(x, y, z));
    }

    int coordAccess(int idx) const {
        return coordMap_m.at(idx);
gsell's avatar
gsell committed
148 149 150
    }

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

    void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
155
                                double &scaleFactor) const override;
gsell's avatar
gsell committed
156 157 158

};

159
#endif