/****************************************************************************** * * AMDiS - Adaptive multidimensional simulations * * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis * * Authors: * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. * * * 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 "solver/SolverMatrix.h" #include "SystemVector.h" namespace AMDiS { using namespace std; struct OperatorPos { int row, col; double *factor, *estFactor; }; /** * \ingroup Problem * * \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 = nullptr); /// Destructor virtual ~ProblemStatSeq(); /// Initialisation of the problem. virtual void initialize(Flag initFlag, ProblemStatSeq *adoptProblem = nullptr, Flag adoptFlag = INIT_NOTHING); /// Used in \ref initialize(). virtual void createMesh(); /// Used in \ref initialize(). virtual void createRefCoarseManager(); /// 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(); /// 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(); /// Implementation of ProblemStatBase::solve(). Deligates the solving /// of problems system to \ref solver. void solve(AdaptInfo *adaptInfo, bool createMatrixData = true, bool storeMatrixData = false) override; /// Implementation of ProblemStatBase::estimate(). Deligates the estimation /// to \ref estimator. void estimate(AdaptInfo *adaptInfo) override; /// Implementation of ProblemStatBase::markElements(). /// Deligated to \ref adapt. Flag markElements(AdaptInfo *adaptInfo) override; /// Implementation of ProblemStatBase::refineMesh(). Deligated to the /// RefinementManager of \ref adapt. Flag refineMesh(AdaptInfo *adaptInfo) override; /// Implementation of ProblemStatBase::coarsenMesh(). Deligated to the /// CoarseningManager of \ref adapt. Flag coarsenMesh(AdaptInfo *adaptInfo) override; /// Implementation of ProblemStatBase::buildBeforeRefine(). /// Does nothing here. void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) override {} /// Implementation of ProblemStatBase::buildBeforeCoarsen(). /// Does nothing here. void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) override {} /// Implementation of ProblemStatBase::buildAfterCoarsen(). /// Assembles \ref A and \ref rhs. With the last two parameters, assembling /// can be restricted to matrices or vectors only. void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag, bool assembleMatrix = true, bool assembleVector = true) override; bool dualMeshTraverseRequired(); void dualAssemble(AdaptInfo *adaptInfo, Flag flag, bool asmMatrix = true, bool asmVector = true); /// 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. Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) override; /// Returns number of managed problems int getNumProblems() override { return 1; } /// Implementation of ProblemStatBase::getNumComponents(), TODO: Wrong!! virtual int getNumComponents() { return nComponents; } /// Returns the problem with the given number. If only one problem /// is managed by this master problem, the number hasn't to be given. ProblemStatBase *getProblem(int number = 0) override { 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 = nullptr, double *estFactor = nullptr); /// Adds an operator to \ref A. void addMatrixOperator(Operator &op, int i, int j, double *factor = nullptr, double *estFactor = nullptr); /// Adds an operator to \ref rhs. void addVectorOperator(Operator *op, int i, double *factor = nullptr, double *estFactor = nullptr); /// Adds an operator to \ref rhs. void addVectorOperator(Operator &op, int i, double *factor = nullptr, double *estFactor = nullptr); /// 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); } /// This function assembles a DOFMatrix and a DOFVector for the case, /// the meshes from row and col FE-space are equal. void assembleOnOneMesh(const 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 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 const 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& getComponentSpaces() { 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 LinearSolver* 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 string getName() override { return name; } /// Returns the name of the problem inline string getComponentName(int comp = 0) { TEST_EXIT(comp < static_cast(componentNames.size()) && comp >= 0) ("invalid component number\n"); return componentNames[comp]; } /// 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(const FiniteElemSpace *feSpace, int comp = 0) { feSpaces[comp] = feSpace; } void setComponentSpace(int comp, const FiniteElemSpace *feSpace) { if (static_cast(componentSpaces.size()) < nComponents) componentSpaces.resize(nComponents); TEST_EXIT(comp >= 0 && comp < nComponents) ("Component number not in feasable range!"); componentSpaces[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(LinearSolver* sol) { solver = sol; } void setSolverMatrix(SolverMatrix >& solverMatrix_) { solverMatrix.setMatrix(*solverMatrix_.getOriginalMat()); } /// 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; } void setMeshes(vector meshes_) { meshes = meshes_; nMeshes = static_cast(meshes.size()); } void setComponentMesh(int comp, Mesh* mesh) { if (static_cast(componentMeshes.size()) < nComponents) componentMeshes.resize(nComponents); TEST_EXIT(comp >= 0 && comp < nComponents) ("Component number not in feasable range!"); componentMeshes[comp] = mesh; } void setRefinementManager(RefinementManager *ref) { refinementManager = ref; } void setCoarseningManager(CoarseningManager *coarse) { coarseningManager = coarse; } /** \} */ /// 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. void serialize(ostream &out) override; /// Function that implements the deserialization procedure. void deserialize(istream &in) override; /// Returns \ref fileWriters. vector& getFileWriterList() { return fileWriters; } protected: /// 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; /// Stores the names for all components. Is used for naming the solution /// vectors, \ref solution. vector componentNames; /// 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; /// 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(). LinearSolver *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; /// 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; /// 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; /// 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; /// 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; /// 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; /// If true, the error is not estimated but computed from the exact solution /// defined by \ref exactSolutionFcts. bool computeExactError; /// 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; }; #ifndef HAVE_PARALLEL_DOMAIN_AMDIS typedef ProblemStatSeq ProblemStat; #endif } #endif