#pragma once #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace AMDiS { // forward declaration template class ProblemInstat; template class ProblemStat : public ProblemStatBase , public StandardProblemIterationAdaptor> { using Self = ProblemStat; friend class ProblemInstat; public: // typedefs and static constants using GlobalBasis = typename Traits::GlobalBasis; using GridView = typename GlobalBasis::GridView; using Grid = AdaptiveGrid_t; using Element = typename GridView::template Codim<0>::Entity; using WorldVector = typename Element::Geometry::GlobalCoordinate; using WorldMatrix = FieldMatrix; /// Dimension of the grid static constexpr int dim = Grid::dimension; /// Dimension of the world static constexpr int dow = Grid::dimensionworld; private: using LinAlgTraits = typename Traits::LinAlgTraits; using Mat = typename LinAlgTraits::template MatrixImpl; using Vec = typename LinAlgTraits::template VectorImpl; using PartitionSet = typename LinAlgTraits::PartitionSet; public: using LinearSolver = LinearSolverInterface; using SystemMatrix = BiLinearForm; using SystemVector = LinearForm; using SolutionVector = DOFVector; public: /** * \brief Constructor. Takes the name of the problem that is used to * access values corresponding to this problem in the parameter file. **/ explicit ProblemStat(std::string const& name) : name_(name) {} /// Constructor taking additionally a grid that is used /// instead of the default created grid, \ref ProblemStat template ProblemStat(std::string const& name, Grid_&& grid) : ProblemStat(name) { adoptGrid(wrap_or_share(FWD(grid))); } /// \brief Constructor taking a grid and basis. /// Wraps both in shared pointers. template , REQUIRES(Concepts::GlobalBasis)> ProblemStat(std::string const& name, Grid_&& grid, Basis_&& globalBasis) : ProblemStat(name, FWD(grid)) { adoptGlobalBasis(wrap_or_share(FWD(globalBasis))); } /// \brief Constructor taking a grid and pre-basis factory to create a global basis /// on the fly. template ::LeafGridView, REQUIRES(Concepts::PreBasisFactory>)> ProblemStat(std::string const& name, Grid_&& grid, PBF_ const& preBasisFactory) : ProblemStat(name, FWD(grid)) { adoptGlobalBasis(makeSharedPtr(GlobalBasis{grid_->leafGridView(), preBasisFactory})); } /// \brief Initialisation of the problem. /** * Parameters read in initialize() * - `[GRID_NAME]->global refinements`: nr of initial global refinements **/ void initialize(Flag initFlag, Self* adoptProblem = nullptr, Flag adoptFlag = INIT_NOTHING); /// \brief Read the grid and solution from backup files and initialize the problem /** * Parameters read in restore() for problem with name 'PROB' * - `[PROB]->restore->grid`: name of the grid backup file * - `[PROB]->restore->solution`: name of the solution backup file **/ void restore(Flag initFlag); /// Add an operator to \ref A. /** @{ */ /// Operator evaluated on the whole element /** * Adds an operator to the list of element operators to be assembled in * quadrature points inside the element. * * \param op A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator * \param row TreePath identifying the sub-basis in the global basis tree * corresponding to the row basis. \see makeTreePath() * \param col TreePath identifying the sub-basis in the global basis tree * corresponding to the column basis. \see makeTreePath() * * Example: * ``` * auto op = makeOperator(tag::test_trial{}, 1.0/tau); * prob.addMatrixOperator(op, _0, _0); * ``` **/ template void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {}) { static constexpr bool isValidTreePath = Concepts::ValidTreePath && Concepts::ValidTreePath; static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addMatrixOperator!"); if constexpr (isValidTreePath) systemMatrix_->addOperator(tag::element_operator{}, op, row, col); } /// Operator evaluated on the boundary of the domain with boundary index `b` /** * Adds an operator to the list of boundary operators to be assembled in * quadrature points on the boundary intersections. * * \param b Boundary identifier where to assemble this operator. Can be * constructed from an integer. \see BoundaryType * \param op A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator * \param row TreePath identifying the sub-basis in the global basis tree * corresponding to the row basis. \see makeTreePath() * \param col TreePath identifying the sub-basis in the global basis tree * corresponding to the column basis. \see makeTreePath() * * Example: * ``` * auto op = makeOperator(tag::test_trial{}, alpha); * prob.addMatrixOperator(BoundaryType{1}, op, _0, _0); * ``` **/ template void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {}) { using I = typename GridView::Intersection; static constexpr bool isValidTreePath = Concepts::ValidTreePath && Concepts::ValidTreePath; static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addMatrixOperator!"); if constexpr (isValidTreePath) systemMatrix_->addOperator(tag::boundary_operator{*boundaryManager_,b}, op, row, col); } /** @} */ /// Add an operator to \ref rhs. /** @{ */ /// Operator evaluated on the whole element /** * Adds an operator to the list of element operators to be assembled in * quadrature points inside the element. * * \param op A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator * \param path TreePath identifying the sub-basis in the global basis tree * corresponding to the row basis. \see makeTreePath() * * Example: * ``` * auto op = makeOperator(tag::test{}, probInstat.getOldSolution(0) / tau); * prob.addVectorOperator(op, _0); * ``` **/ template void addVectorOperator(Operator const& op, TreePath path = {}) { static constexpr bool isValidTreePath = Concepts::ValidTreePath; static_assert(isValidTreePath, "Invalid treepath passed to addVectorOperator!"); if constexpr (isValidTreePath) rhs_->addOperator(tag::element_operator{}, op, path); } /// Operator evaluated on the boundary of the domain with boundary index `b` /** * Adds an operator to the list of boundary operators to be assembled in * quadrature points on the boundary intersections. * * \param b Boundary identifier where to assemble this operator. Can be * constructed from an integer. \see BoundaryType * \param op A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator * \param path TreePath identifying the sub-basis in the global basis tree * corresponding to the row basis. \see makeTreePath() * * Example: * ``` * auto op = makeOperator(tag::test{}, [g](auto const& x) { return g(x); }); * prob.addVectorOperator(BoundaryType{1}, op, _0); * ``` **/ template void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {}) { using I = typename GridView::Intersection; static constexpr bool isValidTreePath = Concepts::ValidTreePath; static_assert(isValidTreePath, "Invalid treepath passed to addVectorOperator!"); if constexpr (isValidTreePath) rhs_->addOperator(tag::boundary_operator{*boundaryManager_,b}, op, path); } /** @} */ /// Add boundary conditions to the system /** @{ */ /// Dirichlet boundary condition /** * Enforce Dirichlet boundary values for the solution vector on boundary * regions identified by the predicate. * * \param predicate Functor `bool(WorldVector)` returning true for all * DOFs on the boundary that should be assigned a value. * \param row TreePath identifying the sub-basis in the global basis tree * corresponding to the row basis. \see makeTreePath() * \param col TreePath identifying the sub-basis in the global basis tree * corresponding to the column basis. \see makeTreePath() * \param values Functor `Range(WorldVector)` or any \ref GridFunction * that is evaluated in the DOFs identified by the predicate. * * Example: * ``` * prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8; }, 0, 0, * [](auto const& x) { return 0.0; }); * ``` **/ template void addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values); template void addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& values); template void addDirichletBC(Identifier&& id, Values&& values) { addDirichletBC(FWD(id), RootTreePath{}, RootTreePath{}, FWD(values)); } /// Add a periodic boundary conditions to the system, by specifying a face transformation /// y = A*x + b of coordinates. We assume, that A is orthonormal. void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b); /** @} */ public: /// Implementation of \ref StandardProblemIteration::oneIteration. Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override { return StandardProblemIteration::oneIteration(adaptInfo, toDo); } /// Implementation of \ref ProblemStatBase::buildAfterCoarse void buildAfterAdapt(AdaptInfo& adaptInfo, Flag flag, bool asmMatrix = true, bool asmVector = true) override; /// \brief Assemble the linear system by calling \ref buildAfterAdapt with /// `asmMatrix` and `asmVector` set to true. void assemble(AdaptInfo& adaptInfo) { buildAfterAdapt(adaptInfo, Flag{0}, true, true); } /// Implementation of \ref ProblemStatBase::solve void solve(AdaptInfo& adaptInfo, bool createMatrixData = true, bool storeMatrixData = false) override; /// Implementation of \ref ProblemStatBase::estimate. void estimate(AdaptInfo& /*adaptInfo*/) override { /* do nothing. */ } /// Implementation of \ref ProblemStatBase::refineMesh. Flag adaptGrid(AdaptInfo& adaptInfo) override; /// Implementation of \ref ProblemStatBase::markElements. Flag markElements(AdaptInfo& adaptInfo) override; /// Uniform global grid coarsening by up to n level Flag globalCoarsen(int n) override; /// Uniform global refinement by n level Flag globalRefine(int n) override; /// Writes output files. If force=true write even if timestep out of write rhythm. void writeFiles(AdaptInfo& adaptInfo, bool force = false); public: // get-methods /// Implementation of \ref ProblemStatBase::name std::string const& name() const override { return name_; } /// Return the \ref grid_ std::shared_ptr grid() { return grid_; } std::shared_ptr grid() const { return grid_; } /// Return the gridView of the basis GridView gridView() const { return globalBasis_->gridView(); } /// Return the boundary manager to identify boundary segments std::shared_ptr> boundaryManager() { return boundaryManager_; } std::shared_ptr const> boundaryManager() const { return boundaryManager_; } /// Return the \ref globalBasis_ std::shared_ptr globalBasis() { return globalBasis_; } std::shared_ptr globalBasis() const { return globalBasis_; } /// Return a reference to the linear solver, \ref linearSolver std::shared_ptr solver() { return linearSolver_; } std::shared_ptr solver() const { return linearSolver_; } /// Returns a reference to system-matrix, \ref systemMatrix_ std::shared_ptr systemMatrix() { return systemMatrix_; } std::shared_ptr systemMatrix() const { return systemMatrix_; } /// Returns a reference to the solution vector, \ref solution_ std::shared_ptr solutionVector() { return solution_; } std::shared_ptr solutionVector() const { return solution_; } /// Return a reference to the rhs system-vector, \ref rhs std::shared_ptr rhsVector() { return rhs_; } std::shared_ptr rhsVector() const { return rhs_; } /// Return a mutable view to a solution component template auto solution(Indices... ii) { assert(bool(solution_) && "You have to call initialize() before."); return valueOf(*solution_, ii...); } /// Return a const view to a solution component template auto solution(Indices... ii) const { assert(bool(solution_) && "You have to call initialize() before."); return valueOf(*solution_, ii...); } public: // set-methods /// Set a new linear solver for the problem template void setSolver(Solver_&& solver) { linearSolver_ = wrap_or_share(FWD(solver)); } /// Set the grid. Stores pointer and initializes feSpaces /// matrices and vectors, as well as markers and file-writers. /// If grid is given as reference, wrap it into a non-destroying shared_ptr template void setGrid(Grid_&& grid) { adoptGrid(wrap_or_share(FWD(grid))); createGlobalBasis(); createMatricesAndVectors(); createMarker(); createFileWriter(); } /// Store the shared_ptr and the name of the marker in the problem /** * Note: multiple markers can be added but must have different names **/ template void addMarker(Marker_&& m) { auto marker = wrap_or_share(FWD(m)); auto it = marker_.emplace(marker->name(), marker); if (marker_.size() > 1) it.first->second->setMaximumMarking(true); } /// Remove a marker with the given name from the problem void removeMarker(std::string name) { std::size_t num = marker_.erase(name); test_warning(num == 1, "A marker with the given name '{}' does not exist.", name); } /// Remove a marker from the problem void removeMarker(Marker const& marker) { removeMarker(marker.name()); } protected: // initialization methods void createGlobalBasis(); void createGrid(); void createMatricesAndVectors(); void createSolver(); void createMarker(); void createFileWriter(); void adoptGlobalBasis(std::shared_ptr globalBasis) { globalBasis_ = std::move(globalBasis); initGlobalBasis(); } void adoptGrid(std::shared_ptr const& grid, std::shared_ptr> const& boundaryManager) { grid_ = grid; boundaryManager_ = boundaryManager; Parameters::get(name_ + "->mesh", gridName_); } void adoptGrid(std::shared_ptr const& grid) { adoptGrid(grid, std::make_shared>(grid)); } void adoptGrid(std::shared_ptr const& hostGrid) { auto grid = std::make_shared(hostGrid); adoptGrid(grid, std::make_shared>(grid)); } private: void createGlobalBasisImpl(std::true_type); void createGlobalBasisImpl(std::false_type); void initGlobalBasis(); private: /// Name of this problem. std::string name_; /// Grid of this problem. std::shared_ptr grid_; /// Name of the grid std::string gridName_ = "mesh"; /// Management of boundary conditions std::shared_ptr> boundaryManager_; /// FE space of this problem. std::shared_ptr globalBasis_; /// A FileWriter object std::list> filewriter_; /// Pointer to the adaptation markers std::map>> marker_; /// Pointer to the estimators for this problem // std::vector estimator; /// An object of the linearSolver Interface std::shared_ptr linearSolver_; /// Matrix that is filled during assembling std::shared_ptr systemMatrix_; /// Vector with the solution components std::shared_ptr solution_; /// Vector (load-vector) corresponding to the right-hand side /// of the equation, filled during assembling std::shared_ptr rhs_; /// A vector with the local element error estimates /// for each node in the basis tree, indexed by [to_string(treePath)][element index] std::map> estimates_; private: // some internal data-structures BoundaryConditions boundaryConditions_; }; namespace Impl { template struct DeducedProblemTraits; template struct DeducedProblemTraits,void> { using type = DefaultProblemTraits>; }; template struct DeducedProblemTraits>>> { using Grid = AdaptiveGrid_t; using GridView = typename Grid::LeafGridView; using Basis = decltype(GlobalBasis{std::declval(),std::declval()}); using type = DefaultProblemTraits; }; template using DeducedProblemTraits_t = typename DeducedProblemTraits::type; } // Deduction guide template ProblemStat(std::string name, Grid&& grid, Basis&& globalBasis) -> ProblemStat,Underlying_t>>; // mark templates as explicitly instantiated in cpp file extern template class ProblemStat>; extern template class ProblemStat>; } // end namespace AMDiS #include "ProblemStat.inc.hpp"