From 79095eeb2af5ebab966ad46061ede7432e96e2b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20M=C3=BCller?= <felix.mueller2@mailbox.tu-dresden.de> Date: Fri, 23 Mar 2018 15:14:52 +0100 Subject: [PATCH] implemented markElements() function --- src/amdis/ProblemStat.hpp | 83 +++++++++++++++++++++++++++- src/amdis/ProblemStat.inc.hpp | 101 +++++++++++++++++++++++++++++++++- 2 files changed, 182 insertions(+), 2 deletions(-) diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index 28ebe217..d42219bf 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -16,9 +16,11 @@ #include <amdis/Assembler.hpp> #include <amdis/CreatorInterface.hpp> #include <amdis/DirichletBC.hpp> +//#include <amdis/Estimator.hpp> #include <amdis/Flag.hpp> #include <amdis/Initfile.hpp> #include <amdis/LinearAlgebra.hpp> +#include <amdis/Marker.cpp> #include <amdis/Mesh.hpp> #include <amdis/ProblemStatBase.hpp> #include <amdis/ProblemStatTraits.hpp> @@ -66,7 +68,7 @@ namespace AMDiS public: /** * \brief Constructor. Takes the name of the problem that is used to - * access values correpsonding to this püroblem in the parameter file. + * access values corresponding to this problem in the parameter file. **/ explicit ProblemStat(std::string name) : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this)) @@ -143,6 +145,18 @@ namespace AMDiS bool storeMatrixData = false) override; + /// Implementation of \ref ProblemStatBase::estimate. + virtual void estimate(AdaptInfo& adaptInfo) override; + + + /// Implementation of \ref ProblemStatBase::markElements. + virtual Flag markElements(AdaptInfo& adaptInfo) override; + + + /// Implementation of \ref ProblemStatBase::refineMesh. + virtual Flag refineMesh(AdaptInfo& adaptInfo) override; + + /// Implementation of \ref ProblemStatBase::buildAfterCoarse virtual void buildAfterCoarsen(AdaptInfo& adaptInfo, Flag flag, @@ -263,6 +277,14 @@ namespace AMDiS systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat"); solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution"); rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs"); + + estimates.resize(leafGridView().indexSet().size(0)); + for (std::size_t i = 0; i < estimates.size(); i++) { + estimates[i].resize(nComponents); + for (std::size_t j = 0; j < estimates[i].size(); j++) { + estimates[i][j] = 0.0; // TODO: Remove when estimate() is implemented + } + } } void createSolver() @@ -276,6 +298,53 @@ namespace AMDiS linearSolver_ = solverCreator->create(name_ + "->solver"); } + void createEstimator() + {/* + for (std::size_t i = 0, i < nComponents, i++) { + std::string estName = ""; + Parameters::get(name + "->estimator->name[" + std::to_string(i) + "]", estName); + auto estimatorCreator + = named(CreatorMap<Estimator>::getCreator(estName, name + "->estimator->name[" + std::to_string(i) + "]")); + if (estimatorCreator) { + estimatorCreator->setName(estName); + estimatorCreator->setRow(i); + estimatorCreator->setSolution(solution->getDOFVector(i)); + estimator[i] = estimatorCreator->create(); + } + + if (estimator[i]) { + for (std::size_t j = 0; j < nComponents; j++) + estimator[i]->addSystem((*systemMatrix)[i][j], + solution->getDOFVector(j), + rhs->getDOFVector(i)); + } + }*/ + } + + + void createMarker() + { + int nMarkersCreated = 0; + + marker.resize(nComponents); + + for (std::size_t i = 0; i < nComponents; i++) { + marker[i] = Marker<Grid>:: + createMarker(name + "->marker[" + std::to_string(i) + "]", i, + estimates, componentGrids[i]); + + if (marker[i]) { + nMarkersCreated++; + + // If there is more than one marker, and all components are defined + // on the same mesh, the maximum marking has to be enabled. + if (nMarkersCreated > 1 && nGrids == 1) + marker[i]->setMaximumMarking(true); + } + } + } + + void createFileWriter(); @@ -321,6 +390,9 @@ namespace AMDiS /// Grid of this problem. std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems + /// Number of grids + int nGrids = 1; + /// Name of the mesh std::string gridName_ = "none"; @@ -331,6 +403,12 @@ namespace AMDiS /// A FileWriter object std::list<std::shared_ptr<FileWriterInterface>> filewriter_; + /// Pointer to the adaptation markers + std::vector<Marker<Grid>* > marker; + + /// Pointer to the estimators for this problem +// std::vector<Estimator*> estimator; + /// An object of the linearSolver Interface std::shared_ptr<LinearSolverType> linearSolver_; @@ -340,6 +418,9 @@ namespace AMDiS /// A block-vector with the solution components std::shared_ptr<SystemVector> solution_; + /// A vector with the local element error estimates + std::vector<std::vector<double> > estimates; + /// A block-vector (load-vector) corresponding to the right.hand side /// of the equation, filled during assembling std::shared_ptr<SystemVector> rhs_; diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 4499009b..c963a668 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -13,6 +13,7 @@ #include <amdis/LocalAssembler.hpp> #include <amdis/GridFunctionOperator.hpp> #include <amdis/common/Loops.hpp> +// #include <amdis/Estimator.hpp> namespace AMDiS { @@ -22,7 +23,7 @@ void ProblemStat<Traits>::initialize( Self* adoptProblem, Flag adoptFlag) { - // create grides + // create grids if (grid_) { warning("grid already created"); } @@ -78,6 +79,7 @@ void ProblemStat<Traits>::initialize( if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) { solution_ = adoptProblem->solution_; + estimates = adoptProblem->estimates; rhs_ = adoptProblem->rhs_; systemMatrix_ = adoptProblem->systemMatrix_; } @@ -101,6 +103,24 @@ void ProblemStat<Traits>::initialize( warning("no solver created\n"); } +/* + // create estimator + estimator.resize(nComponents, NULL); + + if (initFlag.isSet(INIT_ESTIMATOR)) + createEstimator(); + + if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) + estimator = adoptProblem->estimator; +*/ + + // create marker + if (initFlag.isSet(INIT_MARKER)) + createMarker(); + + if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) + marker = adoptProblem->marker; + // create file writer if (initFlag.isSet(INIT_FILEWRITER)) @@ -280,6 +300,84 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) } +template <class Traits> +void ProblemStat<Traits>:: +estimate(AdaptInfo& adaptInfo) +{ + Dune::Timer t; +/* + for (std::size_t i = 0; i < nComponents; i++) { + Estimator *scalEstimator = estimator[i]; + + if (scalEstimator) { + scalEstimator->estimate(adaptInfo->getTimestep()); + + adaptInfo->setEstSum(scalEstimator->getErrorSum(), i); + adaptInfo->setEstMax(scalEstimator->getErrorMax(), i); + + if (adaptInfo->getRosenbrockMode() == false) { + adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i); + adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i); + } + } + } + +//#ifdef HAVE_PARALLEL_DOMAIN_AMDIS +// MPI::COMM_WORLD.Barrier(); +//#endif +*/ + msg("estimation of the error needed ", t.elapsed(), "seconds\n"); +} + + +template <class Traits> +Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo) +{ + Dune::Timer t; + + test_exit(static_cast<unsigned int>(nComponents) == marker.size(), "Wrong number of markers!\n"); + + Flag markFlag = 0; + for (std::size_t i = 0; i < nComponents; i++) + if (marker[i]) + markFlag |= marker[i]->markGrid(adaptInfo); + + msg("markElements needed ", t.elapsed(), " seconds"); + + return markFlag; +} + + +template <class Traits> +Flag ProblemStat<Traits>::refineMesh(AdaptInfo& adaptInfo) +{ + // TODO: data transfer + + grid->preAdapt(); +/* + const auto& gridView = grid->leafGridView(); + const auto& indexSet = gridView.indexSet(); + const auto& idSet = grid−>localIdSet(); + + // Save data in container during adaptation + for (const auto& v : vertices(gridView)) { + persistentContainer[idSet.id(v)] = data[indexSet.index(v)]; + } */ + + Flag flag = grid->adapt(); +/* + // Unpack data + data.resize(gridView.size(dim)); + for (const auto& v : vertices(gridView)) { + data[indexSet.index(v)] = persistentContainer[idSet.id(v)]; + }*/ + + grid->postAdapt(); + + return flag; +} + + template <class Traits> void ProblemStat<Traits>:: buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector) @@ -305,4 +403,5 @@ writeFiles(AdaptInfo& adaptInfo, bool force) msg("writeFiles needed ", t.elapsed(), " seconds"); } + } // end namespace AMDiS -- GitLab