diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 7ae4a826f980642a79f9d4b6c9c1056e33406e4c..79810c0a6c94a89ef6e860318f608f0c85e86434 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 5f7e6307c6b7823e4107807ff636b73bc3d09a3f..38ed3077b35279cf29c0d64ba1e5fc92c46d603d 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 d2bba7c8ae1d818bd11e0908969cbeb332c1b8cf..8f60fc3df707dba29db0f6420efc6a2f043f0462 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 cf7c9ba8a97de93b1d11882cb3c76703daa9bd02..92ddfa4ad6d118cb0599ee057cfe85acbd450e79 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 894c3044684293d6a78f6c7bdf46c0387cd4122e..d58637c107040a8aa3e2ef751eeeba656668f88d 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 165571303cbed4b80a298d46af9b7007f682ed5d..10896d307fc1d477c906a7d92644cee815001212 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 d214f65a5b92af753727c8050b00ed0444bc76c1..69a5f032bf4973d219c54707fad381c001e8ceae 100644
--- a/src/parameterMMDG.ini
+++ b/src/parameterMMDG.ini
@@ -1,4 +1,4 @@
-mu = 1000
+mu0 = 1000
xi = 0.75
d = 1
problem = mmdg2