From 61323dbf235d570ce02057dd45e90f647d1643af Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 2 Sep 2020 16:43:14 +0200
Subject: [PATCH] Move all generic low-level matrix methods to a central place
Several energy density implementations contained code for matrix methods
like sym, skew, transpose etc. This patch moves them all to the file
linearalgebra.hh, to get rid of the redundancy.
---
dune/gfe/cosseratenergystiffness.hh | 45 ++------
dune/gfe/linearalgebra.hh | 127 +++++++++++++++++++-
dune/gfe/nonplanarcosseratshellenergy.hh | 138 +++-------------------
dune/gfe/surfacecosseratenergy.hh | 140 +++--------------------
4 files changed, 158 insertions(+), 292 deletions(-)
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index aafdca5a..4e5fdbc0 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -74,37 +74,6 @@ class CosseratEnergyLocalStiffness
enum {gridDim=GridView::dimension};
enum {dimworld=GridView::dimensionworld};
-
- /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] + A[j][i]);
- return result;
- }
-
- /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] - A[j][i]);
- return result;
- }
-
- /** \brief Return the square of the trace of a matrix */
- template <int N>
- static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
- {
- field_type trace = 0;
- for (int i=0; i<N; i++)
- trace += A[i][i];
- return trace*trace;
- }
-
/** \brief Compute the (row-wise) curl of a matrix R \f$
\param DR The partial derivatives of the matrix R
*/
@@ -171,9 +140,9 @@ public:
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
- return mu_ * sym(UMinus1).frobenius_norm2()
- + mu_c_ * skew(UMinus1).frobenius_norm2()
- + (mu_*lambda_)/(2*mu_ + lambda_) * traceSquared(sym(UMinus1));
+ return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2()
+ + mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+ + (mu_*lambda_)/(2*mu_ + lambda_) * Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1));
}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
@@ -219,7 +188,7 @@ public:
RT detU = U.determinant();
- return mu_ * sym(UMinus1).frobenius_norm2() + mu_c_ * skew(UMinus1).frobenius_norm2()
+ return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2() + mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
}
@@ -275,9 +244,9 @@ public:
for (int k=0; k<3; k++)
RT_DR3[i][j] += R[k][i] * DR[k][2][j];
- return mu_ * sym(RT_DR3).frobenius_norm2()
- + mu_c_ * skew(RT_DR3).frobenius_norm2()
- + mu_*lambda_/(2*mu_+lambda_) * traceSquared(RT_DR3);
+ return mu_ * Dune::GFE::sym(RT_DR3).frobenius_norm2()
+ + mu_c_ * Dune::GFE::skew(RT_DR3).frobenius_norm2()
+ + mu_*lambda_/(2*mu_+lambda_) * Dune::GFE::traceSquared(RT_DR3);
}
/** \brief The shell thickness */
diff --git a/dune/gfe/linearalgebra.hh b/dune/gfe/linearalgebra.hh
index 7febb66e..a751897f 100644
--- a/dune/gfe/linearalgebra.hh
+++ b/dune/gfe/linearalgebra.hh
@@ -1,9 +1,13 @@
-#ifndef LINEAR_ALGEBRA_HH
-#define LINEAR_ALGEBRA_HH
+#ifndef DUNE_GFE_LINEARALGEBRA_HH
+#define DUNE_GFE_LINEARALGEBRA_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
+///////////////////////////////////////////////////////////////////////////////////////////
+// Various vector-space matrix methods
+///////////////////////////////////////////////////////////////////////////////////////////
+
#if ADOLC_ADOUBLE_H
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
template< int m, int n, int p >
@@ -171,4 +175,123 @@ Dune::FieldVector<K,m> operator/ ( const Dune::FieldVector<K, m> &A, const K& s)
}
#endif
+
+///////////////////////////////////////////////////////////////////////////////////////////
+// Various other matrix methods
+///////////////////////////////////////////////////////////////////////////////////////////
+
+namespace Dune {
+
+ namespace GFE {
+
+ /** \brief Return the trace of a matrix */
+ template <class T, int n>
+ static T trace(const FieldMatrix<T,n,n>& A)
+ {
+ T trace = 0;
+ for (int i=0; i<n; i++)
+ trace += A[i][i];
+ return trace;
+ }
+
+ /** \brief Return the square of the trace of a matrix */
+ template <class T, int n>
+ static T traceSquared(const FieldMatrix<T,n,n>& A)
+ {
+ T trace = 0;
+ for (int i=0; i<n; i++)
+ trace += A[i][i];
+ return trace*trace;
+ }
+
+ /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
+ template <class T, int n>
+ static FieldMatrix<T,n,n> sym(const FieldMatrix<T,n,n>& A)
+ {
+ FieldMatrix<T,n,n> result;
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result[i][j] = 0.5 * (A[i][j] + A[j][i]);
+ return result;
+ }
+
+ /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
+ template <class T, int n>
+ static FieldMatrix<T,n,n> skew(const FieldMatrix<T,n,n>& A)
+ {
+ FieldMatrix<T,n,n> result;
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result[i][j] = 0.5 * (A[i][j] - A[j][i]);
+ return result;
+ }
+
+ /** \brief Compute the deviator of a matrix A */
+ template <class T, int n>
+ static FieldMatrix<T,n,n> dev(const FieldMatrix<T,n,n>& A)
+ {
+ FieldMatrix<T,n,n> result = A;
+ auto t = trace(A);
+ for (int i=0; i<n; i++)
+ result[i][i] -= t / n;
+ return result;
+ }
+
+ /** \brief Return the transposed matrix */
+ template <class T, int n>
+ static FieldMatrix<T,n,n> transpose(const FieldMatrix<T,n,n>& A)
+ {
+ FieldMatrix<T,n,n> result;
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result[i][j] = A[j][i];
+
+ return result;
+ }
+
+ /** \brief The Frobenius (i.e., componentwise) product of two matrices */
+ template <class T, int n>
+ static T frobeniusProduct(const FieldMatrix<T,n,n>& A, const FieldMatrix<T,n,n>& B)
+ {
+ T result(0.0);
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result += A[i][j] * B[i][j];
+
+ return result;
+ }
+
+ template <class T, int n>
+ static auto dyadicProduct(const FieldVector<T,n>& A, const FieldVector<T,n>& B)
+ {
+ FieldMatrix<T,n,n> result;
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result[i][j] = A[i]*B[j];
+
+ return result;
+ }
+
+#if ADOLC_ADOUBLE_H
+ template <int n>
+ static auto dyadicProduct(const FieldVector<adouble,n>& A, const FieldVector<double,n>& B)
+ -> FieldMatrix<adouble,n,n>
+ {
+ FieldMatrix<adouble,n,n> result;
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result[i][j] = A[i]*B[j];
+
+ return result;
+ }
+#endif
+
+ }
+}
+
+
#endif
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index 16c97c4f..102115db 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -34,119 +34,6 @@ class NonplanarCosseratShellEnergy
enum {gridDim=GridView::dimension};
enum {dimworld=GridView::dimensionworld};
- /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] + A[j][i]);
- return result;
- }
-
- /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] - A[j][i]);
- return result;
- }
-
- /** \brief Compute the deviator of a matrix A */
- static Dune::FieldMatrix<field_type,dim,dim> dev(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result = A;
- auto t = trace(A);
- for (int i=0; i<dim; i++)
- result[i][i] -= t / dim;
- return result;
- }
-
- /** \brief Return the trace of a matrix */
- template <class T, int N>
- static T trace(const Dune::FieldMatrix<T,N,N>& A)
- {
- T trace = 0;
- for (int i=0; i<N; i++)
- trace += A[i][i];
- return trace;
- }
-
- /** \brief Return the square of the trace of a matrix */
- template <int N>
- static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
- {
- field_type trace = 0;
- for (int i=0; i<N; i++)
- trace += A[i][i];
- return trace*trace;
- }
-
- template <class T, int N>
- static T frobeniusProduct(const Dune::FieldMatrix<T,N,N>& A, const Dune::FieldMatrix<T,N,N>& B)
- {
- T result(0.0);
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result += A[i][j] * B[i][j];
-
- return result;
- }
-
- template <int N>
- static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<adouble,N>& B)
- -> Dune::FieldMatrix<adouble,N,N>
- {
- Dune::FieldMatrix<adouble,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[i]*B[j];
-
- return result;
- }
-
- template <int N>
- static auto dyadicProduct(const Dune::FieldVector<double,N>& A, const Dune::FieldVector<double,N>& B)
- -> Dune::FieldMatrix<double,N,N>
- {
- Dune::FieldMatrix<double,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[i]*B[j];
-
- return result;
- }
-
- template <int N>
- static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<double,N>& B)
- -> Dune::FieldMatrix<adouble,N,N>
- {
- Dune::FieldMatrix<adouble,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[i]*B[j];
-
- return result;
- }
-
- template <class T, int N>
- static Dune::FieldMatrix<T,N,N> transpose(const Dune::FieldMatrix<T,N,N>& A)
- {
- Dune::FieldMatrix<T,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[j][i];
-
- return result;
- }
-
public:
/** \brief Constructor with a set of material parameters
@@ -192,19 +79,20 @@ public:
RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T) const
{
- return mu_ * frobeniusProduct(sym(S), sym(T))
- + mu_c_ * frobeniusProduct(skew(S), skew(T))
- + lambda_ * mu_ / (lambda_ + 2*mu_) * trace(S) * trace(T);
+ return mu_ * Dune::GFE::frobeniusProduct(Dune::GFE::sym(S), Dune::GFE::sym(T))
+ + mu_c_ * Dune::GFE::frobeniusProduct(Dune::GFE::skew(S), Dune::GFE::skew(T))
+ + lambda_ * mu_ / (lambda_ + 2*mu_) * Dune::GFE::trace(S) * Dune::GFE::trace(T);
}
RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S) const
{
- return mu_ * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda_ * 0.5 * traceSquared(S);
+ return mu_ * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda_ * 0.5 * Dune::GFE::traceSquared(S);
}
RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S) const
{
- return mu_ * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
+ return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
+ + b2_ * Dune::GFE::skew(S).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(S));
}
/** \brief The shell thickness */
@@ -297,7 +185,7 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldMatrix<field_type,dim,dim> R;
value.q.matrix(R);
- auto RT = transpose(R);
+ auto RT = Dune::GFE::transpose(R);
Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
@@ -328,7 +216,7 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldMatrix<double,3,3> a(0);
for (int alpha=0; alpha<gridDim; alpha++)
- a += dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
+ a += Dune::GFE::dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
auto a00 = aCovariant[0] * aCovariant[0];
auto a01 = aCovariant[0] * aCovariant[1];
@@ -342,7 +230,7 @@ energy(const typename Basis::LocalView& localView,
for (int alpha=0; alpha<2; alpha++)
for (int beta=0; beta<2; beta++)
- c += sqrt(aScalar) * eps[alpha][beta] * dyadicProduct(aContravariant[alpha], aContravariant[beta]);
+ c += sqrt(aScalar) * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
// Second fundamental form
// The derivative of the normal field
@@ -354,14 +242,14 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldVector<double,3> vec;
for (int i=0; i<3; i++)
vec[i] = normalDerivative[i][alpha];
- b -= dyadicProduct(vec, aContravariant[alpha]);
+ b -= Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
// Gauss curvature
auto K = b.determinant();
// Mean curvatue
- auto H = 0.5 * trace(b);
+ auto H = 0.5 * Dune::GFE::trace(b);
//////////////////////////////////////////////////////////
// Strain tensors
@@ -375,7 +263,7 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldVector<field_type,3> vec;
for (int i=0; i<3; i++)
vec[i] = derivative[i][alpha];
- grad_s_m += dyadicProduct(vec, aContravariant[alpha]);
+ grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
Ee = RT * grad_s_m - a;
@@ -389,7 +277,7 @@ energy(const typename Basis::LocalView& localView,
for (int j=0; j<3; j++)
tmp[i][j] = DR[i][j][alpha];
auto tmp2 = RT * tmp;
- Ke += dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
+ Ke += Dune::GFE::dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
}
//////////////////////////////////////////////////////////
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index 805cd8ae..1aad21af 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -59,121 +59,6 @@ class SurfaceCosseratEnergy
enum {gridDim=GridView::dimension};
static constexpr int boundaryDim = gridDim - 1;
-
- /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] + A[j][i]);
- return result;
- }
-
- /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] - A[j][i]);
- return result;
- }
-
- /** \brief Compute the deviator of a matrix A */
- static Dune::FieldMatrix<field_type,dim,dim> dev(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result = A;
- auto t = trace(A);
- for (int i=0; i<dim; i++)
- result[i][i] -= t / dim;
- return result;
- }
-
- /** \brief Return the trace of a matrix */
- template <class T, int N>
- static T trace(const Dune::FieldMatrix<T,N,N>& A)
- {
- T trace = 0;
- for (int i=0; i<N; i++)
- trace += A[i][i];
- return trace;
- }
-
- /** \brief Return the square of the trace of a matrix */
- template <int N>
- static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
- {
- field_type trace = 0;
- for (int i=0; i<N; i++)
- trace += A[i][i];
- return trace*trace;
- }
-
- template <class T, int N>
- static T frobeniusProduct(const Dune::FieldMatrix<T,N,N>& A, const Dune::FieldMatrix<T,N,N>& B)
- {
- T result(0.0);
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result += A[i][j] * B[i][j];
-
- return result;
- }
-
- template <int N>
- static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<adouble,N>& B)
- -> Dune::FieldMatrix<adouble,N,N>
- {
- Dune::FieldMatrix<adouble,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[i]*B[j];
-
- return result;
- }
-
- template <int N>
- static auto dyadicProduct(const Dune::FieldVector<double,N>& A, const Dune::FieldVector<double,N>& B)
- -> Dune::FieldMatrix<double,N,N>
- {
- Dune::FieldMatrix<double,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[i]*B[j];
-
- return result;
- }
-
- template <int N>
- static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<double,N>& B)
- -> Dune::FieldMatrix<adouble,N,N>
- {
- Dune::FieldMatrix<adouble,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[i]*B[j];
-
- return result;
- }
-
- template <class T, int N>
- static Dune::FieldMatrix<T,N,N> transpose(const Dune::FieldMatrix<T,N,N>& A)
- {
- Dune::FieldMatrix<T,N,N> result;
-
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = A[j][i];
-
- return result;
- }
-
-
/** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
\param value Value of the gfe function at a certain point
\param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
@@ -313,7 +198,7 @@ RT energy(const typename Basis::LocalView& localView,
Dune::FieldMatrix<field_type,dim,dim> R;
value.q.matrix(R);
- auto RT = transpose(R);
+ auto RT = Dune::GFE::transpose(R);
Tensor3<field_type,3,3,boundaryDim> DR;
computeDR(value, derivative2D, DR);
@@ -346,7 +231,7 @@ RT energy(const typename Basis::LocalView& localView,
Dune::FieldMatrix<double,3,3> a(0);
for (int alpha=0; alpha<boundaryDim; alpha++)
- a += dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
+ a += Dune::GFE::dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
auto a00 = aCovariant[0] * aCovariant[0];
auto a01 = aCovariant[0] * aCovariant[1];
@@ -360,7 +245,7 @@ RT energy(const typename Basis::LocalView& localView,
for (int alpha=0; alpha<2; alpha++)
for (int beta=0; beta<2; beta++)
- c += sqrt(aScalar) * eps[alpha][beta] * dyadicProduct(aContravariant[alpha], aContravariant[beta]);
+ c += sqrt(aScalar) * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
// Second fundamental form
// The derivative of the normal field
@@ -372,14 +257,14 @@ RT energy(const typename Basis::LocalView& localView,
Dune::FieldVector<double,3> vec;
for (int i=0; i<3; i++)
vec[i] = normalDerivative[i][alpha];
- b -= dyadicProduct(vec, aContravariant[alpha]);
+ b -= Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
// Gauss curvature
auto K = b.determinant();
// Mean curvatue
- auto H = 0.5 * trace(b);
+ auto H = 0.5 * Dune::GFE::trace(b);
//////////////////////////////////////////////////////////
// Strain tensors
@@ -393,7 +278,7 @@ RT energy(const typename Basis::LocalView& localView,
Dune::FieldVector<field_type,3> vec;
for (int i=0; i<3; i++)
vec[i] = derivative[i][alpha];
- grad_s_m += dyadicProduct(vec, aContravariant[alpha]);
+ grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
Ee = RT * grad_s_m - a;
@@ -407,7 +292,7 @@ RT energy(const typename Basis::LocalView& localView,
for (int j=0; j<3; j++)
tmp[i][j] = DR[i][j][alpha];
auto tmp2 = RT * tmp;
- Ke += dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
+ Ke += Dune::GFE::dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
}
//////////////////////////////////////////////////////////
@@ -440,19 +325,20 @@ RT energy(const typename Basis::LocalView& localView,
RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T, double mu, double lambda) const
{
- return mu * frobeniusProduct(sym(S), sym(T))
- + mu_c_ * frobeniusProduct(skew(S), skew(T))
- + lambda * mu / (lambda + 2*mu) * trace(S) * trace(T);
+ return mu * Dune::GFE::frobeniusProduct(Dune::GFE::sym(S), Dune::GFE::sym(T))
+ + mu_c_ * Dune::GFE::frobeniusProduct(Dune::GFE::skew(S), Dune::GFE::skew(T))
+ + lambda * mu / (lambda + 2*mu) * Dune::GFE::trace(S) * Dune::GFE::trace(T);
}
RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
{
- return mu * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda * 0.5 * traceSquared(S);
+ return mu * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda * 0.5 * Dune::GFE::traceSquared(S);
}
RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S, double mu) const
{
- return mu * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
+ return mu * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
+ + b2_ * Dune::GFE::skew(S).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(S));
}
private:
--
GitLab