Grid.h

#ifndef gridripper_amr1d_Grid_h
#define gridripper_amr1d_Grid_h

#include <gridripper/config.h>
#include <gridripper/amr1d/PDE.h>
#include <gridripper/amr1d/FieldWrapper.h>
#include <gridripper/math/XFArray.h>
#include <gridripper/lang/IllegalArgumentException.h>
#include <gridripper/tvalarray.h>
#include <gridripper/tvector.h>
#include <iostream>
#include <map>

namespace gridripper { namespace amr1d {

class GridIterator;
class Integrator;
class Regridder;
class IntervalList;
class Sigma;

using namespace std;

class Grid
{
    friend class GridIterator;
    friend class Regridder;
    friend class Integrator;
    friend class AMRCore;
    friend class PDE::FiniteIntegral;

public:
    enum {LEFT_BOUND_FLAG = 1, RIGHT_BOUND_FLAG = 2};

    enum {IC_CONSTR_MASK = 3, // Constraint mask in options
          IC_CONSTR_NONE = 0, // Constraints are not applied in the points of the
                           // new subgrid.
          IC_CONSTR_INTERPOLATED_ONLY = 1, // Constraints are applied in the
                                        // interpolated points of the new
                                        // subgrid.
          IC_CONSTR_ALL = 2 // Constraints are applied in all points of the new
                         // subgrid.
    };

protected:
    int flags;

    tvalarray<GReal_t> theData;

    tvalarray<bool> refinedDataExists;

private:
    GReal_t deltaX;

    GReal_t baseMinX;

    GReal_t baseMaxX;

    GReal_t locationX;

    int locationI;

    GReal_t currentTime;

    int timeStepCount;

    int lastErrorCheckTime;

    int numComponents;

    int refinementRatio;

    int refinedMaximumI;

    int dataIndexMultiplier;

    int maximumI;

    tvalarray<GReal_t> energyFlow;

    tvalarray<GReal_t> energyFlowLastUpdateTime;

    tvalarray<int> energyFlowLastUpdateLevel;

    tvalarray<GReal_t> energyLoss;

    Regridder* regridder;

    tvalarray<Grid*> subGrids;

    tvalarray<Grid*> newSubSubgrids;

    int levelInTree;

    Integrator* integrator;

    Grid* parentGrid;

    Grid* nextSibling;

    Grid* offspringGrid;

protected:
    Grid(const Grid& g);

    Grid(GReal_t minx, GReal_t maxx, int maxi, int nc,
         Regridder* refproto, int r, GReal_t dx);

    Grid(const Grid& p, GReal_t xshift);

    Grid(Grid* p, int r1, int loci, int maxi,
         int margin, int vmargin, int minLevel);

public:
    virtual ~Grid();

protected:
    void copy(Grid* orig, map<Grid*, Grid*>& gridmap);

    virtual int getInternalIndex(int i) const =0;

    int getRealIndex(int i) const {
        return i;
    }

    virtual int getSubgridIndexFromIndex(const Grid* g, int i) const =0;

    virtual int getIndexFromSubgridIndex(const Grid* g, int i) const =0;

    virtual Grid* createRefinedGrid(IntervalList* iv,
                                    int margin, int vmargin,
                                    int minLevel, int r1) =0;

    void init(Integrator* proto, const Sigma* sigma, int n);

    void initEnergyLoss(PDE& pde);

    void updateEnergyLoss(const PDE& pde, int i, Grid& baseg, int jbase,
                          GReal_t epsilon);

    XFArray* getDensityData(const PDE::DensityQuantity& q,
                            GReal_t xmin, GReal_t xmax) const;

    virtual void blur(tvalarray<unsigned char>& pf, int i) const =0;

    virtual int getFirstIndexForRefinement(const tvalarray<unsigned char>& pf)
            const =0;

    virtual int findIntervalEndForRefinement(
            const tvalarray<unsigned char>& pf, int i) const =0;

    virtual int intervalRefinedSize(const IntervalList* iv, int r) const =0;

    virtual int gridIndexLeftOf(int i, int dmax) const =0;

    virtual int gridIndexRightOf(int i, int dmax) const =0;

    virtual void fixSubgridNearBounds(IntervalList* l, IntervalList* r,
                                      int margin) const =0;

    virtual int calcNewSubSubgridLocation(int ssloc,
                                          int osloc, int nsloc) =0;

public:
    virtual Grid* cloneGrid() const =0;

    GReal_t getBaseMinX() const {
        return baseMinX;
    }

    GReal_t getBaseMaxX() const {
        return baseMaxX;
    }

    GReal_t getLocationX() const {
        return locationX;
    }

    GridIterator getXIterator() const;

    virtual GReal_t getSizeX() const =0;

    virtual int getLeftmostI() const =0;

    virtual int getRightmostI() const =0;

    virtual bool isPeriodic() const =0;

    virtual GReal_t getLocationX(int i) const =0;

    int getIndexFromX(GReal_t x, double maxerr) const
        throw(IllegalArgumentException&);

    int getLocationI() const {
        return locationI;
    }

    GReal_t getDeltaX() const {
        return deltaX;
    }

    void setTime(GReal_t t) {
        currentTime = t;
    }

    GReal_t getTime() const {
        return currentTime;
    }

    void setTimeStepCount(int tsc) {
        timeStepCount = tsc;
    }

    int getTimeStepCount() const {
        return timeStepCount;
    }

    int getLastErrorCheckTime() const {
        return lastErrorCheckTime;
    }

    int getNumComponents() const {
        return numComponents;
    }

    int getRefinementRatio() const {
        return refinementRatio;
    }

    virtual int getMarginLeft(int ir) const =0;

    virtual int getMarginRight(int ir) const =0;

    int getMaxI() const {
        return maximumI;
    }

    virtual int getNumPointsAboveMaxX() const {
        return 0;
    }

    int getLevel() const {
        return levelInTree;
    }

    Grid* getParent() const {
        return parentGrid;
    }

    Grid* getOffspring() const {
        return offspringGrid;
    }

    Grid* getNextSibling() const {
        return nextSibling;
    }

    bool containsLeftBound() const {
        return (flags & LEFT_BOUND_FLAG) != 0;
    }

    bool containsRightBound() const {
        return (flags & RIGHT_BOUND_FLAG) != 0;
    }

    virtual int distanceFromLeft(int i) const =0;

    virtual int distanceFromRight(int i) const =0;

    int getAbsoluteIndex(int i) const {
        return getAbsoluteIndex(i, 1);
    }

    virtual int getAbsoluteIndex(int i, int mul) const =0;

    void getField(int i, FieldWrapper& f) const {
        i = getInternalIndex(i);
        int off = i*dataIndexMultiplier;
        GReal_t* fdata = f.data();
        unsigned nc = f.size();
        for(unsigned j = 0; j < nc; ++j) {
            fdata[j] = theData[off + j];
        }
    }

    void setEnergyLosses(const tvalarray<GReal_t>& e, PDE* pde);

    GReal_t getEnergyFlow(int i) const {
        return energyFlow[i];
    }

    GReal_t getEnergyLoss(int i) const {
        return energyLoss[i];
    }

    void addTo(int i, GReal_t v[]) const {
        i = getInternalIndex(i)*dataIndexMultiplier;
        for(int j = 0; j < numComponents; ++j) {
            v[j] += theData[i + j];
        }
    }

    void addMultipliedTo(int i, GReal_t c, GReal_t v[]) const {
        i = getInternalIndex(i)*dataIndexMultiplier;
        for(int j = 0; j < numComponents; ++j) {
            v[j] += c*theData[i + j];
        }
    }

    void addMultipliedTo(int i, int exclude, GReal_t c, const GReal_t arr[],
                         GReal_t v[]) const {
        int b = 1;
        i = getInternalIndex(i)*dataIndexMultiplier;
        for(int j = 0; j < numComponents; ++j) {
            if((exclude & b) == 0) {
                v[j] += c*arr[j]*theData[i + j];
            }
            b <<= 1;
        }
    }

    void get(int i, GReal_t v[], int off = 0) const {
        i = getInternalIndex(i);
        int off0 = i*dataIndexMultiplier;
        for(int j = 0; j < numComponents; ++j) {
            v[off + j] = theData[off0 + j];
        }
    }

    void set(int i, GReal_t v[]) {
        int j = getInternalIndex(i);
        int off0 = j*dataIndexMultiplier;
        for(int k = 0; k < numComponents; ++k) {
            theData[off0 + k] = v[k];
        }
    }

    GReal_t getComponent(int i, int k) const {
        i = getInternalIndex(i);
        return theData[i*dataIndexMultiplier + k];
    }

    virtual void initRefinedDataExistenceArray() =0;

    int getRefinedMaxI() const {
        return refinedMaximumI;
    }

    virtual bool hasRefinedData(int i) const =0;

    int getRefinedDataExistence(int i, int n, const PDE& pde) const;

    const Grid* getSubgridAt(int i) const {
        return subGrids.size() == 0? NULL : subGrids[i];
    }

    const Grid* getFinestSubgridAt(GReal_t x) const;

    virtual void getRefinedData(int i, GReal_t yout[]) const =0;

    void setZero();

    Regridder* getRegridder() {
        return regridder;
    }

    void setRegridder(Regridder* r);

    bool hasUncoveredFlaggedPoint(int margin);

    virtual Grid* createTemporaryGrid(int dmargin_L, int dmargin_R,
                                      GReal_t xshift) =0;

    void setTemporaryGridShift(const Grid& p, GReal_t xshift) {
        locationX = p.getLocationX() + xshift;
    }

    virtual void setTemporaryGridMargins(const Grid& p, int dmargin_L,
                                         int dmargin_R) =0;

    Integrator* getIntegrator() const {
        return integrator;
    }

    void injectFineToCoarse(PDE* pde);

    void checkErrors(double critical);

    int generateRefinedGrids(int minLevel, int maxLevel,
                             int numVelocity, int bufzoneradius,
                             int thislevel, int gridi);

    void interpolateNewSubgrids(PDE* pde) throw(IllegalArgumentException&);

    void initConstraints(PDE* pde, FieldWrapper& tmp, int opts);

    bool isEnergyLossStored() const {
        return energyLoss.size() != 0;
    }

    string getDataAsString() const;

    string toString() const;

    string sprintHierarchy() const;

    virtual void interpolateMargins(PDE* pde) =0;

protected:
    virtual void fillRefinedDataArray() =0;

    void initConstraints(PDE* pde, int left, int right,
                         FieldWrapper& tmp, int opts);

private:
    string sprintHierarchy(int tabs, const Grid* gg) const;

    static void fillGridArray(tvalarray<Grid*>& d, int left, int right,
                              Grid* g);

    void initSubsub(Grid* p, int loci, int maxi, int minLevel);
};

} } // namespace gridripper::amr1d

inline std::ostream& operator <<(std::ostream& os,
                                 const gridripper::amr1d::Grid& g)
{
    return os << g.toString();
}

#endif /* gridripper_amr1d_Grid_h */