// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file ProblemVec.h */

#ifndef AMDIS_PROBLEMVEC_H
#define AMDIS_PROBLEMVEC_H

#include "ProblemStatBase.h"
#include "Parameters.h"
#include "Boundary.h"
#include "MatrixVector.h"
#include "StandardProblemIteration.h"
#include "ElementFileWriter.h"
#include <vector>
#include <list>

namespace AMDiS {

  template<typename t1, typename t2> class AbstractFunction;
  class Operator;
  class SystemVector;
  class DOFMatrix;
  class FiniteElemSpace;
  class Estimator;
  class Marker;
  class AdaptStationary;
  class AdaptInfo;
  class FileWriterInterface;
  class CoarseningManager;
  class RefinementManager;
  class Mesh;
  template<typename T> class OEMSolver;
  template<typename T> class Preconditioner;
  template<typename T> class MatVecMultiplier;

  class ProblemVec : public ProblemStatBase,
		     public StandardProblemIteration
  {
  public:
    /** \brief
     * constructor
     */
    ProblemVec(const char* name,
	       ProblemIterationInterface *problemIteration = NULL)
      : StandardProblemIteration(this),
	name_(name),
	nComponents(-1),
	solver_(NULL),
	solution_(NULL),
	rhs_(NULL),
	systemMatrix_(NULL),
	matVec_(NULL),
	leftPrecon_(NULL),
	rightPrecon_(NULL),
	useGetBound_(true),
	info_(10),
	allowFirstRef_(false),
	computeExactError(false)
    {
      GET_PARAMETER(0, name_ + "->components", "%d", &nComponents);
      TEST_EXIT(nComponents > 0)("components not set!\n");    

      estimator_.resize(nComponents, NULL);
      marker.resize(nComponents, NULL);

      assembleMatrixOnlyOnce_.resize(nComponents);
      assembledMatrix_.resize(nComponents);
      for (int i = 0; i < nComponents; i++) {
	assembleMatrixOnlyOnce_[i].resize(nComponents);
	assembledMatrix_[i].resize(nComponents);
	for (int j = 0; j < nComponents; j++) {
	  assembleMatrixOnlyOnce_[i][j] = false;
	  assembledMatrix_[i][j] = false;
	}
      }

      exactSolutionFcts.resize(nComponents);
    };

    /** \brief
     * destructor
     */
    virtual ~ProblemVec() {};

    /** \brief
     * Initialisation of the problem.
     */
    virtual void initialize(Flag initFlag,
			    ProblemVec *adoptProblem = NULL,
			    Flag adoptFlag = INIT_NOTHING);


    /** \brief
     * Used in \ref initialize().
     */
    virtual void createMesh();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void createFESpace();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void createMatricesAndVectors();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void createSolver();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void createEstimator();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void createMarker();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void createFileWriter();

    /** \brief
     * Used in \ref initialize().
     */
    virtual void doOtherStuff();

    /** \brief
     * Implementation of ProblemStatBase::solve(). Deligates the solving
     * of problems system to \ref solver.
     */
    virtual void solve(AdaptInfo *adaptInfo);

    /** \brief
     * Implementation of ProblemStatBase::estimate(). Deligates the estimation
     * to \ref estimator.
     */
    virtual void estimate(AdaptInfo *adaptInfo);

    /** \brief
     * Implementation of ProblemStatBase::markElements().
     * Deligated to \ref adapt.
     */
    virtual Flag markElements(AdaptInfo *adaptInfo);

    /** \brief
     * Implementation of ProblemStatBase::refineMesh(). Deligated to the
     * RefinementManager of \ref adapt.
     */
    virtual Flag refineMesh(AdaptInfo *adaptInfo);

    /** \brief
     * Implementation of ProblemStatBase::coarsenMesh(). Deligated to the
     * CoarseningManager of \ref adapt.
     */
    virtual Flag coarsenMesh(AdaptInfo *adaptInfo);

    /** \brief
     * Implementation of ProblemStatBase::buildBeforeRefine().
     * Does nothing here.
     */
    virtual void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) {};

    /** \brief
     * Implementation of ProblemStatBase::buildBeforeCoarsen().
     * Does nothing here.
     */
    virtual void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) {};

    /** \brief
     * Implementation of ProblemStatBase::buildAfterCoarsen().
     * Assembles \ref A and \ref rhs.
     */
    virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag);

    /** \brief
     * Determines the execution order of the single adaption steps. If adapt is
     * true, mesh adaption will be performed. This allows to avoid mesh adaption,
     * e.g. in timestep adaption loops of timestep adaptive strategies.
     */
    virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);

    /** \brief
     * Returns number of managed problems
     */
    virtual int getNumProblems() { 
      return 1; 
    };

    /** \brief
     * Implementation of ProblemStatBase::getNumComponents()
     */
    virtual int getNumComponents() { 
      return nComponents; 
    }; 

    /** \brief
     * Returns the problem with the given number. If only one problem
     * is managed by this master problem, the number hasn't to be given.
     */
    virtual ProblemStatBase *getProblem(int number = 0) { 
      return this; 
    };

    /** \brief
     * Writes output files.
     */
    void writeFiles(AdaptInfo *adaptInfo, bool force);

    /** \brief
     * Startes parallel output writing;
     */
    void writeDelayedFiles();

    /** \brief
     * Returns true, if there is calculation waiting to be started.
     */
    bool existsDelayedCalculation();

    /** \brief
     * Interpolates fct to \ref solution.
     */
    void interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct);

    /** \brief
     * Adds an Operator to \ref A.
     */
    void addMatrixOperator(Operator *op, int i, int j,
			   double *factor = NULL,
			   double *estFactor = NULL);

    /** \brief
     * Adds an Operator to \ref rhs.
     */
    void addVectorOperator(Operator *op, int i,
			   double *factor = NULL,
			   double *estFactor = NULL);

    /** \brief
     * Adds dirichlet boundary conditions.
     */
    virtual void addDirichletBC(BoundaryType type, int system,
				AbstractFunction<double, WorldVector<double> > *b);

    /** \brief
     * Adds neumann boundary conditions.
     */
    virtual void addNeumannBC(BoundaryType type, int row, int col, 
			      AbstractFunction<double, WorldVector<double> > *n);

    /** \brief
     * Adds robin boundary conditions.
     */
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    AbstractFunction<double, WorldVector<double> > *n,
			    AbstractFunction<double, WorldVector<double> > *r);

    /** \brief
     * Adds periodic boundary conditions.
     */
    virtual void addPeriodicBC(BoundaryType type, int row, int col);


    /** \brief
     * Implementation of ProblemStatBase::allowFirstRefinement().
     */
    inline void allowFirstRefinement() {
      allowFirstRef_ = true;
    }; 
  
    // ===== getting-methods ======================================================

    /** \name getting methods
     * \{ 
     */

    /** \brief
     * Returns \ref solution_.
     */
    inline SystemVector* getSolution() { 
      return solution_; 
    };

    /** \brief
     * Returns \ref rhs_.
     */
    inline SystemVector* getRHS() { 
      return rhs_; 
    };

    /** \brief
     * Returns \ref systemMatrix_.
     */
    inline Matrix<DOFMatrix*> *getSystemMatrix() { 
      return systemMatrix_; 
    };

    /** \brief
     * Returns mesh of given component
     */
    inline Mesh* getMesh(int comp) {
      FUNCNAME("ProblemVec::getMesh()");
      TEST_EXIT(comp < static_cast<int>(componentMeshes_.size()) && comp >= 0)
	("invalid component number\n");
      return componentMeshes_[comp]; 
    };

    /** \brief
     * Returns \ref meshes_
     */
    inline std::vector<Mesh*> getMeshes() {
      return meshes_; 
    };

    /** \brief
     * Returns \ref feSpace_.
     */
    inline FiniteElemSpace* getFESpace(int comp) { 
      FUNCNAME("ProblemVec::getFESpace()");
      TEST_EXIT(comp < static_cast<int>(componentSpaces.size()) && comp >= 0)
	("invalid component number\n");
      return componentSpaces[comp]; 
    };

    /** \brief
     * Returns \ref feSpaces.
     */
    inline std::vector<FiniteElemSpace*> getFESpaces() { 
      return feSpaces; 
    };

    /** \brief
     * Returns \ref componentSpaces;
     */
    inline std::vector<FiniteElemSpace*> getComponentFESpaces() {
      return componentSpaces;
    }

    /** \brief
     * Returns \ref estimator_.
     */
    inline std::vector<Estimator*> getEstimator() { 
      return estimator_; 
    };

    /** \brief
     * Returns \ref estimator_.
     */
    inline Estimator* getEstimator(int comp) { 
      return estimator_[comp]; 
    };

    /** \brief
     * Returns \ref refinementManager_.
     */
    inline RefinementManager* getRefinementManager(int comp) { 
      return refinementManager_; 
    };

    /** \brief
     * Returns \ref refinementManager_.
     */
    inline CoarseningManager* getCoarseningManager(int comp) { 
      return coarseningManager_; 
    };

    /** \brief
     * Returns \ref solver_.
     */
    inline OEMSolver<SystemVector>* getSolver() { 
      return solver_; 
    };

    /** \brief
     * Returns \ref marker.
     */
    inline Marker *getMarker(int comp) { 
      return marker[comp]; 
    };

    /** \brief
     * Returns \ref marker.
     */
    inline std::vector<Marker*> getMarker() { 
      return marker; 
    };

    /** \brief
     * Returns \ref name_.
     */
    inline virtual const std::string& getName() { 
      return name_; 
    };

    /** \brief
     * Returns \ref useGetBound_.
     */
    inline bool getBoundUsed() { 
      return useGetBound_; 
    };

    /** \brief
     * Returns \ref leftPrecon_.
     */
    inline Preconditioner<SystemVector> *getLeftPrecon() { 
      return leftPrecon_; 
    };

    /** \brief
     * Returns \ref rightPrecon_.
     */
    inline Preconditioner<SystemVector> *getRightPrecon() { 
      return rightPrecon_; 
    };

    /** \} */

    // ===== setting-methods ======================================================

    /** \name setting methods
     * \{ 
     */

    /** \brief
     * Sets \ref estimator_.
     */
    inline void setEstimator(std::vector<Estimator*> est) { 
      estimator_ = est; 
    };

    inline void setFESpace(FiniteElemSpace *feSpace, int comp) {
      feSpaces[comp] = feSpace;
    };

    /** \brief
     * Sets \ref estimator_.
     */
    inline void setEstimator(Estimator* est, int comp) { 
      estimator_[comp] = est; 
    };

    inline void setMarker(Marker* mark, int comp) { 
      marker[comp] = mark; 
    };

    /** \brief
     * Sets \ref solver_.
     */
    inline void setSolver(OEMSolver<SystemVector>* sol) { 
      solver_ = sol; 
    };

    /** \brief
     * Sets \ref leftPrecon_.
     */
    inline void setLeftPrecon(Preconditioner<SystemVector> *p) {
      leftPrecon_ = p;
    };

    /** \brief
     * Sets \ref rightPrecon_.
     */
    inline void setRightPrecon(Preconditioner<SystemVector> *p) {
      rightPrecon_ = p;
    };

    inline void setAssembleMatrixOnlyOnce(int i, int j, bool value = true) {
      assembleMatrixOnlyOnce_[i][j] = value;
    }

    void setExactSolutionFct(AbstractFunction<double, WorldVector<double> > *fct,
			     int component) {
      exactSolutionFcts[component] = fct;
    }

    AbstractFunction<double, WorldVector<double> > *getExactSolutionFct(int component) {
      return exactSolutionFcts[component];
    }

    void setComputeExactError(bool v) {
      computeExactError = v;
    }

    /** \} */

    void writeResidualMesh(AdaptInfo *adaptInfo, const std::string name);

    // ===== Serializable implementation =====
  
    /** \brief
     * serialization
     */
    virtual void serialize(std::ostream &out);

    /** \brief
     * deserialization
     */
    virtual void deserialize(std::istream &in);

    std::vector<FileWriterInterface*> getFileWriterList() {
      return fileWriters_;
    };

  protected:
    /** \brief
     * Name of this problem.
     */
    std::string name_;

    /** \brief
     * number of problem components
     */
    int nComponents;

    /** \brief
     * fe spaces of this problem.
     */
    std::vector<FiniteElemSpace*> feSpaces;

    /** \brief
     * meshes of this problem.
     */
    std::vector<Mesh*> meshes_;

    /** \brief
     * Pointer to the fe spaces for the different problem components
     */
    std::vector<FiniteElemSpace*> componentSpaces;

    /** \brief
     * Pointer to the meshes for the different problem components
     */
    std::vector<Mesh*> componentMeshes_;

    /** \brief
     * Responsible for element marking.
     */
    std::vector<Marker*> marker;

    /** \brief
     * Estimator of this problem. Used in \ref estimate().
     */
    std::vector<Estimator*> estimator_;

    /** \brief
     * Linear solver of this problem. Used in \ref solve().
     */
    OEMSolver<SystemVector> *solver_;

    /** \brief
     * system vector  storing the calculated solution of the problem.
     */
    SystemVector *solution_;

    /** \brief
     * system vector for the right hand side 
     */
    SystemVector *rhs_;

    /** \brief
     * System matrix
     */
    Matrix<DOFMatrix*> *systemMatrix_;

    /** \brief
     * Some DOFMatrices of the systemMatrix may be assembled only once (for
     * example if they are independent of the time or older solutions). If
     * [i][j] of this field is set to true, the corresponding DOFMatrix will
     * be assembled only once. All other matrices will be assembled at every
     * time step.
     */
    std::vector< std::vector< bool > > assembleMatrixOnlyOnce_;

    /** \brief
     * If [i][j] of this field is set to true, the corresponding DOFMatrix of
     * the systemMatrix_ has been assembled at least once. This field is used
     * to determine, if assembling of a matrix can be ommitted, if it is set
     * to be assembled only once.
     */
    std::vector< std::vector< bool > > assembledMatrix_;

    /** \brief
     * Matrix-vector multiplication
     */
    MatVecMultiplier<SystemVector> *matVec_;

    /** \brief
     * Left preconditioner. Used in \ref solver.
     */
    Preconditioner<SystemVector> *leftPrecon_;

    /** \brief
     * Right preconditioner. Used in \ref solver.
     */
    Preconditioner<SystemVector> *rightPrecon_;

    /** \brief
     * Determines whether domain boundaries should be considered at assembling.
     */
    bool useGetBound_;

    /** \brief
     * Writes the meshes and solution after the adaption loop.
     */
    std::vector<FileWriterInterface*> fileWriters_;

    /** \brief
     * All actions of mesh refinement are performed by refinementManager.
     * If new refinement algorithms should be realized, one has to override
     * RefinementManager and give one instance of it to AdaptStationary.
     */
    RefinementManager *refinementManager_;

    /** \brief
     * All actions of mesh coarsening are performed by coarseningManager.
     * If new coarsening algorithms should be realized, one has to override
     * CoarseningManager and give one instance of it to AdaptStationary.
     */
    CoarseningManager *coarseningManager_;
  
    /** \brief
     * Info level.
     */
    int info_;

    /** \brief
     * Allows one refinement although the adapt tolerance is reached already.
     */
    bool allowFirstRef_;

    ProblemIterationInterface *problemIteration_;

    /** \brief
     * Used for mesh traversal.
     */
    static ProblemVec *traversePtr_;

    /** \brief
     * This vectors stores pointers to functions defining the exact solution of
     * the problem. This may be used to compute the real error of the computed
     * solution.
     */
    std::vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;

    /** \brief
     * If true, the error is not estimated but computed from the exact solution
     * defined by \ref exactSolutionFcts.
     */
    bool computeExactError;
  };
}

#endif