implemented markElements() function

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_;

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