From 7550c4a71098f4b51dbce2231d31d5bd9f906c50 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Tue, 14 May 2019 17:36:30 +0200
Subject: [PATCH] Add film-on-substrate
Added program to solve a combined minimization problem:
Minmize the sum of an elastic energy (representing the substrate)
and a cosserat energy (representing a thin film on the substrate).
---
dune/gfe/sumcosseratenergy.hh | 61 +++
dune/gfe/surfacecosseratenergy.hh | 599 +++++++++++++++++++++++++++++
src/cantilever-dirichlet-values.py | 17 +
src/film-on-substrate.cc | 482 +++++++++++++++++++++++
src/film-on-substrate.parset | 123 ++++++
5 files changed, 1282 insertions(+)
create mode 100644 dune/gfe/sumcosseratenergy.hh
create mode 100644 dune/gfe/surfacecosseratenergy.hh
create mode 100644 src/cantilever-dirichlet-values.py
create mode 100644 src/film-on-substrate.cc
create mode 100644 src/film-on-substrate.parset
diff --git a/dune/gfe/sumcosseratenergy.hh b/dune/gfe/sumcosseratenergy.hh
new file mode 100644
index 00000000..0f146a84
--- /dev/null
+++ b/dune/gfe/sumcosseratenergy.hh
@@ -0,0 +1,61 @@
+#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
+#define DUNE_GFE_SUMCOSSERATENERGY_HH
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+#include <dune/gfe/localgeodesicfestiffness.hh>
+
+namespace Dune {
+
+namespace GFE {
+template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
+class SumCosseratEnergy
+: public LocalGeodesicFEStiffness<Basis, TargetSpace>
+{
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype ctype;
+ typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
+ typedef typename GridView::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+
+ enum {dim=GridView::dimension};
+
+public:
+
+ /** \brief Constructor with elastic energy and cosserat energy
+ * \param elasticEnergy The elastic energy
+ * \param cosseratEnergy The cosserat energy
+ */
+ SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > , std::vector<GradientRT> > > elasticEnergy,
+ std::shared_ptr<LocalGeodesicFEStiffness<Basis, TargetSpace>> cosseratEnergy)
+
+ : elasticEnergy_(elasticEnergy),
+ cosseratEnergy_(cosseratEnergy)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localSolution) const
+ {
+ auto&& element = localView.element();
+ auto&& localFiniteElement = localView.tree().finiteElement();
+ std::vector<Dune::FieldVector<field_type,dim>> localConfiguration(localSolution.size());
+ for (int i = 0; i < localSolution.size(); i++) {
+ localConfiguration[i] = localSolution[i].r;
+ }
+ return elasticEnergy_->energy(element, localFiniteElement, localConfiguration) + cosseratEnergy_->energy(localView, localSolution);
+ }
+
+private:
+
+ std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> >, std::vector<GradientRT> > > elasticEnergy_;
+
+ std::shared_ptr<LocalGeodesicFEStiffness<Basis, TargetSpace> > cosseratEnergy_;
+};
+
+} // namespace GFE
+
+} // namespace Dune
+
+#endif //#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
new file mode 100644
index 00000000..905ab45d
--- /dev/null
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -0,0 +1,599 @@
+#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
+#define DUNE_GFE_SURFACECOSSERATENERGY_HH
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/gfe/cosseratstrain.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
+#include <dune/gfe/orthogonalmatrix.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/tensor3.hh>
+#include <dune/gfe/vertexnormals.hh>
+
+template <class Basis, std::size_t i>
+class LocalFiniteElementFactory
+{
+public:
+ static auto get(const typename Basis::LocalView& localView,
+ std::integral_constant<std::size_t, i> iType)
+ -> decltype(localView.tree().child(iType).finiteElement())
+ {
+ return localView.tree().child(iType).finiteElement();
+ }
+};
+
+/** \brief Specialize for scalar bases, here we cannot call tree().child() */
+template <class GridView, int order, std::size_t i>
+class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView,order>,i>
+{
+public:
+ static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView,
+ std::integral_constant<std::size_t, i> iType)
+ -> decltype(localView.tree().finiteElement())
+ {
+ return localView.tree().finiteElement();
+ }
+};
+
+template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
+class SurfaceCosseratEnergy
+: public LocalGeodesicFEStiffness<Basis,TargetSpace>
+{
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype ctype;
+ typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
+ typedef typename GridView::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+ typedef typename GridView::template Codim<0>::Entity Entity;
+ typedef typename TargetSpace::template rebind<adouble>::other ATargetSpace;
+
+ enum {dim=GridView::dimension};
+ enum {dimworld=GridView::dimensionworld};
+ enum {gridDim=GridView::dimension};
+ static constexpr int boundaryDim = gridDim - 1;
+
+
+ /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
+ static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
+ {
+ Dune::FieldMatrix<field_type,dim,dim> result;
+ for (int i=0; i<dim; i++)
+ for (int j=0; j<dim; j++)
+ result[i][j] = 0.5 * (A[i][j] + A[j][i]);
+ return result;
+ }
+
+ /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
+ static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
+ {
+ Dune::FieldMatrix<field_type,dim,dim> result;
+ for (int i=0; i<dim; i++)
+ for (int j=0; j<dim; j++)
+ result[i][j] = 0.5 * (A[i][j] - A[j][i]);
+ return result;
+ }
+
+ /** \brief Compute the deviator of a matrix A */
+ static Dune::FieldMatrix<field_type,dim,dim> dev(const Dune::FieldMatrix<field_type,dim,dim>& A)
+ {
+ Dune::FieldMatrix<field_type,dim,dim> result = A;
+ auto t = trace(A);
+ for (int i=0; i<dim; i++)
+ result[i][i] -= t / dim;
+ return result;
+ }
+
+ /** \brief Return the trace of a matrix */
+ template <class T, int N>
+ static T trace(const Dune::FieldMatrix<T,N,N>& A)
+ {
+ T trace = 0;
+ for (int i=0; i<N; i++)
+ trace += A[i][i];
+ return trace;
+ }
+
+ /** \brief Return the square of the trace of a matrix */
+ template <int N>
+ static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
+ {
+ field_type trace = 0;
+ for (int i=0; i<N; i++)
+ trace += A[i][i];
+ return trace*trace;
+ }
+
+ template <class T, int N>
+ static T frobeniusProduct(const Dune::FieldMatrix<T,N,N>& A, const Dune::FieldMatrix<T,N,N>& B)
+ {
+ T result(0.0);
+
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ result += A[i][j] * B[i][j];
+
+ return result;
+ }
+
+ template <int N>
+ static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<adouble,N>& B)
+ -> Dune::FieldMatrix<adouble,N,N>
+ {
+ Dune::FieldMatrix<adouble,N,N> result;
+
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ result[i][j] = A[i]*B[j];
+
+ return result;
+ }
+
+ template <int N>
+ static auto dyadicProduct(const Dune::FieldVector<double,N>& A, const Dune::FieldVector<double,N>& B)
+ -> Dune::FieldMatrix<double,N,N>
+ {
+ Dune::FieldMatrix<double,N,N> result;
+
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ result[i][j] = A[i]*B[j];
+
+ return result;
+ }
+
+ template <int N>
+ static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<double,N>& B)
+ -> Dune::FieldMatrix<adouble,N,N>
+ {
+ Dune::FieldMatrix<adouble,N,N> result;
+
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ result[i][j] = A[i]*B[j];
+
+ return result;
+ }
+
+ template <class T, int N>
+ static Dune::FieldMatrix<T,N,N> transpose(const Dune::FieldMatrix<T,N,N>& A)
+ {
+ Dune::FieldMatrix<T,N,N> result;
+
+ for (int i=0; i<N; i++)
+ for (int j=0; j<N; j++)
+ result[i][j] = A[j][i];
+
+ return result;
+ }
+
+
+ /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
+ \param value Value of the gfe function at a certain point
+ \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
+ \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
+ */
+ static void computeDR(const RigidBodyMotion<field_type,dimworld>& value,
+ const Dune::FieldMatrix<field_type,7,boundaryDim>& derivative,
+ Tensor3<field_type,3,3,boundaryDim>& DR)
+ {
+ // The LocalGFEFunction class gives us the derivatives of the orientation variable,
+ // but as a map into quaternion space. To obtain matrix coordinates we use the
+ // chain rule, which means that we have to multiply the given derivative with
+ // the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
+ // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
+ // However, since the directors of a given unit quaternion are the _columns_ of the
+ // corresponding orthogonal matrix, we need to invert the i and j indices
+ //
+ // So, DR[i][j][k] contains \partial R_ij / \partial k
+ Tensor3<field_type,3 , 3, 4> dd_dq;
+ value.q.getFirstDerivativesOfDirectors(dd_dq);
+
+ DR = field_type(0);
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<boundaryDim; k++)
+ for (int l=0; l<4; l++)
+ DR[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
+ }
+
+
+
+public:
+
+ /** \brief Constructor with a set of material parameters
+ * \param parameters The material parameters
+ */
+ SurfaceCosseratEnergy(const Dune::ParameterTree& parameters, const std::vector<UnitVector<double,3> >& vertexNormals, const BoundaryPatch<GridView>* shellBoundary)
+ : shellBoundary_(shellBoundary),
+ vertexNormals_(vertexNormals)
+ {
+ // 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");
+
+ // Length scale parameter
+ L_c_ = parameters.template get<double>("L_c");
+
+ // Curvature exponent
+ q_ = parameters.template get<double>("q");
+
+ // Shear correction factor
+ kappa_ = parameters.template get<double>("kappa");
+
+ b1_ = parameters.template get<double>("b1");
+ b2_ = parameters.template get<double>("b2");
+ b3_ = parameters.template get<double>("b3");
+ }
+
+ /** \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::GFE::CosseratStrain<field_type,3,dim>& U) const
+ {
+ Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
+ for (int i=0; i<dim; i++)
+ UMinus1[i][i] -= 1;
+
+ return mu_ * sym(UMinus1).frobenius_norm2()
+ + mu_c_ * skew(UMinus1).frobenius_norm2()
+ + (mu_*lambda_)/(2*mu_ + lambda_) * traceSquared(sym(UMinus1));
+ }
+
+ /** \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::GFE::CosseratStrain<field_type,3,dim>& U) const
+ {
+ RT result = 0;
+
+ // shear-stretch energy
+ 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.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
+
+ result += mu_ * sym2x2.frobenius_norm2();
+
+ // first order drill energy
+ 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.matrix()[i][j] - U.matrix()[j][i]);
+
+ result += mu_c_ * skew2x2.frobenius_norm2();
+
+
+ // classical transverse shear energy
+ 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);
+
+ return result;
+ }
+
+ /** \brief Energy for large-deformation problems (private communication by Patrizio Neff)
+ */
+ RT nonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,dim>& U) const
+ {
+ Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
+ for (int i=0; i<dim; i++)
+ UMinus1[i][i] -= 1;
+
+ RT detU = U.determinant();
+
+ return mu_ * sym(UMinus1).frobenius_norm2() + mu_c_ * skew(UMinus1).frobenius_norm2()
+ + (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
+ }
+
+ /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
+ * the second equation of (4.4) in Neff's paper
+ */
+ RT longNonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,dim,dim-1>& U) const
+ {
+ RT result = 0;
+
+ // shear-stretch energy
+ 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.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
+
+ result += mu_ * sym2x2.frobenius_norm2();
+
+ // first order drill energy
+ 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.matrix()[i][j] - U.matrix()[j][i]);
+
+ result += mu_c_ * skew2x2.frobenius_norm2();
+
+
+ // classical transverse shear energy
+ 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
+ RT detU = U.determinant();
+ result += (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
+
+ return result;
+ }
+
+ RT curvatureEnergy(const Tensor3<field_type,3,3,dim-1>& DR) const
+ {
+#ifdef DONT_USE_CURL
+ return mu_ * std::pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0);
+#else
+ return mu_ * std::pow(L_c_ * L_c_ * curl(DR).frobenius_norm2(),q_/2.0);
+#endif
+ }
+
+ RT bendingEnergy(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,dim-1>& DR) const
+ {
+ // left-multiply the derivative of the third director (in DR[][2][]) with R^T
+ Dune::FieldMatrix<field_type,3,3> RT_DR3(0);
+ for (int i=0; i<3; i++)
+ for (int j=0; j<dim; j++)
+ for (int k=0; k<3; k++)
+ RT_DR3[i][j] += R[k][i] * DR[k][2][j];
+
+ return mu_ * sym(RT_DR3).frobenius_norm2()
+ + mu_c_ * skew(RT_DR3).frobenius_norm2()
+ + mu_*lambda_/(2*mu_+lambda_) * traceSquared(RT_DR3);
+ }
+
+RT energy(const typename Basis::LocalView& localView,
+ const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
+{
+ // The element geometry
+ auto element = localView.element();
+ auto geometry = element.geometry();
+
+ // The set of shape functions on this element
+ const auto& localFiniteElement = localView.tree().finiteElement();
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Construct a linear (i.e., non-constant!) normal field on each element
+ ////////////////////////////////////////////////////////////////////////////////////
+ auto gridView = localView.globalBasis().gridView();
+
+ assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
+ std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
+ for (size_t i=0; i<cornerNormals.size(); i++)
+ cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,gridDim)];
+
+ typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
+ typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;
+ //it.type()?
+ P1FiniteElementCache p1FiniteElementCache;
+ const auto& p1LocalFiniteElement = p1FiniteElementCache.get(element.type());
+ Dune::GFE::LocalProjectedFEFunction<gridDim, DT, P1LocalFiniteElement, UnitVector<double,3> > unitNormals(p1LocalFiniteElement, cornerNormals);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Set up the local nonlinear finite element function
+ ////////////////////////////////////////////////////////////////////////////////////
+ typedef LocalGeodesicFEFunction<gridDim, DT, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
+ LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+
+ RT energy = 0;
+
+ for (auto&& it : intersections(shellBoundary_->gridView(), element)) {
+ if (not shellBoundary_->contains(it))
+ continue;
+
+ auto quadOrder = (it.type().isSimplex()) ? localFiniteElement.localBasis().order()
+ : localFiniteElement.localBasis().order() * gridDim;
+
+ const auto& quad = Dune::QuadratureRules<DT, boundaryDim>::rule(it.type(), quadOrder);
+ for (size_t pt=0; pt<quad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());;
+
+ const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
+
+ // The value of the local function
+ RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
+
+ // The derivative of the local function
+ auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+
+ Dune::FieldMatrix<RT,7,2> derivative2D;
+ for (int i = 0; i < 7; i++) {
+ for (int j = 0; j < 2; j++) {
+ derivative2D[i][j] = derivative[i][j];
+ }
+ }
+ //////////////////////////////////////////////////////////
+ // The rotation and its derivative
+ // Note: we need it in matrix coordinates
+ //////////////////////////////////////////////////////////
+
+ Dune::FieldMatrix<field_type,dim,dim> R;
+ value.q.matrix(R);
+ auto RT = transpose(R);
+
+ Tensor3<field_type,3,3,boundaryDim> DR;
+ computeDR(value, derivative2D, DR);
+
+ //////////////////////////////////////////////////////////
+ // Fundamental forms and curvature
+ //////////////////////////////////////////////////////////
+
+ // First fundamental form
+ Dune::FieldMatrix<double,3,3> aCovariant;
+
+ // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
+ // of the element. If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
+ const auto jacobianTransposed = it.geometry().jacobianTransposed(quad[pt].position());
+ // auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
+
+ for (int i=0; i<2; i++)
+ {
+ for (int j=0; j<dimworld; j++)
+ aCovariant[i][j] = jacobianTransposed[i][j];
+ for (int j=dimworld; j<3; i++)
+ aCovariant[i][j] = 0.0;
+ }
+
+ aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
+ aCovariant[2] /= aCovariant[2].two_norm();
+
+ auto aContravariant = aCovariant;
+ aContravariant.invert();
+
+ Dune::FieldMatrix<double,3,3> a(0);
+ for (int alpha=0; alpha<boundaryDim; alpha++)
+ a += dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
+
+ auto a00 = aCovariant[0] * aCovariant[0];
+ auto a01 = aCovariant[0] * aCovariant[1];
+ auto a10 = aCovariant[1] * aCovariant[0];
+ auto a11 = aCovariant[1] * aCovariant[1];
+ auto aScalar = std::sqrt(a00*a11 - a10*a01);
+
+ // Alternator tensor
+ Dune::FieldMatrix<int,2,2> eps = {{0,1},{-1,0}};
+ Dune::FieldMatrix<double,3,3> c(0);
+
+ for (int alpha=0; alpha<2; alpha++)
+ for (int beta=0; beta<2; beta++)
+ c += sqrt(aScalar) * eps[alpha][beta] * dyadicProduct(aContravariant[alpha], aContravariant[beta]);
+
+ // Second fundamental form
+ // The derivative of the normal field
+ auto normalDerivative = unitNormals.evaluateDerivative(quadPos);
+
+ Dune::FieldMatrix<double,3,3> b(0);
+ for (int alpha=0; alpha<boundaryDim; alpha++)
+ {
+ Dune::FieldVector<double,3> vec;
+ for (int i=0; i<3; i++)
+ vec[i] = normalDerivative[i][alpha];
+ b -= dyadicProduct(vec, aContravariant[alpha]);
+ }
+
+ // Gauss curvature
+ auto K = b.determinant();
+
+ // Mean curvatue
+ auto H = 0.5 * trace(b);
+
+ //////////////////////////////////////////////////////////
+ // Strain tensors
+ //////////////////////////////////////////////////////////
+
+ // Elastic shell strain
+ Dune::FieldMatrix<field_type,3,3> Ee(0);
+ Dune::FieldMatrix<field_type,3,3> grad_s_m(0);
+ for (int alpha=0; alpha<boundaryDim; alpha++)
+ {
+ Dune::FieldVector<field_type,3> vec;
+ for (int i=0; i<3; i++)
+ vec[i] = derivative[i][alpha];
+ grad_s_m += dyadicProduct(vec, aContravariant[alpha]);
+ }
+
+ Ee = RT * grad_s_m - a;
+
+ // Elastic shell bending-curvature strain
+ Dune::FieldMatrix<field_type,3,3> Ke(0);
+ for (int alpha=0; alpha<boundaryDim; alpha++)
+ {
+ Dune::FieldMatrix<field_type,3,3> tmp;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ tmp[i][j] = DR[i][j][alpha];
+ auto tmp2 = RT * tmp;
+ Ke += dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
+ }
+
+ //////////////////////////////////////////////////////////
+ // Add the local energy density
+ //////////////////////////////////////////////////////////
+
+ // 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);
+
+ // 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);
+
+ // Add energy density
+ energy += quad[pt].weight() * integrationElement * energyDensity;
+ }
+ }
+
+ return energy;
+}
+
+ RT W_m(const Dune::FieldMatrix<field_type,3,3>& S) const
+ {
+ return W_mixt(S,S);
+ }
+
+ RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T) const
+ {
+ return mu_ * frobeniusProduct(sym(S), sym(T))
+ + mu_c_ * frobeniusProduct(skew(S), skew(T))
+ + lambda_ * mu_ / (lambda_ + 2*mu_) * trace(S) * trace(T);
+ }
+
+ RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S) const
+ {
+ 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
+ {
+ return mu_ * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
+ }
+
+private:
+ /** \brief The Neumann boundary */
+ const BoundaryPatch<GridView>* shellBoundary_;
+
+ /** \brief The normal vectors at the grid vertices. This are used to compute the reference surface curvature. */
+ std::vector<UnitVector<double,3> > vertexNormals_;
+
+ /** \brief The shell thickness */
+ double thickness_;
+
+ /** \brief Lame constants */
+ double mu_, lambda_;
+
+ /** \brief Cosserat couple modulus, preferably 0 */
+ double mu_c_;
+
+ /** \brief Length scale parameter */
+ double L_c_;
+
+ /** \brief Curvature exponent */
+ double q_;
+
+ /** \brief Shear correction factor */
+ double kappa_;
+
+ /** \brief Curvature parameters */
+ double b1_, b2_, b3_;
+
+};
+
+#endif //#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
diff --git a/src/cantilever-dirichlet-values.py b/src/cantilever-dirichlet-values.py
new file mode 100644
index 00000000..4de4795d
--- /dev/null
+++ b/src/cantilever-dirichlet-values.py
@@ -0,0 +1,17 @@
+import math
+
+class DirichletValues:
+ def __init__(self, homotopyParameter):
+ self.homotopyParameter = homotopyParameter
+
+ def deformation(self, x):
+ # Dirichlet b.c. simply clamp the shell in the reference configuration
+ out = x
+
+ return out
+
+
+ def orientation(self, x):
+ rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]
+ return rotation
+
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
new file mode 100644
index 00000000..16dbc219
--- /dev/null
+++ b/src/film-on-substrate.cc
@@ -0,0 +1,482 @@
+#include <iostream>
+#include <fstream>
+
+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
+#include <adolc/taping.h>
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dune/elasticity/materials/exphenckyenergy.hh>
+#include <dune/elasticity/materials/henckyenergy.hh>
+#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
+#include <dune/elasticity/materials/neohookeenergy.hh>
+#include <dune/elasticity/materials/neumannenergy.hh>
+#include <dune/elasticity/materials/sumenergy.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
+#include <dune/fufem/dunepython.hh>
+
+#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/cosseratvtkwriter.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/vertexnormals.hh>
+#include <dune/gfe/surfacecosseratenergy.hh>
+#include <dune/gfe/sumcosseratenergy.hh>
+
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+// grid dimension
+#ifndef WORLD_DIM
+# define WORLD_DIM 3
+#endif
+const int dim = WORLD_DIM;
+const int order = 1;
+
+//differentiation method
+typedef adouble ValueType;
+typedef adouble GradientValueType;
+
+using namespace Dune;
+
+/** \brief A constant vector-valued function, for simple Neumann boundary values */
+struct NeumannFunction
+ : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,dim> >
+{
+ NeumannFunction(const FieldVector<double,dim> values,
+ double homotopyParameter)
+ : values_(values),
+ homotopyParameter_(homotopyParameter)
+ {}
+
+ void evaluate(const FieldVector<double, dim>& x, FieldVector<double,dim>& out) const
+ {
+ out = 0;
+ out.axpy(-homotopyParameter_, values_);
+ }
+
+ FieldVector<double,dim> values_;
+ double homotopyParameter_;
+};
+
+struct VolumeLoad
+ : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
+{
+ VolumeLoad(const FieldVector<double,3> values,
+ double homotopyParameter)
+ : values_(values),
+ homotopyParameter_(homotopyParameter)
+ {}
+
+ void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const {
+ out = 0;
+ out.axpy(homotopyParameter_, values_);
+ }
+
+ FieldVector<double,3> values_;
+ double homotopyParameter_;
+};
+
+
+int main (int argc, char *argv[]) try
+{
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+ // Start Python interpreter
+ Python::start();
+ Python::Reference main = Python::import("__main__");
+ Python::run("import math");
+
+ //feenableexcept(FE_INVALID);
+ Python::runStream()
+ << std::endl << "import sys"
+ << std::endl << "import os"
+ << std::endl << "sys.path.append(os.getcwd() + '/../../src/')"
+ << std::endl;
+
+ typedef BlockVector<FieldVector<double,dim> > SolutionType;
+ typedef RigidBodyMotion<double,dim> TargetSpace;
+ typedef std::vector<TargetSpace> SolutionTypeCosserat;
+
+ const int blocksize = TargetSpace::TangentVector::dimension;
+
+ // parse data file
+ ParameterTree parameterSet;
+
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+ // read solver settings
+ const int numLevels = parameterSet.get<int>("numLevels");
+ int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
+ const double tolerance = parameterSet.get<double>("tolerance");
+ const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
+ const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+ const int multigridIterations = parameterSet.get<int>("numIt");
+ const int nu1 = parameterSet.get<int>("nu1");
+ const int nu2 = parameterSet.get<int>("nu2");
+ const int mu = parameterSet.get<int>("mu");
+ const int baseIterations = parameterSet.get<int>("baseIt");
+ const double mgTolerance = parameterSet.get<double>("mgTolerance");
+ const double baseTolerance = parameterSet.get<double>("baseTolerance");
+ const bool instrumented = parameterSet.get<bool>("instrumented");
+ std::string resultPath = parameterSet.get("resultPath", "");
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ typedef UGGrid<dim> GridType;
+
+ std::shared_ptr<GridType> grid;
+
+ FieldVector<double,dim> lower(0), upper(1);
+
+ if (parameterSet.get<bool>("structuredGrid")) {
+
+ lower = parameterSet.get<FieldVector<double,dim> >("lower");
+ upper = parameterSet.get<FieldVector<double,dim> >("upper");
+
+ std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
+ grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+
+ } else {
+ std::string path = parameterSet.get<std::string>("path");
+ std::string gridFile = parameterSet.get<std::string>("gridFile");
+ grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+ }
+
+ grid->globalRefine(numLevels-1);
+
+ grid->loadBalance();
+
+ if (mpiHelper.rank()==0)
+ std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
+
+ typedef GridType::LeafGridView GridView;
+ GridView gridView = grid->leafGridView();
+
+ // FE basis spanning the FE space that we are working in
+ typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
+ FEBasis feBasis(gridView);
+
+ // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
+ typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
+ FufemFEBasis fufemFEBasis(feBasis);
+
+ // /////////////////////////////////////////
+ // Read Dirichlet values
+ // /////////////////////////////////////////
+
+ BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+ BitSetVector<1> neumannVertices(gridView.size(dim), false);
+ BitSetVector<1> surfaceShellVertices(gridView.size(dim), false);
+
+ const GridView::IndexSet& indexSet = gridView.indexSet();
+
+ // Make Python function that computes which vertices are on the Dirichlet boundary,
+ // based on the vertex positions.
+ std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+
+ // Same for the Neumann boundary
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+
+ // Same for the boundary that will get wrinkled
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("surfaceShellVerticesPredicate", "0") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, bool> pythonSurfaceShellVertices(Python::evaluate(lambda));
+ for (auto&& v : vertices(gridView))
+ {
+ bool isDirichlet;
+ pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
+ dirichletVertices[indexSet.index(v)] = isDirichlet;
+
+ bool isNeumann;
+ pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
+ neumannVertices[indexSet.index(v)] = isNeumann;
+
+ bool isSurfaceShell;
+ pythonSurfaceShellVertices.evaluate(v.geometry().corner(0), isSurfaceShell);
+ surfaceShellVertices[indexSet.index(v)] = isSurfaceShell;
+ }
+
+ BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+ BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
+
+ if (mpiHelper.rank()==0)
+ std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+
+
+ BitSetVector<1> dirichletNodes(feBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+
+ BitSetVector<1> neumannNodes(feBasis.size(), false);
+ constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
+
+ BitSetVector<1> surfaceShellNodes(feBasis.size(), false);
+ constructBoundaryDofs(surfaceShellBoundary,feBasis,surfaceShellNodes);
+
+ BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
+ for (size_t i=0; i<feBasis.size(); i++)
+ for (int j=0; j<blocksize; j++) {
+ if (dirichletNodes[i][0])
+ dirichletDofs[i][j] = true;
+ }
+ for (size_t i=0; i<feBasis.size(); i++)
+ for (int j=3; j<blocksize; j++) {
+ if (not surfaceShellNodes[i][0])
+ dirichletDofs[i][j] = true;
+ }
+
+ // //////////////////////////
+ // Initial iterate
+ // //////////////////////////
+
+ SolutionTypeCosserat x(feBasis.size());
+
+ //Solution in 3D, without the Cosserat directors on the boundary
+ BlockVector<FieldVector<double,3> > v;
+
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
+ ::Functions::interpolate(fufemFEBasis, v, pythonInitialDeformation);
+
+ //Copy over the parts of the boundary
+ for (size_t i=0; i<x.size(); i++) {
+ x[i].r = v[i];
+ }
+
+ ////////////////////////////////////////////////////////
+ // Main homotopy loop
+ ////////////////////////////////////////////////////////
+
+ using namespace Dune::Functions::BasisFactory;
+ auto powerBasis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
+ BlockVector<FieldVector<double,dim>> identity;
+ Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+
+ BlockVector<FieldVector<double,dim> > displacement(v.size());
+ for (int i = 0; i < v.size(); i++) {
+ displacement[i] = v[i];
+ }
+ displacement -= identity;
+
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
+ auto localDisplacementFunction = localFunction(displacementFunction);
+
+ auto neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
+ std::cout << "Neumann values: " << neumannValues << std::endl;
+
+ // We need to subsample, because VTK cannot natively display real second-order functions
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order));
+ vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_0");
+
+ for (int i=0; i<numHomotopySteps; i++)
+ {
+ double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
+
+ // ////////////////////////////////////////////////////////////
+ // Create an assembler for the energy functional
+ // ////////////////////////////////////////////////////////////
+
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+
+ std::shared_ptr<VolumeLoad> volumeLoad;
+
+ std::shared_ptr<NeumannFunction> neumannFunction;
+ if (parameterSet.hasKey("neumannValues"))
+ {
+ neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,dim> >("neumannValues"),
+ homotopyParameter);
+
+ std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,dim> >("neumannValues") << std::endl;
+ }
+
+ if (mpiHelper.rank() == 0)
+ {
+ std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
+ std::cout << "Material parameters:" << std::endl;
+ materialParameters.report();
+ }
+
+ // Assembler using ADOL-C
+ std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
+ std::shared_ptr<LocalFEStiffness<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ std::vector<Dune::FieldVector<ValueType, dim>>,
+ std::vector<GradientValueType> > > elasticEnergy;
+
+ if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
+ elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType> >(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
+ elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType> >(materialParameters);
+
+ if (parameterSet.get<std::string>("energy") == "neohooke")
+ elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType> >(materialParameters);
+
+ if (parameterSet.get<std::string>("energy") == "hencky")
+ elasticEnergy = std::make_shared<HenckyEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType> >(materialParameters);
+
+ if (parameterSet.get<std::string>("energy") == "exphencky")
+ elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType> >(materialParameters);
+
+ if(!elasticEnergy)
+ DUNE_THROW(Exception, "Error: Selected energy not available!");
+
+ auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType> >(&neumannBoundary,neumannFunction.get());
+
+ auto elasticAndNeumann = std::make_shared<SumEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType, GradientValueType>>(elasticEnergy, neumannEnergy);
+
+
+ using LocalEnergyBase = LocalGeodesicFEStiffness<FEBasis,RigidBodyMotion<adouble, dim> >;
+
+ std::shared_ptr<LocalEnergyBase> surfaceCosseratEnergy;
+ std::vector<UnitVector<double,3> > vertexNormals(gridView.size(3));
+ Dune::FieldVector<double,3> vertexNormalRaw = {0,0,1};
+ for (int i = 0; i< vertexNormals.size(); i++) {
+ UnitVector vertexNormal(vertexNormalRaw);
+ vertexNormals[i] = vertexNormal;
+ }
+ surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary);
+
+ std::shared_ptr<LocalEnergyBase> totalEnergy;
+ totalEnergy = std::make_shared<GFE::SumCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>> (elasticAndNeumann, surfaceCosseratEnergy);
+
+
+ LocalGeodesicFEADOLCStiffness<FEBasis,
+ TargetSpace> localGFEADOLCStiffness(totalEnergy.get());
+
+ GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
+
+ // /////////////////////////////////////////////////
+ // Create a Riemannian trust-region solver
+ // /////////////////////////////////////////////////
+
+ RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+ solver.setup(*grid,
+ &assembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+
+ ////////////////////////////////////////////////////////
+ // Set Dirichlet values
+ ////////////////////////////////////////////////////////
+
+ Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
+
+ Python::Callable C = dirichletValuesClass.get("DirichletValues");
+
+ // Call a constructor.
+ Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
+
+ // Extract object member functions as Dune functions
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
+ PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
+
+ std::vector<FieldVector<double,3> > ddV;
+ ::Functions::interpolate(fufemFEBasis, ddV, deformationDirichletValues, dirichletDofs);
+
+ std::vector<FieldMatrix<double,3,3> > dOV;
+ ::Functions::interpolate(fufemFEBasis, dOV, orientationDirichletValues, dirichletDofs);
+
+ for (size_t j=0; j<x.size(); j++)
+ if (dirichletNodes[j][0])
+ {
+ x[j].r = ddV[j];
+ x[j].q.set(dOV[j]);
+ }
+
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ solver.setInitialIterate(x);
+ solver.solve();
+
+ x = solver.getSol();
+
+ /////////////////////////////////
+ // Output result
+ /////////////////////////////////
+ for (size_t i=0; i<x.size(); i++)
+ v[i] = x[i].r;
+
+ // Compute the displacement
+ BlockVector<FieldVector<double,dim> > displacement(v.size());
+ for (int i = 0; i < v.size(); i++) {
+ displacement[i] = v[i];
+ }
+ displacement -= identity;
+
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
+ auto localDisplacementFunction = localFunction(displacementFunction);
+
+ // We need to subsample, because VTK cannot natively display real second-order functions
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order));
+ vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_" + std::to_string(i+1));
+ }
+} catch (Exception& e) {
+ std::cout << e.what() << std::endl;
+}
\ No newline at end of file
diff --git a/src/film-on-substrate.parset b/src/film-on-substrate.parset
new file mode 100644
index 00000000..63b64fc6
--- /dev/null
+++ b/src/film-on-substrate.parset
@@ -0,0 +1,123 @@
+#############################################
+# Grid parameters
+#############################################
+
+structuredGrid = true
+
+# bounding box
+lower = 0 0 0
+upper = 10 10 0.5
+
+elements = 10 10 5
+
+# Number of grid levels
+numLevels = 2
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-3
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 20
+
+trustRegionScaling = 1 1 1 0.01 0.01 0.01
+
+# Initial trust-region radius
+initialTrustRegionRadius = 3.125
+
+# Number of multigrid iterations per trust-region step
+numIt = 400
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 1
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-5
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+# Measure convergence
+instrumented = 0
+
+############################
+# Material parameters
+############################
+
+energy = stvenantkirchhoff
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+# shell thickness
+thickness = 0.6
+
+# Lame parameters
+# corresponds to E = 71240 N/mm^2, nu=0.31
+# However, we use units N/m^2
+mu = 2.7191e+4
+lambda = 4.4364e+4
+
+#Lame parameters for the cosserat material
+mu_cosserat = 2.7191e+4
+lambda_cosserat = 4.4364e+4
+
+# Cosserat couple modulus
+mu_c = 0
+
+# Length scale parameter
+L_c = 0.6
+
+# Curvature exponent
+q = 2
+
+# Shear correction factor
+kappa = 1
+
+b1 = 1
+b2 = 1
+b3 = 1
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = cantilever
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+### Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+surfaceShellVerticesPredicate = "x[2] > 0.49"
+
+### Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 9.99"
+
+## Neumann values
+neumannValues = 1e5 0 0
+
+# Initial deformation
+initialDeformation = "x"
+
+#startFromFile = yes
+#initialIterateFilename = initial_iterate.vtu
--
GitLab