ODE.h

#ifndef gridripper_amr1d_odesolver_ODE_h
#define gridripper_amr1d_odesolver_ODE_h

#include <gridripper/amr1d/FieldArray.h>
#include <gridripper/amr1d/InitCond.h>
#include <gridripper/amr1d/PDE.h>
#include <gridripper/tvalarray.h>
#include <string>

namespace gridripper { namespace odesolver {

using namespace gridripper::amr1d;
using namespace std;

class ODE
{
private:
    int numComponents;

    const string* componentNames;

    const InitCond* thePDEInitCond;

    const PDE* thePDE;

protected:
    ODE(int nc, const string* compnames) {
        numComponents = nc;
        componentNames = compnames;
    }

public: 
    virtual ~ODE() { }

    virtual bool isNumericalConstraintSettingPassedToPDE() const {
        return true;
    }

    virtual bool isConstrainedInInit() const {
        return true;
    }

    virtual void constrain(tvalarray<GReal_t>& mesh) {
    }

    virtual void eval(const GReal_t* f, int offset, GReal_t x, GReal_t* df) =0;

    virtual int getMatrixFieldOrder(int i) const {
        return i;
    }

    virtual bool hasExactSolution() const {
        return false;
    }

    virtual void getExactSolution(tvalarray<GReal_t>& F, GReal_t x) const {
    }

    virtual GReal_t getMinX() const =0;

    virtual GReal_t getMaxX() const =0;

    virtual bool isValid(const GReal_t* F, int nc, GReal_t x) const {
        return true;
    }

    virtual bool dRHS(const GReal_t* f, int offset, GReal_t x, GReal_t* drhs);

    virtual tvalarray<GReal_t> getInitialValuesForRelaxation(int nx) const;

    virtual tvalarray<GReal_t> getInitialX() const;

    virtual FieldArray getInitialValuesForShooting() const;

    virtual GReal_t averageShootingResults(int np, int ip, GReal_t x, int ic,
                                           GReal_t f0, GReal_t f1) const {
        return 0.5*(f0 + f1);
    }

    virtual void getPhysicalField(GReal_t x, GReal_t* f) const {
    }

    const PDE& getPDE() const {
        return *thePDE;
    }

    const InitCond& getPDEInitCond() const {
        return *thePDEInitCond;
    }

    void setPDE(const PDE* pde, const InitCond* init) {
        thePDE = pde;
        thePDEInitCond = init;
    }

    virtual void initPDE(GReal_t x, const GReal_t* fode, FieldWrapper& fpde);

    virtual void initPDE(GReal_t x, const GReal_t* fode,
                         FieldWrapper& fpde, const GReal_t* data, int ndata) {
        initPDE(x, fode, fpde);
    }

    int getNumComponents() const {
        return numComponents;
    }

    string getComponentName(int k) const {
        return componentNames[k];
    }
};

} } // namespace gridripper::odesolver

#endif /* gridripper_amr1d_odesolver_ODE_h */