interface for data transfer during grid adaption

src/amdis/CMakeLists.txt

@@ -24,11 +24,13 @@ install(FILES
     ContextGeometry.hpp
     CreatorInterface.hpp
     CreatorMap.hpp
+    DataTransfer.hpp
     DirichletBC.hpp
     FileWriter.hpp
     Flag.hpp
     GridFunctionOperator.hpp
     GridFunctions.hpp
+    GridTransfer.hpp
     Initfile.hpp
     InitfileParser.hpp
     LinearAlgebra.hpp

src/amdis/DataTransfer.hpp

+#pragma once
+
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/utility/TreeData.hpp>
+#include <amdis/utility/Visitor.hpp>
+
+namespace AMDiS
+{
+  typedef enum {
+    NO_OPERATION = 0,
+    INTERPOLATE  = 1
+  } DataTransferOperation;
+
+  template <class Node, class RangeType>
+  class NodeDataTransfer
+  {
+  public:
+    template <class Basis, class Container>
+    void preAdapt(Basis const& basis, Container const& coeff, bool mightCoarsen) {}
+
+    template <class Basis, class Container>
+    void postAdapt(Basis const& basis, Container& coeff, bool refined) const {}
+  };
+
+
+  template <class Container>
+  class DataTransferInterface
+  {
+  public:
+    /// Virtual destructor
+    virtual ~DataTransferInterface() = default;
+
+    /// Collect data that is needed before grid adaption
+    virtual void preAdapt(Container const& container, bool mightCoarsen) = 0;
+
+    /// Interpolate data to new grid after grid adaption
+    virtual void postAdapt(Container& container, bool refined) const = 0;
+  };
+
+
+  /// Implementation of \ref DataTransferInterface that does not interpolation, but
+  /// just resizes the containers to the dimension of the basis
+  template <class Container>
+  class NoDataTransfer
+      : public DataTransferInterface<Container>
+  {
+  public:
+    virtual void preAdapt(Container const& container, bool) override {}
+
+    virtual void postAdapt(Container& container, bool) const override
+    {
+      container.compress();
+    }
+  };
+
+
+  template <class Container, class Basis>
+  class DataTransfer
+      : public DataTransferInterface<Container>
+  {
+    template <class Node>
+    using NDT = NodeDataTransfer<std::decay_t<Node>, typename Container::value_type>;
+
+  public:
+    DataTransfer(Basis const& basis)
+      : basis_(&basis)
+    {}
+
+    /// Calls \ref NodeDataTransfer::preAdapt() on each basis node
+    virtual void preAdapt(Container const& container, bool mightCoarsen) override
+    {
+      nodeDataTransfer_.init(*basis_);
+      AMDiS::forEachLeafNode_(basis_->localView().tree(), [&](auto const& node, auto const& treePath)
+      {
+        auto subBasis = Dune::Functions::subspaceBasis(*basis_, treePath);
+        nodeDataTransfer_[node].preAdapt(subBasis, container.vector(), mightCoarsen);
+      });
+    }
+
+    /// Calls \ref NodeDataTransfer::postAdapt() on each basis node after compressing the
+    /// Container the dimension of the basis
+    virtual void postAdapt(Container& container, bool refined) const override
+    {
+      container.compress();
+      AMDiS::forEachLeafNode_(basis_->localView().tree(), [&](auto const& node, auto const& treePath)
+      {
+        auto subBasis = Dune::Functions::subspaceBasis(*basis_, treePath);
+        nodeDataTransfer_[node].postAdapt(subBasis, container.vector(), refined);
+      });
+    }
+
+  private:
+    Basis const* basis_;
+    TreeData<Basis, NDT, true> nodeDataTransfer_;
+  };
+
+
+  /// Factory to create DataTransfer objects based on the \ref DataTransferOperation
+  template <class Container>
+  class DataTransferFactory
+  {
+    using Interface = DataTransferInterface<Container>;
+
+  public:
+    template <class Basis>
+    static std::unique_ptr<Interface> create(DataTransferOperation op, Basis const& basis)
+    {
+      switch (op)
+      {
+        case NO_OPERATION:
+          return std::make_unique<NoDataTransfer<Container>>();
+        case INTERPOLATE:
+          return std::make_unique<DataTransfer<Container, Basis>>(basis);
+        default:
+          error_exit("Invalid data transfer\n");
+          return nullptr; // avoid warnings
+      }
+    }
+  };
+
+} // end namespace AMDiS

src/amdis/GridTransfer.hpp

+#pragma once
+
+#include <vector>
+
+#include <amdis/linear_algebra/DOFVectorInterface.hpp>
+
+namespace AMDiS
+{
+  template <class Grid>
+  class GridTransfer
+  {
+    using Self = GridTransfer;
+
+  public:
+    GridTransfer(Grid& grid)
+      : grid_(&grid)
+    {}
+
+    /// Attach a data container to the grid transfer, that gets interpolated during grid change
+    void attach(DOFVectorInterface* vec)
+    {
+      data_.push_back(vec);
+    }
+
+    /// Prepare the grid and the data for the adaption
+    bool preAdapt()
+    {
+      mightCoarsen_ = grid_->preAdapt(); // any element might be coarsened in adapt()
+      for (auto* vec : data_)
+        vec->preAdapt(mightCoarsen_);
+      return mightCoarsen_;
+    }
+
+    /// do the grid adaption
+    bool adapt()
+    {
+      refined_ = grid_->adapt(); // returns true if a least one entity was refined
+      return refined_;
+    }
+
+    // Perform data adaption to the new grid
+    void postAdapt()
+    {
+      if (mightCoarsen_ || refined_) {
+        for (auto* vec : data_)
+          vec->postAdapt(refined_);
+      }
+      grid_->postAdapt();
+    }
+
+  protected:
+    Grid* grid_;
+    std::vector<DOFVectorInterface*> data_;
+    bool mightCoarsen_ = false;
+    bool refined_ = false;
+  };
+
+} // end namespace AMDiS

src/amdis/ProblemStat.hpp

@@ -18,6 +18,7 @@
 #include <amdis/DirichletBC.hpp>
 //#include <amdis/Estimator.hpp>
 #include <amdis/Flag.hpp>
+#include <amdis/GridTransfer.hpp>
 #include <amdis/Initfile.hpp>
 #include <amdis/LinearAlgebra.hpp>
 #include <amdis/LinearSolvers.hpp>
@@ -346,6 +347,11 @@ namespace AMDiS
       setGrid(Dune::stackobject_to_shared_ptr(grid));
     }
 
+    void addAdaptionData(DOFVectorInterface* vec)
+    {
+      assert(bool(gridTransfer_));
+      gridTransfer_->attach(vec);
+    }
 
     void addMarker(std::shared_ptr<Marker<Grid>> const& marker)
     {
@@ -378,6 +384,7 @@ namespace AMDiS
     void adoptGrid(std::shared_ptr<Grid> const& grid)
     {
       grid_ = grid;
+      gridTransfer_ = std::make_shared<GridTransfer<Grid>>(*grid_);
       Parameters::get(name_ + "->mesh", gridName_);
     }
 
@@ -414,6 +421,9 @@ namespace AMDiS
     /// Grid of this problem.
     std::shared_ptr<Grid> grid_;
 
+    /// Handling the adaption of the grid
+    std::shared_ptr<GridTransfer<Grid>> gridTransfer_;
+
     /// Name of the grid
     std::string gridName_ = "mesh";
 

src/amdis/ProblemStat.inc.hpp

@@ -9,8 +9,8 @@
 
 #include <amdis/AdaptInfo.hpp>
 #include <amdis/FileWriter.hpp>
-#include <amdis/LocalAssembler.hpp>
 #include <amdis/GridFunctionOperator.hpp>
+#include <amdis/LocalAssembler.hpp>
 #include <amdis/common/Loops.hpp>
 
 namespace AMDiS {
@@ -36,7 +36,8 @@ void ProblemStat<Traits>::initialize(
         (adoptFlag.isSet(INIT_MESH) ||
         adoptFlag.isSet(INIT_SYSTEM) ||
         adoptFlag.isSet(INIT_FE_SPACE))) {
-      adoptGrid(adoptProblem->grid_);
+      grid_ = adoptProblem->grid_;
+      gridTransfer_ = adoptProblem->gridTransfer_;
     }
   }
 
@@ -121,6 +122,7 @@ void ProblemStat<Traits>::createGrid()
 {
   Parameters::get(name_ + "->mesh", gridName_);
   grid_ = MeshCreator<Grid>::create(gridName_);
+  gridTransfer_ = std::make_shared<GridTransfer<Grid>>(*grid_);
 
   msg("Create grid:");
   msg("#elements = {}"   , grid_->size(0));
@@ -169,8 +171,12 @@ template <class Traits>
 void ProblemStat<Traits>::createMatricesAndVectors()
 {
   systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
-  solution_ = std::make_shared<SystemVector>(*globalBasis_);
-  rhs_ = std::make_shared<SystemVector>(*globalBasis_);
+  solution_ = std::make_shared<SystemVector>(*globalBasis_, INTERPOLATE);
+  rhs_ = std::make_shared<SystemVector>(*globalBasis_, NO_OPERATION);
+
+  assert(bool(gridTransfer_));
+  gridTransfer_->attach(solution_.get());
+  gridTransfer_->attach(rhs_.get());
 
   auto localView = globalBasis_->localView();
   AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
@@ -326,19 +332,14 @@ adaptGrid(AdaptInfo& adaptInfo)
 {
   Dune::Timer t;
 
-  grid_->preAdapt();
-  bool changed = grid_->adapt();
-
-  // update data
-  if (changed) {
+  bool adapted = gridTransfer_->preAdapt();
+  adapted |= gridTransfer_->adapt(); // NOTE: |= does not short-circuit and works with bools as ||=
+  if (adapted)
     globalBasis_->update(gridView());
-    solution_->resize(*globalBasis_);
-  }
-
-  grid_->postAdapt();
+  gridTransfer_->postAdapt();
 
   msg("adaptGrid needed {} seconds", t.elapsed());
-  return changed ? MESH_ADAPTED : Flag(0);
+  return adapted ? MESH_ADAPTED : Flag(0);
 }
 
 
@@ -351,7 +352,6 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
   // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
   systemMatrix_->init(asmMatrix);
   rhs_->init(asmVector);
-  solution_->resize(*globalBasis_);
 
   auto localView = globalBasis_->localView();
   forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
@@ -377,7 +377,6 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
   // 3. finish matrix insertion and apply dirichlet boundary conditions
   systemMatrix_->finish(asmMatrix);
   rhs_->finish(asmVector);
-  (*solution_) = 0;
 
   forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
     forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {

src/amdis/linear_algebra/DOFVectorBase.hpp

@@ -4,6 +4,7 @@
 
 #include <dune/functions/functionspacebases/sizeinfo.hh>
 
+#include <amdis/DataTransfer.hpp>
 #include <amdis/LocalAssemblerList.hpp>
 #include <amdis/common/Math.hpp>
 #include <amdis/common/ScalarTypes.hpp>
@@ -38,10 +39,17 @@ namespace AMDiS
     using BaseVector = typename Backend::BaseVector;
     using ElementVector = Dune::DynamicVector<double>;
 
+    /// Defines an interface to transfer the data during grid adaption
+    using DataTransfer = DataTransferInterface<Self>;
+
+    /// A creator for a concrete data transfer object, depending on \ref DataTransferOperation
+    using DataTransferFactory = AMDiS::DataTransferFactory<Self>;
+
   public:
     /// Constructor. Constructs new BaseVector.
-    DOFVectorBase(BasisType const& basis)
+    DOFVectorBase(BasisType const& basis, DataTransferOperation op = NO_OPERATION)
       : basis_(&basis)
+      , dataTransfer_(DataTransferFactory::create(op, basis))
     {
       compress();
       operators_.init(basis);
@@ -142,6 +150,42 @@ namespace AMDiS
     /// Assemble all vector operators
     void assemble();
 
+    /// Return the associated DataTransfer object
+    DataTransfer const& dataTransfer() const
+    {
+      return *dataTransfer_;
+    }
+
+    /// Return the associated DataTransfer object
+    DataTransfer& dataTransfer()
+    {
+      return *dataTransfer_;
+    }
+
+    /// Create a new DataTransfer object based on the operation type
+    void setDataTransfer(DataTransferOperation op)
+    {
+      dataTransfer_ = DataTransferFactory::create(op, this->basis());
+    }
+
+    /// Assign the DataTransfer object
+    void setDataTransfer(std::shared_ptr<DataTransfer> const& dataTransfer)
+    {
+      dataTransfer_ = dataTransfer;
+    }
+
+    /// Implementation of \ref DOFVectorInterface::preAdapt
+    virtual void preAdapt(bool mightCoarsen) override
+    {
+      dataTransfer_->preAdapt(*this, mightCoarsen);
+    }
+
+    /// Implementation of \ref DOFVectorInterface::postAdapt
+    virtual void postAdapt(bool refined) override
+    {
+      dataTransfer_->postAdapt(*this, refined);
+    }
+
   private:
     /// The finite element space / basis associated with the data vector
     Basis const* basis_;
@@ -154,6 +198,9 @@ namespace AMDiS
 
     /// List of operators associated to nodes
     VectorOperators<Basis> operators_;
+
+    /// Data interpolation when the grid changes
+    std::shared_ptr<DataTransfer> dataTransfer_;
   };
 
 } // end namespace AMDiS

src/amdis/linear_algebra/DOFVectorInterface.hpp

@@ -10,6 +10,19 @@ namespace AMDiS
 
     /// Change dimension of DOFVector to dimension of basis
     virtual void compress() = 0;
+
+    /// \brief Prepare the data interpolation between grids
+    /**
+     * \param mightCoarsen  Indicator that any element in the grid might be
+     *                      coarsened in adapt()
+     **/
+    virtual void preAdapt(bool mightCoarsen) = 0;
+
+    /// \brief Apply the actual interpolation
+    /**
+     * \param refined  Indicator that at least one entity was refined
+     **/
+    virtual void postAdapt(bool refined) = 0;
   };
 
 } // end namespace AMDiS

src/amdis/linear_algebra/eigen/DOFVector.hpp

@@ -86,22 +86,25 @@ namespace AMDiS
 
 
   template <class BasisType, class ValueType = double>
-  struct DOFVector : public DOFVectorBase<BasisType, EigenVector<ValueType>>
+  class DOFVector : public DOFVectorBase<BasisType, EigenVector<ValueType>>
   {
     using Super = DOFVectorBase<BasisType, EigenVector<ValueType>>;
-    using Super::operator=;
-    DOFVector(BasisType const& basis)
-      : Super(basis)
+
+  public:
+    DOFVector(BasisType const& basis, DataTransferOperation op = NO_OPERATION)
+      : Super(basis, op)
     {}
+
+    using Super::operator=;
   };
 
 
   /// Constructor a dofvector from given basis and name
   template <class ValueType = double, class Basis>
   DOFVector<Basis, ValueType>
-  makeDOFVector(Basis const& basis)
+  makeDOFVector(Basis const& basis, DataTransferOperation op = NO_OPERATION)
   {
-    return {basis};
+    return {basis, op};
   }
 
 } // end namespace AMDiS

src/amdis/linear_algebra/istl/DOFVector.hpp

@@ -95,22 +95,25 @@ namespace AMDiS
 
 
   template <class BasisType, class ValueType = double>
-  struct DOFVector : public DOFVectorBase<BasisType, IstlVector<ValueType>>
+  class DOFVector : public DOFVectorBase<BasisType, IstlVector<ValueType>>
   {
     using Super = DOFVectorBase<BasisType, IstlVector<ValueType>>;
-    using Super::operator=;
-    DOFVector(BasisType const& basis)
-      : Super(basis)
+
+  public:
+    DOFVector(BasisType const& basis, DataTransferOperation op = NO_OPERATION)
+      : Super(basis, op)
     {}
+
+    using Super::operator=;
   };
 
 
   /// Constructor a dofvector from given basis and name
   template <class ValueType = double, class Basis>
   DOFVector<Basis, ValueType>
-  makeDOFVector(Basis const& basis)
+  makeDOFVector(Basis const& basis, DataTransferOperation op = NO_OPERATION)
   {
-    return {basis};
+    return {basis, op};
   }
 
 } // end namespace AMDiS

src/amdis/linear_algebra/mtl/DOFVector.hpp

@@ -86,22 +86,25 @@ namespace AMDiS
 
 
   template <class BasisType, class ValueType = double>
-  struct DOFVector : public DOFVectorBase<BasisType, MtlVector<ValueType>>
+  class DOFVector : public DOFVectorBase<BasisType, MtlVector<ValueType>>
   {
     using Super = DOFVectorBase<BasisType, MtlVector<ValueType>>;
-    using Super::operator=;
-    DOFVector(BasisType const& basis)
-      : Super(basis)
+
+  public:
+    DOFVector(BasisType const& basis, DataTransferOperation op = NO_OPERATION)
+      : Super(basis, op)
     {}
+
+    using Super::operator=;
   };
 
 
   /// Constructor a dofvector from given basis and name
   template <class ValueType = double, class Basis>
   DOFVector<Basis, ValueType>
-  makeDOFVector(Basis const& basis)
+  makeDOFVector(Basis const& basis, DataTransferOperation op = NO_OPERATION)
   {
-    return {basis};
+    return {basis, op};
   }
 
 } // end namespace AMDiS