ProblemStat.h

// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



/** \file ProblemStat.h */

#ifndef AMDIS_PROBLEM_STAT_H
#define AMDIS_PROBLEM_STAT_H

#include <vector>
#include <list>
#include "AMDiS_fwd.h"
#include "ProblemStatBase.h"
#include "Initfile.h"
#include "Boundary.h"
#include "MatrixVector.h"
#include "StandardProblemIteration.h"
#include "io/ElementFileWriter.h"
#include "ComponentTraverseInfo.h"
#include "AbstractFunction.h"
#include "SolverMatrix.h"
#include "SystemVector.h"

namespace AMDiS {

  using namespace std;

  struct OperatorPos 
  {
    int row, col;
    double *factor, *estFactor;
    Flag operatorType;
  };

  /** \brief
   * This class defines the stationary problem definition in sequential
   * computations. For parallel computations, see 
   * \ref ParallelProblemStatBase.
   */
  class ProblemStatSeq : public ProblemStatBase,
			 public StandardProblemIteration
  {
  protected:
    // Defines a mapping type from dof indices to world coordinates.
    typedef map<DegreeOfFreedom, WorldVector<double> > DofToCoord;

    // Defines a mapping type from dof indices to world coordinates.
    typedef map<WorldVector<double>, DegreeOfFreedom> CoordToDof;

  public:
    /// Constructor
    ProblemStatSeq(string nameStr,
	       ProblemIterationInterface *problemIteration = NULL)
      : StandardProblemIteration(this),
	name(nameStr),
	nComponents(-1),
	nMeshes(0),
	traverseInfo(0),
	solver(NULL),
	solution(NULL),
	rhs(NULL),
	systemMatrix(NULL),
	useGetBound(true),
	info(10),
	deserialized(false),
	computeExactError(false),
	boundaryConditionSet(false),
	writeAsmInfo(false)
    {
      Parameters::get(name + "->components", 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);
    }

    /// Destructor
    virtual ~ProblemStatSeq() {}

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

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

    /// Used in \ref initialize().
    virtual void createFeSpace(DOFAdmin *admin);

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

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

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

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

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

    /** \brief
     * Used in \ref initialize(). This function is deprecated and should not be used
     * anymore. There is no guarantee that it will work in the future.
     */
    virtual void doOtherStuff();

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

    /** \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. With the last two parameters, assembling
     * can be restricted to matrices or vectors only.
     */
    virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				   bool assembleMatrix = true,
				   bool assembleVector = true);

    void buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag);

    bool dualMeshTraverseRequired();

    void dualAssemble(AdaptInfo *adaptInfo, Flag flag, 
		      bool asmMatrix = true, bool asmVector = true);

    void createPrecon();


    /** \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);

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

    /// 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; 
    }

    /// Writes output files. TODO: Make obsolete.
    void writeFiles(AdaptInfo *adaptInfo, bool force);

    /// Writes output files.
    void writeFiles(AdaptInfo &adaptInfo, bool force);

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

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

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

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

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

    /// Adds a Dirichlet boundary condition, where the rhs is given by an 
    /// abstract function.
    virtual void addDirichletBC(BoundaryType type, int row, int col,
				AbstractFunction<double, WorldVector<double> > *b);

    /// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
    /// vector.
    virtual void addDirichletBC(BoundaryType type, int row, int col,
				DOFVector<double> *vec);

    /// Adds a Neumann boundary condition, where the rhs is given by an
    /// abstract function.
    virtual void addNeumannBC(BoundaryType type, int row, int col, 
			      AbstractFunction<double, WorldVector<double> > *n);

    /// Adds a Neumann boundary condition, where the rhs is given by an DOF
    /// vector.
    virtual void addNeumannBC(BoundaryType type, int row, int col, 
			      DOFVector<double> *n);

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

    /// Adds Robin boundary condition.
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    DOFVector<double> *n,
			    DOFVector<double> *r);

    /// Adds Robin boundary condition.
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    Operator *n,
			    Operator *r);

    /// Adds a periodic boundary condition.
    virtual void addPeriodicBC(BoundaryType type, int row, int col);

    /** \brief
     * This function assembles a DOFMatrix and a DOFVector for the case,
     * the meshes from row and col FE-space are equal.
     */
    void assembleOnOneMesh(FiniteElemSpace *feSpace, 
			   Flag assembleFlag,
			   DOFMatrix *matrix, DOFVector<double> *vector);


    ///
    void assembleBoundaryConditions(DOFVector<double> *rhs,
				    DOFVector<double> *solution,
				    Mesh *mesh,
				    Flag assembleFlag);
  
    /** \name getting methods
     * \{ 
     */

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

    inline DOFVector<double>* getSolution(int i)
    {
      return solution->getDOFVector(i);
    }

    /// Returns \ref rhs.
    inline SystemVector* getRhs() 
    { 
      return rhs; 
    }

    inline const DOFVector<double>* getRhsVector(int i = 0)
    {
      return rhs->getDOFVector(i);
    }

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

    /// Returns a pointer to the corresponding DOFMatrix.
    inline DOFMatrix* getSystemMatrix(int row, int col) 
    {
      return (*systemMatrix)[row][col];
    }

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

    /// Returns \ref meshes
    inline vector<Mesh*> getMeshes() 
    {
      return meshes; 
    }

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

    /// Returns \ref feSpaces.
    inline vector<FiniteElemSpace*> getFeSpaces() 
    { 
      return feSpaces; 
    }

    /// Returns \ref componentSpaces;
    inline vector<FiniteElemSpace*> getComponentFeSpaces() 
    {
      return componentSpaces;
    }

    /// Returns \ref estimator.
    inline vector<Estimator*> getEstimators() 
    { 
      return estimator; 
    }

    /// Returns \ref estimator.
    inline Estimator* getEstimator(int comp = 0) 
    { 
      return estimator[comp]; 
    }

    /// Returns \ref refinementManager.
    inline RefinementManager* getRefinementManager(int comp = 0) 
    { 
      return refinementManager; 
    }

    /// Returns \ref refinementManager.
    inline CoarseningManager* getCoarseningManager(int comp = 0) 
    { 
      return coarseningManager; 
    }

    /// Returns \ref solver.
    inline OEMSolver* getSolver() 
    { 
      return solver; 
    }

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

    /// Returns \ref marker.
    inline vector<Marker*> getMarkers() 
    { 
      return marker; 
    }

    /// Returns the name of the problem
    inline virtual string getName() 
    { 
      return name; 
    }

    /// Returns \ref useGetBound.
    inline bool getBoundUsed() 
    { 
      return useGetBound; 
    }

    /// Returns \ref info.
    int getInfo()
    {
      return info;
    }

    /// Returns \ref deserialized;
    bool isDeserialized()
    {
      return deserialized;
    }

    /** \} */

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

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

    /// Sets the FE space for the given component.
    inline void setFeSpace(FiniteElemSpace *feSpace, int comp = 0) 
    {
      feSpaces[comp] = feSpace;
    }

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

    /// Sets \ref marker.
    inline void setMarker(Marker* mark, int comp = 0) 
    { 
      marker[comp] = mark; 
    }

    /// Sets \ref solver.
    inline void setSolver(OEMSolver* sol) 
    { 
      solver = sol; 
    }

    ///
    inline void setAssembleMatrixOnlyOnce(int i = 0, int j = 0, 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 i = 0) 
    {
      return exactSolutionFcts[i];
    }

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

    /// Sets \ref writeAsmInfo;
    void setWriteAsmInfo(bool b)
    {
      writeAsmInfo = b;
    }

    /** \} */

    /** \brief
     * Outputs the mesh of the given component, but the values are taken from the 
     * residual error estimator. 
     */
    void writeResidualMesh(int comp, AdaptInfo *adaptInfo, string name);

    /// Function that implements the serialization procedure.
    virtual void serialize(ostream &out);

    /// Function that implements the deserialization procedure.
    virtual void deserialize(istream &in);


    /// Returns \ref fileWriters.
    vector<FileWriterInterface*>& getFileWriterList() 
    {
      return fileWriters;
    }

  protected:
    /** \brief
     * If the exact solution is known, the problem can compute the exact
     * error instead of the error estimation. This is done in this function.
     */
    void computeError(AdaptInfo *adaptInfo);

  protected:
    
    /// Name of this problem.
    string name;

    /// Number of problem components
    int nComponents;

    /** \brief
     * Number of problem meshes. If all components are defined on the same mesh,
     * this number is 1. Otherwise, this variable is the number of different meshes
     * within the problem.
     */
    int nMeshes;

    /// FE spaces of this problem.
    vector<FiniteElemSpace*> feSpaces;

    /// Meshes of this problem.
    vector<Mesh*> meshes;

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

    /// Pointer to the meshes for the different problem components
    vector<Mesh*> componentMeshes;

    /** \brief
     * Stores information about which meshes must be traversed to assemble the
     * specific components. I.e., it was implemented to make use of different
     * meshes for different components.
     */
    ComponentTraverseInfo traverseInfo;
    
    /// Responsible for element marking.
    vector<Marker*> marker;

    /// Estimator of this problem. Used in \ref estimate().
    vector<Estimator*> estimator;

    /// Linear solver of this problem. Used in \ref solve().
    OEMSolver *solver;

    /// System vector  storing the calculated solution of the problem.
    SystemVector *solution;

    /// System vector for the right hand side 
    SystemVector *rhs;

    /// System matrix
    Matrix<DOFMatrix*> *systemMatrix;

    /// Composed system matrix
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;

    /** \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.
     */
    vector<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.
     */
    vector<vector<bool> > assembledMatrix;

    /// Determines whether domain boundaries should be considered at assembling.
    bool useGetBound;

    /// Writes the meshes and solution after the adaption loop.
    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;
  
    /// Info level.
    int info;

    /// If true, the stationary problem was deserialized from some serialization file.
    bool deserialized;

    /** \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.
     */
    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;
    
    /** \brief
     * If at least on boundary condition is set, this variable is true. It is used
     * to ensure that no operators are added after boundary condition were set. If
     * this would happen, boundary conditions could set wrong on off diagonal matrices.
     */
    bool boundaryConditionSet;

    /// If true, AMDiS prints information about the assembling procedure to the screen.
    bool writeAsmInfo;

    map<Operator*, vector<OperatorPos> > operators;

    map<Operator*, Flag> opFlags;
  };

#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
  typedef ProblemStatSeq ProblemStat;
#endif
}

#endif