// // Class BoxCornerDomain // :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 . // #ifndef BOXCORNER_DOMAIN_H #define BOXCORNER_DOMAIN_H #include #include #include #include "Solvers/RegularDomain.h" /* A and B are the half aperture of the box / (A,B) / / / L1 / ------------ --------------+ (-A,B) | L2 | | C| | | |------| | / ..... | / (0,0)---.......-----------------+ / ..... | / z | / | | / --------------------------------+/ (-A,-B) Length_m Test in which of the 3 parts of the geometry we are in. if((z < L1) || (z > (L1 + L2))) b = B; else b = B-C; A = getXRangeMax() B = getYRangeMax() L1 = getZRangeMin() L2 = getZRangeMax() - getZRangeMin */ class BoxCornerDomain : public RegularDomain { public: /** * \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 L1, double L2, IntVector_t nr, Vector_t hr, std::string interpl); ~BoxCornerDomain(); /// as a function of z, determine the height (B) of the geometry inline double getB(double z) const { if((z < getZRangeMin()) || (z > getZRangeMax())) return getYRangeMax(); else return getYRangeMax() - C_m; } /// queries if a given (x,y,z) coordinate lies inside the domain inline bool isInside(int x, int y, int z) const { 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]); } void compute(Vector_t hr, NDIndex<3> localId); 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, 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; /// height of the corner double C_m; inline double getXIntersection(double cx, int /*z*/) const { return (cx < 0) ? getXRangeMin() : getXRangeMax(); } inline double getYIntersection(double cy, int z) const { return (cy < 0) ? getYRangeMin() : getB(z * hr_m[2]); } /// conversion from (x,y,z) to index in xyz plane inline int toCoordIdx(int x, int y, int z) const { return (z * nr_m[1] + y) * nr_m[0] + x; } /// conversion from (x,y,z) to index on the 3D grid 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); } /// different interpolation methods for boundary points void linearInterpolation(int x, int y, int z, StencilValue_t& value, double &scaleFactor) const override; void quadraticInterpolation(int x, int y, int z, StencilValue_t& value, double &scaleFactor) const override; }; #endif