From d2987209089f5c38331d2f71e3c09106652d97df Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 27 Dec 2015 17:45:50 +0100
Subject: [PATCH] Add more matrix operations, with an eye on adouble types
This patch adds several more implementations of operator+ etc for FieldMatrices.
In particular, you can now combine matrices of adoubles (the 'active' type from
the ADOL-C automatic differentiation library) with regular double types.
These new operators are needed for the NonplanarCosseratShellEnergy class,
which is coming in the next patch.
---
dune/gfe/linearalgebra.hh | 126 +++++++++++++++++++++++++++++++++++---
1 file changed, 117 insertions(+), 9 deletions(-)
diff --git a/dune/gfe/linearalgebra.hh b/dune/gfe/linearalgebra.hh
index 8e06338f..12ef869d 100644
--- a/dune/gfe/linearalgebra.hh
+++ b/dune/gfe/linearalgebra.hh
@@ -3,17 +3,17 @@
#include <dune/common/fmatrix.hh>
-//! calculates ret = A * B
-template< class K, int m, int n, int p >
-Dune::FieldMatrix< K, m, p > operator* ( const Dune::FieldMatrix< K, m, n > &A, const Dune::FieldMatrix< K, n, p > &B)
+template< int m, int n, int p >
+auto operator* ( const Dune::FieldMatrix< adouble, m, n > &A, const Dune::FieldMatrix< adouble, n, p > &B)
+ -> Dune::FieldMatrix<adouble, m, p>
{
- typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
- Dune::FieldMatrix< K, m, p > ret;
+ typedef typename Dune::FieldMatrix< adouble, m, p > :: size_type size_type;
+ Dune::FieldMatrix< adouble, m, p > ret;
for( size_type i = 0; i < m; ++i ) {
for( size_type j = 0; j < p; ++j ) {
- ret[ i ][ j ] = K( 0 );
+ ret[ i ][ j ] = 0.0;
for( size_type k = 0; k < n; ++k )
ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
}
@@ -21,20 +21,115 @@ Dune::FieldMatrix< K, m, p > operator* ( const Dune::FieldMatrix< K, m, n > &A,
return ret;
}
-//! calculates ret = A - B
+template< int m, int n, int p >
+auto operator* ( const Dune::FieldMatrix< adouble, m, n > &A, const Dune::FieldMatrix< double, n, p > &B)
+ -> Dune::FieldMatrix<adouble, m, p>
+{
+ typedef typename Dune::FieldMatrix< adouble, m, p > :: size_type size_type;
+ Dune::FieldMatrix< adouble, m, p > ret;
+
+ for( size_type i = 0; i < m; ++i ) {
+
+ for( size_type j = 0; j < p; ++j ) {
+ ret[ i ][ j ] = 0.0;
+ for( size_type k = 0; k < n; ++k )
+ ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
+ }
+ }
+ return ret;
+}
+
+template< int m, int n, int p >
+auto operator* ( const Dune::FieldMatrix< double, m, n > &A, const Dune::FieldMatrix< adouble, n, p > &B)
+ -> Dune::FieldMatrix<adouble, m, p>
+{
+ typedef typename Dune::FieldMatrix< adouble, m, p > :: size_type size_type;
+ Dune::FieldMatrix< adouble, m, p > ret;
+
+ for( size_type i = 0; i < m; ++i ) {
+
+ for( size_type j = 0; j < p; ++j ) {
+ ret[ i ][ j ] = 0.0;
+ for( size_type k = 0; k < n; ++k )
+ ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
+ }
+ }
+ return ret;
+}
+
+//! calculates ret = A + B
template< class K, int m, int n>
-Dune::FieldMatrix<K,m,n> operator- ( const Dune::FieldMatrix<K, m, n> &A, const Dune::FieldMatrix<K,m,n> &B)
+Dune::FieldMatrix<K,m,n> operator+ ( const Dune::FieldMatrix<K, m, n> &A, const Dune::FieldMatrix<K,m,n> &B)
{
typedef typename Dune::FieldMatrix<K,m,n> :: size_type size_type;
Dune::FieldMatrix<K,m,n> ret;
for( size_type i = 0; i < m; ++i )
for( size_type j = 0; j < n; ++j )
- ret[i][j] = A[i][j] - B[i][j];
+ ret[i][j] = A[i][j] + B[i][j];
return ret;
}
+//! calculates ret = A - B
+template <int m, int n>
+auto operator- ( const Dune::FieldMatrix< adouble, m, n > &A, const Dune::FieldMatrix< double, m, n > &B)
+ -> Dune::FieldMatrix<adouble, m, n>
+{
+ Dune::FieldMatrix<adouble,m,n> result;
+ typedef typename decltype(result)::size_type size_type;
+
+ for( size_type i = 0; i < m; ++i )
+ for( size_type j = 0; j < n; ++j )
+ result[i][j] = A[i][j] - B[i][j];
+
+ return result;
+}
+
+//! calculates ret = A - B
+template <int m, int n>
+auto operator- ( const Dune::FieldMatrix< adouble, m, n > &A, const Dune::FieldMatrix< adouble, m, n > &B)
+ -> Dune::FieldMatrix<adouble, m, n>
+{
+ Dune::FieldMatrix<adouble,m,n> result;
+ typedef typename decltype(result)::size_type size_type;
+
+ for( size_type i = 0; i < m; ++i )
+ for( size_type j = 0; j < n; ++j )
+ result[i][j] = A[i][j] - B[i][j];
+
+ return result;
+}
+
+//! calculates ret = s*A
+template< int m, int n>
+auto operator* ( const double& s, const Dune::FieldMatrix<adouble, m, n> &A)
+ -> Dune::FieldMatrix<adouble,m,n>
+{
+ typedef typename Dune::FieldMatrix<adouble,m,n> :: size_type size_type;
+ Dune::FieldMatrix<adouble,m,n> ret;
+
+ for( size_type i = 0; i < m; ++i )
+ for( size_type j = 0; j < n; ++j )
+ ret[i][j] = s * A[i][j];
+
+ return ret;
+}
+
+//! calculates ret = s*A
+template< int m, int n>
+auto operator* ( const double& s, const Dune::FieldMatrix<double, m, n> &A)
+ -> Dune::FieldMatrix<double,m,n>
+{
+ typedef typename Dune::FieldMatrix<double,m,n> :: size_type size_type;
+ Dune::FieldMatrix<double,m,n> ret;
+
+ for( size_type i = 0; i < m; ++i )
+ for( size_type j = 0; j < n; ++j )
+ ret[i][j] = s * A[i][j];
+
+ return ret;
+}
//! calculates ret = A/s
template< class K, int m, int n>
@@ -50,4 +145,17 @@ Dune::FieldMatrix<K,m,n> operator/ ( const Dune::FieldMatrix<K, m, n> &A, const
return ret;
}
+//! calculates ret = A/s
+template< class K, int m>
+Dune::FieldVector<K,m> operator/ ( const Dune::FieldVector<K, m> &A, const K& s)
+{
+ typedef typename Dune::FieldVector<K,m> :: size_type size_type;
+ Dune::FieldVector<K,m> result;
+
+ for( size_type i = 0; i < m; ++i )
+ result[i] = A[i] / s;
+
+ return result;
+}
+
#endif
--
GitLab