ODERelax.h

#ifndef gridripper_odesolver_relax_ODERelax_h
#define gridripper_odesolver_relax_ODERelax_h

#include <gridripper/config.h>
#include <gridripper/lang/Exception.h>
#include <gridripper/util/GrrProperties.h>
#include <gridripper/tvalarray.h>
#include <string>

namespace gridripper { namespace odesolver {

class ODE;

namespace relax {

using namespace gridripper;
using namespace gridripper::util;
using namespace std;

class ODERelax
{
public:
    class Method
    {
    public:
        virtual ~Method() { }

        virtual void set(tvalarray<GReal_t>& fieldValues,
                tvalarray<GReal_t>& matrix, tvalarray<GReal_t>& rhsvec,
                const tvalarray<GReal_t>& deltaField) =0;

        virtual int getNumSubsteps() const =0;

        int getRowLength() {
            int substeps = getNumSubsteps();
            return 2*substeps + 1;
        }
    };

    class RelaxationException: public Exception
    {
    public:
        RelaxationException(const RelaxationException& other):
            Exception(other) { }

        RelaxationException(const string& msg, const string& where):
            Exception(msg, where) { }

        RelaxationException(const Exception& orig, const string& where):
            Exception(orig, where) { }

        RelaxationException(const Exception& orig, const string& msg,
                             const string& where):
            Exception(orig, msg, where) { }
    };

private:
    string dataFormat;
    ODE* theODE;
    int maxi;
    int order;
    tvalarray<GReal_t> matrix;
    tvalarray<GReal_t> rhsvec;
    tvalarray<GReal_t> fieldValues;
    tvalarray<GReal_t> deltaFieldValues;
    Method* method;
    int numIterations;
    double relativeError;

public:
    ODERelax(const GrrProperties& props) throw(RelaxationException&);

    ODERelax(const GrrProperties& props, int maxi)
        throw(RelaxationException&);

    ODERelax(ODE* ode, tvalarray<GReal_t>& data, int order)
            throw(RelaxationException&);

    ~ODERelax();

    ODE& getODE() const {
        return *theODE;
    }

    int getMaxI() {
        return maxi;
    }

    string getDataFormat() {
        return dataFormat;
    }

    double getRelativeError() {
        return relativeError;
    }

    int getNumIterations() {
        return numIterations;
    }

    void setNumIterations(int n) {
        numIterations = n;
    }

    void solve(tvalarray<GReal_t>& sol) throw(RelaxationException&);

    double nextIteration() throw(RelaxationException&);

private:
    void init(ODE* ode, tvalarray<GReal_t>& data, int order)
            throw(RelaxationException&);

    void diagonalize();

    void unite() throw(RelaxationException&);

    string getOrderedComponentName(int j) const;
};

} } } // namespace gridripper::odesolver::relax

#endif /* gridripper_odesolver_relax_ODERelax_h */