EllipticDomain.h

Go to the documentation of this file.

//
// Class EllipticDomain
//   This class provides an elliptic beam pipe. The mesh adapts to the bunch size
//   in the longitudinal direction. At the intersection of the mesh with the
//   beam pipe, three stencil interpolation methods are available.
//
// 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 ELLIPTICAL_DOMAIN_H
#define ELLIPTICAL_DOMAIN_H
 
#include <vector>
#include <map>
#include <string>
#include <cmath>
#include "Solvers/RegularDomain.h"
#include "Structure/BoundaryGeometry.h"
#include "Utilities/OpalException.h"
 
class EllipticDomain : public RegularDomain {
 
public:
    EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr,
                   Vector_t hr, std::string interpl);
 
    ~EllipticDomain();
 
    bool isInside(int x, int y, int z) const override {
        double xx = getXRangeMin() + hr_m[0] * (x + 0.5);
        double yy = getYRangeMin() + hr_m[1] * (y + 0.5);
 
        bool isInsideEllipse = (xx * xx / (getXRangeMax() * getXRangeMax()) +
                                yy * yy / (getYRangeMax() * getYRangeMax()) < 1);
 
        return (isInsideEllipse && z >= 0 && z < nr_m[2]);
    }
 
    void compute(Vector_t hr, NDIndex<3> localId) override;
 
private:
 
    typedef std::multimap<int, double> EllipticPointList_t;
 
    EllipticPointList_t intersectXDir_m;
 
    EllipticPointList_t intersectYDir_m;
 
    int toCoordIdx(int x, int y) const { return y * nr_m[0] + x; }
 
    int indexAccess(int x, int y, int z) const override {
        return idxMap_m.at(toCoordIdx(x, y)) + z * getNumXY();
    }
 
    int coordAccess(int idx) const override {
        int ixy = idx % getNumXY();
        return coordMap_m.at(ixy);
    }
 
    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 //#ifdef ELLIPTICAL_DOMAIN_H