From 5d525c212cb8b90002f5380dc85b4254770c46d5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Tue, 7 Jan 2020 18:19:43 +0100
Subject: [PATCH] [WIP] begin implementation of dg
---
dune/mmdg/clockhelper.hh | 20 +++++
src/CMakeLists.txt | 2 +
src/dg.cc | 170 ++++++++++++++++++++++++++++++++++-----
src/grids/cube1d.dgf | 10 +++
src/grids/cube2d.dgf | 10 +++
src/grids/cube3d.dgf | 10 +++
src/parameterDG.ini | 2 +
src/parameterMMDG.ini | 0
8 files changed, 206 insertions(+), 18 deletions(-)
create mode 100644 dune/mmdg/clockhelper.hh
create mode 100644 src/grids/cube1d.dgf
create mode 100644 src/grids/cube2d.dgf
create mode 100644 src/grids/cube3d.dgf
create mode 100644 src/parameterDG.ini
create mode 100644 src/parameterMMDG.ini
diff --git a/dune/mmdg/clockhelper.hh b/dune/mmdg/clockhelper.hh
new file mode 100644
index 0000000..05c3ddc
--- /dev/null
+++ b/dune/mmdg/clockhelper.hh
@@ -0,0 +1,20 @@
+#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
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a8348fd..5088d17 100644
--- a/src/CMakeLists.txt
+++ b/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)
diff --git a/src/dg.cc b/src/dg.cc
index 61ea9c4..559b344 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -1,30 +1,164 @@
-// -*- 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;
}
}
diff --git a/src/grids/cube1d.dgf b/src/grids/cube1d.dgf
new file mode 100644
index 0000000..10e02bb
--- /dev/null
+++ b/src/grids/cube1d.dgf
@@ -0,0 +1,10 @@
+DGF
+
+Interval
+0
+1
+5
+#
+
+Simplex
+#
diff --git a/src/grids/cube2d.dgf b/src/grids/cube2d.dgf
new file mode 100644
index 0000000..f1c490a
--- /dev/null
+++ b/src/grids/cube2d.dgf
@@ -0,0 +1,10 @@
+DGF
+
+Interval
+0 0
+1 1
+2 2
+#
+
+Simplex
+#
diff --git a/src/grids/cube3d.dgf b/src/grids/cube3d.dgf
new file mode 100644
index 0000000..d201d0a
--- /dev/null
+++ b/src/grids/cube3d.dgf
@@ -0,0 +1,10 @@
+DGF
+
+Interval
+0 0 0
+1 1 1
+20 20 20
+#
+
+Simplex
+#
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
new file mode 100644
index 0000000..4ec0b2d
--- /dev/null
+++ b/src/parameterDG.ini
@@ -0,0 +1,2 @@
+K = 1
+mu = 3
diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini
new file mode 100644
index 0000000..e69de29
--
GitLab