diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7431b7d0fb4d1f93c021fb256f55d2ecfcf5f029..dd48f7bd05b9f1db1ada9208d075213efda0ebf4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
# Master
+- In the `BulkCosseratDensity` class: Make the type of curvature
+ tensor controllable at run-time (previously, ugly preprocessor switches
+ where used). Among its parameters `BulkCosseratDensity` now expects
+ a string `curvatureType`, which has to take one of the values `norm`,
+ `curl`, or `wryness`.
+
- Building `dune-gfe` now requires Dune 2.9 or newer. Older configurations
have not gotten CI-tested for quite a while, anyway.
diff --git a/dune/gfe/densities/bulkcosseratdensity.hh b/dune/gfe/densities/bulkcosseratdensity.hh
index c37a7d23d8b0496e10b73c18fced82a6a02e9704..6237dd9a0d9a00a665ff80fd06181abe955066ea 100644
--- a/dune/gfe/densities/bulkcosseratdensity.hh
+++ b/dune/gfe/densities/bulkcosseratdensity.hh
@@ -1,9 +1,6 @@
#ifndef DUNE_GFE_DENSITIES_BULKCOSSERATDENSITY_HH
#define DUNE_GFE_DENSITIES_BULKCOSSERATDENSITY_HH
-#define DONT_USE_CURL
-#define CURVATURE_WITH_WRYNESS
-
#include <dune/common/fmatrix.hh>
#include <dune/gfe/cosseratstrain.hh>
@@ -80,6 +77,9 @@ namespace Dune::GFE {
public:
+ /** \brief The different types of curvature tensors */
+ enum CurvatureType {Norm, Curl, Wryness};
+
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
@@ -94,6 +94,17 @@ namespace Dune::GFE {
// Length scale parameter
L_c_ = parameters.template get<double>("L_c");
+ // Curvature type
+ const auto& curvatureString = parameters.template get<std::string>("curvatureType");
+ if (curvatureString=="norm")
+ curvatureType_ = CurvatureType::Norm;
+ else if (curvatureString=="curl")
+ curvatureType_ = CurvatureType::Curl;
+ else if (curvatureString=="wryness")
+ curvatureType_ = CurvatureType::Wryness;
+ else
+ DUNE_THROW(NotImplemented, "Curvature type '" << curvatureString << "' is not implemented!");
+
// Curvature exponent
q_ = parameters.template get<double>("q");
@@ -105,8 +116,8 @@ namespace Dune::GFE {
/** \brief Constructor with a set of material parameters
*/
- BulkCosseratDensity(double mu, double lambda, double mu_c, double L_c, double q, const std::array<double,3>& b)
- : mu_(mu), lambda_(lambda), mu_c_(mu_c), L_c_(L_c), q_(q), b1_(b[0]), b2_(b[1]), b3_(b[2])
+ BulkCosseratDensity(double mu, double lambda, double mu_c, double L_c, CurvatureType curvatureType, double q, const std::array<double,3>& b)
+ : mu_(mu), lambda_(lambda), mu_c_(mu_c), L_c_(L_c), curvatureType_(curvatureType), q_(q), b1_(b[0]), b2_(b[1]), b3_(b[2])
{}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
@@ -159,30 +170,54 @@ namespace Dune::GFE {
////////////////////////////////////////////////////////////////////
Tensor3<field_type,3,3,gridDim> DR = rotationValue.quaternionTangentToMatrixTangent(orientationDerivative);
- strainEnergyDensity = quadraticEnergy(U);
-
- #ifdef CURVATURE_WITH_WRYNESS
- auto curvatureTensor = wryness(R,DR);
-
- // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
- auto norm = (b1_ * GFE::dev(GFE::sym(curvatureTensor)).frobenius_norm2()
- + b2_ * GFE::skew(curvatureTensor).frobenius_norm2()
- + b3_ * GFE::traceSquared(curvatureTensor));
-
- strainEnergyDensity += mu_ * L_c_ * L_c_ * norm;
- #else
- #ifdef DONT_USE_CURL
- auto argument = DR;
- #else
- auto argument = curl(DR);
- #endif
+ // Compute the bending energy density
+ strainEnergyDensity = quadraticEnergy(U);
- auto norm = (b1_ * GFE::dev(GFE::sym(argument)).frobenius_norm2()
- + b2_ * GFE::skew(argument).frobenius_norm2() + b3_ * GFE::traceSquared(argument));
+ // Compute the curvature energy density
+ field_type norm = 0;
+
+ switch (curvatureType_)
+ {
+ case CurvatureType::Norm :
+ {
+ /* This is a simplified version of the curvature energy that appears in
+ *
+ * P. Neff: "A geometrically exact Cosserat shell-model including size effects,
+ * avoiding degeneracy in the thin shell limit. Part I: Formal dimensional reduction
+ * for elastic plates and existence of minimizers for positive Cosserat couple modulus"
+ * Continuum Mech. Thermodyn. (2004) 16: 577–628
+ * DOI: 10.1007/s00161-004-0182-4
+ *
+ * The curvature energy given there is much more complicated, but it
+ * does use the full norm of DR to define the curvature.
+ */
+ norm = DR.frobenius_norm2();
+ break;
+ }
+ case CurvatureType::Curl :
+ {
+ auto curvatureTensor = curl(DR);
+
+ // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
+ norm = (b1_ * GFE::dev(GFE::sym(curvatureTensor)).frobenius_norm2()
+ + b2_ * GFE::skew(curvatureTensor).frobenius_norm2()
+ + b3_ * GFE::traceSquared(curvatureTensor));
+ break;
+ }
+ case CurvatureType::Wryness :
+ {
+ auto curvatureTensor = wryness(R,DR);
+
+ // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
+ norm = (b1_ * GFE::dev(GFE::sym(curvatureTensor)).frobenius_norm2()
+ + b2_ * GFE::skew(curvatureTensor).frobenius_norm2()
+ + b3_ * GFE::traceSquared(curvatureTensor));
+ break;
+ }
+ }
strainEnergyDensity += mu_ * std::pow(L_c_ * L_c_ * norm,q_/2.0);
- #endif
return strainEnergyDensity;
}
@@ -190,7 +225,11 @@ namespace Dune::GFE {
// Construct a copy of this density but using 'adouble' as the number type
virtual std::unique_ptr<LocalDensity<Position,ATargetSpace> > makeActiveDensity() const
{
- return std::make_unique<BulkCosseratDensity<Position,adouble> >(mu_, lambda_, mu_c_, L_c_, q_, std::array<double,3>{b1_, b2_, b3_});
+ // curvatureType_ is a local enum type, and therefore its type changes
+ // together with the type of the surrounding class.
+ auto activeCurvatureType = static_cast<typename BulkCosseratDensity<Position,adouble>::CurvatureType>(curvatureType_);
+
+ return std::make_unique<BulkCosseratDensity<Position,adouble> >(mu_, lambda_, mu_c_, L_c_, activeCurvatureType, q_, std::array<double,3>{b1_, b2_, b3_});
}
/** \brief The density depends on the microrotation value, but not on the deformation
@@ -216,6 +255,9 @@ namespace Dune::GFE {
/** \brief Length scale parameter */
double L_c_;
+ /** \brief The precise way to compute the curvature tensor */
+ enum CurvatureType curvatureType_;
+
/** \brief Curvature exponent */
double q_;
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index bcc99d9ac31295d44987116dd40799d425f5c210..f09fef8116d13bf166c06b892754ebc7bbd97164 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -188,6 +188,7 @@ int main (int argc, char *argv[])
parameters["lambda"] = "1";
parameters["mu_c"] = "1";
parameters["L_c"] = "0.1";
+ parameters["curvatureType"] = "wryness";
parameters["q"] = "2";
parameters["b1"] = "1";
parameters["b2"] = "1";
diff --git a/test/localintegralstiffnesstest.cc b/test/localintegralstiffnesstest.cc
index 57ca582cc8c6a61e8cd269bed1135989371a8e40..156c3d0d44c054a6a75689d0ae8489798c3d2cde 100644
--- a/test/localintegralstiffnesstest.cc
+++ b/test/localintegralstiffnesstest.cc
@@ -288,6 +288,7 @@ int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
parameters["lambda"] = "1";
parameters["mu_c"] = "1";
parameters["L_c"] = "0.1";
+ parameters["curvatureType"] = "wryness";
parameters["q"] = "2";
parameters["kappa"] = "1";
parameters["b1"] = "1";