diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index aafdca5a6983154dd90f452eb894f506281dca92..4e5fdbc0ecc069ef52a5b62d7e7b03580f621f91 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 7febb66e833af845024d8d7bfb11d6a5df28a9e6..a751897f8cca8a0fdab4fbd1ad63fcac7ad4e952 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 16c97c4f0cd595ab06f7bd550aa61622a4a02ab4..102115dbc0a371ce6eb83a33e2e5dd739bd031d3 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 805cd8ae3eb67b1d176cef37f0cc0749879d45dd..1aad21afc16e0769232786189914aa4772a8305a 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: