// ============================================================================ // == == // == 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 #include #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 > DofToCoord; // Defines a mapping type from dof indices to world coordinates. typedef map, 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); /** \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 >*> *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 > *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 *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 > *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 *n); /// Adds Robin boundary condition. virtual void addRobinBC(BoundaryType type, int row, int col, AbstractFunction > *n, AbstractFunction > *r); /// Adds Robin boundary condition. virtual void addRobinBC(BoundaryType type, int row, int col, DOFVector *n, DOFVector *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); /// add boundary operator to matrix side virtual void addBoundaryMatrixOperator(BoundaryType type, Operator *op, int row, int col); virtual void addBoundaryMatrixOperator(BoundaryType type, Operator &op, int row, int col) { addBoundaryMatrixOperator(type, &op, row, col); } /// add boundary operator to rhs (vector) side virtual void addBoundaryVectorOperator(BoundaryType type, Operator *op, int row); virtual void addBoundaryVectorOperator(BoundaryType type, Operator &op, int row) { addBoundaryVectorOperator(type, &op, row); } /** \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 *vector); /// void assembleBoundaryConditions(DOFVector *rhs, DOFVector *solution, Mesh *mesh, Flag assembleFlag); /** \name getting methods * \{ */ /// Returns \ref solution. inline SystemVector* getSolution() { return solution; } inline DOFVector* getSolution(int i) { return solution->getDOFVector(i); } /// Returns \ref rhs. inline SystemVector* getRhs() { return rhs; } inline const DOFVector* getRhsVector(int i = 0) { return rhs->getDOFVector(i); } /// Returns \ref systemMatrix. inline Matrix *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(componentMeshes.size()) && comp >= 0) ("invalid component number\n"); return componentMeshes[comp]; } /// Returns \ref meshes inline vector getMeshes() { return meshes; } /// Returns \ref feSpace_. inline FiniteElemSpace* getFeSpace(int comp = 0) { FUNCNAME("ProblemStatSeq::getFeSpace()"); TEST_EXIT(comp < static_cast(componentSpaces.size()) && comp >= 0) ("invalid component number\n"); return componentSpaces[comp]; } /// Returns \ref feSpaces. inline vector getFeSpaces() { return feSpaces; } /// Returns \ref componentSpaces; inline vector getComponentFeSpaces() { return componentSpaces; } /// Returns \ref estimator. inline vector 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 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 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 > *fct, int component) { exactSolutionFcts[component] = fct; } /// AbstractFunction > *getExactSolutionFct(int i = 0) { return exactSolutionFcts[i]; } /// vector< AbstractFunction >* > getExactSolutionFcts() { return exactSolutionFcts; } /// 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& 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 feSpaces; /// Meshes of this problem. vector meshes; /// Pointer to the fe spaces for the different problem components vector componentSpaces; /// Pointer to the meshes for the different problem components vector 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; /// Estimator of this problem. Used in \ref estimate(). vector 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 *systemMatrix; /// Composed system matrix SolverMatrix > 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 > 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 > assembledMatrix; /// Determines whether domain boundaries should be considered at assembling. bool useGetBound; /// Writes the meshes and solution after the adaption loop. vector 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 >*> 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 > operators; map opFlags; }; #ifndef HAVE_PARALLEL_DOMAIN_AMDIS typedef ProblemStatSeq ProblemStat; #endif } #endif