diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 35acb2579cd8accc343403988b1fb2996dc4cb6a..83676e7aab71bd914fabbcf736bcde04119aaf6c 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -6,6 +6,7 @@
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
+#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/geometry/quadraturerules.hh>
@@ -16,10 +17,12 @@
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
-template<class GridView, class Mapper, class Problem>
+template<class Grid, class GridView, class Mapper, class Problem>
class DG
-{
+{ //TODO: Ordnung (private / public)
public:
+ //using GridView = typename Grid::LeafGridView;
+ //using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
@@ -35,21 +38,23 @@ public:
static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
//constructor
- DG (const GridView& gridView, const Mapper& mapper,
- const Problem& problem) :
- gridView_(gridView), mapper_(mapper), problem_(problem),
- dof((1 + dim) * gridView.size(0))
+ DG (const Grid& grid, const GridView& gridView, const Mapper mapper, const Problem& problem) :
+ //NOTE: why does this not work?
+ //grid_(grid), gridView_(grid.leafGridView()), mapper_(*(new Mapper(gridView_, Dune::mcmgElementLayout()))), problem_(problem), dof_((1 + dim) * gridView_.size(0))
+ grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem),
+ dof_((1 + dim) * gridView.size(0))
{
//initialize stiffness matrix A, load vector b and solution vector d
- A = std::make_shared<Matrix>(dof, dof, 10, 1, Matrix::implicit); //TODO
- b = Vector(dof);
- d = Vector(dof);
+ A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+ b = Vector(dof_);
+ d = Vector(dof_);
}
const void operator() (const Scalar mu)
{
//assemble stiffness matrix A and load vector b
assembleSLE(mu);
+ A->compress();
Dune::InverseOperatorResult result;
Dune::UMFPack<Matrix> solver(*A);
@@ -59,13 +64,25 @@ public:
writeVTKoutput();
}
+ const Grid& grid_;
const GridView& gridView_;
const Mapper& mapper_; //element mapper for gridView_
const Problem& problem_; //the DG problem to be solved
- const int dof; //degrees of freedom
+ const int dof_;
+
+protected:
+ DG (const Grid& grid, const GridView& gridView, const Mapper mapper,
+ const Problem& problem, const int dof) :
+ grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem),
+ dof_(dof)
+ {
+ //initialize stiffness matrix A, load vector b and solution vector d
+ A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+ b = Vector(dof_);
+ d = Vector(dof_);
+ }
-private:
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar mu)
{
@@ -145,6 +162,14 @@ private:
//iterate over all intersection with the boundary of elem
for (const auto& intersct : intersections(gridView_, elem))
{
+ if constexpr (dim != 1)
+ {
+ if (grid_.isInterface(intersct))
+ {
+ continue;
+ }
+ }
+
const auto& normal = intersct.centerUnitOuterNormal();
const auto& intersctGeo = intersct.geometry();
const double intersctVol = intersctGeo.volume();
@@ -377,7 +402,7 @@ private:
// int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
//with K*grad(phi_elem,i) = K[:,i]
for (int i = 0; i < dim; i++)
- {
+ {
b[elemIdxSLE + i + 1] +=
weight * integrationElem * (mu * qpGlobal[i] -
(KEvaluation[i] * normal)) * boundaryEvaluation;
@@ -423,8 +448,6 @@ private:
}
}
}
-
- A->compress();
}
//writes the solution to a vtk file
@@ -433,8 +456,8 @@ private:
//storage for pressure data, for each element we store the
//pressure at the corners of the element, the VTKFunction will
//create the output using a linear interpolation for each element
- Dune::DynamicVector<Scalar> pressure(dof, 0.0);
- Dune::DynamicVector<Scalar> exactPressure(dof, 0.0);
+ Dune::DynamicVector<Scalar> pressure(dof_, 0.0);
+ Dune::DynamicVector<Scalar> exactPressure(dof_, 0.0);
for (const auto& elem : elements(gridView_))
{
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index e5f79b72ba5f852fc346544895e6236b590bea62..be1b24ac60acd88c0f047352fb0cc2b701049a16 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -1,6 +1,177 @@
#ifndef DUNE_MMDG_HH
#define DUNE_MMDG_HH
-// add your classes here
+#include <dune/mmdg/dg.hh>
-#endif // DUNE_MMDG_HH
+template<class Grid, class GridView, class Mapper, class InterfaceGridView,
+ class InterfaceMapper, class Problem>
+class MMDG : public DG<Grid, GridView, Mapper, Problem>
+{
+public:
+ static constexpr int dim = GridView::dimension;
+
+ using Base = DG<Grid, GridView, Mapper, Problem>;
+ using Scalar = typename Base::Scalar;
+ using Matrix = typename Base::Matrix;
+ using Vector = typename Base::Vector;
+ using Coordinate = Dune::FieldVector<Scalar, dim>;
+ using QRuleInterface = Base::QRuleEdge;
+ using QRulesInterface = Base::QRuleEdges;
+ static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);
+
+ //constructor
+ MMDG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
+ const InterfaceGridView& iGridView, const InterfaceMapper& iMapper,
+ const Problem& problem) :
+ Base(grid, gridView, mapper, problem,
+ (1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
+ iGridView_(iGridView), iMapper_(iMapper) {}
+
+ const void operator() (const Scalar mu)
+ {
+ assembleSLE(mu);
+ Base::A->compress();
+ }
+
+ const InterfaceGridView& iGridView_;
+ const InterfaceMapper& iMapper_;
+
+private:
+ //assemble stiffness matrix A and load vector b
+ void assembleSLE (const Scalar mu)
+ {
+ //quadrature rules
+ const QRuleInterface& rule_qInterface =
+ QRulesInterface::rule(iSimplex, problem_.quadratureOrder_qInterface());
+
+ //we use the bulk basis
+ // phi_elem,0 (x) = indicator(elem);
+ // phi_elem,i (x) = x[i]*indicator(elem);
+ //for elem in elements(gridView) and i = 1,...,dim
+
+ //we use the interface basis
+ // phi_iElem,0 (x) = indicator(elem);
+ // phi_iElem,i (x) = (x * tau_i)*indicator(elem);
+ //for iElem in elements(iGridView) and i = 1,...,dim-1
+ //where {tau_i | i=1,...,dim-1} is the interface coordinate system given
+ //by the method interfaceFrame(normal)
+
+ //in total we use the basis given by the union of
+ // { (phi_elem,i , 0) | elem in elements(gridView), i = 0,...,dim }
+ //and
+ // { (0 , phi_iElem,i ) | iElem in elements(iGridView), i = 0,...,dim-1 }
+
+ //thus the stiffness matrix A and the load vector b exhibit the following
+ //block structure
+ // A = [ bulk, coupling b = [ bulk,
+ // coupling, interface], interface ]
+
+ //call base class method for the assignment to fill stiffness matrix A
+ //and load vector b with bulk entries
+ Base::assembleSLE(mu);
+
+ for (const auto& iElem : elements(iGridView_))
+ {
+ const int iElemIdx = iMapper_.index(iElem);
+ const auto& iGeo = iElem.geometry();
+
+ //get basis of the interface coordinate system for evaluation of
+ //basis function phi_iElem,i with i=1,...,dim-1
+ const std::array<Coordinate, dim-1> iFrame = interfaceFrame(
+ Base::grid_.asIntersection(iElem).centerUnitOuterNormal());
+
+ //in the total system of linear equations (SLE) Ad = b,
+ //the index iElemIdxSLE refers to the basis function (0, phi_iElem,0)
+ //and the indices iElemIdxSLE + i + 1, i = 1,...,dim-1, refer to the
+ //basis functions (0, phi_iElem,i)
+ const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
+
+ //fill load vector b with interface entries using a quadrature rule
+ for (const auto& ip : rule_qInterface)
+ {
+ //quadrature point in local and global coordinates
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = iGeo.global(qpLocal);
+
+ const Scalar weight(ip.weight());
+ const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+ const Scalar qInterfaceEvaluation(problem_.qInterface(qpGlobal));
+
+ //quadrature for int_elem q*phi_elem,0 dV
+ b[iElemIdxSLE] += weight * integrationElem * qEvaluation;
+
+ //quadrature for int_elem q*phi_elem,i dV
+ for (int i = 0; i < dim - 1; i++)
+ {
+ b[iElemIdxSLE + i + 1] +=
+ weight * integrationElem * (qpGlobal * iFrame[i]) * qEvaluation;
+ }
+ }
+
+ //TODO: fill stiffness matrix A with interface entries
+
+ //fill stiffness matrix A with coupling entries
+ for (const auto& elem : elements(Base::gridView_))
+ {
+ //TODO
+ }
+ }
+ }
+
+ //returns a basis of the interface coordinate system
+ //2d: tau = [ normal[1],
+ // -normal[0] ]
+ //3d: tau_1 ~ [ normal[1] + normal[2], tau_2 ~ [ -normal[1],
+ // -normal[0], normal[0] + normal[2]
+ // -normal[0] ], -normal[1] ]
+ std::array<Coordinate, dim-1> interfaceFrame(const Coordinate& normal)
+ {
+ for (int i = 0; i < dim - 1; i++)
+ {
+ Coordinate tau(0.0);
+ for (int j = 0; j < dim - 1; j++)
+ {
+ if (j == i)
+ {
+ continue;
+ }
+
+ tau[i] += normal[j];
+ tau[j] = -normal[i];
+ }
+ tau /= tau.two_norm();
+
+ interfaceFrame[i] = tau;
+ }
+/* //NOTE: which version is better?
+ if constexpr (dim == 2)
+ {
+ Coordinate tau;
+ tau[0] = normal[1];
+ tau[1] = -normal[0];
+
+ interfaceFrame = {tau};
+ }
+ else if constexpr (dim == 3)
+ {
+ Coordinate tau1;
+ tau1[0] = normal[1] + normal[2];
+ tau1[1] = -normal[0];
+ tau1[2] = -normal[0];
+ Coordinate tau2;
+ tau2[0] = -normal[1];
+ tau2[1] = normal[0] + normal[2];
+ tau2[2] = -normal[1];
+ interfaceFrame = {tau1, tau2};
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid dimension: dim = " + std::to_string(dim));
+ }*/
+ }
+
+
+};
+
+#endif
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index 64cc67fd5b9c163db702e1e184b8ccdbf0d13138..982015eaf089143a0fd010dba1fb30f72efa75d5 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/dune/mmdg/problems/dgproblem.hh
@@ -3,13 +3,11 @@
//abstract class providing an interface for a problem that is to be solved with
//a discontinuous Galerkin scheme
-template<class Coordinate, class ctype = double>
+template<class Vector, class Scalar = double>
class DGProblem
{
public:
- using Scalar = ctype;
- using Vector = Coordinate;
- static constexpr int dimension = Coordinate::dimension;
+ static constexpr int dimension = Vector::dimension;
using Matrix = Dune::FieldMatrix<Scalar, dimension, dimension>;
//the exact solution at position pos
diff --git a/src/dg.cc b/src/dg.cc
index 8b160fec3ca62e53607e1a442eb1fd2c8e92d804..9f7f07c5cecb2f2a0ccacd6cf148c512350d824c 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -41,6 +41,8 @@ int main(int argc, char** argv)
using GridView = typename Grid::LeafGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Coordinate = Dune::FieldVector<double, dim>;
+ using Problem = DGProblem<Coordinate>;
+ using DG = DG<Grid, GridView, Mapper, Problem>;
//get parameters from file "parameterDG.ini"
Dune::ParameterTree pt;
@@ -50,15 +52,14 @@ int main(int argc, char** argv)
const double mu = std::stod(pt["mu"]);
//create a grid from .dgf file
- GridFactory gridFactory( "grids/cube" + std::to_string(dim) +
- "d.dgf" );
+ 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());
- DGProblem<Coordinate>* problem; //problem to be solved
+ Problem* problem; //problem to be solved
//determine problem type
if (pt["problem"] == "poisson")
@@ -72,7 +73,8 @@ int main(int argc, char** argv)
"parameterDG.ini (use 'poisson' or TODO)");
}
- const DG dg(gridView, mapper, *problem);
+ //DG<GridView, Mapper, Problem> dg(gridView, mapper, *problem);
+ DG dg(grid, gridView, mapper, *problem);
dg(mu);
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 61ea9c46f1b6dd98d48505af1f71279d9678ee9e..aea7f9a8660c145c718a1f1a56addde93fccb578 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -1,30 +1,102 @@
-// -*- 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 <sys/stat.h>
+#include <sys/types.h>
+#include <type_traits>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.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/mmesh/mmesh.hh>
+
+#include <dune/mmdg/mmdg.hh>
+#include <dune/mmdg/clockhelper.hh>
+#include <dune/mmdg/problems/dgproblem.hh>
+#include <dune/mmdg/problems/poisson.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;
+ 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::MovingMesh<dim>;
+ using GridFactory = Dune::DGFGridFactory<Grid>;
+ using GridView = typename Grid::LeafGridView;
+ using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+ using IGrid = typename Grid::InterfaceGrid;
+ using IGridView = typename IGrid::LeafGridView;
+ using IMapper =
+ typename Dune::MultipleCodimMultipleGeomTypeMapper<IGridView>;
+ using Coordinate = Dune::FieldVector<double, dim>;
+ using Problem = DGProblem<Coordinate>;
+ using MMDG = MMDG<Grid, GridView, Mapper, IGridView, IMapper, Problem>;
+
+ //get parameters from file "parameterDG.ini"
+ Dune::ParameterTree pt;
+ Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
+
+ //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();
+ const IGrid& iGrid = grid.interfaceGrid();
+ const IGridView& iGridView = iGrid.leafGridView();
+
+ //element mapper for the leaf grid views
+ const Mapper mapper(gridView, Dune::mcmgElementLayout());
+ const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
+
+ Problem* problem; //problem to be solved
+
+ //determine problem type
+ if (pt["problem"] == "poisson")
+ {
+ problem = new Poisson<dim>();
+ }
else
- std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
- <<" processes!"<<std::endl;
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid problem type or parameter 'problem' not specified in file "
+ "parameterDG.ini (use 'poisson' or TODO)");
+ }
+
+ MMDG mmdg(grid, gridView, mapper, iGridView, iMapper, *problem);
+ mmdg(mu);
+
+ //print total run time
+ const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
+ std::cout << "Finished with a total run time of " << timeTotal <<
+ " Seconds.\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/parameterDG.ini b/src/parameterDG.ini
index 808b309c788e7fb01abf6fe8127d74a91f0a991e..165571303cbed4b80a298d46af9b7007f682ed5d 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,3 +1,2 @@
-K = 1
mu = 1000
problem = poisson
diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..165571303cbed4b80a298d46af9b7007f682ed5d 100644
--- a/src/parameterMMDG.ini
+++ b/src/parameterMMDG.ini
@@ -0,0 +1,2 @@
+mu = 1000
+problem = poisson