BoxCornerDomain.h 5.15 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 31 32 33 34 35 36 37 38

#include <map>
#include <string>
#include <cmath>
#include <iostream>  // Neeeded for stream I/O
#include <fstream>   // Needed for file I/O
#include "IrregularDomain.h"


/*

39
    A and B are the half apperture of the box
gsell's avatar
gsell committed
40 41


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

58
            Length_m
gsell's avatar
gsell committed
59 60 61

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

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

67 68 69 70 71

A  = max_m(0)
B  = max_m(1)
L1 = min_m(2)
L2 = max_m(2) - min_m(2)
gsell's avatar
gsell committed
72 73 74 75 76
*/

class BoxCornerDomain : public IrregularDomain {

public:
77 78
    using IrregularDomain::StencilIndex_t;
    using IrregularDomain::StencilValue_t;
gsell's avatar
gsell committed
79

80 81 82 83 84 85 86 87 88 89
    /**
     * \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
     */
    BoxCornerDomain(double A, double B, double C, double length,
                    double L1, double L2, Vector_t nr, Vector_t hr,
gsell's avatar
gsell committed
90 91 92 93 94
                    std::string interpl);
    ~BoxCornerDomain();

    /// as a function of z, determine the hight (B) of the geometry
    inline double getB(double z) {
95 96
      if((z < min_m(2)) || (z > max_m(2)))
            return max_m(1);
gsell's avatar
gsell committed
97
        else
98
            return max_m(1) - C_m;
gsell's avatar
gsell committed
99 100 101 102
    }

    /// queries if a given (x,y,z) coordinate lies inside the domain
    inline bool isInside(int x, int y, int z) {
frey_m's avatar
frey_m committed
103 104 105 106
        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]);
        return (xx < getXRangeMax() && yy < b && z != 0 && z != nr_m[2] - 1);
gsell's avatar
gsell committed
107 108
    }

109
    void compute(Vector_t hr, NDIndex<3> localId);
gsell's avatar
gsell committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130

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;

131 132
    /// length of the structure
    double length_m;
gsell's avatar
gsell committed
133

134 135
    /// height of the corner
    double C_m;
gsell's avatar
gsell committed
136 137 138



139
    inline double getXIntersection(double cx, int /*z*/) {
140
        return (cx < 0) ? min_m(0) : max_m(0);
gsell's avatar
gsell committed
141 142 143
    }

    inline double getYIntersection(double cy, int z) {
frey_m's avatar
frey_m committed
144
        return (cy < 0) ? min_m(1) : getB(z * hr_m[2]);
gsell's avatar
gsell committed
145 146 147 148
    }

    /// conversion from (x,y,z) to index in xyz plane
    inline int toCoordIdx(int x, int y, int z) {
frey_m's avatar
frey_m committed
149
        return (z * nr_m[1] + y) * nr_m[0] + x;
gsell's avatar
gsell committed
150 151 152
    }

    /// conversion from (x,y,z) to index on the 3D grid
frey_m's avatar
frey_m committed
153 154
    int indexAccess(int x, int y, int z) {
        return idxMap_m[toCoordIdx(x, y, z)];
gsell's avatar
gsell committed
155 156 157
    }

    /// different interpolation methods for boundary points
158
    void constantInterpolation(int x, int y, int z, StencilValue_t& value,
159
                               double &scaleFactor) override;
160 161

    void linearInterpolation(int x, int y, int z, StencilValue_t& value,
162
                             double &scaleFactor) override;
163 164

    void quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
165
                                double &scaleFactor) override;
gsell's avatar
gsell committed
166 167 168

};

169 170 171 172 173 174 175 176 177
#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: