From bc5052be6489b6a4b558691afe8af323152aab88 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sat, 28 Mar 2020 21:06:11 +0100
Subject: [PATCH] implement penalty paramter for dg
---
dune/mmdg/dg.hh | 51 +++++++++++++++++++++++++++++++--
dune/mmdg/problems/dgproblem.hh | 1 +
src/dg.cc | 11 +++----
src/dgAnalysis.py | 2 +-
src/mmdg.cc | 4 +--
src/parameterDG.ini | 2 +-
src/parameterMMDG.ini | 2 +-
7 files changed, 58 insertions(+), 15 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 7ae4a82..79810c0 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -26,6 +26,7 @@ public:
static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
using Scalar = typename GridView::ctype;
+ using Coordinate = typename Problem::CoordinateType;
// using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
// using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
@@ -56,10 +57,10 @@ public:
// d = Vector(dof_);
}
- void operator() (const Scalar mu)
+ void operator() (const Scalar mu0)
{
//assemble stiffness matrix A and load vector b
- assembleSLE(mu);
+ assembleSLE(mu0);
//A->compress();
(*A).solve(d, b);
@@ -143,7 +144,7 @@ protected:
}
//assemble stiffness matrix A and load vector b
- void assembleSLE (const Scalar mu)
+ void assembleSLE (const Scalar mu0)
{
//quadrature rules
const QRuleElem& rule_q =
@@ -166,6 +167,7 @@ protected:
{
const int elemIdx = mapper_.index(elem);
const auto& geo = elem.geometry();
+ const auto& center = geo.center();
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
@@ -233,6 +235,9 @@ protected:
const double intersctVol = intersctGeo.volume();
const auto& intersctCenter = intersctGeo.center();
+ //determine penalty parameter
+ const Scalar mu = getPenaltyParameter(mu0, intersct, center);
+
//create storage for multiply used integral values
//linearIntegrals[i] = int_intersct x[i] ds
Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
@@ -572,6 +577,46 @@ private:
writeVTKoutput(dof_);
}
+ //returns the penalty parameter mu for a facet given as intersection of the
+ //grid where center is the center of intersection.inside()
+ const Scalar getPenaltyParameter (const Scalar mu0, const auto& intersection,
+ const auto& center) const
+ {
+ Scalar mu = getLargestEigenvalue(problem_.K(center))
+ / getDiameter(intersection.inside());
+
+ if (intersection.neighbor())
+ {
+ const auto& neighbor = intersection.outside();
+ mu = std::max(mu,
+ getLargestEigenvalue(problem_.K(neighbor.geometry().center()))
+ / getDiameter(neighbor));
+ }
+
+ return mu * mu0 * 2.0 * (1.0 + dim);
+ }
+
+ //returns the diameter of a grid element elem
+ const Scalar getDiameter (const auto& elem) const
+ {
+ Scalar diameter(0.0);
+
+ for (const auto& intersct : intersections(gridView_, elem))
+ {
+ diameter = std::max(diameter, intersct.geometry().volume());
+ }
+
+ return diameter;
+ }
+
+ //determines the largest absolute eigenvalue of a matrix
+ const Scalar getLargestEigenvalue (const auto& matrix) const
+ {
+ Coordinate eigenvalues;
+ Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+ return eigenvalues.infinity_norm();
+ }
+
};
#endif
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index 5f7e630..38ed307 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/dune/mmdg/problems/dgproblem.hh
@@ -9,6 +9,7 @@ class DGProblem
public:
static constexpr int dim = Coordinate::dimension;
using Matrix = Dune::FieldMatrix<Scalar, dim, dim>;
+ using CoordinateType = Coordinate;
//the exact solution at position pos
virtual Scalar exactSolution (const Coordinate& pos) const
diff --git a/src/dg.cc b/src/dg.cc
index d2bba7c..8f60fc3 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -39,9 +39,7 @@ int main(int argc, char** argv)
//use YaspGrid if dim = 1 and MMesh otherwise
using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
Dune::MovingMesh<dim>>::type;
- //using GridFactory = Dune::DGFGridFactory<Grid>;
- using GridFactory = Dune::GmshGridFactory<Grid, /*implicit=*/false>;
-
+ using GridFactory = Dune::DGFGridFactory<Grid>;
using GridView = typename Grid::LeafGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using Coordinate = Dune::FieldVector<double, dim>;
@@ -53,11 +51,10 @@ int main(int argc, char** argv)
Dune::ParameterTreeParser::readINITree("parameterDG.ini", pt);
//assume constant penalty parameter
- const double mu = std::stod(pt["mu"]);
+ const double mu0 = std::stod(pt["mu0"]);
//create a grid from .dgf file
- //GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
- GridFactory gridFactory( "grids/mmdg2.msh" );
+ GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
const Grid& grid = *gridFactory.grid();
const GridView& gridView = grid.leafGridView();
@@ -89,7 +86,7 @@ int main(int argc, char** argv)
DG dg(gridView, mapper, *problem);
- dg(mu);
+ dg(mu0);
//error and total run time
const double error = dg.computeL2error();
diff --git a/src/dgAnalysis.py b/src/dgAnalysis.py
index cf7c9ba..92ddfa4 100644
--- a/src/dgAnalysis.py
+++ b/src/dgAnalysis.py
@@ -8,7 +8,7 @@ from math import log
# writes a dgf file for a grid with N elements per directions and dimension dim
def writeDGF(N, dim):
# open dgf file
- file = open(join("grids", "cube" + str(dim) + "d.dgf"), "rw")
+ file = open(join("grids", "cube" + str(dim) + "d.dgf"), "w")
file.write("DGF\n\nInterval\n")
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 894c304..d58637c 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -54,7 +54,7 @@ int main(int argc, char** argv)
Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
//assume constant penalty parameter //TODO
- const double mu = std::stod(pt["mu"]);
+ const double mu0 = std::stod(pt["mu0"]);
const double xi = std::stod(pt["xi"]); //coupling parameter
const double aperture = std::stod(pt["d"]);
@@ -112,7 +112,7 @@ int main(int argc, char** argv)
const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
- mmdg(mu, xi);
+ mmdg(mu0, xi);
//error and total run time
const double error = mmdg.computeL2error();
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index 1655713..10896d3 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,2 +1,2 @@
-mu = 1000
+mu0 = 1000
problem = poisson
diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini
index d214f65..69a5f03 100644
--- a/src/parameterMMDG.ini
+++ b/src/parameterMMDG.ini
@@ -1,4 +1,4 @@
-mu = 1000
+mu0 = 1000
xi = 0.75
d = 1
problem = mmdg2
--
GitLab