diff --git a/dune/gfe/assemblers/discretekirchhoffbendingenergy.hh b/dune/gfe/assemblers/discretekirchhoffbendingenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..542cca43532a1b4c07b64d92e97ea0a9013dc145
--- /dev/null
+++ b/dune/gfe/assemblers/discretekirchhoffbendingenergy.hh
@@ -0,0 +1,94 @@
+#ifndef DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
+#define DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/assemblers/localenergy.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+#include <dune/gfe/tensor3.hh>
+
+#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
+
+/**
+ * \brief Assemble the discrete Kirchhoff bending energy for a single element.
+ *
+ * The energy is essentially the harmonic energy of a
+ * discrete Jacobian operator which is represented as a
+ * linear combination in a P2 finite element space.
+ * The gradient of that linear combination serves as an
+ * discrete Hessian in the energy.
+ */
+namespace Dune::GFE
+{
+ template<class Basis, class DiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
+ class DiscreteKirchhoffBendingEnergy
+ : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+ {
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+ typedef typename DiscreteKirchhoffFunction::DerivativeType DerivativeType;
+
+ // Grid dimension
+ constexpr static int gridDim = GridView::dimension;
+
+ // P2-LocalFiniteElement used to represent the discrete Jacobian on the current element.
+ typedef typename Dune::LagrangeSimplexLocalFiniteElement<DT, double, gridDim, 2> P2LocalFiniteElement;
+
+ public:
+ DiscreteKirchhoffBendingEnergy(DiscreteKirchhoffFunction &localFunction)
+ : localFunction_(localFunction)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+ {
+ DUNE_THROW(NotImplemented, "!");
+ }
+
+
+ /** \brief Assemble the energy for a single element */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localConfiguration) const override
+ {
+
+ const auto element = localView.element();
+
+ int quadOrder = 3;
+ const auto &quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
+
+ /**
+ bind the underlying Hermite basis to the current element.
+ Note: The DiscreteKirchhoffFunction object stores a vector of global
+ and local basis coefficients that are accessed for evaluation
+ methods (linear combinations of hermite basis function and local coefficients).
+ Therefore the local configuration is passed to the bind method as well since
+ these coefficient vectors need to be updated with the new local configuration
+ previous to evaluation.
+ */
+ localFunction_.bind(element,localConfiguration);
+
+ RT bendingEnergy = 0;
+
+ for (auto&& quadPoint : quad)
+ {
+ auto discreteHessian= localFunction_.evaluateDiscreteHessian(quadPoint.position());
+
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ bendingEnergy += weight * discreteHessian.frobenius_norm2();
+ }
+
+ return 0.5 * bendingEnergy;
+ }
+
+ private:
+ DiscreteKirchhoffFunction& localFunction_;
+
+ P2LocalFiniteElement lagrangeLFE_;
+ };
+
+} // namespace Dune::GFE
+#endif
diff --git a/dune/gfe/assemblers/forceenergy.hh b/dune/gfe/assemblers/forceenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..916068c9b676a6a020fc287e35a0f5bcc4103183
--- /dev/null
+++ b/dune/gfe/assemblers/forceenergy.hh
@@ -0,0 +1,97 @@
+#ifndef DUNE_GFE_FORCEENERGY_HH
+#define DUNE_GFE_FORCEENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/assemblers/localenergy.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+
+
+/**
+ * \brief Assemble the energy of a force term (force times deformation) for a single element.
+ */
+namespace Dune::GFE
+{
+ template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
+ class ForceEnergy
+ : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+ {
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+
+ // Grid dimension
+ constexpr static int gridDim = GridView::dimension;
+
+ public:
+ ForceEnergy(LocalDiscreteKirchhoffFunction &localFunction,
+ LocalForce &localForce)
+ : localFunction_(localFunction),
+ localForce_(localForce)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+ {
+ DUNE_THROW(NotImplemented, "!");
+ }
+
+
+ /** \brief Assemble the energy for a single element */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localConfiguration) const override
+ {
+ RT energy = 0;
+
+ int quadOrder = 3;
+
+ const auto &quad = QuadratureRules<double, gridDim>::rule(GeometryTypes::simplex(gridDim), quadOrder);
+
+ const auto element = localView.element();
+
+ /** bind the local force to the current element */
+ localForce_.bind(element);
+
+ /**
+ bind the underlying Hermite basis to the current element.
+ Note: The DiscreteKirchhoffFunction object stores a vector of global
+ and local basis coefficients that are accessed for evaluation
+ methods (linear comibnations of hermite basis function and local coefficients).
+ Therefore the local configuration is passed to the bind method as well since
+ these coefficient vectors need to be updated with the new local configuration
+ previous to evaluation.
+ */
+ localFunction_.bind(element,localConfiguration);
+
+ /**
+ * @brief Compute contribution of the force Term
+ *
+ * Integrate the scalar product of the force
+ * with the local deformation function 'localFunction_'
+ * over the current element.
+ */
+ RT forceTermEnergy = 0;
+
+ for (auto&& quadPoint : quad)
+ {
+ auto deformationValue = localFunction_.evaluate(quadPoint.position());
+ auto forceValue = localForce_(quadPoint.position());
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
+ }
+
+ return forceTermEnergy;
+ }
+
+ private:
+ LocalDiscreteKirchhoffFunction& localFunction_;
+
+ LocalForce& localForce_;
+
+ };
+
+} // namespace Dune::GFE
+#endif
diff --git a/dune/gfe/bendingisometryhelper.hh b/dune/gfe/bendingisometryhelper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c3075959de53fb93281317b0e5702c6a055506c7
--- /dev/null
+++ b/dune/gfe/bendingisometryhelper.hh
@@ -0,0 +1,399 @@
+#ifndef DUNE_GFE_BENDINGISOMETRYHELPER_HH
+#define DUNE_GFE_BENDINGISOMETRYHELPER_HH
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/functions/functionspacebases/cubichermitebasis.hh>
+#include <dune/functions/common/differentiablefunction.hh>
+#include <dune/functions/common/differentiablefunctionfromcallables.hh>
+#include <dune/functions/common/functionconcepts.hh>
+
+
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+#include <dune/gfe/functions/embeddedglobalgfefunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
+
+#include <dune/matrix-vector/crossproduct.hh>
+
+/** \file
+ * \brief Various helper functions related to bending isometries.
+ */
+
+
+namespace Dune::GFE::Impl
+{
+
+
+ /** \brief The identity function in R^d, and its derivative
+ *
+ * This is needed to compute the displacement from the deformation (for visualization).
+ * Since we want to represent it in a Hermite finite element space, we need
+ * the derivative as well.
+ *
+ * Is there no shorter way to implement this?
+ */
+ template<class K>
+ class IdentityGridEmbedding
+ {
+ public:
+ //! Evaluate identity
+ Dune::FieldVector<K,3> operator() (const Dune::FieldVector<K,2>& x) const
+ {
+ return {x[0], x[1], 0.0};
+ }
+
+ /**
+ * \brief Obtain derivative of identity function
+ */
+ friend auto derivative(const IdentityGridEmbedding& p)
+ {
+ return [](const Dune::FieldVector<K,2>& x) { return Dune::FieldMatrix<K,3,2>({{1,0}, {0,1}, {0,0}}); };
+ }
+ };
+
+
+ /**
+ Define some Generic lambda expressions (C++ 20)
+ */
+ auto cancelLastMatrixColumn = []<typename RT>(Dune::FieldMatrix<RT,3,3> matrix) -> Dune::FieldMatrix<RT,3,2>
+ {
+ return Dune::FieldMatrix<RT,3,2>{{matrix[0][0],matrix[0][1]},
+ {matrix[1][0],matrix[1][1]},
+ {matrix[2][0],matrix[2][1]}
+ };
+ };
+
+ auto rotationExtraction = []<typename RT>(Dune::FieldVector<RT,7> vec) -> Rotation<RT,3>
+ {
+ std::array<RT,4> quaternion {vec[3],vec[4],vec[5],vec[6]};
+ return Rotation<RT,3>(quaternion);
+ };
+
+ auto deformationExtraction = []<typename RT>(Dune::FieldVector<RT,7> vec) -> Dune::FieldVector<RT,3>
+ {
+ return Dune::FieldVector<RT,3>{vec[0],vec[1],vec[2]};
+ };
+
+
+
+
+
+
+
+ /**
+ * @brief Data structure conversion from 'VectorSpaceCoefficients' to 'IsometryCoefficients'.
+ */
+ template<class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients, class IsometryCoefficients>
+ void vectorToIsometryCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
+ const CoefficientBasis& coefficientBasis,
+ const Coefficients& vectorCoefficients,
+ IsometryCoefficients& isometryCoefficients)
+ {
+ isometryCoefficients.resize(coefficientBasis.size());
+ auto localView = discreteKirchhoffBasis.localView();
+ auto localViewCoefficients = coefficientBasis.localView();
+
+ using RT = typename IsometryCoefficients::value_type::ctype;
+
+ for (const auto &element : elements(coefficientBasis.gridView()))
+ {
+ localView.bind(element);
+ localViewCoefficients.bind(element);
+
+ /**
+ Create data structures to store deformation values and the deformation gradient (Dofs)
+ on the current element. We store the values in an std::vector where each entry
+ contains the deformation vector (resp. x,y derivative vectors) at one node of the element.
+ */
+ std::vector<Dune::FieldVector<RT,3> > deformationDofs(3);
+ std::vector<Dune::FieldVector<RT,3> > xDerivativeDofs(3);
+ std::vector<Dune::FieldVector<RT,3> > yDerivativeDofs(3);
+
+ // Loop over components of power basis
+ for(std::size_t k=0; k<3 ; k++)
+ {
+ // This could be placed out of the loop if we assume a power basis.
+ const auto& hermiteLFE = localView.tree().child(k).finiteElement();
+
+ for(std::size_t i=0; i<hermiteLFE.size(); i++)
+ {
+ auto localIdx = localView.tree().child(k).localIndex(i);
+ auto globalIdx = localView.index(localIdx);
+ // Get current node index (Caution: This assumes that the hermite LocalFiniteElement only contains vertex dofs.)
+ auto nodeIdx = hermiteLFE.localCoefficients().localKey(i).subEntity();
+
+ /**
+ We use the functionalDescriptor of the (discreteKirchhoffBasis) hermite LocalFiniteELement
+ to determine whether a current DOF corresponds to a function evaluation or partial derivative w.r.t
+ the first or second variable.
+ (Alternatively to the - experimental - FunctionalDescriptor: use localKey(i).index(): 0~value, 1~xDerivative, 2~yDerivative.)
+ */
+ if(std::array<unsigned int,2> val {0,0} ; hermiteLFE.localInterpolation().functionalDescriptor(i).partialDerivativeOrder() == val )
+ {
+ deformationDofs[nodeIdx][k] = vectorCoefficients[globalIdx[0]][globalIdx[1]];
+ }
+ if(std::array<unsigned int,2> val {1,0} ; hermiteLFE.localInterpolation().functionalDescriptor(i).partialDerivativeOrder() == val )
+ {
+ xDerivativeDofs[nodeIdx][k] = vectorCoefficients[globalIdx[0]][globalIdx[1]];
+ }
+ if(std::array<unsigned int,2> val {0,1} ; hermiteLFE.localInterpolation().functionalDescriptor(i).partialDerivativeOrder() == val )
+ {
+ yDerivativeDofs[nodeIdx][k] = vectorCoefficients[globalIdx[0]][globalIdx[1]];
+ }
+ }
+ }
+
+ // Loop over vertices of the current element
+ for(std::size_t c=0; c<3 ; c++)
+ {
+ /** Normalize derivative vectors (this is necessary due to the h-dependent scaling used in the hermite-Basis). */
+ xDerivativeDofs[c] /= xDerivativeDofs[c].two_norm();
+ yDerivativeDofs[c] /= yDerivativeDofs[c].two_norm();
+
+ /**
+ cross product to get last column of rotation matrix.
+ */
+ Dune::FieldVector<RT,3> cross = Dune::MatrixVector::crossProduct(xDerivativeDofs[c],yDerivativeDofs[c]);
+
+ Dune::FieldMatrix<RT,3,3> rotMatrix(0);
+
+ for(std::size_t k=0; k<3; k++)
+ {
+ rotMatrix[k][0] = xDerivativeDofs[c][k];
+ rotMatrix[k][1] = yDerivativeDofs[c][k];
+ rotMatrix[k][2] = cross[k];
+ }
+
+ /**
+ Check if derivative vectors are orthonormal.
+ */
+ Dune::FieldMatrix<double,3,3> I = Dune::ScaledIdentityMatrix<double,3>(1);
+
+ if(((rotMatrix.transposed()*rotMatrix) - I).frobenius_norm()>1e-8)
+ DUNE_THROW(Dune::Exception, "Error: input rotation is not orthonormal!"<< ((rotMatrix.transposed()*rotMatrix) - I).frobenius_norm());
+
+ size_t localIdxOut = localViewCoefficients.tree().localIndex(c);
+ size_t globalIdxOut = localViewCoefficients.index(localIdxOut);
+
+ using namespace Dune::Indices;
+ isometryCoefficients[globalIdxOut][_1].set(rotMatrix);
+ isometryCoefficients[globalIdxOut][_0] = (RealTuple<RT, 3>)deformationDofs[c];
+ }
+ }
+ }
+
+
+ /**
+ * @brief Data structure conversion from 'IsometryCoefficients' to 'VectorSpaceCoefficients'.
+ */
+ template<class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients, class IsometryCoefficients>
+ void isometryToVectorCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
+ const CoefficientBasis& coefficientBasis,
+ Coefficients& vectorCoefficients,
+ const IsometryCoefficients& isometryCoefficients)
+ {
+ auto localView = discreteKirchhoffBasis.localView();
+ auto localViewCoefficients = coefficientBasis.localView();
+
+ using GridView = typename DiscreteKirchhoffBasis::GridView;
+ auto gridView = discreteKirchhoffBasis.gridView();
+
+ /** Create an EmbeddedGlobalGFEFunction that serves as a GridViewFunction later on.
+ We need a function to be passed into the interpolation method of the hermite basis
+ (although the values in between grid nodes do not matter).
+
+ The range type of EmbeddedGlobalFEFunction is a FieldVector<double,7>
+ where the first three entries correspond to the deformation and the
+ last four entries correspond to a (quaternion) rotation.
+ */
+ typedef Dune::GFE::LocalProjectedFEFunction<GridView::dimension, double, typename CoefficientBasis::LocalView::Tree::FiniteElement, Dune::GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > LocalInterpolationRule;
+ Dune::GFE::EmbeddedGlobalGFEFunction<CoefficientBasis, LocalInterpolationRule, Dune::GFE::ProductManifold<RealTuple<double,3>, Rotation<double,3> > > embeddedGlobalFunction(coefficientBasis, isometryCoefficients);
+
+ auto productSpaceGridViewFunction = Dune::Functions::makeAnalyticGridViewFunction(embeddedGlobalFunction, gridView);
+
+ /** Extract the deformation and rotation into separate functions. */
+ auto deformationGlobalFunction = Dune::Functions::makeComposedGridFunction(deformationExtraction,productSpaceGridViewFunction);
+
+
+
+ auto rotationGlobalFunction = Dune::Functions::makeComposedGridFunction(cancelLastMatrixColumn,
+ Dune::Functions::makeComposedGridFunction(Rotation<double,3>::quaternionToMatrix,
+ Dune::Functions::makeComposedGridFunction(rotationExtraction,productSpaceGridViewFunction)));
+
+ using Domain = const Dune::FieldVector<double,2>&;
+ using Range = Dune::FieldVector<double,3>;
+ auto globalIsometryFunction = makeDifferentiableFunctionFromCallables(Dune::Functions::SignatureTag<Range(Domain)>(),deformationGlobalFunction,rotationGlobalFunction);
+
+ /** Interpolate into the hermite basis to get the coefficient vector. */
+ interpolate(discreteKirchhoffBasis, vectorCoefficients, globalIsometryFunction);
+ }
+
+
+ template<class LocalDiscreteKirchhoffFunction, class NormalBasis>
+ auto computeDiscreteSurfaceNormal(LocalDiscreteKirchhoffFunction& deformationFunction,
+ const NormalBasis& normalBasis)
+ {
+ std::vector<Dune::FieldVector<double,3> > discreteNormalVectorCoefficients(normalBasis.size());
+
+ auto localView = normalBasis.localView();
+ for (auto&& element : elements(normalBasis.gridView()))
+ {
+ auto geometry = element.geometry();
+
+ localView.bind(element);
+ const auto nSf = localView.tree().child(0).finiteElement().localBasis().size();
+
+ deformationFunction.bind(element);
+
+ for (int i=0; i<geometry.corners(); i++)
+ {
+ auto numDir = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(i)));
+
+ /**
+ Get the normal vector via cross-product.
+ */
+ Dune::FieldVector<double,3> normal = {numDir[1][0]*numDir[2][1] - numDir[1][1]*numDir[2][0],
+ numDir[2][0]*numDir[0][1] - numDir[0][0]*numDir[2][1],
+ numDir[0][0]*numDir[1][1] - numDir[1][0]*numDir[0][1]};
+ // printvector(std::cout, normal, "discrete normal:" , "--");
+ for (size_t k=0; k < 3; k++)
+ for (size_t i=0; i < nSf; i++)
+ {
+ size_t localIdx = localView.tree().child(k).localIndex(i); // hier i:leafIdx
+ size_t globalIdx = localView.index(localIdx);
+ discreteNormalVectorCoefficients[globalIdx] = normal[k];
+ }
+ }
+ }
+
+ return discreteNormalVectorCoefficients;
+ }
+
+
+ /**
+ * @brief Compute Coefficients of the local discrete gradient operator.
+ *
+ * The discrete Gradient is a special linear combination represented in a [P2]^2 space (locally by a representation local finite element)
+ * The coefficients of this linear combination correspond to certain linear combinations of the Gradients of localfunction_ .
+ * The coefficients are stored in the form [Basisfunctions x components x gridDim]
+ * in a BlockVector<FieldMatrix<RT, 3, gridDim>> .
+ */
+ template<class Coefficient,class LocalP2Element, class DeformationFunction >
+ auto discreteGradientCoefficients(Coefficient& discreteGradientCoefficients,
+ LocalP2Element& representationLFE,
+ DeformationFunction& deformationFunction,
+ auto& element)
+ {
+ using RT = typename DeformationFunction::RT;
+
+ auto geometry = element.geometry();
+ constexpr static int gridDim = DeformationFunction::gridDim;
+ auto gridView = deformationFunction.gridView();
+ const auto &indexSet = gridView.indexSet();
+ discreteGradientCoefficients.resize(representationLFE.size());
+ /**
+ * @brief On the current element we need to assign the coefficients of the discrete Gradient to the right P2-Lagrange-Basisfunctions.
+ * this is done by iterating over all P2-Lagrange-Basisfunctions and determine the corresponding coefficients.
+ * We distinguish two situations:
+ * - [1] If the Lagrange-basisfunction is assigned to a node: the corresponding coefficient is just the evaluation of the deformation gradient at that specific node.
+ * - [2] If the Lagrange-basisfunction is assigned to an edge center: the corresponding coefficient is decomposed into a normal-& tangential part,
+ * where the normal-part: is given by the affine combination of the deformation gradient values at the edge-vertices
+ * and the tangential-part: is given by the tangential part of the deformation gradient at the edge-center.
+ *
+ * * We do this for all (k=3) components of the deformation-function simulatneously.
+ */
+ for(std::size_t i=0; i<representationLFE.size(); i++)
+ {
+ // [1] Determine coefficients for basis functions assigned to a vertex
+ if (representationLFE.localCoefficients().localKey(i).codim() == 2)
+ {
+ size_t nodeIndex = representationLFE.localCoefficients().localKey(i).subEntity(); // This gives index on reference Element ?
+ // Alternative: use subEntity to get (local) node position
+
+ typename DeformationFunction::DerivativeType whJacobianValue = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(nodeIndex)));
+ discreteGradientCoefficients[i] = whJacobianValue;
+ }
+
+ // [2] Determine coefficients for basis functions assigned to an edge
+ if (representationLFE.localCoefficients().localKey(i).codim() == 1)
+ {
+ /**
+ * @brief On each edge the (componentwise) value of the discete jacobian on the edge midpoint
+ * is decomposed into a normal and tangent component.
+ * The normal component at the edge center is given as the affine combination of the
+ * normal component of the jacobian of localfunctions_ at the two vertices on each edge.
+ * The tangential component is the tangent component of the jacobian of localfunctions_ on
+ * the edge center.
+ */
+ size_t edgeIndex = representationLFE.localCoefficients().localKey(i).subEntity();
+ auto edgeGeometry = element.template subEntity<1>(edgeIndex).geometry();
+ auto edgeEntity = element.template subEntity<1>(edgeIndex);
+
+ /**
+ * @brief Get the right orientation for tangent and normal vectors on current edge.
+ */
+ const auto &refElement = referenceElement(element);
+ // Local vertex indices within the element
+ auto localV0 = refElement.subEntity(edgeIndex, gridDim-1, 0, gridDim);
+ auto localV1 = refElement.subEntity(edgeIndex, gridDim-1, 1, gridDim);
+
+ // Global vertex indices within the grid
+ auto globalV0 = indexSet.subIndex(element, localV0, gridDim);
+ auto globalV1 = indexSet.subIndex(element, localV1, gridDim);
+
+ double edgeOrientation = 1.0;
+
+ if ((localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1))
+ edgeOrientation = -1.0;
+
+ // Get tangent vector on current edge
+ Dune::FieldVector<RT,2> tmp = (geometry.corner(localV1) - geometry.corner(localV0)) / abs(edgeGeometry.volume());
+
+ //convert to FieldMatrix
+ Dune::FieldMatrix<RT,2,1> edgeTangent = {tmp[0], tmp[1]};
+ edgeTangent *= edgeOrientation;
+
+ /**
+ * @brief Get edge normal by rotating the edge-tangent by by pi/2 clockwise
+ */
+ Dune::FieldMatrix<RT,2,1> edgeNormal = {edgeTangent[1], -1.0*edgeTangent[0]};
+
+ /**
+ * @brief Evaluate Jacobian of localfunctions_ at edge-vertices and edge-midpoints.
+ */
+ typename DeformationFunction::DerivativeType whJacobianValueCenter = deformationFunction.evaluateDerivative(geometry.local(edgeGeometry.center())); // Q: Are these transformations correct?
+ typename DeformationFunction::DerivativeType whJacobianValueFirstVertex = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(localV0)));
+ typename DeformationFunction::DerivativeType whJacobianValueSecondVertex = deformationFunction.evaluateDerivative(geometry.local(geometry.corner(localV1)));
+
+ // @OS: Is there a elementwise/hadamard matrix multiplication in Dune?
+ Dune::FieldMatrix<RT,3,gridDim> normalComponent(0);
+ Dune::FieldMatrix<RT,3,gridDim> tangentialComponent(0);
+
+ auto normalComponentFactor = 0.5 * (whJacobianValueFirstVertex + whJacobianValueSecondVertex) * edgeNormal;
+ auto tangentialComponentFactor = whJacobianValueCenter * edgeTangent;
+
+ /**
+ * @brief Compute normal and tangential component on current edge.
+ */
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ {
+ normalComponent[k][l] = normalComponentFactor[k][0]*edgeNormal[l];
+ tangentialComponent[k][l] = tangentialComponentFactor[k][0]*edgeTangent[l];
+ }
+ discreteGradientCoefficients[i] = normalComponent + tangentialComponent;
+
+ }
+ } //end of P2-indices
+ }
+
+
+
+} // end namespace Dune::GFE::Impl
+
+#endif
diff --git a/dune/gfe/functions/discretekirchhoffbendingisometry.hh b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a7296ae5869fe5b81df2b661b0a2d5b4b8b314e6
--- /dev/null
+++ b/dune/gfe/functions/discretekirchhoffbendingisometry.hh
@@ -0,0 +1,315 @@
+#ifndef DUNE_GFE_FUNCTIONS_DISCRETEKIRCHHOFFBENDINGISOMETRY_HH
+#define DUNE_GFE_FUNCTIONS_DISCRETEKIRCHHOFFBENDINGISOMETRY_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+#include <dune/gfe/functions/localisometrycomponentfunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+#include <dune/gfe/tensor3.hh>
+
+#include <dune/functions/functionspacebases/cubichermitebasis.hh>
+
+namespace Dune::GFE
+{
+ /** \brief Interpolate a Discrete Kirchhoff finite element function from a set of RealTuple x Rotation coefficients
+ *
+ * The Discrete Kirchhoff finite element is a Hermite-type element.
+ * The Rotation component of the coefficient set represents the deformation gradient.
+ *
+ * The class stores both a global coefficient vector 'coefficients_'
+ * as well as a local coefficient vector 'localHermiteCoefficients_'
+ * that is used for local evaluation after binding to an element.
+ *
+ *
+ * \tparam DiscreteKirchhoffBasis The basis used to compute function values
+ * \tparam CoefficientBasis The basis used to index the coefficient vector.
+ * \tparam Coefficients The container of global coefficients
+ */
+ template <class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients>
+ class DiscreteKirchhoffBendingIsometry
+ {
+ using GridView = typename DiscreteKirchhoffBasis::GridView;
+ using ctype = typename GridView::ctype;
+
+ public:
+ constexpr static int gridDim = GridView::dimension;
+ using Coefficient = typename Coefficients::value_type;
+ using RT = typename Coefficient::ctype;
+
+ static constexpr int components_ = 3;
+
+ typedef Dune::GFE::LocalProjectedFEFunction<GridView::dimension, ctype, typename CoefficientBasis::LocalView::Tree::FiniteElement, Dune::GFE::ProductManifold<RealTuple<RT,components_>, Rotation<RT,components_> > > LocalInterpolationRule;
+
+ /** \brief The type used for derivatives */
+ typedef Dune::FieldMatrix<RT, components_, gridDim> DerivativeType;
+
+ /** \brief The type used for discrete Hessian */
+ typedef Tensor3<RT, components_, gridDim, gridDim> discreteHessianType;
+
+
+ /** \brief Constructor
+ * \param basis An object of Hermite basis type
+ * \param coefficientBasis An object of type LagrangeBasis<GridView,1>
+ * \param globalIsometryCoefficients Values and derivatives of the function at the Lagrange points
+ */
+ DiscreteKirchhoffBendingIsometry(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
+ const CoefficientBasis& coefficientBasis,
+ Coefficients& globalIsometryCoefficients)
+ : basis_(discreteKirchhoffBasis),
+ coefficientBasis_(coefficientBasis),
+ localView_(basis_),
+ localViewCoefficient_(coefficientBasis_),
+ globalIsometryCoefficients_(globalIsometryCoefficients)
+ {}
+
+ void bind(const typename GridView::template Codim<0>::Entity& element)
+ {
+ localView_.bind(element);
+
+ /** extract the local coefficients from the global coefficient vector */
+ std::vector<Coefficient> localIsometryCoefficients;
+ getLocalIsometryCoefficients(localIsometryCoefficients,element);
+ updateLocalHermiteCoefficients(localIsometryCoefficients,element);
+
+ /**
+ * @brief Get Coefficients of the discrete Gradient
+ *
+ * The discrete Gradient is a special linear combination represented in a [P2]^2 - (Lagrange) space
+ * The coefficients of this linear combination correspond to certain linear combinations of the Gradients of
+ * the deformation (hermite) basis. The coefficients are stored in the form [Basisfunctions x components x gridDim]
+ * in a BlockVector<FieldMatrix<RT, 3, gridDim>>.
+ */
+ Impl::discreteGradientCoefficients(discreteJacobianCoefficients_,rotationLFE_,*this,element);
+ }
+
+ /** bind to an element and update the coefficient member variables for a set of new local coefficients */
+ void bind(const typename GridView::template Codim<0>::Entity& element,const Coefficients& newLocalCoefficients)
+ {
+ this->bind(element);
+
+ // Update the global isometry coefficients.
+ for (unsigned int vertex=0; vertex<3; ++vertex)
+ {
+ size_t localIdx = localViewCoefficient_.tree().localIndex(vertex);
+ size_t globalIdx = localViewCoefficient_.index(localIdx);
+ globalIsometryCoefficients_[globalIdx] = newLocalCoefficients[vertex];
+ }
+
+ // Update the local hermite coefficients.
+ updateLocalHermiteCoefficients(newLocalCoefficients,element);
+
+ // After setting a new local configurarion, the coefficients of the discrete Gradient need to be updated as well.
+ Impl::discreteGradientCoefficients(discreteJacobianCoefficients_,rotationLFE_,*this,element);
+ }
+
+ /** \brief The type of the element we are bound to */
+ Dune::GeometryType type() const
+ {
+ return localView_.element().type();
+ }
+
+
+ /** \brief Evaluate the function */
+ auto evaluate(const Dune::FieldVector<ctype, gridDim>& local) const
+ {
+ const auto &localFiniteElement = localView_.tree().child(0).finiteElement();
+
+ // Evaluate the shape functions
+ std::vector<FieldVector<double, 1> > values;
+ localFiniteElement.localBasis().evaluateFunction(local, values);
+
+ FieldVector<RT,components_> result(0);
+
+ for(size_t i=0; i<localFiniteElement.size(); i++)
+ result.axpy(values[i][0], localHermiteCoefficients_[i]);
+
+ return (RealTuple<RT, components_>)result;
+ }
+
+ /** \brief Evaluate the derivative of the function */
+ DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, gridDim>& local) const
+ {
+ const auto &localFiniteElement = localView_.tree().child(0).finiteElement();
+
+ const auto jacobianInverse = localView_.element().geometry().jacobianInverse(local);
+ std::vector<FieldMatrix<double,1,gridDim> > jacobianValues;
+ localFiniteElement.localBasis().evaluateJacobian(local, jacobianValues);
+ for (size_t i = 0; i<jacobianValues.size(); i++)
+ jacobianValues[i] = jacobianValues[i] * jacobianInverse;
+
+ DerivativeType result(0);
+
+ for(size_t i=0; i<localFiniteElement.size(); i++)
+ for (unsigned int k = 0; k<components_; k++)
+ for (unsigned int j = 0; j<gridDim; ++j)
+ result[k][j] += (jacobianValues[i][0][j] * localHermiteCoefficients_[i][k]);
+
+ return result;
+ }
+
+ /** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
+ * \param local Local coordinates in the reference element where to evaluate the derivative
+ * \param q Value of the local gfe function at 'local'. If you provide something wrong here the result will be wrong, too!
+ */
+ DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, gridDim>& local,
+ const Coefficient& q) const
+ {
+ return evaluateDerivative(local);
+ }
+
+
+
+ /** \brief Evaluate the discrete gradient operator (This quantity lives in the P2-rotation space) */
+ DerivativeType evaluateDiscreteGradient(const Dune::FieldVector<ctype, gridDim>& local) const
+ {
+ //Get values of the P2-Basis functions on local point
+ std::vector<FieldVector<double,1> > basisValues;
+ rotationLFE_.localBasis().evaluateFunction(local, basisValues);
+
+ DerivativeType discreteGradient(0);
+
+ for (int k=0; k<components_; k++)
+ for (int l=0; l<gridDim; l++)
+ for (std::size_t i=0; i<rotationLFE_.size(); i++)
+ discreteGradient[k][l] += discreteJacobianCoefficients_[i][k][l]*basisValues[i];
+
+ return discreteGradient;
+ }
+
+
+ /** \brief Evaluate the discrete hessian operator (This quantity lives in the P2-rotation space) */
+ discreteHessianType evaluateDiscreteHessian(const Dune::FieldVector<ctype, gridDim>& local) const
+ {
+ // Get gradient values of the P2-Basis functions on local point
+ const auto geometryJacobianInverse = localView_.element().geometry().jacobianInverse(local);
+ std::vector<FieldMatrix<double,1,gridDim> > gradients;
+ rotationLFE_.localBasis().evaluateJacobian(local, gradients);
+
+ for (size_t i=0; i< gradients.size(); i++)
+ gradients[i] = gradients[i] * geometryJacobianInverse;
+
+ discreteHessianType discreteHessian(0);
+
+ for (int k=0; k<components_; k++)
+ for (int l=0; l<gridDim; l++)
+ for (std::size_t i=0; i<rotationLFE_.size(); i++)
+ {
+ discreteHessian[k][l][0] += discreteJacobianCoefficients_[i][k][l]*gradients[i][0][0];
+ discreteHessian[k][l][1] += discreteJacobianCoefficients_[i][k][l]*gradients[i][0][1];
+ }
+
+ return discreteHessian;
+ }
+
+
+ //! Obtain the grid view that the basis is defined on
+ const GridView &gridView() const
+ {
+ return basis_.gridView();
+ }
+
+ void updateglobalIsometryCoefficients(const Coefficients& newCoefficients)
+ {
+ globalIsometryCoefficients_ = newCoefficients;
+ }
+
+ auto getDiscreteJacobianCoefficients()
+ {
+ return discreteJacobianCoefficients_;
+ }
+
+ /**
+ Update the local hermite basis coefficient vector
+ */
+ void updateLocalHermiteCoefficients(const Coefficients& localIsometryCoefficients,const typename GridView::template Codim<0>::Entity& element)
+ {
+ localViewCoefficient_.bind(element);
+
+ //Create a LocalProjected-Finite element from the local coefficients used for interpolation.
+ auto P1LagrangeLFE = localViewCoefficient_.tree().finiteElement();
+ LocalInterpolationRule localPBfunction(P1LagrangeLFE,localIsometryCoefficients);
+
+ /**
+ Interpolate into the local hermite space for each component.
+ */
+ const auto &localHermiteFiniteElement = localView_.tree().child(0).finiteElement();
+
+ for(size_t k=0; k<components_; k++)
+ {
+ std::vector<RT> hermiteComponentCoefficients;
+ Dune::GFE::Impl::LocalIsometryComponentFunction<double,LocalInterpolationRule> localIsometryComponentFunction(localPBfunction,k);
+
+ localHermiteFiniteElement.localInterpolation().interpolate(localIsometryComponentFunction,hermiteComponentCoefficients);
+
+ for(size_t i=0; i<localHermiteFiniteElement.size(); i++)
+ localHermiteCoefficients_[i][k] = hermiteComponentCoefficients[i];
+ }
+ }
+
+ /** Extract the local isometry coefficients. */
+ void getLocalIsometryCoefficients(std::vector<Coefficient>& in,const typename GridView::template Codim<0>::Entity& element)
+ {
+ localViewCoefficient_.bind(element);
+
+ in.resize(3);
+ for (unsigned int vertex = 0; vertex<3; ++vertex)
+ {
+ size_t localIdx = localViewCoefficient_.tree().localIndex(vertex);
+ size_t globalIdx = localViewCoefficient_.index(localIdx);
+ in[vertex] = globalIsometryCoefficients_[globalIdx];
+ }
+ }
+
+ /** \brief Return the representation LocalFiniteElement for the rotation space */
+ auto const &getRotationFiniteElement() const
+ {
+ return rotationLFE_;
+ }
+
+ /** \brief Return the discrete Jacobian coefficients */
+ BlockVector<FieldMatrix<RT, components_, gridDim> > getDiscreteJacobianCoefficients() const
+ {
+ return discreteJacobianCoefficients_;
+ }
+
+ /** \brief Get the i'th base coefficient. */
+ Coefficients coefficient(int i) const
+ {
+ return globalIsometryCoefficients_[i];
+ }
+
+ private:
+
+ /** \brief The Lagrange basis used to assign manifold-valued degrees of freedom */
+ const CoefficientBasis& coefficientBasis_;
+
+ /** \brief The hermite basis used to represent deformations */
+ const DiscreteKirchhoffBasis& basis_;
+
+ /** \brief LocalFiniteElement (P2-Lagrange) used to represent the discrete Jacobian (rotation) on the current element. */
+ Dune::LagrangeSimplexLocalFiniteElement<ctype, double, gridDim, 2> rotationLFE_;
+
+ /** \brief The coefficients of the discrete Gradient/Hessian operator */
+ BlockVector<FieldMatrix<RT, components_, gridDim> > discreteJacobianCoefficients_;
+
+ mutable typename DiscreteKirchhoffBasis::LocalView localView_;
+ mutable typename CoefficientBasis::LocalView localViewCoefficient_;
+
+ /** \brief The global coefficient vector */
+ Coefficients& globalIsometryCoefficients_;
+
+ static constexpr int localHermiteBasisSize_ = Dune::Functions::Impl::CubicHermiteLocalCoefficients<gridDim,true>::size();
+
+ /** \brief The local coefficient vector of the hermite basis. */
+ std::array<Dune::FieldVector<RT,components_>,localHermiteBasisSize_> localHermiteCoefficients_;
+
+ };
+
+} // end namespace Dune::GFE
+#endif // DUNE_GFE_FUNCTIONS_DISCRETEKIRCHHOFFBENDINGISOMETRY_HH
diff --git a/dune/gfe/functions/localisometrycomponentfunction.hh b/dune/gfe/functions/localisometrycomponentfunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..558bbab6d6cef2d84ddc533767c7ecaeb4eaba9a
--- /dev/null
+++ b/dune/gfe/functions/localisometrycomponentfunction.hh
@@ -0,0 +1,48 @@
+#ifndef DUNE_GFE_FUNCTIONS_LOCALISOMETRYCOMPONENTFUNCTION_HH
+#define DUNE_GFE_FUNCTIONS_LOCALISOMETRYCOMPONENTFUNCTION_HH
+
+#include <dune/gfe/bendingisometryhelper.hh>
+
+namespace Dune::GFE::Impl
+{
+ /**
+ * \brief Local isometry component function that can be used as a Dune::Functions::DifferentiableFunction
+ *
+ * \tparam ctype Type used for coordinates on the reference element
+ * \tparam LocalInterpolation Local function mapping into a product manifold of type RealTuple x Rotations
+ */
+ template<class ctype, class LocalInterpolation>
+ class LocalIsometryComponentFunction
+ {
+ public:
+ LocalIsometryComponentFunction(LocalInterpolation& localInterpolation,
+ const int component)
+ : localInterpolation_(localInterpolation),
+ component_(component)
+ {}
+
+ auto operator() (const Dune::FieldVector<ctype,2>& x) const
+ {
+ return localInterpolation_.evaluate(x)[Dune::Indices::_0].globalCoordinates()[component_];
+ }
+
+ /**
+ * \brief Obtain derivative function
+ */
+ friend auto derivative(const LocalIsometryComponentFunction& p)
+ {
+ return [&p](const Dune::FieldVector<ctype,2>& x)
+ {
+ auto derivativeQuaternion = p.localInterpolation_.evaluate(x)[Dune::Indices::_1];
+ auto derivativeMatrix = Rotation<ctype,3>::quaternionToMatrix(derivativeQuaternion);
+ auto reducedDerivativeMatrix = cancelLastMatrixColumn(derivativeMatrix);
+ return reducedDerivativeMatrix[p.component_];
+ };
+ }
+ LocalInterpolation& localInterpolation_;
+ const int component_;
+ };
+
+} // end namespace Dune::GFE::Impl
+
+#endif //#ifndef DUNE_GFE_FUNCTIONS_LOCALISOMETRYCOMPONENTFUNCTION_HH
diff --git a/problems/L-clamped-Plate.py b/problems/L-clamped-Plate.py
new file mode 100644
index 0000000000000000000000000000000000000000..52d9bf820f710f5f0b71d091833ce8d8f6ae5bf9
--- /dev/null
+++ b/problems/L-clamped-Plate.py
@@ -0,0 +1,96 @@
+import math
+
+class ParameterSet(dict):
+ def __init__(self, *args, **kwargs):
+ super(ParameterSet, self).__init__(*args, **kwargs)
+ self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+#############################################
+# Paths
+#############################################
+parameterSet.resultPath = '/home/klaus/Desktop/Dune-Hermite/dune-gfe/outputs_L-clamped-Plate'
+parameterSet.instrumentedPath = '/home/klaus/Desktop/Dune-Hermite/dune-gfe/outputs_L-clamped-Plate/instrumented'
+parameterSet.executablePath = '/home/klaus/Desktop/Dune-Hermite/dune-gfe/build-cmake/src'
+parameterSet.baseName= 'L-clamped-Plate'
+
+#############################################
+# Grid parameters
+#############################################
+nX=1
+nY=1
+
+parameterSet.structuredGrid = 'simplex'
+parameterSet.lower = '0 0'
+parameterSet.upper = '4 4'
+parameterSet.elements = str(nX)+' '+ str(nY)
+
+parameterSet.macroGridLevel = 3
+#############################################
+# Solver parameters
+#############################################
+# Choose solver: "RNHM" (Default:Riemannian Newton with Hessian modification), "RiemannianTR" (Riemannain Trust-region method)
+parameterSet.Solver = "RNHM"
+# Tolerance of the multigrid solver
+parameterSet.tolerance = 1e-12
+# Maximum number of multigrid iterations
+parameterSet.maxProximalNewtonSteps = 100
+# Initial regularization
+parameterSet.initialRegularization = 1
+# Measure convergence
+parameterSet.instrumented = 0
+
+############################
+# Problem specifications
+############################
+# Dimension of the domain (only used for numerical-test python files)
+parameterSet.dim = 2
+
+# Dimension of the target space
+parameterSet.targetDim = 3
+
+parameterSet.targetSpace = 'BendingIsometry'
+#############################################
+# Options
+#############################################
+# Write discrete solution as .vtk-File
+parameterSet.writeVTK = 1
+
+# Write Dof-Vector to .txt-file
+parameterSet.writeDOFvector = 0
+
+#############################################
+# Dirichlet boundary indicator
+#############################################
+def dirichlet_indicator(x) :
+ if( (x[0] <= 0.001) or (x[1]<=0.001)):
+ return True
+ else:
+ return False
+
+
+#############################################
+# Initial iterate function
+#############################################
+def f(x):
+ return [x[0], x[1], 0]
+
+
+def df(x):
+ return ((1,0),
+ (0,1),
+ (0,0))
+
+
+fdf = (f, df)
+
+#############################################
+# Force
+############################################
+def force(x):
+ return [0, 0, 0.025]
+
+
+
+
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 658ef54c6c2c21aa88a712195bef103147ea1fa9..68deff572ca4f35264eba4b8d19b385b7c9c945b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -44,6 +44,11 @@ add_executable("cosserat-continuum-3d-in-3d-2-1" cosserat-continuum.cc)
set_property(TARGET "cosserat-continuum-3d-in-3d-2-1" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=3; WORLD_DIM=3; LFE_ORDER=2; GFE_ORDER=1")
set(programs cosserat-continuum-3d-in-3d-2-1 ${programs})
+if(dune-functions_VERSION VERSION_GREATER_EQUAL 2.11.0)
+ add_executable("bending-plate" bending-plate.cc)
+ set(programs bending-plate ${programs})
+endif()
+
foreach(_program ${programs})
target_link_dune_default_libraries(${_program})
add_dune_pythonlibs_flags(${_program})
diff --git a/src/bending-plate.cc b/src/bending-plate.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6edfea3ada4f781068ac87e8da0248b4907b83da
--- /dev/null
+++ b/src/bending-plate.cc
@@ -0,0 +1,337 @@
+#include <config.h>
+#include <memory>
+#include <array>
+#include <math.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.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/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/composedgridfunction.hh>
+#include <dune/functions/functionspacebases/cubichermitebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+#include <dune/fufem/dunepython.hh>
+#include <dune/fufem/discretizationerror.hh>
+
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+#include <dune/gfe/functions/discretekirchhoffbendingisometry.hh>
+#include <dune/gfe/functions/embeddedglobalgfefunction.hh>
+#include <dune/gfe/functions/localprojectedfefunction.hh>
+#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/assemblers/geodesicfeassembler.hh>
+#include <dune/gfe/assemblers/discretekirchhoffbendingenergy.hh>
+#include <dune/gfe/assemblers/forceenergy.hh>
+#include <dune/gfe/assemblers/sumenergy.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+#include <dune/gfe/riemannianpnsolver.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+
+#include <dune/gmsh4/gmsh4reader.hh>
+#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
+
+#include <dune/vtk/vtkwriter.hh>
+#include <dune/vtk/writers/unstructuredgridwriter.hh>
+#include <dune/vtk/datacollectors/continuousdatacollector.hh>
+#include <dune/vtk/datacollectors/discontinuousdatacollector.hh>
+#include <dune/vtk/datacollectors/quadraticdatacollector.hh>
+#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
+
+const int dim = 2;
+
+using namespace Dune;
+
+/**
+ Source file for the computation of bending isometries.
+
+ Usage:
+ ./bending-isometries <python path> <python module without extension>
+ */
+
+int main(int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+ Dune::Timer globalTimer;
+
+ if (argc < 3)
+ DUNE_THROW(Exception, "Usage: ./bending-isometries <python path> <python module without extension>");
+
+ // Start Python interpreter
+ Python::start();
+ auto pyMain = Python::main();
+ pyMain.runStream()
+ << std::endl << "import math"
+ << std::endl << "import sys"
+ << std::endl << "sys.path.append('" << argv[1] << "')" << std::endl;
+ auto pyModule = pyMain.import(argv[2]);
+
+ std::cout << "Current path is " << std::filesystem::current_path() << '\n';
+ std::filesystem::path file_path = (__FILE__);
+ std::cout<< "File path: " << file_path<<std::endl;
+ std::cout << "dir_path: " << file_path.parent_path() << std::endl;
+ std::string dir_path = file_path.parent_path();
+
+ // parse data file
+ ParameterTree parameterSet;
+ pyModule.get("parameterSet").toC(parameterSet);
+
+ // read possible further parameters from the command line
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+ // Print all parameters, to make them appear in the log file
+ std::cout << "Executable: bending-isometries, with parameters:" << std::endl;
+ parameterSet.report();
+
+ bool PRINT_DEBUG = parameterSet.get<bool>("print_debug", 0);
+
+ /////////////////////////////////////////
+ // Create the grid
+ /////////////////////////////////////////
+ using GridType = UGGrid<dim>;
+
+ std::shared_ptr<GridType> grid;
+ FieldVector<double,dim> lower(0), upper(1);
+ std::array<unsigned int,dim> elementsArray;
+
+ std::string structuredGridType = parameterSet["structuredGrid"];
+ if (structuredGridType != "false" )
+ {
+ lower = parameterSet.get<FieldVector<double,dim> >("lower");
+ upper = parameterSet.get<FieldVector<double,dim> >("upper");
+ elementsArray = parameterSet.get<std::array<unsigned int,dim> >("elements");
+
+ if (structuredGridType == "simplex")
+ grid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elementsArray);
+ else if (structuredGridType == "cube")
+ grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elementsArray);
+ else
+ DUNE_THROW(Exception, "Unknown structured grid type '" << structuredGridType << "' found!");
+ } else {
+ std::cout << "Read GMSH grid." << std::endl;
+ std::string gridPath = parameterSet.get<std::string>("gridPath");
+ std::string gridFile = parameterSet.get<std::string>("gridFile");
+ GridFactory<GridType> factory;
+ Gmsh4::LagrangeGridCreator creator{factory};
+ Gmsh4Reader reader{creator};
+ reader.read(gridPath + "/" + gridFile);
+ grid = factory.createGrid();
+ }
+
+ const int macroGridLevel = parameterSet.get<int>("macroGridLevel");
+ grid->globalRefine(macroGridLevel-1);
+
+ using GridView = typename GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ ///////////////////////////////////////////////////////
+ // Set up the function spaces
+ ///////////////////////////////////////////////////////
+ // General coefficient vector of a Discrete Kirchhoff deformation function
+ using VectorSpaceCoefficients = BlockVector<FieldVector<double,3> >;
+
+ // Coefficient vector of a Discrete Kirchhoff deformation function that is constrained to be an isometry.
+ // we need both 'double' and 'adouble' versions.
+ using Coefficient = GFE::ProductManifold<GFE::RealTuple<double,3>, GFE::Rotation<double,3> >;
+ using IsometryCoefficients = std::vector<Coefficient>;
+ using ACoefficient = typename Coefficient::template rebind<adouble>::other;
+ using AIsometryCoefficients = std::vector<ACoefficient>;
+
+ using namespace Functions::BasisFactory;
+ auto deformationBasis = makeBasis(gridView,
+ power<3>(reducedCubicHermite(),
+ blockedInterleaved()));
+ using DeformationBasis = decltype(deformationBasis);
+
+
+ // The next basis is used to assign (nonlinear) degrees of freedom to the grid vertices.
+ // The actual basis function values are never used.
+ using CoefficientBasis = Functions::LagrangeBasis<GridView, 1>;
+ CoefficientBasis coefficientBasis(gridView);
+
+ // A basis for the tangent space (used to set DirichletNodes)
+ auto tangentBasis = makeBasis(gridView,
+ power<Coefficient::TangentVector::dimension>(
+ lagrange<1>(),
+ blockedInterleaved()));
+
+
+ // Print some information on the grid and degrees of freedom.
+ std::cout << "Coefficient::TangentVector::dimension: " << Coefficient::TangentVector::dimension<< std::endl;
+ std::cout << "Number of Elements in the grid: " << gridView.size(0)<< std::endl;
+ std::cout << "Number of Nodes in the grid: " << gridView.size(dim)<< std::endl;
+ std::cout << "deformationBasis.size(): " << deformationBasis.size() << std::endl;
+ std::cout << "deformationBasis.dimension(): " << deformationBasis.dimension() << std::endl;
+ std::cout << "Degrees of Freedom: " << deformationBasis.dimension() << std::endl;
+
+ ///////////////////////////////////////////
+ // Read Dirichlet values
+ ///////////////////////////////////////////
+ BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+ const typename GridView::IndexSet &indexSet = gridView.indexSet();
+ BitSetVector<Coefficient::TangentVector::dimension> dirichletNodes(tangentBasis.size(), false); //tangentBasis.size()=coefficientBasis.size()
+
+ // Make Python function that computes which vertices are on the Dirichlet boundary,
+ // based on the vertex positions.
+ auto dirichletIndicatorFunction = Python::make_function<bool>(pyModule.get("dirichlet_indicator"));
+
+ // If we want to clamp DOFs inside the domain, we cannot use 'BoundaryPatch'
+ // and 'constructBoundaryDofs'. This is a workaround for now.
+ for (auto &&vertex : vertices(gridView))
+ {
+ dirichletVertices[indexSet.index(vertex)] = dirichletIndicatorFunction(vertex.geometry().corner(0));
+
+ if(dirichletIndicatorFunction(vertex.geometry().corner(0)))
+ {
+ dirichletNodes[indexSet.index(vertex)] = true;
+ if(PRINT_DEBUG)
+ std::cout << "Dirichlet Vertex with coordinates:" << vertex.geometry().corner(0) << std::endl;
+ }
+ }
+
+ ///////////////////////////////////////////
+ // Get initial Iterate
+ ///////////////////////////////////////////
+ auto pythonInitialIterate = Python::makeDifferentiableFunction<FieldVector<double,3>(FieldVector<double,2>)>(pyModule.get("f"), pyModule.get("df"));
+ VectorSpaceCoefficients x(deformationBasis.size());
+ interpolate(deformationBasis, x, pythonInitialIterate);
+
+
+
+ // We need to setup DiscreteKirchhoffBendingIsometry with a coefficient
+ // vector of ctype 'adouble' while the solver gets a coefficient vector
+ // of ctype 'double'.
+ IsometryCoefficients isometryCoefficients(coefficientBasis.size());
+ AIsometryCoefficients isometryCoefficients_adouble(coefficientBasis.size());
+
+ using namespace Dune::GFE::Impl;
+ // Copy the current iterate into a data type that encapsulates the isometry constraint
+ // i.e. convert coefficient data structure from 'VectorSpaceCoefficients' to 'IsometryCoefficients'
+ vectorToIsometryCoefficientMap(deformationBasis,coefficientBasis,x,isometryCoefficients);
+ vectorToIsometryCoefficientMap(deformationBasis,coefficientBasis,x,isometryCoefficients_adouble);
+
+ // Create a DiscreteKirchhoffBendingIsometry.
+ // This serves as the deformation function.
+ using LocalDKFunction = GFE::DiscreteKirchhoffBendingIsometry<DeformationBasis, CoefficientBasis, AIsometryCoefficients>;
+ LocalDKFunction localDKFunction(deformationBasis, coefficientBasis, isometryCoefficients_adouble);
+
+ // Read the force term.
+ auto pythonForce = Python::make_function<FieldVector<double,3> >(pyModule.get("force"));
+ auto forceGVF = Dune::Functions::makeGridViewFunction(pythonForce, gridView);
+ auto localForce = localFunction(forceGVF);
+
+ ////////////////////////////////////////////////////////////////////////
+ // Create an assembler for the Discrete Kirchhoff Energy Functional (using ADOL-C)
+ ////////////////////////////////////////////////////////////////////////
+
+ // Setup nonconforming energy and assembler - only option for this minimal example.
+ auto forceEnergy = std::make_shared<GFE::ForceEnergy<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient> >(localDKFunction, localForce);
+ auto localEnergy_nonconforming = std::make_shared<GFE::DiscreteKirchhoffBendingEnergy<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient> >(localDKFunction);
+
+ auto sumEnergy = std::make_shared<GFE::SumEnergy<CoefficientBasis, GFE::RealTuple<adouble,3>, GFE::Rotation<adouble,3> > >();
+ sumEnergy->addLocalEnergy(localEnergy_nonconforming);
+ sumEnergy->addLocalEnergy(forceEnergy);
+
+ auto localGFEADOLCStiffness_nonconforming= std::make_shared<Dune::GFE::LocalGeodesicFEADOLCStiffness<CoefficientBasis, Coefficient> >(sumEnergy);
+ std::shared_ptr<Dune::GFE::GeodesicFEAssembler<CoefficientBasis, Coefficient> > assembler_nonconforming;
+ assembler_nonconforming = std::make_shared<Dune::GFE::GeodesicFEAssembler<CoefficientBasis, Coefficient> >(coefficientBasis, localGFEADOLCStiffness_nonconforming);
+
+ // Create a solver:
+ // * Riemannian Newton with Hessian modification
+ // * Riemannian Trust-region
+ Dune::GFE::RiemannianProximalNewtonSolver<CoefficientBasis, Coefficient> RNHMsolver;
+ Dune::GFE::RiemannianTrustRegionSolver<CoefficientBasis, Coefficient> RTRsolver;
+
+ std::string Solver_name = parameterSet.get<std::string>("Solver", "RNHM");
+ double numerical_energy; //final discrete energy
+
+ if(Solver_name == "RNHM")
+ {
+ RNHMsolver.setup(*grid,
+ &(*assembler_nonconforming),
+ isometryCoefficients,
+ dirichletNodes,
+ parameterSet);
+
+ RNHMsolver.solve();
+ isometryCoefficients = RNHMsolver.getSol();
+ numerical_energy = RNHMsolver.getStatistics().finalEnergy;
+ } else if (Solver_name =="RiemannianTR")
+ {
+ std::cout << "Using Riemannian Trust-region method for energy minimization." << std::endl;
+ RTRsolver.setup(*grid,
+ &(*assembler_nonconforming),
+ isometryCoefficients,
+ dirichletNodes,
+ parameterSet);
+
+ RTRsolver.solve();
+ isometryCoefficients = RTRsolver.getSol();
+ numerical_energy = RTRsolver.getStatistics().finalEnergy;
+ } else
+ DUNE_THROW(Dune::Exception, "Unknown Solver type for bending isometries.");
+
+
+ /** Convert coefficient data structure from 'IsometryCoefficients' back to 'VectorSpaceCoefficients' */
+ VectorSpaceCoefficients x_out(deformationBasis.size());
+ isometryToVectorCoefficientMap(deformationBasis,coefficientBasis,x_out,isometryCoefficients);
+
+ ////////////////////////////////
+ // Output result
+ ////////////////////////////////
+ std::string baseNameDefault = "bending-isometries-";
+ std::string baseName = parameterSet.get("baseName", baseNameDefault);
+ std::string resultFileName = parameterSet.get("resultPath", "")
+ + "/" + baseName
+ + "-level" + std::to_string(parameterSet.get<int>("macroGridLevel"));
+
+
+ if (parameterSet.get<bool>("writeVTK", 1))
+ {
+ std::cout << "write VTK to Filename: " << resultFileName << std::endl;
+
+ // Compute the displacement from the deformation.
+ // interpolate the identity and subtract from the coefficient vector.
+ VectorSpaceCoefficients identity(deformationBasis.size());
+ IdentityGridEmbedding<double> identityGridEmbedding;
+ interpolate(deformationBasis, identity, identityGridEmbedding);
+
+ // Compute the displacement
+ auto displacement = x_out;
+ displacement -= identity;
+
+ auto deformationFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(deformationBasis, x_out);
+ auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(deformationBasis, displacement);
+
+ // Use VTK writer from dune-vtk.
+ // Setup a DataCollector of order 3.
+ Dune::Vtk::LagrangeDataCollector<GridView, 3> lagrangeDataCollector(gridView);
+ Dune::Vtk::UnstructuredGridWriter duneVTKwriter(lagrangeDataCollector, Dune::Vtk::FormatTypes::ASCII, Vtk::DataTypes::FLOAT32);
+
+ // Write discrete displacement
+ duneVTKwriter.addPointData(displacementFunction, Dune::Vtk::FieldInfo{"Displacement",3, Vtk::RangeTypes::VECTOR});
+
+ duneVTKwriter.write(resultFileName);
+ }
+
+ std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
+ return 0;
+}