Restructuring of Matrix and Vector and linear-algebra backends.

Added init-finalize also for vectors.
Require gather/scatter operations instead of direct vector/matrix access.
Rewritten the interpolation.
Added ParallelIndexSet and DOFMapping to be used also in other backends.

bin/amdisproject

@@ -677,7 +677,7 @@ cat> "$PROJECT/init/$NAME.dat" <<EOF
 mesh->global refinements: 0
 
 $NAME->mesh: mesh
-$NAME->solver->name: default
+$NAME->solver: default
 $NAME->solver->relative tolerance: 1.e-6
 $NAME->solver->info: 1
 

examples/ellipt.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-#include <iostream>
-
 #include <amdis/AMDiS.hpp>
 #include <amdis/Integrate.hpp>
 #include <amdis/LocalOperators.hpp>
 #include <amdis/ProblemStat.hpp>
 #include <amdis/common/Literals.hpp>
 
+#include <dune/grid/yaspgrid.hh>
+using Grid = Dune::YaspGrid<GRIDDIM>;
+
 using namespace AMDiS;
 using namespace Dune::Indices;
 
-// 1 component with polynomial degree 1
-using Param   = YaspGridBasis<GRIDDIM, 2>;
-using ElliptProblem = ProblemStat<Param>;
-
 int main(int argc, char** argv)
 {
   Environment env(argc, argv);
+  Dune::Timer t;
 
-  int numLevels = GRIDDIM == 2 ? 8 : 5;
+  int numLevels = GRIDDIM == 2 ? 6 : 4;
   if (argc > 2)
     numLevels = std::atoi(argv[2]);
 
@@ -34,7 +29,8 @@ int main(int argc, char** argv)
     return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
   };
 
-  ElliptProblem prob("ellipt");
+  using Param = LagrangeBasis<Grid, 2>;
+  ProblemStat<Param> prob("ellipt");
   prob.initialize(INIT_ALL);
 
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
@@ -44,40 +40,44 @@ int main(int argc, char** argv)
   prob.addVectorOperator(opForce, 0);
 
   // set boundary condition
-  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
+  prob.addDirichletBC(1, 0, 0, g);
 
   AdaptInfo adaptInfo("adapt");
 
   std::vector<double> errL2; errL2.reserve(numLevels);
   std::vector<double> errH1; errH1.reserve(numLevels);
   std::vector<double> widths; widths.reserve(numLevels);
-  for (int i = 0; i < numLevels; ++i) {
-    prob.grid()->globalRefine(1);
-    auto gridView = prob.gridView();
+  for (int l = 0; l < numLevels; ++l) {
+    prob.globalRefine(1);
 
     double h = 0;
-    for (auto const& e : edges(gridView))
-      h = std::max(h, e.geometry().volume());
+    for (auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) {
+      auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type());
+      auto geo = e.geometry();
+      for (int i = 0; i < refElem.size(GRIDDIM-1); ++i) { // edges
+        auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM));
+        auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM));
+        h = std::max(h, (v0 - v1).two_norm());
+      }
+    }
     widths.push_back(h);
 
-    prob.globalBasis()->update(gridView);
-    prob.solutionVector()->compress();
     prob.assemble(adaptInfo);
     prob.solve(adaptInfo);
 
-    double errorL2 = integrate(sqr(g - prob.solution(0)), gridView, 6);
+    double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
     errL2.push_back(std::sqrt(errorL2));
-    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), gridView, 6);
+    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), prob.gridView(), 6);
     errH1.push_back(std::sqrt(errorH1));
 
 #if WRITE_FILES
     Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
     vtkWriter.addVertexData(prob.solution(0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
-    vtkWriter.write("u_" + std::to_string(i));
+    vtkWriter.write("u_" + std::to_string(l));
 #endif
   }
 
-  std::cout << "\n";
+  msg("");
   msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
       "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
   msg_("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
@@ -92,5 +92,6 @@ int main(int argc, char** argv)
     msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
         i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
 
+  msg("elapsed time: {} seconds", t.elapsed());
   return 0;
 }

examples/init/boundary.dat.2d

-elliptMesh->global refinements: 3
+mesh->global refinements: 3
 
-ellipt->mesh:                   elliptMesh
+ellipt->mesh:                   mesh
 
-ellipt->solver->name:           bicgstab
+ellipt->solver:                 bicgstab
 ellipt->solver->max iteration:  10000
 ellipt->solver->relative tolerance: 1.e-10
-ellipt->solver->info:           10
-ellipt->solver->left precon:    ilu
-ellipt->solver->right precon:   no
-ellipt->solver->precon->name:   ilu
+ellipt->solver->info:           1
+ellipt->solver->precon:         ilu
 
 ellipt->output[0]->filename:    boundary.2d
 ellipt->output[0]->name:        u

examples/init/cahn_hilliard.dat.2d

@@ -12,9 +12,9 @@ refinement->interface: 12
 refinement->bulk: 2
 
 ch->mesh:                       chMesh
-ch->solver->name:               direct
+ch->solver:                     direct
 ch->solver->max iteration:      1000
-ch->solver->absolute tolerance: 1e-6
+ch->solver->relative tolerance: 1e-6
 ch->solver->break if tolerance not reached: 1
 ch->solver->info: 2
 

examples/init/ellipt.dat.2d

 elliptMesh->global refinements: 0
+elliptMesh->num cells: 8 8
+elliptMesh->overlap: 1
 
-ellipt->mesh:                   elliptMesh
+ellipt->mesh: elliptMesh
+ellipt->symmetry: spd
 
-ellipt->solver->name:           bicgstab
-ellipt->solver->max iteration:  10000
+% ISTL backend parameters
+% =======================
+ellipt->solver: pcg
+ellipt->solver->info: 1
+ellipt->solver->max iteration: 10000
 ellipt->solver->relative tolerance: 1.e-10
-ellipt->solver->info:           10
-ellipt->solver->left precon:    ilu
-ellipt->solver->right precon:   no
-ellipt->solver->precon->name:   ilu
+ellipt->solver->precon: ilu
 
-ellipt->output[0]->filename:    ellipt.2d
-ellipt->output[0]->name:        u
+ellipt->output[0]->filename: ellipt.2d
+ellipt->output[0]->name: u
 ellipt->output[0]->output directory: ./output

examples/init/ellipt.dat.3d

 elliptMesh->global refinements: 0
+elliptMesh->num cells: 8 8 8
+elliptMesh->overlap: 1
 
-ellipt->mesh:                   elliptMesh
+ellipt->mesh: elliptMesh
+ellipt->symmetry: spd
 
-ellipt->solver->name:           bicgstab
-ellipt->solver->max iteration:  1000
-ellipt->solver->relative tolerance: 1.e-9
-ellipt->solver->info:           10
-ellipt->solver->left precon:    diag
-ellipt->solver->right precon:   no
+% ISTL backend parameters
+% =======================
+ellipt->solver: pcg
+ellipt->solver->info: 1
+ellipt->solver->max iteration: 10000
+ellipt->solver->relative tolerance: 1.e-10
+ellipt->solver->precon: ilu
 
-ellipt->output[0]->filename:    ellipt.3d
-ellipt->output[0]->name:        u
+ellipt->output[0]->filename: ellipt.3d
+ellipt->output[0]->name: u
 ellipt->output[0]->output directory: ./output

examples/init/ellipt.ini.2d

-dimension of world = 2
-
-[elliptMesh]
-macro file name = ./macro/macro.stand.2d
-global refinements = 5
-
-[ellipt]
-mesh = elliptMesh
-
-[ellipt.solver]
-name = cg
-max iteration = 1000
-absolute tolerance = 1e-6
-info = 1
-left precon = diag
-
-[ellipt.output[0]]
-filename = ellipt.2d
-name = u
-output directory = ./output

examples/init/heat.dat.2d

@@ -3,14 +3,14 @@ heatMesh->min corner: 0 0
 heatMesh->max corner: 2 2
 heatMesh->num cells: 8 8
 
-heat->mesh:  heatMesh
+heat->mesh: heatMesh
 heat->names: u
-heat->solver->name: cg
-heat->solver->max iteration:  1000
+
+heat->solver: pcg
+heat->solver->max iteration: 1000
 heat->solver->absolute tolerance: 1e-6
-heat->solver->break if tolerance not reached: 1
 heat->solver->info: 1
-heat->solver->left precon:    diag
+heat->solver->precon: ilu
 
 heat->output[0]->filename: heat.2d
 heat->output[0]->output directory: output

examples/init/marker.dat.2d

-dimension of world: 2
-
-test->mesh: testMesh
-testMesh->global refinements: 0

examples/init/stokes.dat.2d

 stokesMesh->global refinements: 0
 stokesMesh->max corner: 1.0 1.0
-stokesMesh->num cells: 4 4
+stokesMesh->num cells: 8 8
 
 stokes->mesh: stokesMesh
 
-stokes->solver->name: bicgstag
-stokes->solver->ell: 3
+stokes->solver: bicgstab
 stokes->solver->max iteration: 1000
 stokes->solver->relative tolerance: 1e-8
-stokes->solver->restart: 50
-stokes->solver->precon->name: ilu
+stokes->solver->precon: ilu
 stokes->solver->info:	2
 
 stokes->output[0]->filename: stokes_u.2d

examples/init/surface.dat.2d

@@ -3,11 +3,11 @@ surfaceMesh->global refinements: 1
 
 surface->mesh:               surfaceMesh
 
-surface->solver->name:           bicgstab
+surface->solver:                 bicgstab
 surface->solver->max iteration:  10000
 surface->solver->relative tolerance: 1.e-6
-surface->solver->info:           10
-surface->solver->left precon:    ilu
+surface->solver->info:           1
+surface->solver->precon:         ilu
 
 surface->output[0]->filename:    surface.2d
 surface->output[0]->name:        u

src/amdis/linearalgebra/DOFMatrixBase.hpp → src/amdis/BiLinearForm.hpp

@@ -2,12 +2,11 @@
 
 #include <cmath>
 
-#include <dune/common/timer.hh>
-
+#include <amdis/LinearAlgebra.hpp>
 #include <amdis/OperatorList.hpp>
 #include <amdis/common/FlatMatrix.hpp>
-#include <amdis/common/Math.hpp>
-#include <amdis/typetree/MultiIndex.hpp>
+#include <amdis/common/TypeTraits.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 #include <amdis/typetree/TreePath.hpp>
 
 namespace AMDiS
@@ -19,11 +18,17 @@ namespace AMDiS
    *
    * \tparam RB  Basis of the matrix rows
    * \tparam CB  Basis of matrix columns
-   * \tparam Backend  A linear-algebra backend for the matrix storage
+   * \tparam T   Coefficient type to store in the matrix
    **/
-  template <class RB, class CB, class Backend>
-  class DOFMatrixBase
+  template <class RB, class CB, class T = double>
+  class BiLinearForm
+      : public MatrixBase<RB,CB,MatrixBackend<BackendTraits<RB,T>>>
   {
+    using Self = BiLinearForm;
+    using Traits = BackendTraits<RB,T>;
+    using Backend = MatrixBackend<Traits>;
+    using Super = MatrixBase<RB,CB,Backend>;
+
   public:
     /// The type of the finite element space / basis of the row
     using RowBasis = RB;
@@ -33,88 +38,29 @@ namespace AMDiS
     using ColBasis = CB;
     using ColLocalView = typename ColBasis::LocalView;
 
-    using Element = typename RowLocalView::Element;
-    using Geometry = typename Element::Geometry;
-
-    /// The index/size - type
-    using size_type = typename RowBasis::size_type;
-    using value_type = typename Backend::value_type;
-
-    /// The type of the data matrix used in the backend
-    using BaseMatrix = typename Backend::BaseMatrix;
+    /// The type of the elements of the DOFVector
+    using CoefficientType = typename Traits::CoefficientType;
 
     /// The type of the matrix filled on an element with local contributions
-    using ElementMatrix = FlatMatrix<value_type>;
+    using ElementMatrix = FlatMatrix<CoefficientType>;
 
   public:
-    /// Constructor. Stores the shared_ptr to the bases.
-    DOFMatrixBase(std::shared_ptr<RowBasis> rowBasis, std::shared_ptr<ColBasis> colBasis)
-      : rowBasis_(rowBasis)
-      , colBasis_(colBasis)
-    {
-      operators_.init(*rowBasis_, *colBasis_);
-    }
-
     /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into a new shared_ptr.
-    template <class RB_, class CB_>
-    DOFMatrixBase(RB_&& rowBasis, CB_&& colBasis)
-      : DOFMatrixBase(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
-    {}
-
-    /// Return the row-basis \ref rowBasis of the matrix
-    std::shared_ptr<RowBasis const> rowBasis() const
-    {
-      return rowBasis_;
-    }
-
-    /// Return the col-basis \ref colBasis of the matrix
-    std::shared_ptr<ColBasis const> colBasis() const
-    {
-      return colBasis_;
-    }
-
-    /// Return the data-matrix
-    BaseMatrix const& matrix() const
-    {
-      return backend_.matrix();
-    }
-
-    /// Return the data-matrix
-    BaseMatrix& matrix()
+    template <class RB_, class CB_, class Comm_>
+    BiLinearForm(RB_&& rowBasis, CB_&& colBasis, Comm_&& comm)
+      : Super(FWD(rowBasis), FWD(colBasis), FWD(comm))
     {
-      return backend_.matrix();
+      operators_.init(*Super::rowBasis(), *Super::colBasis());
     }
 
-    /// Return the size of the \ref rowBasis_
-    size_type rows() const
+    /// Prepare the matrix for insertion. Clears all the entries.
+    void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
     {
-      return rowBasis_->dimension();
-    }
+      Super::init(symmetry);
 
-    /// Return the size of the \ref colBasis_
-    size_type cols() const
-    {
-      return colBasis_->dimension();
-    }
-
-    /// Initialize the matrix for insertion, e.g. allocate the non-zero pattern
-    /// If \p setToZero is true, the matrix is set to 0
-    void init(bool asmMatrix);
-
-    /// Finish the matrix insertion, e.g. cleanup or final insertion
-    void finish(bool asmMatrix);
-
-    /// Insert a block of values into the matrix (add to existing values)
-    /// The global matrix indices are determined by the corresponding localviews.
-    void insert(RowLocalView const& rowLocalView,
-                ColLocalView const& colLocalView,
-                ElementMatrix const& elementMatrix);
-
-    /// Insert a single value into the matrix (add to existing value)
-    template <class RowIndex, class ColIndex>
-    void insert(RowIndex row, ColIndex col, typename Backend::value_type const& value)
-    {
-      backend_.insert(flatMultiIndex(row), flatMultiIndex(col), value);
+      auto const rowSize = this->rowBasis()->localView().maxSize();
+      auto const colSize = this->colBasis()->localView().maxSize();
+      elementMatrix_.resize(rowSize, colSize);
     }
 
     /// \brief Associate a local operator with this DOFMatrix
@@ -146,23 +92,9 @@ namespace AMDiS
                   ColLocalView const& colLocalView);
 
     /// Assemble all matrix operators, TODO: incooperate boundary conditions
-    void assemble();
-
-    /// Number of nonzeros in the matrix
-    size_type nnz() const
-    {
-      return backend_.nnz();
-    }
+    void assemble(SymmetryStructure symmetry = SymmetryStructure::unknown);
 
   protected:
-    /// The finite element space / basis associated with the rows
-    std::shared_ptr<RowBasis> rowBasis_;
-
-    /// The finite element space / basis associated with the columns
-    std::shared_ptr<ColBasis> colBasis_;
-
-    /// Data backend
-    Backend backend_;
 
     /// Dense matrix to store coefficients during \ref assemble()
     ElementMatrix elementMatrix_;
@@ -173,4 +105,4 @@ namespace AMDiS
 
 } // end namespace AMDiS
 
-#include "DOFMatrixBase.inc.hpp"
+#include <amdis/BiLinearForm.inc.hpp>

src/amdis/linearalgebra/DOFMatrixBase.inc.hpp → src/amdis/BiLinearForm.inc.hpp

 #pragma once
 
+#include <utility>
+
 #include <amdis/Assembler.hpp>
 #include <amdis/LocalOperator.hpp>
 #include <amdis/typetree/Traversal.hpp>
@@ -7,47 +9,9 @@
 
 namespace AMDiS {
 
-template <class RB, class CB, class B>
-void DOFMatrixBase<RB,CB,B>::
-init(bool asmMatrix)
-{
-  backend_.init(*rowBasis_, *colBasis_, asmMatrix);
-
-  auto const rowSize = rowBasis_->localView().maxSize();
-  auto const colSize = colBasis_->localView().maxSize();
-  elementMatrix_.resize(rowSize, colSize);
-}
-
-
-template <class RB, class CB, class B>
-void DOFMatrixBase<RB,CB,B>::
-finish(bool /*asmMatrix*/)
-{
-  backend_.finish();
-}
-
-
-template <class RB, class CB, class B>
-void DOFMatrixBase<RB,CB,B>::
-insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
-       ElementMatrix const& elementMatrix)
-{
-  using std::abs;
-  for (size_type i = 0; i < rowLocalView.size(); ++i) {
-    size_type const row = flatMultiIndex(rowLocalView.index(i));
-    for (size_type j = 0; j < colLocalView.size(); ++j) {
-      if (abs(elementMatrix[i][j]) > threshold<double>) {
-        size_type const col = flatMultiIndex(colLocalView.index(j));
-        backend_.insert(row, col, elementMatrix[i][j]);
-      }
-    }
-  }
-}
-
-
-template <class RB, class CB, class B>
+template <class RB, class CB, class T>
   template <class ContextTag, class Expr, class RowTreePath, class ColTreePath>
-void DOFMatrixBase<RB,CB,B>::
+void BiLinearForm<RB,CB,T>::
 addOperator(ContextTag contextTag, Expr const& expr,
             RowTreePath row, ColTreePath col)
 {
@@ -56,24 +20,26 @@ addOperator(ContextTag contextTag, Expr const& expr,
   static_assert( Concepts::PreTreePath<ColTreePath>,
       "col must be a valid treepath, or an integer/index-constant");
 
-  auto i = child(rowBasis_->localView().tree(), makeTreePath(row));
-  auto j = child(colBasis_->localView().tree(), makeTreePath(col));
+  auto i = child(this->rowBasis()->localView().tree(), makeTreePath(row));
+  auto j = child(this->colBasis()->localView().tree(), makeTreePath(col));
 
   using LocalContext = typename ContextTag::type;
   using Traits = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
-  auto op = makeLocalOperator<LocalContext>(expr, rowBasis_->gridView());
+  auto op = makeLocalOperator<LocalContext>(expr, this->rowBasis()->gridView());
   auto localAssembler = makeUniquePtr(makeAssembler<Traits>(std::move(op), i, j));
 
   operators_[i][j].push(contextTag, std::move(localAssembler));
 }
 
 
-template <class RB, class CB, class B>
-void DOFMatrixBase<RB,CB,B>::
+template <class RB, class CB, class T>
+void BiLinearForm<RB,CB,T>::
 assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
 {
+  elementMatrix_.resize(rowLocalView.size(), colLocalView.size());
   elementMatrix_ = 0;
-  auto const& gv = rowBasis_->gridView();
+
+  auto const& gv = this->rowBasis()->gridView();
   auto const& element = rowLocalView.element();
   auto geometry = element.geometry();
 
@@ -88,30 +54,30 @@ assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
     });
   });
 
-  insert(rowLocalView, colLocalView, elementMatrix_);
+  this->scatter(rowLocalView, colLocalView, elementMatrix_);
 }
 
 
-template <class RB, class CB, class B>
-void DOFMatrixBase<RB,CB,B>::
-assemble()
+template <class RB, class CB, class T>
+void BiLinearForm<RB,CB,T>::
+assemble(SymmetryStructure symmetry)
 {
-  auto rowLocalView = rowBasis_->localView();
-  auto colLocalView = colBasis_->localView();
+  auto rowLocalView = this->rowBasis()->localView();
+  auto colLocalView = this->colBasis()->localView();
 
-  init(true);
-  for (auto const& element : elements(rowBasis_->gridView())) {
+  this->init(symmetry);
+  for (auto const& element : elements(this->rowBasis()->gridView(), typename Traits::PartitionSet{})) {
     rowLocalView.bind(element);
-    if (rowBasis_ == colBasis_)
-      assemble(rowLocalView, rowLocalView);
+    if (this->rowBasis() == this->colBasis())
+      this->assemble(rowLocalView, rowLocalView);
     else {
       colLocalView.bind(element);
-      assemble(rowLocalView, colLocalView);
+      this->assemble(rowLocalView, colLocalView);
       colLocalView.unbind();
     }
     rowLocalView.unbind(element);
   }
-  finish(true);
+  this->finish();
 }
 
 

src/amdis/CMakeLists.txt

@@ -23,6 +23,8 @@ install(FILES
     Assembler.hpp
     AssemblerInterface.hpp
     BackupRestore.hpp
+    BiLinearForm.hpp
+    BiLinearForm.inc.hpp
     Boundary.hpp
     BoundaryCondition.hpp
     BoundaryManager.hpp
@@ -33,6 +35,9 @@ install(FILES
     DataTransfer.inc.hpp
     DirichletBC.hpp
     DirichletBC.inc.hpp
+    DOFVector.hpp
+    DOFVector.inc.hpp
+    DOFVectorInterface.hpp
     Environment.hpp
     FileWriter.hpp
     FileWriterInterface.hpp
@@ -45,6 +50,8 @@ install(FILES
     InitfileParser.hpp
     Integrate.hpp
     LinearAlgebra.hpp
+    LinearForm.hpp
+    LinearForm.inc.hpp
     LocalOperator.hpp
     LocalOperators.hpp
     Marker.hpp

src/amdis/linearalgebra/DOFVectorBase.hpp → src/amdis/DOFVector.hpp

 #pragma once