diff --git a/CHANGELOG.md b/CHANGELOG.md
index b7a71f58d5c52cc4112fddb1bc74a5f80b2965f9..5f2f6dd8e3609ae8df829bd6d8633349c52af710 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,10 @@
# Master
+- Do not scale the density functions with the thickness to avoid confusion
+ since some densities need to be scaled and some do not need to be scaled
+ with the thickness depending on the dimension of the grid, their direction
+ and their kind (Neumann or volume load).
+
- Fix bug in the `RealTuple::log` method: Calling `log(a,b)`returned `a-b`
instead of `b-a`.
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 6f308de6207222da4808764ee1089b9b0489d9a8..2f7c2523b039144999a72dc051a80e406d4eb80c 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -99,7 +99,8 @@ public:
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
: neumannBoundary_(neumannBoundary),
- neumannFunction_(neumannFunction)
+ neumannFunction_(neumannFunction),
+ volumeLoad_(volumeLoad)
{
// The shell thickness
thickness_ = parameters.template get<double>("thickness");
@@ -371,7 +372,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
- energy += thickness_ * (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
}
@@ -405,7 +406,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
- energy += thickness_ * (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
}
@@ -513,6 +514,20 @@ energy(const typename Basis::LocalView& localView,
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
+ ///////////////////////////////////////////////////////////
+ // Volume load contribution
+ ///////////////////////////////////////////////////////////
+ if (not volumeLoad_)
+ continue;
+
+ // Value of the volume load density at the current position
+ auto volumeLoadDensity = volumeLoad_(element.geometry().global(quad[pt].position()));
+
+
+ // Only translational dofs are affected by the volume load
+ for (size_t i=0; i<volumeLoadDensity.size(); i++)
+ energy += (volumeLoadDensity[i] * deformationValue[i]) * quad[pt].weight() * integrationElement;
+
}
//////////////////////////////////////////////////////////////////////////////
@@ -544,7 +559,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
- energy += thickness_ * (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
+ energy += (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index e8a3678eed7a6c7880bea11b2a400e9ab451bb53..ee86f785219561c18e862e4ba0d0601b3e46c165 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -80,11 +80,33 @@ public:
b1_ = parameters.template get<double>("b1");
b2_ = parameters.template get<double>("b2");
b3_ = parameters.template get<double>("b3");
+
+ // Indicator to use the alternative energy W_Coss from Birsan 2021:
+ useAlternativeEnergyWCoss_ = parameters.template get<bool>("useAlternativeEnergyWCoss", false);
}
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const;
+ /* Sources:
+ Birsan 2019: Derivation of a refined six-parameter shell model, equation (111)
+ Birsan 2021: Alternative derivation of the higher-order constitudtive model for six-parameter elastic shells, equations (119) and (126)
+ */
+ RT W_Coss(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<double,3,3>& a, const Dune::FieldVector<double,3>& n0) const
+ {
+ return W_Coss_mixt(S,S,a,n0);
+ }
+
+ RT W_Coss_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T, const Dune::FieldMatrix<double,3,3>& a, const Dune::FieldVector<double,3>& n0) const
+ {
+ auto planarPart = W_mixt(a*S,a*T);
+ Dune::FieldVector<field_type, 3> n0S;
+ Dune::FieldVector<field_type, 3> n0T;
+ S.mtv(n0, n0S);
+ T.mtv(n0, n0T);
+ field_type normalPart = 2*mu_*mu_c_* n0S * n0T /(mu_ + mu_c_);
+ return planarPart + normalPart;
+ }
RT W_m(const Dune::FieldMatrix<field_type,3,3>& S) const
{
@@ -124,6 +146,10 @@ public:
/** \brief Curvature parameters */
double b1_, b2_, b3_;
+ /** \brief Indicator to use the alternative energy W_Coss from Birsan 2021:
+ Alternative derivation of the higher-order constitudtive model for six-parameter elastic shells, equations (119) and (126). */
+
+ bool useAlternativeEnergyWCoss_;
/** \brief The geometry used for assembling */
const StressFreeStateGridFunction* stressFreeStateGridFunction_;
@@ -307,10 +333,18 @@ energy(const typename Basis::LocalView& localView,
//////////////////////////////////////////////////////////
// Add the membrane energy density
- auto energyDensity = (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_m(Ee)
- + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_m(Ee*b + c*Ke)
- + Dune::power(thickness_,3) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
- + Dune::power(thickness_,5) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
+ field_type energyDensity = 0;
+ if (useAlternativeEnergyWCoss_) {
+ energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_Coss(Ee, a, aContravariant[2])
+ + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_Coss(Ee*b + c*Ke, a, aContravariant[2])
+ + Dune::power(thickness_,3) / 6.0 * W_Coss_mixt(Ee, c*Ke*b - 2*H*c*Ke, a, aContravariant[2])
+ + Dune::power(thickness_,5) / 80.0 * W_Coss( (Ee*b + c*Ke)*b, a, aContravariant[2]);
+ } else {
+ energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_m(Ee)
+ + (Dune::power(thickness_,3) / 12.0 - K * Dune::power(thickness_,5) / 80.0)*W_m(Ee*b + c*Ke)
+ + Dune::power(thickness_,3) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
+ + Dune::power(thickness_,5) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
+ }
// Add the bending energy density
energyDensity += (thickness_ - K*Dune::power(thickness_,3) / 12.0) * W_curv(Ke)
@@ -332,7 +366,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
- energy += thickness_ * (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
}
@@ -365,7 +399,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
- energy += thickness_ * (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
}
}
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index cc0ea99152d42c3817531764e1f633f279721d74..3e5666ae9a5fc498aa0cb19eab33654eb2476db7 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -199,11 +199,11 @@ int main (int argc, char *argv[]) try
using namespace Dune::Functions::BasisFactory;
- const int dimRotation = Rotation<double,dim>::embeddedDim;
+ const int dimRotation = Rotation<double,3>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
composite(
- power<dim>(
+ power<3>(
lagrange<displacementOrder>()
),
power<dimRotation>(
@@ -403,15 +403,14 @@ int main (int argc, char *argv[]) try
return nV;
};
- FieldVector<double,3> volumeLoadValues {0,0,0};
- if (parameterSet.hasKey("volumeLoad"))
- volumeLoadValues = parameterSet.get<FieldVector<double,3> >("volumeLoad");
+ Python::Reference volumeLoadClass = Python::import(parameterSet.get<std::string>("volumeLoadPythonFunction", "zero-volume-load"));
+ Python::Callable volumeLoadCallable = volumeLoadClass.get("VolumeLoad");
- auto volumeLoad = [&]( FieldVector<double,dimworld>) {
- auto vL = volumeLoadValues;
- vL *= (-homotopyParameter);
- return vL;
- };
+ // Call a constructor
+ Python::Reference volumeLoadPythonObject = volumeLoadCallable(homotopyParameter);
+
+ // Extract object member functions as Dune functions
+ auto volumeLoad = Python::make_function<FieldVector<double,3> > (volumeLoadPythonObject.get("volumeLoad"));
if (mpiHelper.rank() == 0) {
std::cout << "Material parameters:" << std::endl;
diff --git a/test/adolctest-scalar-and-vector-mode.cc b/test/adolctest-scalar-and-vector-mode.cc
index 2f4fffce180161186a8c908ef764d0dd940b61a4..8990660ec830ff719bc9e135e6214f51ff7607dc 100644
--- a/test/adolctest-scalar-and-vector-mode.cc
+++ b/test/adolctest-scalar-and-vector-mode.cc
@@ -82,17 +82,6 @@ int main (int argc, char *argv[])
using namespace Functions::BasisFactory;
- auto compositeBasisMixed = makeBasis(
- gridView,
- composite(
- power<dim>(
- lagrange<2>()
- ),
- power<dim>(
- lagrange<1>()
- )
- ));
-
auto compositeBasis = makeBasis(
gridView,
composite(