From f3689993d56a66c1198c230b5635d4df4a14e274 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 8 Mar 2012 10:38:01 +0000
Subject: [PATCH] A class for orthogonal matrices
[[Imported from SVN: r8554]]
---
dune/gfe/Makefile.am | 1 +
dune/gfe/orthogonalmatrix.hh | 93 ++++++++++++++++++++++++++++++++++++
2 files changed, 94 insertions(+)
create mode 100644 dune/gfe/orthogonalmatrix.hh
diff --git a/dune/gfe/Makefile.am b/dune/gfe/Makefile.am
index 70642e14..85557597 100644
--- a/dune/gfe/Makefile.am
+++ b/dune/gfe/Makefile.am
@@ -20,6 +20,7 @@ srcinclude_HEADERS = averagedistanceassembler.hh \
localgeodesicfestiffness.hh \
localgfetestfunctionbasis.hh \
maxnormtrustregion.hh \
+ orthogonalmatrix.hh \
pktop1mgtransfer.hh \
quaternion.hh \
realtuple.hh \
diff --git a/dune/gfe/orthogonalmatrix.hh b/dune/gfe/orthogonalmatrix.hh
new file mode 100644
index 00000000..62c29480
--- /dev/null
+++ b/dune/gfe/orthogonalmatrix.hh
@@ -0,0 +1,93 @@
+#ifndef ORTHOGONAL_MATRIX_HH
+#define ORTHOGONAL_MATRIX_HH
+
+#include <dune/common/fmatrix.hh>
+
+/** \brief An orthogonal \f$ n \times n \f$ matrix
+ * \tparam T Type of the matrix entries
+ * \tparam N Space dimension
+ */
+template <class T, int N>
+class OrthogonalMatrix
+{
+public:
+
+ /** \brief Constructor from a general matrix
+ *
+ * The input matrix gets normalized to det = 1. Since computing
+ * the determinant is expensive you may not want to use this method
+ * if you know your input matrix to be orthogonal anyways.
+ */
+ explicit OrthogonalMatrix(const Dune::FieldMatrix<T,N,N>& matrix)
+ : data_(matrix)
+ {
+ data_ /= matrix.determinant();
+ }
+
+ /** \brief Const access to the matrix entries */
+ const Dune::FieldMatrix<T,N,N>& data() const
+ {
+ return data_;
+ }
+
+ /** \brief Project the input matrix onto the tangent space at "this"
+ *
+ * See Absil, Mahoney, Sepulche, Example 5.3.2 for the formula
+ * \f[ P_X \xi = (I - XX^T) \xi + X \operatorname{skew} (X^T \xi) \f]
+ */
+ Dune::FieldMatrix<T,N,N> projectOntoTangentSpace(const Dune::FieldMatrix<T,N,N>& matrix) const
+ {
+ // rename to get the same notation as Absil & Co
+ const Dune::FieldMatrix<T,N,N>& X = data_;
+
+ // I - XX^T
+ Dune::FieldMatrix<T,N,N> IdMinusXXT;
+
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++) {
+ IdMinusXXT = (i==j);
+ for (int k=0; k<N; k++)
+ IdMinusXXT[i][j] -= X[i][k] * X[j][k];
+ }
+
+ // (I - XX^T) \xi
+ Dune::FieldMatrix<T,N,N> result(0);
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ for (int k=0; k<N; k++)
+ result[i][j] += IdMinusXXT[i][k] * matrix[k][j];
+
+ // X^T \xi
+ Dune::FieldMatrix<T,N,N> XTxi(0);
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ for (int k=0; k<N; k++)
+ XTxi[i][j] += X[k][i] * matrix[k][j];
+
+ // skew (X^T \xi)
+ Dune::FieldMatrix<T,N,N> skewXTxi(0);
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ skewXTxi[i][j] = 0.5 * ( XTxi[i][j] - XTxi[j][i] );
+
+ // X skew (X^T \xi)
+ Dune::FieldMatrix<T,N,N> XskewXTxi(0);
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ for (int k=0; k<N; k++)
+ XskewXTxi[i][j] += X[i][k] * skewXTxi[k][j];
+
+ //
+ result += XskewXTxi;
+
+ return result;
+ }
+
+private:
+
+ /** \brief The actual data */
+ Dune::FieldMatrix<T,N,N> data_;
+
+};
+
+#endif // ORTHOGONAL_MATRIX_HH
\ No newline at end of file
--
GitLab