A shell energy for a geometrically exact Cosserat shell with a nonplanar reference configuration

This model was given to me by Patrizio Neff in private communication.

The code is a prototype.  It appears to work, but it has not been tested very
thoroughly, and quite likely there are some bugs left.  Also, the code contains
the shape of the reference configuration hard-wired (a sphere).

dune/gfe/nonplanarcosseratshellenergy.hh

+#ifndef DUNE_GFE_NONPLANARCOSSERATSHELLENERGY_HH
+#define DUNE_GFE_NONPLANARCOSSERATSHELLENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/tensor3.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+
+
+template<class Basis, int dim, class field_type=double>
+class NonplanarCosseratShellEnergy
+  : public LocalGeodesicFEStiffness<Basis,RigidBodyMotion<field_type,dim> >
+{
+  // grid types
+  typedef typename Basis::GridView GridView;
+  typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
+  typedef typename GridView::ctype DT;
+  typedef RigidBodyMotion<field_type,dim> TargetSpace;
+  typedef typename TargetSpace::ctype RT;
+  typedef typename GridView::template Codim<0>::Entity Entity;
+
+  // some other sizes
+  enum {gridDim=GridView::dimension};
+  enum {dimworld=GridView::dimensionworld};
+
+  /** \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,3>& value,
+                        const Dune::FieldMatrix<field_type,7,gridDim>& derivative,
+                        Tensor3<field_type,3,3,gridDim>& 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<gridDim; 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
+   */
+  NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
+                               const BoundaryPatch<GridView>* neumannBoundary,
+                               const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
+                               const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
+  : neumannBoundary_(neumannBoundary),
+    neumannFunction_(neumannFunction),
+    volumeLoad_(volumeLoad)
+  {
+    // The shell thickness
+    thickness_ = parameters.template get<double>("thickness");
+
+    // Lame constants
+    mu_ = parameters.template get<double>("mu");
+    lambda_ = parameters.template get<double>("lambda");
+
+    // Cosserat couple modulus
+    mu_c_ = parameters.template get<double>("mu_c");
+
+    // Length scale parameter
+    L_c_ = parameters.template get<double>("L_c");
+
+    // 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 */
+  RT energy (const Entity& e,
+             const LocalFiniteElement& localFiniteElement,
+             const std::vector<TargetSpace>& localSolution) const;
+
+  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));
+  }
+
+  /** \brief The shell thickness */
+  double thickness_;
+
+  /** \brief Lame constants */
+  double mu_, lambda_;
+
+  /** \brief Cosserat couple modulus */
+  double mu_c_;
+
+  /** \brief Length scale parameter */
+  double L_c_;
+
+  /** \brief Curvature parameters */
+  double b1_, b2_, b3_;
+
+  /** \brief The Neumann boundary */
+  const BoundaryPatch<GridView>* neumannBoundary_;
+
+  /** \brief The function implementing the Neumann data */
+  const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction_;
+
+  /** \brief The function implementing a volume load */
+  const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
+};
+
+template <class Basis, int dim, class field_type>
+typename NonplanarCosseratShellEnergy<Basis,dim,field_type>::RT
+NonplanarCosseratShellEnergy<Basis,dim,field_type>::
+energy(const Entity& element,
+       const typename Basis::LocalView::Tree::FiniteElement& localFiniteElement,
+       const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
+{
+  assert(element.type() == localFiniteElement.type());
+
+  // The element geometry
+  auto geometry = element.geometry();
+
+  ////////////////////////////////////////////////////////////////////////////////////
+  //  Construct a linear (i.e., non-constant!) normal field on the surface
+  ////////////////////////////////////////////////////////////////////////////////////
+
+  std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
+  for (size_t i=0; i<cornerNormals.size(); i++)
+  {
+    // TODO: WARNING WARNING WARNING:  Shape is hard-coded here!
+    // The sphere:
+    cornerNormals[i] = geometry.corner(i) / geometry.corner(i).two_norm();
+  }
+
+  typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
+  typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;
+  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;
+
+  auto quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                                : localFiniteElement.localBasis().order() * gridDim;
+
+  const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+  for (size_t pt=0; pt<quad.size(); pt++)
+  {
+    // Local position of the quadrature point
+    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+    const DT integrationElement = geometry.integrationElement(quadPos);
+
+    DT weight = quad[pt].weight() * integrationElement;
+
+    // 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);
+
+    //////////////////////////////////////////////////////////
+    //  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,gridDim> DR;
+    computeDR(value, derivative, DR);
+
+    //////////////////////////////////////////////////////////
+    //  Fundamental forms and curvature
+    //////////////////////////////////////////////////////////
+
+    // First fundamental form
+    Dune::FieldMatrix<double,3,3> aCovariant;
+
+    aCovariant[0] = geometry.jacobianTransposed(quadPos)[0];
+    aCovariant[1] = geometry.jacobianTransposed(quadPos)[1];
+    aCovariant[2] = Arithmetic::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<gridDim; 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<gridDim; 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<gridDim; 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<gridDim; 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
+    energy += weight * (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_m(Ee)
+            + weight * (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_m(Ee*b + c*Ke)
+            + weight * Dune::Power<3>::eval(thickness_) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
+            + weight * Dune::Power<5>::eval(thickness_) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
+
+    // Add the bending energy density
+    energy += weight * (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_curv(Ke)
+            + weight * (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_curv(Ke*b)
+            + weight * Dune::Power<5>::eval(thickness_) / 80.0 * W_curv(Ke*b*b);
+
+    ///////////////////////////////////////////////////////////
+    // Volume load contribution
+    ///////////////////////////////////////////////////////////
+
+    if (not volumeLoad_)
+      continue;
+
+    // Value of the volume load density at the current position
+    Dune::FieldVector<double,3> volumeLoadDensity;
+
+    if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
+      std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
+    else
+      volumeLoad_->evaluate(geometry.global(quad[pt].position()), volumeLoadDensity);
+
+    // 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;
+  }
+
+
+  //////////////////////////////////////////////////////////////////////////////
+  //   Assemble boundary contributions
+  //////////////////////////////////////////////////////////////////////////////
+
+  if (not neumannFunction_)
+    return energy;
+
+  for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
+  {
+    if (not neumannBoundary_ or not neumannBoundary_->contains(it))
+      continue;
+
+    const auto& quad = Dune::QuadratureRules<DT, gridDim-1>::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);
+
+      // Value of the Neumann data at the current position
+      Dune::FieldVector<double,3> neumannValue;
+
+      if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
+        std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
+      else
+        neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+
+      // 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;
+    }
+
+  }
+
+  return energy;
+}
+
+#endif   //#ifndef DUNE_GFE_NONPLANARCOSSERATSHELLENERGY_HH
+