//
// Class RectangularDomain
//   :FIXME: add brief description
//
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// 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)
//
// 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/>.
//
#ifndef RECTANGULAR_DOMAIN_H
#define RECTANGULAR_DOMAIN_H

#include <vector>
#include <string>
#include "IrregularDomain.h"

class RectangularDomain : public IrregularDomain {

public:
    using IrregularDomain::StencilIndex_t;
    using IrregularDomain::StencilValue_t;

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

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

    /// 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, StencilValue_t& value,
                            double &scaleFactor);

    /// returns index of neighbours at (x,y,z)
    using IrregularDomain::getNeighbours;
    /// queries if a given (x,y,z) coordinate lies inside the domain
    inline bool isInside(int x, int y, int /*z*/) {
        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; }
    double getZRangeMin() { return getMinZ(); }
    double getZRangeMax() { return getMaxZ(); }


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

};

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