[WIP] begin implementation of dg

dune/mmdg/clockhelper.hh

+#ifndef __DUNE_MMDG_CLOCKHELPER_HH__
+#define __DUNE_MMDG_CLOCKHELPER_HH__
+
+#include <chrono>
+
+class ClockHelper
+{
+  public:
+    using Clock = std::chrono::steady_clock;
+
+    //returns the elapsed run time since the time point begin in seconds with a
+    //precision of a microsecond
+    static double elapsedRuntime(const Clock::time_point begin)
+    {
+      return 1e-6*std::chrono::duration_cast<std::chrono::microseconds>
+        (Clock::now() - begin).count();;
+    }
+};
+
+#endif

src/CMakeLists.txt

@@ -9,3 +9,5 @@ target_link_dune_default_libraries("dg-2d")
 add_executable("mmdg-2d" mmdg.cc)
 target_compile_definitions("mmdg-2d" PRIVATE GRIDDIM=2)
 target_link_dune_default_libraries("mmdg-2d")
+
+dune_symlink_to_source_files(FILES grids parameterMMDG.ini parameterDG.ini)

src/dg.cc

-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
 #ifdef HAVE_CONFIG_H
 # include "config.h"
 #endif
-#include <iostream>
-#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
-#include <dune/common/exceptions.hh> // We use exceptions
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/mmdg/clockhelper.hh>
 
 int main(int argc, char** argv)
 {
-  try{
-    // Maybe initialize MPI
-    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
-    std::cout << "Hello World! This is dune-mmdg." << std::endl;
-    if(Dune::MPIHelper::isFake)
-      std::cout<< "This is a sequential program." << std::endl;
-    else
-      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
-        <<" processes!"<<std::endl;
+  try
+  {
+    using Clock = std::chrono::steady_clock;
+
+    //start run time measurement
+    const Clock::time_point timeInitial = Clock::now();
+
+    static constexpr int dim = GRIDDIM;
+
+    using Grid = Dune::YaspGrid<dim>;
+    using GridFactory = Dune::DGFGridFactory<Grid>;
+    using GridView = typename Grid::LeafGridView;
+    using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+    using Coordinate = Dune::FieldVector<double, dim>;
+
+    //Message Passing Interface (MPI) for parallel computation
+    Dune::MPIHelper::instance(argc, argv);
+
+    //get parameters from file "parameterDG.ini"
+    Dune::ParameterTree pt;
+    Dune::ParameterTreeParser::readINITree("parameterDG.ini", pt);
+
+    //assume isotropic medium such that permeability tensor reduces to a scalar
+    const double K = std::stod(pt["K"]);
+
+    //assume constant penalty parameter
+    const double mu = std::stod(pt["mu"]);
+
+    //create a grid from .dgf file
+    GridFactory gridFactory( "grids/cube" + std::to_string(dim) +
+      "d.dgf" );
+    const Grid& grid = *gridFactory.grid();
+    const GridView& gridView = grid.leafGridView();
+
+    //element mapper for the leaf grid view
+    const Mapper mapper(gridView, Dune::mcmgElementLayout());
+
+    const int dof = (1 + dim)*gridView.size(0); //degrees of freedom
+    using Matrix = Dune::DynamicMatrix<double>; //Dune::FieldMatrix<double, dof, dof>;
+    using Vector = Dune::DynamicVector<double>; //Dune::FieldVector<double, dof>;
+    Matrix A(dof, dof, 0.0); //stiffness matrix
+    Vector b(dof, 0.0);
+    Vector d;
+
+    const auto q = [](const Coordinate& x) {return x.two_norm();}; //source term
+
+    const int order = 3; //order of the quadrature rule
+    for (const auto& elem : elements(gridView))
+    {
+      const auto& geo = elem.geometry();
+      const int elemIdx = mapper.index(elem);
+      //TODO: can be done outside of the loop?
+      const Dune::QuadratureRule<double,dim>& rule =
+        Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
+
+      for (const auto& ip : rule)
+      {
+        const auto qp = ip.position();
+        //NOTE: is the volume of the element taken into account automatically?
+        const auto weight = ip.weight();
+
+        //basis function phi_{elem,0}(x) = indicator(elem);
+        b[(1 + dim)*elemIdx] += weight * q(qp);
+        for (int i = 0; i < dim; i++)
+        {
+          //basis function phi_{elem,i}(x) = x[i]*indicator(elem);
+          b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp);
+        }
+      }
+
+      for (int i = 0; i < dim; i++)
+      {
+        //basis function phi_{elem,i}(x) = x[i]*indicator(elem);
+        A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + i + 1] =
+          K * geo.volume();  // \int_elem K\nabla\phi_i\cdot\nabla\phi_i dV
+      }
+
+      //NOTE: Is there a neat way to evaluate the jump and average operators?
+      // (loop over facets?)
+      //NOTE: jump / average operators disappear as basis function are supported
+      //on a single element?
+
+      //iterate over all intersection with the boundary of elem
+      for (const auto& intersct : intersections(gridView, elem))
+      {
+        const auto& outerNormal = intersct.centerUnitOuterNormal();
+        const auto& intersctGeo = intersct.geometry();
+
+        const Dune::QuadratureRule<double,dim-1>& rule =
+          Dune::QuadratureRules<double,dim-1>::rule(intersct.type(),order);
+
+        A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] += mu * intersctGeo.volume();
+
+        for (const auto& ip : rule)
+        {
+          const auto qp = ip.position();
+          const auto weight = ip.weight();
+
+          //TODO
+        }
+
+
+
+        if (intersct.neighbor()) //intersct has neighboring element
+        {
+          //index of the neighboring element
+          const int neighborIndex = mapper.index(intersct.outside());
+
+          // A[][] = 0.5* ...; //note that interior facets are considered twice
+
+          //TODO
+        }
+        else //boundary element
+        {
+            //TODO
+        }
+      }
+    }
+
+    A.solve(d, b);
+
+    //NOTE: what is a suitable way to create the output?
+
+    //print total run time
+    const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
+    std::cout << "Finished with a total run time of " << timeTotal << ".\n";
+
     return 0;
   }
-  catch (Dune::Exception &e){
-    std::cerr << "Dune reported error: " << e << std::endl;
+  catch (Dune::Exception& e)
+  {
+    std::cerr << "Dune reported error: " << e.what() << std::endl;
+  }
+  catch (std::exception& e)
+  {
+    std::cerr << "Error: " << e.what() << std::endl;
   }
-  catch (...){
+  catch (...)
+  {
     std::cerr << "Unknown exception thrown!" << std::endl;
   }
 }

src/grids/cube1d.dgf

+DGF
+
+Interval
+0	
+1	
+5	
+#
+
+Simplex
+#

src/grids/cube2d.dgf

+DGF
+
+Interval
+0	0	
+1	1	
+2	2	
+#
+
+Simplex
+#

src/grids/cube3d.dgf

+DGF
+
+Interval
+0	0	0	
+1	1	1	
+20	20	20
+#
+
+Simplex
+#

src/parameterDG.ini

+K = 1
+mu = 3