RectangularDomain.h 3.36 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4
//
// Class RectangularDomain
//   :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 9
// All rights reserved
//
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 25
//
// 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
26 27 28 29
#ifndef RECTANGULAR_DOMAIN_H
#define RECTANGULAR_DOMAIN_H

#include <vector>
30
#include <string>
gsell's avatar
gsell committed
31 32 33 34 35
#include "IrregularDomain.h"

class RectangularDomain : public IrregularDomain {

public:
36 37
    using IrregularDomain::StencilIndex_t;
    using IrregularDomain::StencilValue_t;
gsell's avatar
gsell committed
38 39 40 41 42 43

    /// constructor
    RectangularDomain(Vector_t nr, Vector_t hr);
    /// constructor
    RectangularDomain(double a, double b, Vector_t nr, Vector_t hr);

44
    /// calculates intersection with the beam pipe
45
    void compute(Vector_t hr, NDIndex<3> /*localId*/);
46

gsell's avatar
gsell committed
47 48 49
    /// returns number of nodes in xy plane (here independent of z coordinate)
    int getNumXY(int z);
    /// returns discretization at (x,y,z)
50 51 52
    void getBoundaryStencil(int x, int y, int z, StencilValue_t& value,
                            double &scaleFactor);

gsell's avatar
gsell committed
53
    /// returns index of neighbours at (x,y,z)
54
    using IrregularDomain::getNeighbours;
gsell's avatar
gsell committed
55
    /// queries if a given (x,y,z) coordinate lies inside the domain
56
    inline bool isInside(int x, int y, int /*z*/) {
gsell's avatar
gsell committed
57 58 59 60 61 62 63 64 65 66 67 68
        double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
        double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
        return (xx <= a_m && yy < b_m);
    }

    void setB_m(double b) {b_m = b;}
    void setA_m(double a) {a_m = a;}

    double getXRangeMin() { return -a_m; }
    double getXRangeMax() { return a_m; }
    double getYRangeMin() { return -b_m; }
    double getYRangeMax() { return b_m; }
69 70 71
    double getZRangeMin() { return getMinZ(); }
    double getZRangeMax() { return getMaxZ(); }

gsell's avatar
gsell committed
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101

    int getStartIdx() {return 0;}

private:

    /// longer side a of the rectangles
    double a_m;
    /// shorter side b of the rectangles
    double b_m;
    /// number of nodes in the xy plane (for this case: independent of the z coordinate)
    int nxy_m;

    /// conversion from (x,y,z) to index on the 3D grid
    inline int getIdx(int x, int y, int z) {
        if(isInside(x, y, z) && x >= 0 && y >= 0 && z >= 0)
            return y * nr[0] + x + z * nxy_m;
        else
            return -1;
    }
    /// 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;
        int inr = nr[0];
        x = ixy % inr;
        y = (ixy - x) / nr[0];
        z = (idx - ixy) / nxy_m;
    }

};

102 103 104 105 106 107 108 109 110 111
#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: