diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 8a00b9bc87c1205c1abb56a5b8e237ad5f8ef674..a2aa5715c1e90e12491dd7df1854c08e2fe895ad 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -17,7 +17,6 @@
#include <dune/istl/umfpack.hh>
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
-#include <dune/mmdg/eigenvaluehelper.hh>
#include <dune/mmdg/solvertype.hh>
template<class GridView, class Mapper, class Problem>
@@ -550,14 +549,6 @@ protected:
return diameter;
}
- //determines the largest absolute eigenvalue of a matrix
- const Scalar getLargestEigenvalue (const auto& matrix) const
- {
- Coordinate eigenvalues;
- EigenvalueHelper::eigenValues(matrix, eigenvalues);
- return eigenvalues.infinity_norm();
- }
-
//apply the quadrature rule "rule" on the given geometry as specified in the
//lambda function "lambda", that takes the quadrature point as global
//coordinate and the quadrature weight (resp. weight * integrationElement) as
@@ -595,15 +586,12 @@ private:
const Scalar getPenaltyParameter (const Scalar mu0, const auto& intersection,
const auto& center) const
{
- Scalar mu = getLargestEigenvalue(problem_.K(center))
- / getDiameter(intersection.inside());
+ Scalar mu = 1.0 / getDiameter(intersection.inside());
if (intersection.neighbor())
{
const auto& neighbor = intersection.outside();
- mu = std::max(mu,
- getLargestEigenvalue(problem_.K(neighbor.geometry().center()))
- / getDiameter(neighbor));
+ mu = std::max(mu, 1.0 / getDiameter(neighbor));
}
return mu * mu0 * 2.0 * (1.0 + dim);
diff --git a/dune/mmdg/eigenvaluehelper.hh b/dune/mmdg/eigenvaluehelper.hh
deleted file mode 100644
index 7fe56886e13b5457466458d04b91826281c81655..0000000000000000000000000000000000000000
--- a/dune/mmdg/eigenvaluehelper.hh
+++ /dev/null
@@ -1,130 +0,0 @@
-#ifndef DUNE_MMDG_EIGENVALUEHELPER_HH
-#define DUNE_MMDG_EIGENVALUEHELPER_HH
-
-#include <cmath>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/math.hh>
-
-/** \file
- * \brief Eigenvalue computations for the Dune::FieldMatrix class
- * NOTE: taken from DUNE-common (fmatrixev.hh) with minor changes
- */
-class EigenvalueHelper {
-
-public:
- /** \brief calculates the eigenvalues of a symmetric field matrix
- \param[in] matrix matrix eigenvalues are calculated for
- \param[out] eigenvalues Dune::FieldVector that contains eigenvalues in
- ascending order
- */
- template <typename K>
- static void eigenValues(const Dune::FieldMatrix<K, 1, 1>& matrix,
- Dune::FieldVector<K, 1>& eigenvalues)
- {
- eigenvalues[0] = matrix[0][0];
- }
-
- /** \brief calculates the eigenvalues of a symmetric field matrix
- \param[in] matrix matrix eigenvalues are calculated for
- \param[out] eigenvalues Dune::FieldVector that contains eigenvalues in
- ascending order
- */
- template <typename K>
- static void eigenValues(const Dune::FieldMatrix<K, 2, 2>& matrix,
- Dune::FieldVector<K, 2>& eigenvalues)
- {
- if ( std::abs(matrix[1][0]) < 1e-15 && std::abs(matrix[0][1]) < 1e-15 )
- { //diagonal matrix
- eigenvalues[0] = std::min(matrix[0][0], matrix[1][1]);
- eigenvalues[1] = std::max(matrix[0][0], matrix[1][1]);
- return;
- }
-
- using std::sqrt;
- const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
- const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
- K q = p * p - detM;
- if( q < 0 && q > -1e-13/*-1e-14*/ ) q = 0;
- if (q < 0)
- {
- std::cout << std::numeric_limits<K>::epsilon() << "\n" << q << "\n" << matrix << std::endl;
- // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
- DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
- }
-
- // get square root
- q = sqrt(q);
-
- // store eigenvalues in ascending order
- eigenvalues[0] = p - q;
- eigenvalues[1] = p + q;
- }
-
- /** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
- \param[in] matrix matrix eigenvalues are calculated for
- \param[out] eigenvalues Eigenvalues in ascending order
-
- \note If the input matrix is not symmetric the behavior of this method is undefined.
-
- This implementation was adapted from the pseudo-code (Python?) implementation found on
- http://en.wikipedia.org/wiki/Eigenvalue_algorithm (retrieved late August 2014).
- Wikipedia claims to have taken it from
- Smith, Oliver K. (April 1961), Eigenvalues of a symmetric 3 × 3 matrix.,
- Communications of the ACM 4 (4): 168, doi:10.1145/355578.366316
- */
- template <typename K>
- static void eigenValues(const Dune::FieldMatrix<K, 3, 3>& matrix,
- Dune::FieldVector<K, 3>& eigenvalues)
- {
- using std::sqrt;
- using std::acos;
- using real_type = typename Dune::FieldTraits<K>::real_type;
- const K pi = Dune::MathematicalConstants<K>::pi();
- K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
-
- if (p1 <= 1e-8)
- {
- // A is diagonal.
- eigenvalues[0] = matrix[0][0];
- eigenvalues[1] = matrix[1][1];
- eigenvalues[2] = matrix[2][2];
- }
- else
- {
- // q = trace(A)/3
- K q = 0;
- for (int i=0; i<3; i++)
- q += matrix[i][i]/3.0;
-
- K p2 = (matrix[0][0] - q)*(matrix[0][0] - q) + (matrix[1][1] - q)*(matrix[1][1] - q) + (matrix[2][2] - q)*(matrix[2][2] - q) + 2 * p1;
- K p = sqrt(p2 / 6);
- // B = (1 / p) * (A - q * I); // I is the identity matrix
- Dune::FieldMatrix<K,3,3> B;
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- B[i][j] = (real_type(1.0)/p) * (matrix[i][j] - q*(i==j));
-
- K r = B.determinant() / 2.0;
-
- // In exact arithmetic for a symmetric matrix -1 <= r <= 1
- // but computation error can leave it slightly outside this range.
- K phi;
- if (r <= -1)
- phi = pi / 3.0;
- else if (r >= 1)
- phi = 0;
- else
- phi = acos(r) / 3;
-
- // the eigenvalues satisfy eig[2] <= eig[1] <= eig[0]
- eigenvalues[2] = q + 2 * p * cos(phi);
- eigenvalues[0] = q + 2 * p * cos(phi + (2*pi/3));
- eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2]; // since trace(matrix) = eig1 + eig2 + eig3
- }
- }
-};
-
-#endif
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 0889dafd7bd1519bd9a1b3c4de3be1d53bebee0b..7aec00703505c1292ba99b465ae1e1ee0e3b39a4 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -976,18 +976,15 @@ private:
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());
+ Scalar mu = 1.0 / 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));
+ mu = std::max(mu, 1.0 / Base::getDiameter(neighbor));
}
- return mu * mu0 * 2.0 * (1.0 + dim);
+ return mu * mu0 * 2.0 * dim;
}
//returns the outer product v * v^T of the vector v