From 33c3fd6cdc998a21c707f750204f5074131dec9b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sun, 29 Mar 2020 15:45:15 +0200
Subject: [PATCH] add penalty parameter to mmdg, [bugfix] getDiameter now runs
over subentities with codim=1
---
dune/mmdg/dg.hh | 44 ++++++++++++++++++++++----------------------
dune/mmdg/mmdg.hh | 35 ++++++++++++++++++++++++++++++-----
src/dg.cc | 2 +-
src/dgAnalysis.py | 4 ++--
src/mmdg.cc | 2 +-
src/mmdgAnalysis.py | 4 ++--
6 files changed, 58 insertions(+), 33 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 79810c0..d77c032 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -32,7 +32,6 @@ public:
using Matrix = Dune::DynamicMatrix<Scalar>;
using Vector = Dune::DynamicVector<Scalar>;
-
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<Scalar>>;
@@ -566,6 +565,28 @@ protected:
vtkWriter.pwrite("pressureData", "", "data");
}
+ //returns the diameter of a grid element elem
+ const Scalar getDiameter (const auto& elem) const
+ {
+ Scalar diameter(0.0);
+
+ for (int i = 0; i < ((dim + 1) * dim) / 2; i++)
+ {
+ const auto& edge = elem.template subEntity<1>(i);
+ diameter = std::max(diameter, edge.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();
+ }
+
std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
@@ -596,27 +617,6 @@ private:
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/mmdg.hh b/dune/mmdg/mmdg.hh
index 4d58fd3..489b352 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -31,9 +31,9 @@ public:
(1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
iGridView_(iGridView), iMapper_(iMapper) {}
- const void operator() (const Scalar mu, const Scalar xi)
+ const void operator() (const Scalar mu0, const Scalar xi)
{
- assembleSLE(mu, xi);
+ assembleSLE(mu0, xi);
//Base::A->compress();
(*Base::A).solve(Base::d, Base::b);
@@ -100,7 +100,7 @@ public:
private:
//assemble stiffness matrix A and load vector b
- void assembleSLE (const Scalar mu, const Scalar xi)
+ void assembleSLE (const Scalar mu0, const Scalar xi)
{
//modell parameter
const auto beta = [this, &xi](const Coordinate& pos) -> Scalar
@@ -142,12 +142,13 @@ private:
//call base class method for the assignment to fill stiffness matrix A
//and load vector b with bulk entries
- Base::assembleSLE(mu);
+ Base::assembleSLE(mu0);
for (const auto& iElem : elements(iGridView_))
{
const int iElemIdx = iMapper_.index(iElem);
const auto& iGeo = iElem.geometry();
+ const auto& iGeoCenter = iGeo.center();
const auto& intersct = Base::gridView_.grid().asIntersection(iElem);
//get basis of the interface coordinate system for evaluation of
@@ -175,7 +176,8 @@ private:
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar betaEvaluation = beta(qpGlobal);
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] += quadratureFactor * betaEvaluation;
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
+ quadratureFactor * betaEvaluation;
for (int i = 0; i < dim - 1; i++)
{
@@ -478,6 +480,9 @@ private:
const auto& normal = intersct.centerUnitOuterNormal();
Coordinate KparDotTau_i;
+ //determine penalty parameter
+ const Scalar mu = getPenaltyParameter(mu0, intersct, iGeoCenter);
+
const Scalar dEvaluation = Base::problem_.aperture(center);
const Scalar dEvaluation2 = dEvaluation * dEvaluation;
@@ -762,6 +767,26 @@ private:
vtkWriter.pwrite("interfacePressureData", "", "data");
}
+ //returns the interface penalty parameter mu^Gamma for an interface facet
+ //given as intersection of the interface grid where center is the center of
+ //intersection.inside()
+ const Scalar getPenaltyParameter (const Scalar mu0, const auto& intersection,
+ const auto& center) const
+ {
+ Scalar mu = Base::getLargestEigenvalue(Base::problem_.Kparallel(center))
+ / Base::getDiameter(intersection.inside());
+
+ if (intersection.neighbor())
+ {
+ const auto& neighbor = intersection.outside();
+ mu = std::max(mu, Base::getLargestEigenvalue(
+ Base::problem_.Kparallel(neighbor.geometry().center()))
+ / Base::getDiameter(neighbor));
+ }
+
+ return mu * mu0 * 2.0 * (1.0 + dim);
+ }
+
};
#endif
diff --git a/src/dg.cc b/src/dg.cc
index 8f60fc3..dcc4659 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -114,7 +114,7 @@ int main(int argc, char** argv)
//print total run time and error
std::cout << "Finished with a total run time of " << timeTotal <<
- " Seconds.\nL2-error: " << error << ".\n";
+ " Seconds.\nL2-error: " << error << ".\n\n";
return 0;
}
diff --git a/src/dgAnalysis.py b/src/dgAnalysis.py
index 92ddfa4..d8a7262 100644
--- a/src/dgAnalysis.py
+++ b/src/dgAnalysis.py
@@ -97,11 +97,11 @@ for dim in dimensions:
EOC[i] = -1.0*float('inf')
# end for i
- print("\n\n-------------------------------------------------------------\n")
# print EOCs
+ print("\n-------------------------------------------------------------\n")
print("EOCs (" + str(dim) + "d):\t")
print(EOC)
- print("\n\n")
+ print("\n=============================================================\n\n")
# write data to file
np.savetxt(join("plots", "dgAnalysis_" + str(dim) + "d"), data)
diff --git a/src/mmdg.cc b/src/mmdg.cc
index d58637c..726d33f 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -140,7 +140,7 @@ int main(int argc, char** argv)
//print total run time and error
std::cout << "Finished with a total run time of " << timeTotal <<
- " Seconds.\nL2-error: " << error << ".\n";
+ " Seconds.\nL2-error: " << error << ".\n\n";
return 0;
}
diff --git a/src/mmdgAnalysis.py b/src/mmdgAnalysis.py
index 89c1f6e..e8312ed 100644
--- a/src/mmdgAnalysis.py
+++ b/src/mmdgAnalysis.py
@@ -102,11 +102,11 @@ for i in range(len(gridfiles) - 1):
EOC[i] = -1.0*float('inf')
# end for i
-print("\n\n-------------------------------------------------------------\n")
# print EOCs
+print("\n-------------------------------------------------------------\n")
print("EOCs (" + str(dim) + "d):\t")
print(EOC)
-print("\n\n")
+print("\n=============================================================\n\n")
# write data to file
np.savetxt(join("plots", "mmdgAnalysis_" + str(dim) + "d"), data)
--
GitLab