From df922afd6194367612bdeb5719960c6b09135393 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 22 May 2014 19:33:20 +0000
Subject: [PATCH] Use the new CosseratStrain object, but without any
optimizations yet
The CosseratStrain object is a 3x3 matrix, where the last column
is (0,0,1)^T. This special structure allows to use specially tuned
matrix methods (e.g., the determinant), which I intent to do in the
future.
[[Imported from SVN: r9752]]
---
dune/gfe/cosseratenergystiffness.hh | 41 ++++------------
dune/gfe/cosseratstrain.hh | 76 +++++++++++++++++++++++++++++
2 files changed, 86 insertions(+), 31 deletions(-)
create mode 100644 dune/gfe/cosseratstrain.hh
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 6998469b..a51ccd85 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -13,6 +13,7 @@
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/orthogonalmatrix.hh>
+#include <dune/gfe/cosseratstrain.hh>
#define DONT_USE_CURL
@@ -149,7 +150,7 @@ public:
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the first equation of (4.4) in Neff's paper
*/
- RT quadraticMembraneEnergy(const Dune::FieldMatrix<field_type,3,3>& U) const
+ RT quadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
Dune::FieldMatrix<field_type,3,3> UMinus1 = U;
for (int i=0; i<dim; i++)
@@ -163,7 +164,7 @@ public:
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the second equation of (4.4) in Neff's paper
*/
- RT longQuadraticMembraneEnergy(const Dune::FieldMatrix<field_type,3,3>& U) const
+ RT longQuadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
RT result = 0;
@@ -171,7 +172,7 @@ public:
Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
for (int i=0; i<dim-1; i++)
for (int j=0; j<dim-1; j++)
- sym2x2[i][j] = 0.5 * (U[i][j] + U[j][i]) - (i==j);
+ sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
result += mu_ * sym2x2.frobenius_norm2();
@@ -179,13 +180,13 @@ public:
Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
for (int i=0; i<dim-1; i++)
for (int j=0; j<dim-1; j++)
- skew2x2[i][j] = 0.5 * (U[i][j] - U[j][i]);
+ skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
result += mu_c_ * skew2x2.frobenius_norm2();
// classical transverse shear energy
- result += kappa_ * (mu_ + mu_c_)/2 * (U[2][0]*U[2][0] + U[2][1]*U[2][1]);
+ result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
// elongational stretch energy
result += mu_*lambda_ / (2*mu_ + lambda_) * traceSquared(sym2x2);
@@ -195,9 +196,9 @@ public:
/** \brief Energy for large-deformation problems (private communication by Patrizio Neff)
*/
- RT nonquadraticMembraneEnergy(const Dune::FieldMatrix<field_type,3,3>& U) const
+ RT nonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
- Dune::FieldMatrix<field_type,3,3> UMinus1 = U;
+ Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
@@ -312,29 +313,7 @@ energy(const Entity& element,
Dune::FieldMatrix<field_type,dim,dim> R;
value.q.matrix(R);
- /* Compute F, the deformation gradient.
- In the continuum case this is
- \f$ \hat{F} = \nabla m \f$
- In the case of a shell it is
- \f$ \hat{F} = (\nabla m | \overline{R}_3) \f$
- */
- Dune::FieldMatrix<field_type,dim,dim> F;
- for (int i=0; i<dim; i++)
- for (int j=0; j<gridDim; j++)
- F[i][j] = derivative[i][j];
-
- for (int i=0; i<dim; i++)
- for (int j=gridDim; j<dim; j++)
- F[i][j] = R[i][j];
-
- // U = R^T F
- Dune::FieldMatrix<field_type,dim,dim> U;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++) {
- U[i][j] = 0;
- for (int k=0; k<dim; k++)
- U[i][j] += R[k][i] * F[k][j];
- }
+ Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(derivative,R);
//////////////////////////////////////////////////////////
// Compute the derivative of the rotation
@@ -347,7 +326,7 @@ energy(const Entity& element,
// Add the local energy density
if (gridDim==2) {
#ifdef QUADRATIC_MEMBRANE_ENERGY
- //energy += weight * thickness_ * quadraticMembraneEnergy(U);
+ //energy += weight * thickness_ * quadraticMembraneEnergy(U.matrix());
energy += weight * thickness_ * longQuadraticMembraneEnergy(U);
#else
energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
diff --git a/dune/gfe/cosseratstrain.hh b/dune/gfe/cosseratstrain.hh
new file mode 100644
index 00000000..cc72745d
--- /dev/null
+++ b/dune/gfe/cosseratstrain.hh
@@ -0,0 +1,76 @@
+#ifndef DUNE_GFE_COSSERATSTRAIN_HH
+#define DUNE_GFE_COSSERATSTRAIN_HH
+
+#include <dune/common/fmatrix.hh>
+
+namespace Dune {
+
+ namespace GFE {
+
+ /** \brief Strain tensor of a Cosserat material
+ *
+ * Mathematically, this object is a dimworld x dimworld matrix.
+ * However, the last dimworld-dim columns consists of the corresponding
+ * columns of the identity matrix. Knowing this allows more efficient
+ * implementations of several frequently occurring operations.
+ *
+ * \tparam T Type used for real numbers
+ * \tparam dimworld Dimension of the world space
+ * \tparam dim Dimension of the domain
+ */
+ template <class T, int dimworld, int dim>
+ class CosseratStrain
+ {
+ public:
+ /** \brief Compute strain from the deformation gradient and the microrotation
+ *
+ * \param R The microrotation
+ */
+ CosseratStrain(const FieldMatrix<T,dimworld+4,dim>& deformationGradient,
+ const FieldMatrix<T,dimworld,dimworld>& R)
+ {
+ /* Compute F, the deformation gradient.
+ In the continuum case (domain dimension == world dimension) this is
+ \f$ \hat{F} = \nabla m \f$
+ In the case of a shell it is
+ \f$ \hat{F} = (\nabla m | \overline{R}_3) \f$
+ */
+ FieldMatrix<T,dimworld,dimworld> F;
+ for (int i=0; i<dimworld; i++)
+ for (int j=0; j<dim; j++)
+ F[i][j] = deformationGradient[i][j];
+
+ for (int i=0; i<dimworld; i++)
+ for (int j=dim; j<dimworld; j++)
+ F[i][j] = R[i][j];
+
+ // U = R^T F
+ for (int i=0; i<dimworld; i++)
+ for (int j=0; j<dimworld; j++) {
+ data_[i][j] = 0;
+ for (int k=0; k<dimworld; k++)
+ data_[i][j] += R[k][i] * F[k][j];
+ }
+
+ }
+
+ T determinant() const
+ {
+ return data_.determinant();
+ }
+
+ /** \brief Temporary: get the raw matrix data */
+ const FieldMatrix<T,dimworld,dimworld>& matrix() const
+ {
+ return data_;
+ }
+
+ private:
+ FieldMatrix<T,dimworld,dimworld> data_;
+ };
+
+ } // namespace GFE
+
+} // namespace Dune
+
+#endif // DUNE_GFE_COSSERATSTRAIN_HH
\ No newline at end of file
--
GitLab