ellipt demo with convergence study added

.gitmodules

+[submodule "externals/fmt"]
+	path = externals/fmt
+	url = https://github.com/fmtlib/fmt.git

CMakeLists.txt

@@ -19,6 +19,7 @@ add_subdirectory("test")
 add_subdirectory("examples" EXCLUDE_FROM_ALL)
 add_subdirectory("doc")
 add_subdirectory("cmake/modules")
+add_subdirectory("externals" EXCLUDE_FROM_ALL)
 
 target_include_directories(amdis PUBLIC
   $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src>)

examples/CMakeLists.txt

@@ -6,7 +6,7 @@ foreach(project ${projects2d})
     add_executable(${project}.2d ${project}.cc)
     add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project}.2d)
     target_link_dune_default_libraries(${project}.2d)
-    target_link_libraries(${project}.2d "amdis")
+    target_link_libraries(${project}.2d amdis fmt)
     target_compile_definitions(${project}.2d PRIVATE AMDIS_DIM=2 AMDIS_DOW=2)
     add_dependencies(examples ${project}.2d)
 endforeach()
@@ -17,7 +17,7 @@ foreach(project ${projects3d})
     add_executable(${project}.3d ${project}.cc)
     add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 ${project}.3d)
     target_link_dune_default_libraries(${project}.3d)
-    target_link_libraries(${project}.3d "amdis")
+    target_link_libraries(${project}.3d amdis fmt)
     target_compile_definitions(${project}.3d PRIVATE AMDIS_DIM=3 AMDIS_DOW=3)
     add_dependencies(examples ${project}.3d)
 endforeach()
\ No newline at end of file

examples/ellipt.cc

@@ -3,55 +3,95 @@
 
 #include <iostream>
 
+#include <fmt/core.h>
+
 #include <amdis/AMDiS.hpp>
 #include <amdis/ProblemStat.hpp>
 #include <amdis/Operators.hpp>
 #include <amdis/common/Literals.hpp>
+#include <amdis/gridfunctions/Integrate.hpp>
 
 using namespace AMDiS;
+using namespace Dune::Indices;
 
 // 1 component with polynomial degree 1
-//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
-using ElliptParam   = YaspGridBasis<AMDIS_DIM, 1>;
-using ElliptProblem = ProblemStat<ElliptParam>;
+using Param   = YaspGridBasis<AMDIS_DIM, 1>;
+using ElliptProblem = ProblemStat<Param>;
 
 int main(int argc, char** argv)
 {
   AMDiS::init(argc, argv);
 
-  using namespace Dune::Indices;
+  int numLevels = AMDIS_DIM == 2 ? 8 : 5;
+  if (argc > 2)
+    numLevels = std::atoi(argv[2]);
+
+  auto f = [](auto const& x) {
+    double r2 = dot(x,x);
+    double ux = std::exp(-10.0 * r2);
+    return -(400.0 * r2 - 20.0 * x.size()) * ux;
+  };
+  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
+  auto grad_g = [](auto const& x) -> FieldMatrix<double,1,AMDIS_DIM> {
+    return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
+  };
 
   ElliptProblem prob("ellipt");
   prob.initialize(INIT_ALL);
 
-  AdaptInfo adaptInfo("adapt");
-
   auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
   prob.addMatrixOperator(opL, _0, _0);
 
-  auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
+  auto opForce = makeOperator(tag::test{}, f, 6);
   prob.addVectorOperator(opForce, _0);
 
-  auto opForce2 = makeOperator(tag::test{}, [](auto const& x) { return -2.0; }, 0);
-  prob.addVectorOperator(BoundaryType{0}, opForce2, _0);
-
   // set boundary condition
-  auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
-  auto dbcValues = [](auto const& x){ return 0.0; }; // set value
-  prob.addDirichletBC(predicate, _0, _0, dbcValues);
+  auto boundary = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };
+  prob.addDirichletBC(boundary, _0, _0, g);
+
+  AdaptInfo adaptInfo("adapt");
 
-  prob.buildAfterCoarsen(adaptInfo, Flag(0));
+  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.getGrid()->globalRefine(1);
+    auto gridView = prob.getGlobalBasis()->gridView();
+
+    double h = 0;
+    for (auto const& e : edges(gridView))
+      h = std::max(h, e.geometry().volume());
+    widths.push_back(h);
+
+    prob.buildAfterCoarsen(adaptInfo, Flag(0));
+    prob.solve(adaptInfo);
+
+    double errorL2 = integrate(sqr(g - prob.getSolution(_0)), gridView, 6);
+    errL2.push_back(std::sqrt(errorL2));
+    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.getSolution(_0))), gridView, 6);
+    errH1.push_back(std::sqrt(errorH1));
+
+#if WRITE_FILES
+    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(prob.getGlobalBasis()->gridView());
+    vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
+    vtkWriter.write("u_" + std::to_string(i));
+#endif
+  }
 
-  // write matrix to file
-  if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 5) {
-    mtl::io::matrix_market_ostream out("matrix.mtx");
-    out << prob.getSystemMatrix()->getMatrix();
+  std::cout << "\n";
+  fmt::print("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}\n",
+             "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
+  fmt::print("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
 
-    std::cout << prob.getSystemMatrix()->getMatrix() << '\n';
+  std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
+  for (int i = 1; i < numLevels; ++i) {
+    eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]);
+    eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]);
   }
 
-  prob.solve(adaptInfo);
-  prob.writeFiles(adaptInfo, true);
+  for (int i = 0; i < numLevels; ++i)
+    fmt::print("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}\n",
+                i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
 
   AMDiS::finalize();
   return 0;

examples/init/ellipt.dat.2d

-dimension of world:             2
-
-elliptMesh->macro file name:    ./macro/macro.stand.2d
-elliptMesh->global refinements: 5
+elliptMesh->global refinements: 0
 
 ellipt->mesh:                   elliptMesh
 
-ellipt->solver->name:           cg
+ellipt->solver->name:           bicgstab
 ellipt->solver->max iteration:  1000
-ellipt->solver->tolerance:      1.e-8
+ellipt->solver->relative tolerance: 1.e-8
 ellipt->solver->info:           10
 ellipt->solver->left precon:    diag
 ellipt->solver->right precon:   no

examples/init/ellipt.dat.3d

+elliptMesh->global refinements: 0
+
+ellipt->mesh:                   elliptMesh
+
+ellipt->solver->name:           bicgstab
+ellipt->solver->max iteration:  1000
+ellipt->solver->relative tolerance: 1.e-8
+ellipt->solver->info:           10
+ellipt->solver->left precon:    diag
+ellipt->solver->right precon:   no
+
+ellipt->output[0]->filename:    ellipt.3d
+ellipt->output[0]->name:        u
+ellipt->output[0]->output directory: ./output

externals/CMakeLists.txt

+add_subdirectory(fmt)
\ No newline at end of file

fmt @ c6d9730d

+Subproject commit c6d9730ddbfaa5dbefd5d2b48ef047ad362157f4

src/amdis/DirichletBC.hpp

@@ -80,6 +80,7 @@ namespace AMDiS
       // subtract columns of dirichlet nodes from rhs
       // for (auto const& triplet : columns)
       //   rhs[triplet.row] -= triplet.value * solution[triplet.col];
+      initialized_ = false;
     }
 
   protected:

src/amdis/GridFunctionOperator.hpp

@@ -61,24 +61,9 @@ namespace AMDiS
      **/
     GridFunctionOperatorBase(GridFunction const& gridFct, int termOrder)
       : gridFct_(gridFct)
-      , localFct_(localFunction(gridFct_))
       , termOrder_(termOrder)
     {}
 
-    // Copy constructor
-    GridFunctionOperatorBase(GridFunctionOperatorBase const& that)
-      : gridFct_(that.gridFct_)
-      , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
-      , termOrder_(that.termOrder_)
-    {}
-
-    // Move constructor
-    GridFunctionOperatorBase(GridFunctionOperatorBase&& that)
-      : gridFct_(std::move(that.gridFct_))
-      , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
-      , termOrder_(std::move(that.termOrder_))
-    {}
-
     /// \brief Binds operator to `element` and `geometry`.
     /**
      * Binding an operator to the currently visited element in grid traversal.
@@ -91,21 +76,24 @@ namespace AMDiS
     void bind_impl(Element const& element, Geometry const& geometry)
     {
       assert( bool(quadFactory_) );
-      localFct_.bind(element);
-      quadFactory_->bind(localFct_);
+      localFct_.emplace(localFunction(gridFct_));
+      localFct_->bind(element);
+      quadFactory_->bind(localFct_.value());
     }
 
     /// Unbinds operator from element.
     void unbind_impl()
     {
-      localFct_.unbind();
+      localFct_->unbind();
+      localFct_.reset();
     }
 
     /// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
     template <class PreQuadFactory>
     void setQuadFactory(PreQuadFactory const& pre)
     {
-      quadFactory_ = makeQuadratureFactoryPtr<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>(pre);
+      using ctype = typename Geometry::ctype;
+      quadFactory_.emplace(makeQuadratureFactory<ctype, LocalContext::mydimension, LocalFunction>(pre));
     }
 
   protected:
@@ -113,7 +101,7 @@ namespace AMDiS
     auto coefficient(LocalCoordinate const& local) const
     {
       assert( this->bound_ );
-      return localFct_(local);
+      return (*localFct_)(local);
     }
 
     /// Create a quadrature rule using the \ref QuadratureFactory by calculating the
@@ -131,10 +119,10 @@ namespace AMDiS
     GridFunction gridFct_;
 
     /// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
-    LocalFunction localFct_;
+    Dune::Std::optional<LocalFunction> localFct_;
 
     /// Assign each element type a quadrature rule
-    std::shared_ptr<QuadFactory> quadFactory_ = std::make_shared<QuadFactory>();
+    Dune::Std::optional<QuadFactory> quadFactory_;
 
     /// the derivative order of this operator (in {0, 1, 2})
     int termOrder_ = 0;

src/amdis/common/Math.hpp

@@ -118,7 +118,7 @@ namespace AMDiS
   }
 
   template <class T>
-  constexpr T threshold = T(1.e-16); //Math::sqr(std::numeric_limits<T>::epsilon());
+  constexpr T threshold = T(1.e-18L); //Math::sqr(std::numeric_limits<T>::epsilon());
 
 
   /// Calculates factorial of i

src/amdis/linear_algebra/LinearSolverInterface.hpp

@@ -45,6 +45,7 @@ namespace AMDiS
     void solve(Matrix const& A, VectorX& x, VectorB const& b,
                SolverInfo& solverInfo)
     {
+      assert( num_cols(A) == size(x) && num_rows(A) == size(b) );
       solveImpl(A, x, b, solverInfo);
     }
 

src/amdis/linear_algebra/mtl/DOFMatrix.hpp

@@ -146,7 +146,9 @@ namespace AMDiS
       test_exit(!insertionMode && !inserter, "Matrix already in insertion mode!");
 
       calculateSlotSize();
-      matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
+      matrix->change_dim(rowFeSpace.dimension(), colFeSpace.dimension());
+      test_exit(num_rows(*matrix) == rowFeSpace.dimension() && num_cols(*matrix) == colFeSpace.dimension(),
+        "Wrong dimensions in matrix");
       if (prepareForInsertion) {
         set_to_zero(*matrix);
         inserter = new Inserter(*matrix, slot_size);
@@ -281,7 +283,7 @@ namespace AMDiS
 
   private:
     /// The minimal number of nonzeros per row
-    static constexpr int MIN_NNZ_PER_ROW = 20;
+    static constexpr int MIN_NNZ_PER_ROW = 50;
 
     /// The finite element space / basis of the row
     RowFeSpace const& rowFeSpace;
@@ -302,7 +304,7 @@ namespace AMDiS
     bool insertionMode = false;
 
     /// The size of the slots used to insert values per row
-    int slot_size = 20;
+    int slot_size = MIN_NNZ_PER_ROW;
 
     // friend class declarations
     template <class, class> friend class SystemMatrix;

src/amdis/linear_algebra/mtl/DOFVector.hpp

@@ -45,7 +45,6 @@ namespace AMDiS
       , vector(ClonablePtr<BaseVector>::make())
     {
       compress();
-      *vector = ValueType(0);
     }
 
     /// Constructor. Takes reference to existing BaseVector
@@ -89,8 +88,10 @@ namespace AMDiS
     /// Resize the \ref vector to the size of the \ref feSpace.
     void compress()
     {
-      if (num_rows(*vector) != getSize())
+      if (num_rows(*vector) != getSize()) {
         resize(getSize());
+        *vector = ValueType(0);
+      }
     }
 
     template <class SizeInfo>