Skip to content
Snippets Groups Projects
Commit 2f5d0bd3 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

remove permeability eigenvalue from penalty parameter

parent 983196fd
No related branches found
No related tags found
No related merge requests found
......@@ -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);
......
#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
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment