diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 2f1eb8b27af6b78dcc48e9a7f4dd14843df4965b..df4749eb9be221a20381c136f34bfa3218e99013 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -21,7 +21,13 @@
#define DONT_USE_CURL
-//#define QUADRATIC_MEMBRANE_ENERGY
+#define CURVATURE_WITH_WRYNESS //only relevant for gridDim == 3
+
+//#define QUADRATIC_2006 //only relevant for gridDim == 3
+// Use the formula for the quadratic membrane energy (first equation of (4.4)) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
+// instead of the equation (2.27) of 2019: Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature
+
+//#define QUADRATIC_MEMBRANE_ENERGY //only relevant for gridDim == 2
/** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
*
@@ -102,7 +108,7 @@ public:
neumannFunction_(neumannFunction),
volumeLoad_(volumeLoad)
{
- // The shell thickness
+ // The shell thickness // only relevant for dim == 2
thickness_ = parameters.template get<double>("thickness");
// Lame constants
@@ -118,8 +124,13 @@ public:
// Curvature exponent
q_ = parameters.template get<double>("q");
- // Shear correction factor
+ // Shear correction factor // only relevant for dim == 2
kappa_ = parameters.template get<double>("kappa");
+
+ // Curvature parameters
+ b1_ = parameters.template get<double>("b1");
+ b2_ = parameters.template get<double>("b2");
+ b3_ = parameters.template get<double>("b3");
}
/** \brief Assemble the energy for a single element */
@@ -133,6 +144,7 @@ public:
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the first equation of (4.4) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
+ * OR: the equation (2.27) of 2019: Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature
*/
RT quadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
@@ -142,7 +154,11 @@ public:
return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2()
+ mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+#ifdef QUADRATIC_2006
+ (mu_*lambda_)/(2*mu_ + lambda_) * Dune::GFE::traceSquared(UMinus1); // Dune::GFE::traceSquared(UMinus1) = Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1))
+#else
+ + lambda_/2 * Dune::GFE::traceSquared(UMinus1); // Dune::GFE::traceSquared(UMinus1) = Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1))
+#endif
}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
@@ -236,6 +252,36 @@ public:
#endif
}
+ RT curvatureWithWryness(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
+ {
+ // construct Wryness tensor \Gamma as in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
+ Dune::FieldMatrix<field_type,3,3> dRx1(0); //Derivative of R with respect to x1
+ Dune::FieldMatrix<field_type,3,3> dRx2(0); //Derivative of R with respect to x2
+ Dune::FieldMatrix<field_type,3,3> dRx3(0); //Derivative of R with respect to x3, this of course only exists if gridDim = 3 (then dim = 3 also, because dim >= gridDim, always)
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++) {
+ dRx1[i][j] = DR[i][j][0];
+ dRx2[i][j] = DR[i][j][1];
+ dRx3[i][j] = gridDim == 3 ? DR[i][j][2] : 0;
+ }
+
+ Dune::FieldMatrix<field_type,3,3> wryness(0);
+
+ auto axialVectorx1 = SkewMatrix<field_type,3>(transpose(R)*dRx1).axial();
+ auto axialVectorx2 = SkewMatrix<field_type,3>(transpose(R)*dRx2).axial();
+ auto axialVectorx3 = SkewMatrix<field_type,3>(transpose(R)*dRx3).axial();
+ for (int i=0; i<3; i++) {
+ wryness[i][0] = axialVectorx1[i];
+ wryness[i][1] = axialVectorx2[i];
+ wryness[i][2] = gridDim == 3 ? axialVectorx3[i] : 0;
+ }
+
+ // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
+ return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(wryness)).frobenius_norm2()
+ + b2_ * Dune::GFE::skew(wryness).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(wryness));
+ }
+
+
RT bendingEnergy(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
{
// left-multiply the derivative of the third director (in DR[][2][]) with R^T
@@ -268,6 +314,9 @@ public:
/** \brief Shear correction factor */
double kappa_;
+ /** \brief Curvature parameters */
+ double b1_, b2_, b3_;
+
/** \brief The Neumann boundary */
const BoundaryPatch<GridView>* neumannBoundary_;
@@ -357,7 +406,11 @@ energy(const typename Basis::LocalView& localView,
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
} else if (gridDim==3) {
energy += weight * quadraticMembraneEnergy(U);
+#ifdef CURVATURE_WITH_WRYNESS
+ energy += weight * curvatureWithWryness(R,DR);
+#else
energy += weight * curvatureEnergy(DR);
+#endif
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
@@ -511,7 +564,11 @@ energy(const typename Basis::LocalView& localView,
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
} else if (gridDim==3) {
energy += weight * quadraticMembraneEnergy(U);
+#ifdef CURVATURE_WITH_WRYNESS
+ energy += weight * curvatureWithWryness(R,DR);
+#else
energy += weight * curvatureEnergy(DR);
+#endif
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
diff --git a/test/adolctest-scalar-and-vector-mode.cc b/test/adolctest-scalar-and-vector-mode.cc
index 8990660ec830ff719bc9e135e6214f51ff7607dc..036ce5de75b07503afa64f379ec3a535dd5badce 100644
--- a/test/adolctest-scalar-and-vector-mode.cc
+++ b/test/adolctest-scalar-and-vector-mode.cc
@@ -108,6 +108,9 @@ int main (int argc, char *argv[])
parameters["L_c"] = "0.01";
parameters["q"] = "2";
parameters["kappa"] = "1";
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
//Mixed space
CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergyMixed(parameters,
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 2ae115f3157b50917d5e14b27f8c85676da3f665..9a408d0441756b551aef2b1fd5c27b8e47ce9912 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -490,6 +490,9 @@ int main (int argc, char *argv[]) try
materialParameters["L_c"] = "1";
materialParameters["q"] = "2";
materialParameters["kappa"] = "1";
+ materialParameters["b1"] = "1";
+ materialParameters["b2"] = "1";
+ materialParameters["b3"] = "1";
///////////////////////////////////////////////////////////////////////
// Assemblers for the Euclidean derivatives in an embedding space
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 9f4b5bc40e667f1355bf4531f6c25b5d889ecd44..1d41eff2f9fddc017348dc3b04261745ab9faa4f 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -61,6 +61,9 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
materialParameters["L_c"] = "0.1";
materialParameters["q"] = "2.5";
materialParameters["kappa"] = "0.1";
+ materialParameters["b1"] = "1";
+ materialParameters["b2"] = "1";
+ materialParameters["b3"] = "1";
typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> FEBasis;
FEBasis feBasis(grid->leafGridView());
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index 8030dd8365598d5ff9d91b2686eaa7c10c749612..173a5d7c72df493b6b97e882c19cf035145a4d39 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -129,6 +129,9 @@ int main (int argc, char *argv[])
parameters["L_c"] = "0.01";
parameters["q"] = "2";
parameters["kappa"] = "1";
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
FieldVector<double,dim> values_ = {3e4,2e4,1e4};