robust implementation of eigenvalue computation, prevents complex eigenvalues...

robust implementation of eigenvalue computation, prevents complex eigenvalues due to round-off errors

dune/tectonic/spatial-solving/tnnmg/functional.hh

@@ -16,6 +16,7 @@
 
 #include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
 
+#include "../../utils/maxeigenvalue.hh"
 #include "../../utils/debugutils.hh"
 
 template<class M, class V, class N, class L, class U, class R>
@@ -423,18 +424,19 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
   {
      for (size_t i=0; i<this->quadraticPart().N(); ++i) {
        const auto& Aii = this->quadraticPart()[i][i];
+       maxEigenvalues_[i] = maxEigenvalue(Aii);
 
        // due to numerical imprecision, Aii might not be symmetric leading to complex eigenvalues
        // consider Toeplitz decomposition of Aii and take only symmetric part
 
-       auto symBlock = Aii;
+       /*auto symBlock = Aii;
 
-       /*for (size_t j=0; j<symBlock.N(); j++) {
+       for (size_t j=0; j<symBlock.N(); j++) {
            for (size_t k=0; k<symBlock.M(); k++) {
                if (symBlock[j][k]/symBlock[j][j] < 1e-14)
                 symBlock[j][k] = 0;
            }
-       }*/
+       }
 
        try {
        typename Vector::block_type eigenvalues;
@@ -442,8 +444,11 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
        maxEigenvalues_[i] =
            *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
        } catch (Dune::Exception &e) {
+            print(symBlock, "symBlock");
+            typename Vector::block_type eigenvalues;
+            Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
             std::cout << "complex eig in dof: " << i << std::endl;
-       }
+       } */
      }
   }
 

dune/tectonic/utils/CMakeLists.txt

@@ -4,6 +4,7 @@ add_custom_target(tectonic_dune_utils SOURCES
   diameter.hh
   geocoordinate.hh
   index-in-sorted-range.hh
+  maxeigenvalue.hh
   tobool.hh
 )
 
@@ -14,5 +15,6 @@ install(FILES
   diameter.hh
   geocoordinate.hh
   index-in-sorted-range.hh
+  maxeigenvalue.hh
   tobool.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/utils/maxeigenvalue.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+
+/** \file
+ * \brief max eigenvalue computations for the FieldMatrix class
+ */
+
+#include <iostream>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
+    return matrix[0][0];
+}
+
+/** \brief calculates the eigenvalues of a symmetric field matrix
+    \param[in]  matrix matrix eigenvalues are calculated for
+    \param[out] max eigenvalue
+*/
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
+    const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
+    const K p2 = p - matrix[1][1];
+    K q = p2 * p2 + matrix[1][0] * matrix[0][1];
+    if( q < 0 && q > -1e-14 ) q = 0;
+    if (q < 0) {
+        std::cout << 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 = std::sqrt(q);
+
+    return p + q;
+}
+
+/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
+    \param[in]  matrix eigenvalues are calculated for
+    \param[out] max eigenvalue
+
+    \note If the input matrix is not symmetric the behavior of this method is undefined.
+ */
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
+    Dune::FieldVector<K, 3> eigenvalues;
+    Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+    return eigenvalues[0];
+}
+
+  /** @} end documentation */
+
+#endif