diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index e6a0bb1aec297aa3090c13d648915e58eb055374..101c37b9986b73c96746603e2e8014d9bf130d7d 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -213,18 +213,15 @@ public:
SurfaceCosseratEnergy(const Dune::ParameterTree& parameters,
const std::vector<UnitVector<double,3> >& vertexNormals,
const BoundaryPatch<GridView>* shellBoundary,
- const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dim-1, dim>>& geometriesOnShellBoundary)
+ const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dim-1, dim>>& geometriesOnShellBoundary,
+ const std::function<double(Dune::FieldVector<double,dim>)> thicknessF,
+ const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF)
: shellBoundary_(shellBoundary),
vertexNormals_(vertexNormals),
- geometriesOnShellBoundary_(geometriesOnShellBoundary)
+ geometriesOnShellBoundary_(geometriesOnShellBoundary),
+ thicknessF_(thicknessF),
+ lameF_(lameF)
{
- // The shell thickness
- thickness_ = parameters.template get<double>("thickness");
-
- // Lame constants
- mu_ = parameters.template get<double>("mu_cosserat");
- lambda_ = parameters.template get<double>("lambda_cosserat");
-
// Cosserat couple modulus
mu_c_ = parameters.template get<double>("mu_c");
@@ -293,6 +290,13 @@ RT energy(const typename Basis::LocalView& localView,
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());;
+ // Global position of the quadrature point
+ auto quadPosGlobal = it.geometry().global(quad[pt].position());
+ double thickness = thicknessF_(quadPosGlobal);
+ auto lameConstants = lameF_(quadPosGlobal);
+ double mu = lameConstants[0];
+ double lambda = lameConstants[1];
+
const DT integrationElement = boundaryGeometry.integrationElement(quad[pt].position());
// The value of the local function
@@ -416,15 +420,15 @@ RT energy(const typename Basis::LocalView& localView,
//////////////////////////////////////////////////////////
// Add the membrane energy density
- auto energyDensity = (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_m(Ee);
- energyDensity += (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_m(Ee*b + c*Ke);
- energyDensity += Dune::Power<3>::eval(thickness_) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke);
- energyDensity += Dune::Power<5>::eval(thickness_) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
+ auto energyDensity = (thickness - K*Dune::Power<3>::eval(thickness) / 12.0) * W_m(Ee, mu, lambda);
+ energyDensity += (Dune::Power<3>::eval(thickness) / 12.0 - K * Dune::Power<5>::eval(thickness) / 80.0)*W_m(Ee*b + c*Ke, mu, lambda);
+ energyDensity += Dune::Power<3>::eval(thickness) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke, mu, lambda);
+ energyDensity += Dune::Power<5>::eval(thickness) / 80.0 * W_mp( (Ee*b + c*Ke)*b, mu, lambda);
// Add the bending energy density
- energyDensity += (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_curv(Ke)
- + (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_curv(Ke*b)
- + Dune::Power<5>::eval(thickness_) / 80.0 * W_curv(Ke*b*b);
+ energyDensity += (thickness - K*Dune::Power<3>::eval(thickness) / 12.0) * W_curv(Ke, mu)
+ + (Dune::Power<3>::eval(thickness) / 12.0 - K * Dune::Power<5>::eval(thickness) / 80.0)*W_curv(Ke*b, mu)
+ + Dune::Power<5>::eval(thickness) / 80.0 * W_curv(Ke*b*b, mu);
// Add energy density
energy += quad[pt].weight() * integrationElement * energyDensity;
@@ -434,26 +438,26 @@ RT energy(const typename Basis::LocalView& localView,
return energy;
}
- RT W_m(const Dune::FieldMatrix<field_type,3,3>& S) const
+ RT W_m(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
{
- return W_mixt(S,S);
+ return W_mixt(S,S, mu, lambda);
}
- RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T) const
+ 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))
+ return mu * frobeniusProduct(sym(S), sym(T))
+ mu_c_ * frobeniusProduct(skew(S), skew(T))
- + lambda_ * mu_ / (lambda_ + 2*mu_) * trace(S) * trace(T);
+ + lambda * mu / (lambda + 2*mu) * trace(S) * trace(T);
}
- RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S) const
+ 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 * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda * 0.5 * traceSquared(S);
}
- RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S) const
+ 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_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
}
private:
@@ -463,11 +467,14 @@ private:
/** \brief Stress-free geometries of the shell elements*/
const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dim-1, dim>> geometriesOnShellBoundary_;
- /** \brief The normal vectors at the grid vertices. This are used to compute the reference surface curvature. */
+ /** \brief The normal vectors at the grid vertices. They are used to compute the reference surface curvature. */
std::vector<UnitVector<double,3> > vertexNormals_;
- /** \brief The shell thickness */
- double thickness_;
+ /** \brief The shell thickness as a function*/
+ std::function<double(Dune::FieldVector<double,dim>)> thicknessF_;
+
+ /** \brief The Lamé-parameters as a function*/
+ std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF_;
/** \brief Lame constants */
double mu_, lambda_;
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 5e9d376a1c8f4f26ed1639f871097e732c6b8bd0..471ea478564f912ead8f141f3acc0c03caf597d9 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -434,6 +434,14 @@ int main (int argc, char *argv[]) try
}
}
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+ Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
+ Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
+
+ Python::Reference pythonObject = surfaceShellCallable();
+ PythonFunction<Dune::FieldVector<double, dim>, double> fThickness(pythonObject.get("thickness"));
+ PythonFunction<Dune::FieldVector<double, dim>, Dune::FieldVector<double, 2>> fLame(pythonObject.get("lame"));
+
for (int i=0; i<numHomotopySteps; i++)
{
double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
@@ -442,7 +450,6 @@ int main (int argc, char *argv[]) try
// Create an assembler for the energy functional
// ////////////////////////////////////////////////////////////
- const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
std::shared_ptr<NeumannFunction> neumannFunction;
@@ -547,7 +554,7 @@ int main (int argc, char *argv[]) try
vertexNormals[i] = vertexNormal;
}
- surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary));
+ surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
std::shared_ptr<LocalEnergyBase> totalEnergy;
totalEnergy = std::make_shared<GFE::SumCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble>> (elasticAndNeumann, surfaceCosseratEnergy);
diff --git a/src/film-on-substrate.parset b/src/film-on-substrate.parset
index ccfdc893004623560cf8cdb141aabce01e6c3949..645c408bf9efca144da5cdc2c0dfe1b101901442 100644
--- a/src/film-on-substrate.parset
+++ b/src/film-on-substrate.parset
@@ -97,16 +97,12 @@ energy = stvenantkirchhoff
## For the Wriggers L-shape example
[materialParameters]
-# shell thickness
-thickness = 0.6
+surfaceShellParameters = surface-shell-parameters-1-3
## Lame parameters for stvenantkirchhoff, E = mu(3*lambda + 2*mu)/(lambda + mu)
# mu = 2.7191e+4
# lambda = 4.4364e+4
-#Lame parameters for the cosserat material
-mu_cosserat = 7.64e+6 #Parameters for E = 20.4 MPa, nu=0.335
-lambda_cosserat = 1.55e+7
# Cosserat couple modulus
mu_c = 0
@@ -144,5 +140,4 @@ mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute th
# log: Generalized Rivlin model or polynomial hyperelastic model, using 0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
# square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
-[]
-
+[]
\ No newline at end of file
diff --git a/src/surface-shell-parameters-1-3.py b/src/surface-shell-parameters-1-3.py
new file mode 100644
index 0000000000000000000000000000000000000000..17ebcee453e862d3fcfd9fd9401f0404dd68e4a0
--- /dev/null
+++ b/src/surface-shell-parameters-1-3.py
@@ -0,0 +1,18 @@
+import math
+
+class SurfaceShellParameters:
+ def thickness(self, x):
+ if (x[1] > 50.0):
+ return 1.0
+ else:
+ return 3.0
+
+ #Lame parameters for the surface shell material:
+ #Parameters for E = 20.4 MPa, nu=0.335
+ #mu_cosserat = 7.64e+6
+ #lambda_cosserat = 1.55e+7
+ def lame(self, x):
+ if (x[1] > 50.0):
+ return [7.64e+6,1.55e+7]
+ else:
+ return [7.64e+6,1.55e+7]