RectangularDomain.h 3.82 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 36 37 38 39 40 41
#include "IrregularDomain.h"

class RectangularDomain : public IrregularDomain {

public:

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

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

gsell's avatar
gsell committed
45 46 47 48 49 50 51
    /// returns number of nodes in xy plane (here independent of z coordinate)
    int getNumXY(int z);
    /// returns discretization at (x,y,z)
    void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
    /// returns discretization at 3D index
    void getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
    /// returns index of neighbours at (x,y,z)
52
    using IrregularDomain::getNeighbours;
gsell's avatar
gsell committed
53 54 55 56
    void getNeighbours(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B);
    /// returns index of neighbours at 3D index
    void getNeighbours(int idx, double &W, double &E, double &S, double &N, double &F, double &B);
    /// returns type of boundary condition
57
    std::string getType() {return "Rectangular";}
gsell's avatar
gsell committed
58
    /// queries if a given (x,y,z) coordinate lies inside the domain
59
    inline bool isInside(int x, int y, int /*z*/) {
gsell's avatar
gsell committed
60 61 62 63 64 65 66 67 68 69 70 71
        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; }
72 73 74
    double getZRangeMin() { return getMinZ(); }
    double getZRangeMax() { return getMaxZ(); }

gsell's avatar
gsell committed
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 102 103 104

    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;
    }

};

105 106 107 108 109 110 111 112 113 114
#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: