diff --git a/dune/microstructure/energies/CMakeLists.txt b/dune/microstructure/energies/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8c71e5f4c8702a8129b3b10a3423164bbf7f2b71
--- /dev/null
+++ b/dune/microstructure/energies/CMakeLists.txt
@@ -0,0 +1,3 @@
+#install headers
+install(FILES discretekirchhoffbendingenergyprestrained.hh
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/energies)
diff --git a/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d6e46988f5cb1f46813831d270d6dff8e56f54bc
--- /dev/null
+++ b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
@@ -0,0 +1,673 @@
+#ifndef DUNE_MICROSTRUCTURE_ENERGIES_DISCRETEKIRCHHOFFBENDINGENERGYPRESTRAINED_HH
+#define DUNE_MICROSTRUCTURE_ENERGIES_DISCRETEKIRCHHOFFBENDINGENERGYPRESTRAINED_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>
+
+#include <cmath>
+/**
+ * \brief Assemble the discrete Kirchhoff bending energy for a single element.
+ *
+ * The Kirchhoff bending energy consists of two parts:
+ *
+ * 1. The energy of Qhom(II - Beff) where II is the second fundamental form of the surface parametrized by the deformation function
+ * and Beff is a (given) effective prestrain tensor.
+ *
+ * This contribution is split intro three parts which corresponds to the binomial formula (a+b)^2 = a^2 + 2ab + b^2
+ * each term is discretized separately.
+ *
+ * 2. An integral over the scalar product of a forcing term and the discrete deformation
+ * (i.e. the evaluation of localFunction_ ).
+ */
+namespace Dune::GFE
+{
+ template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
+ class DiscreteKirchhoffBendingEnergyPrestrained
+ : 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 LocalDiscreteKirchhoffFunction::DerivativeType DerivativeType;
+ typedef typename Dune::FieldMatrix<RT,2,2> PrestrainType;
+ typedef typename Dune::FieldMatrix<RT,3,3> MatrixType;
+ // 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:
+ DiscreteKirchhoffBendingEnergyPrestrained(LocalDiscreteKirchhoffFunction &localFunction,
+ LocalForce &localForce,
+ const Dune::ParameterTree& parameterSet,
+ Python::Module pyModule)
+ : localFunction_(localFunction),
+ localForce_(localForce),
+ parameterSet_(parameterSet),
+ module_(pyModule)
+ {
+ setupEffectiveQuantities();
+ }
+
+
+
+
+ // Dune::FieldMatrix<RT,3,3> Qhom_(0);
+
+ void setupEffectiveQuantities()
+ {
+
+
+ Dune::FieldMatrix<double,3,3> Qhom_double(0);
+ module_.get("effectiveQuadraticForm").toC<Dune::FieldMatrix<double,3,3>>(Qhom_double);
+ convertFieldMatrix<RT,decltype(Qhom_double)>(Qhom_double,Qhom_);
+ printmatrix(std::cout, Qhom_, "Setup effective Quadratic form (Qhom) as ", "--");
+
+ Dune::FieldMatrix<double,2,2> Beff_double(0);
+ module_.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff_double);
+ convertFieldMatrix<RT,decltype(Beff_double)>(Beff_double,Beff_);
+ printmatrix(std::cout, Beff_, "Setup effective prestrain (Beff) as ", "--");
+
+
+ }
+
+ void resetEffectiveQuantities()
+ {
+
+
+ }
+
+
+
+ // auto applyQhom(const Dune::FieldMatrix<RT,1,3>& A, const Dune::FieldMatrix<RT,1,3>& B) const
+ // {
+ // std::cout << "Apply Qhom" << std::endl;
+
+ // return Qhom_[0][0]*A[0][0]*B[0][0] + Qhom_[0][1]*A[0][1]*B[0][0] + Qhom_[0][2]*A[0][2]*B[0][0]
+ // + Qhom_[1][0]*A[0][0]*B[0][1] + Qhom_[1][1]*A[0][1]*B[0][1] + Qhom_[1][2]*A[0][2]*B[0][1]
+ // + Qhom_[2][0]*A[0][0]*B[0][2] + Qhom_[2][1]*A[0][1]*B[0][2] + Qhom_[2][2]*A[0][2]*B[0][2];
+ // }
+
+
+ /**
+ * @brief
+ *
+ * @param A coefficient vector of first input
+ * @param B coefficient vector of second input
+ * @return auto template<class VectorType>
+ */
+ auto applyQhom(const Dune::FieldVector<RT,3>& A, const Dune::FieldVector<RT,3>& B) const
+ {
+ // std::cout << "Apply Qhom" << std::endl;
+
+ return Qhom_[0][0]*A[0]*B[0] + Qhom_[0][1]*A[1]*B[0] + Qhom_[0][2]*A[2]*B[0]
+ + Qhom_[1][0]*A[0]*B[1] + Qhom_[1][1]*A[1]*B[1] + Qhom_[1][2]*A[2]*B[1]
+ + Qhom_[2][0]*A[0]*B[2] + Qhom_[2][1]*A[1]*B[2] + Qhom_[2][2]*A[2]*B[2];
+ }
+ // template<class Vector>
+ // auto applyQhom(const Vector& A, const Vector& B) const
+ // {
+ // // std::cout << "Apply Qhom" << std::endl;
+
+
+ // return Qhom_[0][0] + Qhom_[0][1]*A[1]*B[0] + Qhom_[0][2]*A[2]*B[0]
+ // + Qhom_[1][0]*A[0]*B[1] + Qhom_[1][1]*A[1]*B[1] + Qhom_[1][2]*A[2]*B[1]
+ // + Qhom_[2][0]*A[0]*B[2] + Qhom_[2][1]*A[1]*B[2] + Qhom_[2][2]*A[2]*B[2];
+ // }
+
+// // Template specialization for Matrix:
+// template<>
+// string toCustomString<std::string>(std::vector<std::string> vec) {
+// // strings
+// }
+
+// template<>
+// auto applyQhom<Dune:>(const Vector& A, const Vector& B) const
+// {
+// // std::cout << "Apply Qhom" << std::endl;
+
+// return Qhom_[0][0] + Qhom_[0][1]*A[1]*B[0] + Qhom_[0][2]*A[2]*B[0]
+// + Qhom_[1][0]*A[0]*B[1] + Qhom_[1][1]*A[1]*B[1] + Qhom_[1][2]*A[2]*B[1]
+// + Qhom_[2][0]*A[0]*B[2] + Qhom_[2][1]*A[1]*B[2] + Qhom_[2][2]*A[2]*B[2];
+// }
+
+
+
+ /** \brief Assemble the energy for a single element */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localSolution) const
+ {
+ RT energy = 0;
+
+
+ // if (parameterSet_.get<bool>("macroscopically_constant_microstructure", 1))
+ /**
+ * @brief Get effective prestrain Tensor
+ */
+ // Dune::FieldVector<RT,3> Beffvec = parameterSet_.get<Dune::FieldVector<RT,3>>("effectivePrestrain", {0.0, 0.0, 0.0});
+ // PrestrainType Beff = {{Beffvec[0], Beffvec[2]}, {Beffvec[2], Beffvec[1]} };
+
+
+
+
+ // printmatrix(std::cout, Beff, "effective prestrain (Beff): ", "--");
+
+ /**
+ * @brief Get effective quadratic form Qhom
+ *
+ * input-vector: [q_1,q_2,q_3,q_12,q_13,q_23]
+ * is assembled into a matrix where the off-diagonal entries are divided by 2 (compare with definition in the paper)
+ * ( q_1 , 0.5*q_12 , 0.5*q_13 )
+ * Q = ( 0.5*q_12 , q_2 , 0.5*q_23 )
+ * ( 0.5*q_13 , 0.5*q_23 , q_3 )
+ */
+ // Dune::FieldVector<RT,6> Qhomvec = parameterSet_.get<Dune::FieldVector<RT,6>>("effectiveQuadraticForm", {1.0, 1.0, 1.0, 0.0, 0.0, 0.0});
+ // MatrixType Qhom = {{Qhomvec[0] , 0.5*Qhomvec[3], 0.5*Qhomvec[4]},
+ // {0.5*Qhomvec[3], Qhomvec[1], 0.5*Qhomvec[5]},
+ // {0.5*Qhomvec[4], 0.5*Qhomvec[5], Qhomvec[2]}};
+
+ // Dune::FieldMatrix<RT,3,3> TestMatrix(0);
+
+ // Dune::FieldMatrix<double,3,3> TestMatrix = parameterSet_.get<Dune::FieldMatrix<double,3,3>>("TestMatrixString");
+ // module_.get("TestMatrix").toC<Dune::FieldMatrix<RT,3,3>>(TestMatrix);
+ // printmatrix(std::cout, TestMatrix, "TestMatrix from pyModule: ", "--");
+
+
+ // Qhom_ = {{Qhomvec[0] , 0.5*Qhomvec[3], 0.5*Qhomvec[4]},
+ // {0.5*Qhomvec[3], Qhomvec[1], 0.5*Qhomvec[5]},
+ // {0.5*Qhomvec[4], 0.5*Qhomvec[5], Qhomvec[2]}};
+ // printmatrix(std::cout, Qhom, "effective quadratic form (Qhom): ", "--");
+
+
+ // print some quantities for debugging purposes.
+ bool PRINT_DEBUG = parameterSet_.get<bool>("print_debug", 0);
+
+
+#if 1
+ // TODO: More serious order selection
+ int quadOrder = 3;
+ // int quadOrder = 6;
+#else
+ if (localFiniteElement.localBasis().order()==1)
+ quadOrder = 2;
+ else
+ quadOrder = (localFiniteElement.type().isSimplex()) ? (localFiniteElement.localBasis().order() - 1) * 2
+ : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
+#endif
+
+ const auto element = localView.element();
+ auto geometry = element.geometry();
+ auto gridView = localFunction_.gridView();
+ const auto &indexSet = gridView.indexSet();
+
+ // bind to current element.
+ localFunction_.bind(element);
+ localForce_.bind(element);
+ /**
+ * @brief We need to update the Coefficients of localFunction_
+ * on the current element.
+ */
+ localFunction_.updateLocalCoefficients(localSolution,element);
+
+
+
+ if(PRINT_DEBUG)
+ {
+ std::cout << " -------------------------------------------- " << std::endl;
+ std::cout << "ELEMENT NUMBER( indexSet.index(element)) ):" << indexSet.index(element) << std::endl;
+ std::cout << "geometry.corner(0): " << geometry.corner(0) << std::endl;
+ std::cout << "geometry.corner(1)): " << geometry.corner(1) << std::endl;
+ std::cout << "geometry.corner(2): " << geometry.corner(2) << std::endl;
+ }
+
+
+ /**
+ * @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 localfunction_ .
+ * The coefficients are stored in the form [Basisfunctions x components x gridDim]
+ * in a BlockVector<FieldMatrix<RT, 3, gridDim>> .
+ */
+ BlockVector<FieldMatrix<RT, 3, gridDim>> discreteJacobianCoefficients;
+ discreteGradientCoefficients(discreteJacobianCoefficients,lagrangeLFE_,localFunction_,element);
+
+
+
+
+
+ /**
+ * @brief Compute harmonic energy contribution:
+ */
+ RT harmonicEnergy = 0;
+
+
+ /**
+ * @brief Mixed version using the binomal formula (a-b)^2 = a^2 - 2ab + b^2
+ * for | II - Beff |^2 and discretize every term individually.
+ */
+ // Gauss-Quadrature:
+ const auto &quadRule = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
+ // Trapezoidal-rule:
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ // Trapezoidal-rule + evaluation on edge centers (only for testing purposes):
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{0.5,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,0.5}, 1.0/6.0}, {{0.5,0.5}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ for (auto&& quadPoint : quadRule)
+ {
+ const auto& quadPos = quadPoint.position();
+ const auto integrationElement = geometry.integrationElement(quadPos);
+
+ //Get values of the P2-Basis functions on current quadrature point
+ std::vector<FieldVector<double,1>> basisValues;
+ lagrangeLFE_.localBasis().evaluateFunction(quadPos, basisValues);
+
+ // Get Jacobians of the P2-Basis functions on current quadrature point
+ // std::vector<FieldMatrix<double, 1, gridDim>> referenceGradients;
+ // lagrangeLFE_.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+ // const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+ // std::vector<FieldVector<RT, gridDim>> gradients(referenceGradients.size());
+ // for (size_t i = 0; i<gradients.size(); i++)
+ // jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+ // New compact way to write Gradient:
+ const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
+ std::vector<FieldMatrix<double,1,gridDim> > gradients;
+ lagrangeLFE_.localBasis().evaluateJacobian(quadPos, gradients);
+
+ for (size_t i=0; i< gradients.size(); i++)
+ gradients[i] = gradients[i] * geometryJacobianInverse;
+
+
+ Tensor3<RT,3,gridDim,gridDim> discreteHessian(0);
+ FieldMatrix<RT, 3, gridDim> discreteGradient(0);
+
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ for (std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ {
+ discreteGradient[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues[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];
+ }
+
+ if(PRINT_DEBUG)
+ {
+ /**
+ * @brief print the (three) 2x2 slices of the discrete Hessian and check for symmetry
+ */
+ FieldMatrix<RT, gridDim, gridDim> discreteHessian_slice1(0);
+ FieldMatrix<RT, gridDim, gridDim> discreteHessian_slice2(0);
+ FieldMatrix<RT, gridDim, gridDim> discreteHessian_slice3(0);
+
+ for (int i=0; i<gridDim; i++)
+ for (int l=0; l<gridDim; l++)
+ {
+ discreteHessian_slice1[i][l] = discreteHessian[0][i][l];
+ discreteHessian_slice2[i][l] = discreteHessian[1][i][l];
+ discreteHessian_slice3[i][l] = discreteHessian[2][i][l];
+ }
+ printmatrix(std::cout, discreteHessian_slice1, "discreteHessian_slice1: ", "--");
+ printmatrix(std::cout, discreteHessian_slice2, "discreteHessian_slice2: ", "--");
+ printmatrix(std::cout, discreteHessian_slice3, "discreteHessian_slice3: ", "--");
+ printmatrix(std::cout, sym(discreteHessian_slice1), "SYM discreteHessian_slice1: ", "--");
+ printmatrix(std::cout, sym(discreteHessian_slice2), "SYM discreteHessian_slice2: ", "--");
+ printmatrix(std::cout, sym(discreteHessian_slice3), "SYM discreteHessian_slice3: ", "--");
+ }
+
+
+
+ /**
+ * @brief Compute the normal vector of the surface parametrized by the
+ * discrete deformation. This is given by the
+ * cross product of partial derivatives
+ */
+ // FieldVector<RT,3> surfaceNormal(0);
+ FieldVector<RT,3> partialX = {discreteGradient[0][0], discreteGradient[1][0], discreteGradient[2][0]};
+ FieldVector<RT,3> partialY = {discreteGradient[0][1], discreteGradient[1][1], discreteGradient[2][1]};
+ // Normalize if requested
+ // auto normX = partialX.two_norm();
+ // auto normY = partialY.two_norm();
+ // partialX /= normX;
+ // partialY /= normY;
+ // surfaceNormal[0] = partialX[1]*partialY[2] - partialX[2]*partialY[1];
+ // surfaceNormal[1] = partialX[2]*partialY[0] - partialX[0]*partialY[2];
+ // surfaceNormal[2] = partialX[0]*partialY[1] - partialX[1]*partialY[0];
+ FieldVector<RT,3> surfaceNormal = {partialX[1]*partialY[2] - partialX[2]*partialY[1],
+ partialX[2]*partialY[0] - partialX[0]*partialY[2],
+ partialX[0]*partialY[1] - partialX[1]*partialY[0]};
+
+
+
+ /**
+ * @brief Compute the second fundamental form of the surface parametrized by the
+ * discrete deformation.
+ * We have to multiply by (-1.0) to get the right sign/normal/orientation.
+ */
+ FieldMatrix<RT,gridDim,gridDim> secondFF(0);
+
+ for(size_t m=0; m<2; m++)
+ for(size_t n=0; n<2; n++)
+ for(size_t k=0; k<3; k++)
+ secondFF[m][n] += -1.0*surfaceNormal[k]*discreteHessian[k][n][m];
+ // To get the right sign:
+ // secondFF = -1.0*secondFF;
+
+
+ //TODO : Add method "applyQhom"
+
+ /**
+ * @brief Compute the first term: Qhom(II)
+ */
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPos);
+
+ RT term1 = 0.0;
+
+ // FieldMatrix<RT,3,3> X_old = tensorToSymCoefficients(discreteHessian);
+ std::vector<FieldVector<RT,3>> X(3);
+ tensorToSymCoefficients(discreteHessian, X);
+
+
+ // printmatrix(std::cout, X_old, "X_old", "--");
+ // for(int k=0; k<3; k++)
+ // {
+ // printvector(std::cout, X[k], "X[k]", "--");
+
+ // }
+
+
+
+ RT term1_old = 0.0;
+
+ for(int k=0; k<3; k++)
+ {
+ // term1 += Qhom[0][0]*pow((2.0*discreteHessian[k][0][0]),2)+ Qhom[0][1]* ... more efficient (TODO) Qhom is symmetric
+ term1_old += Qhom_[0][0]*pow(X[k][0],2) + Qhom_[0][1]*X[k][1]*X[k][0] + Qhom_[0][2]*X[k][2]*X[k][0]
+ + Qhom_[1][0]*X[k][0]*X[k][1] + Qhom_[1][1]*pow(X[k][1],2) + Qhom_[1][2]*X[k][2]*X[k][1]
+ + Qhom_[2][0]*X[k][0]*X[k][2] + Qhom_[2][1]*X[k][1]*X[k][2] + Qhom_[2][2]*pow(X[k][2],2);
+
+ // apply Qhom to each slice of the Hessian...
+ term1 += applyQhom(X[k], X[k]);
+ }
+ // term1 = term1 * 0.5 * weight;
+
+ // std::cout << "term1: " << term1 << std::endl;
+ // std::cout << "term1_old: " << term1_old << std::endl;
+
+ term1 = term1 * weight;
+ term1_old = term1_old * weight;
+
+
+ if(abs(term1-term1_old) > 1e-6)
+ std::cout << "DIFFERENCE term1-term1_old " << term1-term1_old << std::endl;
+
+
+
+ // /**
+ // * @brief Compute the second (mixed) term
+ // */
+ FieldVector<RT,3> secondFFCoefficients = matrixToSymCoefficients(secondFF);
+ FieldVector<RT,3> BeffCoefficients = matrixToSymCoefficients(Beff_);
+
+ RT term2 = 0.0;
+ RT term2_old = 0.0;
+ term2_old += Qhom_[0][0]*BeffCoefficients[0]*secondFFCoefficients[0] + Qhom_[0][1]*BeffCoefficients[1]*secondFFCoefficients[0] + Qhom_[0][2]*BeffCoefficients[2]*secondFFCoefficients[0]
+ + Qhom_[1][0]*BeffCoefficients[0]*secondFFCoefficients[1] + Qhom_[1][1]*BeffCoefficients[1]*secondFFCoefficients[1] + Qhom_[1][2]*BeffCoefficients[2]*secondFFCoefficients[1]
+ + Qhom_[2][0]*BeffCoefficients[0]*secondFFCoefficients[2] + Qhom_[2][1]*BeffCoefficients[1]*secondFFCoefficients[2] + Qhom_[2][2]*BeffCoefficients[2]*secondFFCoefficients[2];
+
+
+
+ // term2 += applyQhom(secondFFCoefficients,BeffCoefficients);
+ term2 += applyQhom(BeffCoefficients,secondFFCoefficients);
+
+ // std::cout << "term2: " << term2 << std::endl;
+ // std::cout << "term2_old: " << term2_old << std::endl;
+
+ term2 = 2.0 * term2 * weight;
+ term2_old = 2.0 * term2_old * weight;
+
+ if(abs(term2-term2_old) > 1e-6)
+ std::cout << "DIFFERENCE term2-term2_old " << term2-term2_old << std::endl;
+
+
+
+ // auto term3 = 0.5 * weight* Beff.frobenius_norm2();
+
+ RT term3 = 0.0;
+ RT term3_old = 0.0;
+
+ term3_old += Qhom_[0][0]*BeffCoefficients[0]*BeffCoefficients[0] + Qhom_[0][1]*BeffCoefficients[1]*BeffCoefficients[0] + Qhom_[0][2]*BeffCoefficients[2]*BeffCoefficients[0]
+ + Qhom_[1][0]*BeffCoefficients[0]*BeffCoefficients[1] + Qhom_[1][1]*BeffCoefficients[1]*BeffCoefficients[1] + Qhom_[1][2]*BeffCoefficients[2]*BeffCoefficients[1]
+ + Qhom_[2][0]*BeffCoefficients[0]*BeffCoefficients[2] + Qhom_[2][1]*BeffCoefficients[1]*BeffCoefficients[2] + Qhom_[2][2]*BeffCoefficients[2]*BeffCoefficients[2];
+ // // term3 = term3 * 0.5 * weight ;
+ term3 += applyQhom(BeffCoefficients,BeffCoefficients);
+
+ // std::cout << "term3: " << term3 << std::endl;
+ // std::cout << "term3_old: " << term3_old << std::endl;
+
+
+ term3 = term3 * weight ;
+ term3_old = term3_old * weight ;
+
+ if(abs(term3-term3_old) > 1e-6)
+ std::cout << "DIFFERENCE term3-term3_old " << term3-term3_old << std::endl;
+
+
+
+
+
+ harmonicEnergy += term1 - term2 + term3;
+ // harmonicEnergy += term1_old - term2_old + term3_old;
+
+ // std::cout << "harmonicEnergy: " << harmonicEnergy << std::endl;
+ }
+ // #endif
+
+
+ #if 0
+ /**
+ * @brief VERSION - 3 : Mixed version using the binomal formula (a-b)^2 = a^2 - 2ab + b^2
+ * for | II - Beff |^2 and discretize every term individually.
+ * *But with |II|^2 instead of |D^2y|^2
+ */
+ // Gauss-Quadrature:
+ const auto &quadRule = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
+ // Trapezoidal-rule:
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ // Trapezoidal-rule + evaluation on edge centers (only for testing purposes):
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{0.5,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,0.5}, 1.0/6.0}, {{0.5,0.5}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ for (auto&& quadPoint : quadRule)
+ {
+
+ //Get values of the P2-Basis functions on current quadrature point
+ std::vector<FieldVector<double,1>> basisValues;
+ lagrangeLFE_.localBasis().evaluateFunction(quadPoint.position(), basisValues);
+
+
+ // Get Jacobians of the P2-Basis functions on current quadrature point
+ std::vector<FieldMatrix<double, 1, gridDim>> referenceGradients;
+
+ lagrangeLFE_.localBasis().evaluateJacobian(quadPoint.position(), referenceGradients);
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+ const auto integrationElement = geometry.integrationElement(quadPoint.position());
+
+ std::vector<FieldVector<RT, gridDim>> gradients(referenceGradients.size());
+ for (size_t i = 0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+ Tensor3<RT,3,2,2> discreteHessian(0);
+ FieldMatrix<RT, 3, gridDim> discreteGradient(0);
+
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ {
+ discreteGradient[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues[i];
+ discreteHessian[k][l][0] += discreteJacobianCoefficients[i][k][l]*gradients[i][0];
+ discreteHessian[k][l][1] += discreteJacobianCoefficients[i][k][l]*gradients[i][1];
+
+ }
+
+
+ //compute normal vector (cross product of partial derivatives)
+ FieldVector<RT,3> surfaceNormal(0);
+
+ FieldVector<RT,3> partialX = {discreteGradient[0][0], discreteGradient[1][0], discreteGradient[2][0]};
+ FieldVector<RT,3> partialY = {discreteGradient[0][1], discreteGradient[1][1], discreteGradient[2][1]};
+
+ auto normX = partialX.two_norm();
+ auto normY = partialY.two_norm();
+
+ // Normalize
+ // partialX /= normX;
+ // partialY /= normY;
+
+
+ surfaceNormal[0] = partialX[1]*partialY[2] - partialX[2]*partialY[1];
+ surfaceNormal[1] = partialX[2]*partialY[0] - partialX[0]*partialY[2];
+ surfaceNormal[2] = partialX[0]*partialY[1] - partialX[1]*partialY[0];
+ // printvector(std::cout, surfaceNormal, "surfaceNormal", "--");
+ // printvector(std::cout, surfaceNormal.two_norm(), "surfaceNormal.two_norm()", "--");
+
+ // RT surfaceNormalMag = sqrt(surfaceNormal[0]*surfaceNormal[0] + surfaceNormal[1]*surfaceNormal[1] + surfaceNormal[2]*surfaceNormal[2]);
+ // std::cout << "surfaceNormalMag: " << surfaceNormalMag << std::endl;
+
+
+ FieldMatrix<RT,2,2> secondFF(0);
+
+ for(size_t m=0; m<2; m++)
+ for(size_t n=0; n<2; n++)
+ for(size_t k=0; k<3; k++)
+ {
+ secondFF[m][n] += surfaceNormal[k]*discreteHessian[k][n][m];
+ //Test: transposed version
+ // secondFF[m][n] += discreteHessian[k][m][n]*surfaceNormal[k];
+ // secondFF[0][0] += surfaceNormal[k]*discreteHessian[k][0][0];
+ // secondFF[0][1] += surfaceNormal[k]*discreteHessian[k][1][0];
+ // secondFF[1][0] += surfaceNormal[k]*discreteHessian[k][0][1];
+ // secondFF[1][1] += surfaceNormal[k]*discreteHessian[k][1][1];
+ }
+
+ // substract effective prestrain
+ // G = G + Beff;
+
+ //Take symmetric part?
+ // auto symG = 0.5*(G.transposed() + G);
+
+
+
+
+ // printmatrix(std::cout, G , "matrix G: ", "--");
+ // printmatrix(std::cout, symG , "matrix symG: ", "--");
+
+ // std::cout << " discreteHessian.frobenius_norm2(): " << discreteHessian.frobenius_norm2() << std::endl;
+ // std::cout << " G.frobenius_norm2() : " << G.frobenius_norm2() << std::endl;
+ // std::cout << " DIFFERENCE : " << G.frobenius_norm2()- discreteHessian.frobenius_norm2() << std::endl;
+
+ /**
+ * @brief Compute squared Frobenius norm of the discrete Hessian
+ */
+ // auto quadResult = discreteHessian.frobenius_norm2();
+ // std::cout << "quadResult: " << weight*quadResult << std::endl;
+ // auto weight = geometry.volume()/3.0 * element.geometry().integrationElement(geometry.local(geometry.corner(c)));
+
+ // auto weight = geometry.volume()/3.0;
+
+ // std::cout << "geometry.volume()/3.0 : " << geometry.volume()/3.0 << std::endl;
+
+ // we need to use the surface area of the reference element (which is 1/2)
+ // std::cout << "geometry.volume(): " << geometry.volume() << std::endl;
+ // auto weight = (1.0/6.0) * element.geometry().integrationElement(geometry.local(geometry.corner(c)));
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ // std::cout << "weight: " << weight << std::endl;
+
+ // harmonicEnergy += weight * discreteHessian.frobenius_norm2();
+ // harmonicEnergy += weight * symG.frobenius_norm2();
+
+ // auto term1 = 0.5 * weight * discreteHessian.frobenius_norm2();
+ auto term1 = 0.5 * weight * secondFF.frobenius_norm2();
+
+ RT term2 = 0.0;
+
+ for(size_t i=0; i<2; i++)
+ for(size_t j=0; j<2; j++)
+ {
+ term2 += weight * secondFF[i][j] * Beff[i][j];
+ }
+
+ auto term3 = 0.5 * weight* Beff.frobenius_norm2();
+
+ // std::cout << "term1: " << term1 << std::endl;
+ // std::cout << "term2: " << term2 << std::endl;
+ // std::cout << "term3: " << term3 << std::endl;
+ // std::cout << "term1 - term2 + term3: " << term1 - term2 + term3 << std::endl;
+
+ // harmonicEnergy += weight * G.frobenius_norm2();
+ harmonicEnergy += term1 - term2 + term3; //eigentlich muss man term2 abziehen bei | II - Z |
+
+ // std::cout << "harmonicEnergy: " << harmonicEnergy << std::endl;
+ }
+ #endif
+
+
+ /**
+ * @brief Compute contribution of the force Term
+ *
+ * Integrate the scalar product of the force 'f'
+ * with the local deformation function 'localFunction_'
+ * over the current element.
+ */
+ RT forceTermEnergy = 0;
+
+ if (parameterSet_.get<bool>("assemble_force_term", 1))
+ {
+ for (auto&& quadPoint : quadRule)
+ {
+ auto deformationValue = localFunction_.evaluate(quadPoint.position());
+ auto forceValue = localForce_(quadPoint.position());
+ // printvector(std::cout, deformationValue.globalCoordinates(), "deformationValue", "--");
+ // printvector(std::cout, forceValue, "forceValue", "--");
+ // const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
+ // std::cout << "forceTermEnergy:" << forceTermEnergy << std::endl;
+ }
+ }
+
+ return harmonicEnergy - forceTermEnergy;
+ }
+
+ protected:
+ LocalDiscreteKirchhoffFunction& localFunction_;
+
+ LocalForce& localForce_;
+
+ P2LocalFiniteElement lagrangeLFE_;
+
+ const Dune::ParameterTree parameterSet_;
+
+ // PrestrainType Beff_;
+ // Dune::FieldMatrix<RT,3,3> Qhom_;
+ MatrixType Qhom_;
+ PrestrainType Beff_;
+
+ Python::Module module_;
+
+ // const PrestrainType& prestrain_;
+ // const MatrixType& quadraticForm_;
+ };
+
+
+
+} // namespace Dune::GFE
+#endif
\ No newline at end of file
diff --git a/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained_debug.hh b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained_debug.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c6c5091307d64b0c7902b242a0732bdae9dfd733
--- /dev/null
+++ b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained_debug.hh
@@ -0,0 +1,950 @@
+#ifndef DUNE_GFE_ENERGIES_DISCRETEKIRCHHOFFBENDINGENERGYPRESTRAINED_HH
+#define DUNE_GFE_ENERGIES_DISCRETEKIRCHHOFFBENDINGENERGYPRESTRAINED_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>
+
+#include <cmath>
+/**
+ * \brief Assemble the discrete Kirchhoff bending energy for a single element.
+ *
+ * The Kirchhoff bending energy consists of two parts:
+ *
+ * 1. The energy of Qhom(II - Beff) where II is the second fundamental form of the surface parametrized by the deformation function
+ * and Beff is a (given) effective prestrain tensor.
+ *
+ * This contribution is split intro three parts which corresponds to the binomial formula (a+b)^2 = a^2 + 2ab + b^2
+ * each term is discretized separately.
+ *
+ * 2. An integral over the scalar product of a forcing term and the discrete deformation
+ * (i.e. the evaluation of localFunction_ ).
+ */
+namespace Dune::GFE
+{
+ template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
+ class DiscreteKirchhoffBendingEnergyPrestrained
+ : 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 LocalDiscreteKirchhoffFunction::DerivativeType DerivativeType;
+ typedef typename Dune::FieldMatrix<RT,2,2> PrestrainType;
+ typedef typename Dune::FieldMatrix<RT,3,3> MatrixType;
+ // 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:
+ DiscreteKirchhoffBendingEnergyPrestrained(LocalDiscreteKirchhoffFunction &localFunction,
+ LocalForce &localForce,
+ const Dune::ParameterTree& parameterSet)
+ : localFunction_(localFunction),
+ localForce_(localForce),
+ parameterSet_(parameterSet)
+ {}
+
+
+ /** \brief Assemble the energy for a single element */
+ virtual RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localSolution) const
+ {
+ RT energy = 0;
+
+ /**
+ * @brief Get Prestrain Tensor
+ */
+ // Dune::FieldMatrix<RT,2,2> Beff = prestrain_;
+
+ /**
+ * @brief Get effective quadratic form
+ */
+ // MatrixType Qhom = quadraticForm_;
+
+
+ /**
+ * @brief Get effective prestrain Tensor
+ */
+ Dune::FieldVector<RT,3> Beffvec = parameterSet_.get<Dune::FieldVector<RT,3>>("effectivePrestrain", {0.0, 0.0, 0.0});
+ Dune::FieldMatrix<RT,2,2> Beff = {{Beffvec[0], Beffvec[2]}, {Beffvec[2], Beffvec[1]} };
+ // printmatrix(std::cout, Beff, "effective prestrain (Beff): ", "--");
+
+
+ /**
+ * @brief Get effective quadratic form Qhom
+ *
+ * input-vector: [q_1,q_2,q_3,q_12,q_13,q_23]
+ * is assembled into a matrix where the off-diagonal entries are divided by 2 (compare with definition in the paper)
+ * ( q_1 , 0.5*q_12 , 0.5*q_13 )
+ * Q = ( 0.5*q_12 , q_2 , 0.5*q_23 )
+ * ( 0.5*q_13 , 0.5*q_23 , q_3 )
+ */
+ Dune::FieldVector<RT,6> Qhomvec = parameterSet_.get<Dune::FieldVector<RT,6>>("effectiveQuadraticForm", {1.0, 1.0, 1.0, 0.0, 0.0, 0.0});
+ MatrixType Qhom = {{Qhomvec[0] , 0.5*Qhomvec[3], 0.5*Qhomvec[4]},
+ {0.5*Qhomvec[3], Qhomvec[1], 0.5*Qhomvec[5]},
+ {0.5*Qhomvec[4], 0.5*Qhomvec[5], Qhomvec[2]}};
+ // printmatrix(std::cout, Qhom, "effective quadratic form (Qhom): ", "--");
+
+
+
+
+ bool PRINT_DEBUG = 0;
+
+
+#if 1
+ // TODO: More serious order selection
+ int quadOrder = 3;
+ // int quadOrder = 6;
+#else
+ if (localFiniteElement.localBasis().order()==1)
+ quadOrder = 2;
+ else
+ quadOrder = (localFiniteElement.type().isSimplex()) ? (localFiniteElement.localBasis().order() - 1) * 2
+ : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
+#endif
+
+ const auto element = localView.element();
+ auto geometry = element.geometry();
+ auto gridView = localFunction_.gridView();
+ const auto &indexSet = gridView.indexSet();
+
+ // bind to current element.
+ localFunction_.bind(element);
+ localForce_.bind(element);
+ /**
+ * @brief We need to update the Coefficients of localFunction_
+ * on the current element.
+ */
+ localFunction_.updateLocalCoefficients(localSolution,element);
+
+
+
+ // std::cout << " -------------------------------------------- " << std::endl;
+ // std::cout <<"ELEMENT NUMBER( indexSet.index(element)) ):" << indexSet.index(element) << std::endl;
+ if(PRINT_DEBUG)
+ {
+ std::cout << "geometry.corner(0): " << geometry.corner(0) << std::endl;
+ std::cout << "geometry.corner(1)): " << geometry.corner(1) << std::endl;
+ std::cout << "geometry.corner(2): " << geometry.corner(2) << std::endl;
+ }
+
+
+ /**
+ * @brief Get Coefficients of the discrete Jacobian
+ *
+ * The discrete Jacobian 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 Jacobians of localfunction_ .
+ * The coefficients are stored in the form [Basisfunctions x components x gridDim]
+ * in a BlockVector<FieldMatrix<RT, 3, gridDim>> .
+ */
+ BlockVector<FieldMatrix<RT, 3, gridDim>> discreteJacobianCoefficients;
+ discreteGradientCoefficients(discreteJacobianCoefficients,lagrangeLFE_,localFunction_,element);
+
+
+
+
+
+ /**
+ * @brief Compute harmonic energy contribution:
+ */
+ RT harmonicEnergy = 0;
+
+
+
+ //VERSION - 1
+
+ #if 0
+ /**
+ * @brief Setup Quadrature rule.
+ *
+ */
+ // Gauss-Quadrature:
+ const auto &quadRule = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
+ // Trapezoidal-rule:
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ // Trapezoidal-rule + evaluation on edge centers (only for testing purposes):
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{0.5,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,0.5}, 1.0/6.0}, {{0.5,0.5}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ for (auto&& quadPoint : quadRule)
+ {
+ //Get values of the P2-Basis functions on current quadrature point
+ std::vector<FieldVector<double,1>> basisValues;
+ lagrangeLFE_.localBasis().evaluateFunction(quadPoint.position(), basisValues);
+
+ // if(PRINT_DEBUG)
+ // {
+ // std::cout << "- * - * - * - * - * - *" << std::endl;
+ // std::cout << "quadPoint.position():" << quadPoint.position() << std::endl;
+ // std::cout << "geometry.global(quadPoint.position()): " << geometry.global(quadPoint.position()) << std::endl;
+ // }
+
+ // Get Jacobians of the P2-Basis functions on current quadrature point
+ std::vector<FieldMatrix<double, 1, gridDim>> referenceGradients;
+
+ lagrangeLFE_.localBasis().evaluateJacobian(quadPoint.position(), referenceGradients);
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+ const auto integrationElement = geometry.integrationElement(quadPoint.position());
+ // printmatrix(std::cout, jacobian , "jacobian : ", "--");
+
+ std::vector<FieldVector<RT, gridDim>> gradients(referenceGradients.size());
+ for (size_t i = 0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+ //TEST: Do not transform P2 basis
+ // std::vector<FieldMatrix<double, 1, gridDim>> gradients;
+ // lagrangeLFE_.localBasis().evaluateJacobian(quadPoint.position(), gradients);
+
+
+ Tensor3<RT,3,2,2> discreteHessian(0);
+ FieldMatrix<RT, 3, gridDim> discreteGradient(0);
+
+
+ /**
+ * @brief Compute the discrete Hessian and discrete Gradient as a linear combination of
+ * function and gradient evaluations of the P2-basis functions with the coefficients given
+ * by 'discreteJacobianCoefficients'.
+ * Since the deformation function has k=3 components we get
+ * - discreteGradient : 3x2 matrix
+ * - discreteHessian: 3x2x2 tensor
+ */
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ {
+ // Compute discrete gradient
+ discreteGradient[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues[i];
+ // Compute discrete Hessian
+ discreteHessian[k][l][0] += discreteJacobianCoefficients[i][k][l]*gradients[i][0];
+ discreteHessian[k][l][1] += discreteJacobianCoefficients[i][k][l]*gradients[i][1];
+ // discreteHessian[k][l][0] += discreteJacobianCoefficients[i][k][l]*gradients[i][0][0];
+ // discreteHessian[k][l][1] += discreteJacobianCoefficients[i][k][l]*gradients[i][0][1];
+ }
+
+
+ // DerivativeType whJacobianValue = localFunction_.evaluateDerivative(quadPoint.position());
+
+ // Check difference in discrete gradients:
+ // printmatrix(std::cout, localFunction_.evaluateDerivative(quadPoint.position()) , "localFunction_.evaluateDerivative(quadPoint.position()): ", "--");
+ // printmatrix(std::cout, discreteGradient, "discreteGradient: ", "--");
+
+ // printmatrix(std::cout, discreteGradient-localFunction_.evaluateDerivative(quadPoint.position()) , "gradient difference: ", "--");
+
+
+ //
+
+
+
+
+ /**
+ * @brief Check if isometry constraint is satisfied (at the current quadrature point).
+ * This should at least be satisfied on the nodes of the triangulation.
+ */
+ // FieldMatrix<RT, 2,2> isometryConstraint(0);
+ // for(size_t k=0; k<3; k++)
+ // for(size_t l=0; l<2; l++)
+ // for(size_t i=0; i<2; i++)
+ // {
+ // isometryConstraint[l][i] += discreteGradient[k][l]*discreteGradient[k][i];
+ // }
+ // printmatrix(std::cout, isometryConstraint, "isometryConstraint: ", "--");
+
+
+
+
+
+ // if(PRINT_DEBUG)
+ // {
+ /**
+ * @brief print the (two) different 3x2 slices of the discrete Hessian
+ */
+ FieldMatrix<RT, 3, 2> discreteHessian_slice1(0);
+ FieldMatrix<RT, 3, 2> discreteHessian_slice2(0);
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ {
+ discreteHessian_slice1[k][l] = discreteHessian[k][l][0];
+ discreteHessian_slice2[k][l] = discreteHessian[k][l][1];
+ }
+ // printmatrix(std::cout, discreteHessian_slice1, "discreteHessian_slice1: ", "--");
+ // printmatrix(std::cout, discreteHessian_slice2, "discreteHessian_slice2: ", "--");
+ // }
+
+
+ /**
+ * @brief Test: evaluate the deformation function and gradient at the current quadrature point.
+ */
+ // auto test1= localFunction_.evaluateDerivative(quadPoint.position());
+ // auto test2= localFunction_.evaluate(quadPoint.position()).globalCoordinates();
+
+
+
+ /**
+ * @brief Compute the surface normal given by the cross-product of the two different partial derivatives
+ * of the discrete gradient.
+ * This is needed to compute the second fundamental form.
+ */
+ FieldVector<RT,3> surfaceNormal(0);
+ FieldVector<RT,3> partialX = {discreteGradient[0][0], discreteGradient[1][0], discreteGradient[2][0]};
+ FieldVector<RT,3> partialY = {discreteGradient[0][1], discreteGradient[1][1], discreteGradient[2][1]};
+
+ auto normX = partialX.two_norm();
+ auto normY = partialY.two_norm();
+
+ /**
+ * @brief Test: normalize
+ */
+ // partialX /= normX;
+ // partialY /= normY;
+
+ surfaceNormal[0] = partialX[1]*partialY[2] - partialX[2]*partialY[1];
+ surfaceNormal[1] = partialX[2]*partialY[0] - partialX[0]*partialY[2];
+ surfaceNormal[2] = partialX[0]*partialY[1] - partialX[1]*partialY[0];
+ // printvector(std::cout, surfaceNormal, "surfaceNormal", "--");
+ // std::cout << "surfaceNormal.two_norm() : " << surfaceNormal.two_norm() << std::endl;
+
+ // surfaceNormal *= -1.0;
+
+ /**
+ * @brief Test: normalize
+ */
+ // auto norm = surfaceNormal.two_norm();
+ // surfaceNormal /= norm;
+ // printvector(std::cout, surfaceNormal, "surfaceNormal", "--");
+
+
+ /**
+ * @brief Compute the second fundamental form. This is a 2x2 matrix.
+ */
+ FieldMatrix<RT,2,2> secondFF(0);
+
+ for(size_t m=0; m<2; m++)
+ for(size_t n=0; n<2; n++)
+ for(size_t k=0; k<3; k++)
+ {
+ secondFF[m][n] += surfaceNormal[k]*discreteHessian[k][n][m];
+
+
+ //Test: transposed version
+ // secondFF[m][n] += discreteHessian[k][m][n]*surfaceNormal[k];
+ // secondFF[0][0] += surfaceNormal[k]*discreteHessian[k][0][0];
+ // secondFF[0][1] += surfaceNormal[k]*discreteHessian[k][1][0];
+ // secondFF[1][0] += surfaceNormal[k]*discreteHessian[k][0][1];
+ // secondFF[1][1] += surfaceNormal[k]*discreteHessian[k][1][1];
+ }
+ // printmatrix(std::cout, secondFF , "secondFF: ", "--");
+
+
+ /**
+ * @brief Check orthogonality between discrete Hessian and discrete gradient
+ * at nodes.
+ */
+ // std::cout << "---- ORTHOGONALITY CHECK 1: ----" << std::endl;
+ // for(int i=0; i<2; i++)
+ // for(int j=0; j<2; j++)
+ // for(int l=0; l<2; l++)
+ // {
+ // RT test = 0.0;
+ // for(int k=0; k<3; k++)
+ // {
+ // test += discreteHessian[k][j][i] * discreteGradient[k][l];
+ // // test += discreteHessian[k][i][j] * discreteGradient[k][l];
+ // }
+ // std::cout << "Partial_" << i << j << "U * Partial_" << l << "U = " << test << std::endl;
+ // }
+
+
+ /**
+ * @brief Print the (discrete) second partial derivatives.
+ */
+ // for(int i=0; i<2; i++)
+ // for(int j=0; j<2; j++)
+ // {
+ // FieldVector<RT,3> partialTMP(0);
+ // for(int k=0; k<3; k++)
+ // {
+ // partialTMP[k] = discreteHessian[k][j][i];
+ // }
+ // std::cout << "Partial_" << i << j << "U = " << std::endl;
+
+ // printvector(std::cout , partialTMP, "PartialTMP:", "--");
+ // }
+
+ /**
+ * @brief Check if Norm-identity holds
+ */
+ // std::cout << "---- ORTHOGONALITY CHECK 2: ----" << std::endl;
+ // for(int i=0; i<2; i++)
+ // for(int j=0; j<2; j++)
+ // {
+ // RT secondPartialDerNorm = 0.0;
+ // RT secondPartialDerNorm_normalPart = 0.0;
+
+ // for(int k=0; k<3; k++)
+ // {
+ // secondPartialDerNorm += discreteHessian[k][j][i] * discreteHessian[k][j][i];
+ // secondPartialDerNorm_normalPart += (discreteHessian[k][j][i] * surfaceNormal[k]) * (discreteHessian[k][j][i] * surfaceNormal[k]);
+ // // secondPartialDerNorm += discreteHessian[k][i][j] * discreteHessian[k][i][j];
+ // // secondPartialDerNorm_normalPart += (discreteHessian[k][i][j] * surfaceNormal[k]) * (discreteHessian[k][i][j] * surfaceNormal[k]);
+ // }
+ // std::cout << "secondPartialDerNorm (" << i << j << ") = " << secondPartialDerNorm << std::endl;
+ // std::cout << "secondPartialDerNorm_normalPart (" << i << j << ") = " << secondPartialDerNorm_normalPart << std::endl;
+ // }
+
+
+ /**
+ * @brief Subtract the effective prestrain
+ */
+ auto G = secondFF - Beff;
+
+
+ // auto symG = 0.5*(G.transposed() + G);
+ // printmatrix(std::cout, G , "matrix G: ", "--");
+ // printmatrix(std::cout, symG , "matrix symG: ", "--");
+
+
+
+ // auto weight = (1.0/6.0) * element.geometry().integrationElement(quadPoint.position());
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+
+ /**
+ * @brief TEST: Check the difference between the squared Frobenius norm of the
+ * second fundamental form and discrete Hessian.
+ * (Note that for isometries, these should be equal.)
+ */
+ // auto DIFFERENCE = secondFF.frobenius_norm2() - discreteHessian.frobenius_norm2();
+ // if (abs(DIFFERENCE)>1e-1)
+ // {
+ // std::cout << " discreteHessian.frobenius_norm2(): " << discreteHessian.frobenius_norm2() << std::endl;
+ // std::cout << " secondFF.frobenius_norm2() : " << secondFF.frobenius_norm2() << std::endl;
+ // std::cout << " DIFFERENCE : " << secondFF.frobenius_norm2() - discreteHessian.frobenius_norm2() << std::endl;
+ // printmatrix(std::cout, secondFF , "secondFF: ", "--");
+ // std::cout << "G.frobenius_norm2(): " << G.frobenius_norm2() << std::endl;
+ // printmatrix(std::cout, isometryConstraint, "isometryConstraint: ", "--");
+
+ // printvector(std::cout, surfaceNormal, "surfaceNormal", "--");
+ // printmatrix(std::cout, discreteHessian_slice1, "discreteHessian_slice1: ", "--");
+ // printmatrix(std::cout, discreteHessian_slice2, "discreteHessian_slice2: ", "--");
+
+ // //compare with discretization that used binomial formula (a-b)^2 = a^2 - 2ab + b^2
+ // auto term1 = 0.5 * weight * discreteHessian.frobenius_norm2();
+ // RT term2 = 0.0;
+ // for(size_t i=0; i<2; i++)
+ // for(size_t j=0; j<2; j++)
+ // {
+ // term2 += weight * secondFF[i][j] * Beff[i][j];
+ // }
+ // auto term3 = 0.5 * weight * Beff.frobenius_norm2();
+ // std::cout << "term1-term2+term3 :" << term1-term2+term3 << std::endl;
+ // }
+
+ /**
+ * @brief Compute harmonic energy contribution of the current quadrature point.
+ */
+ harmonicEnergy += 0.5 * weight * G.frobenius_norm2();
+ // std::cout << "harmonicEnergy: " << harmonicEnergy << std::endl;
+ }
+ #endif
+
+
+ /**
+ * @brief Debugging Tests
+ */
+ // std::cout << "Auswertung auf Kantenmitte: " << std::endl;
+ // auto edge1Geometry = element.template subEntity<1>(0).geometry();
+ // auto edge1Entity = element.template subEntity<1>(0);
+ // auto edge2Geometry = element.template subEntity<1>(1).geometry();
+ // auto edge2Entity = element.template subEntity<1>(1);
+ // auto edge3Geometry = element.template subEntity<1>(2).geometry();
+ // auto edge3Entity = element.template subEntity<1>(2);
+
+ // auto edge1Center = edge1Geometry.center();
+ // auto edge2Center = edge2Geometry.center();
+ // auto edge3Center = edge3Geometry.center();
+
+ // std::cout << "edge1Center: " << edge1Center << std::endl;
+ // std::cout << "edge2Center: " << edge2Center << std::endl;
+ // std::cout << "edge3Center: " << edge3Center << std::endl;
+
+ // std::cout << "geometry.local(edge1Center):" << geometry.local(edge1Center) << std::endl;
+ // std::cout << "geometry.local(edge2Center):" << geometry.local(edge2Center) << std::endl;
+ // std::cout << "geometry.local(edge3Center):" << geometry.local(edge3Center) << std::endl;
+ // // std::cout << "geometry.global({0.5,0.0}): "<< geometry.global({0.5,0.0}) << std::endl;
+ // // std::cout << "geometry.global({0.5,0.5}): "<< geometry.global({0.5,0.5}) << std::endl;
+ // // std::cout << "geometry.global({0.0,0.5}): "<< geometry.global({0.0,0.5}) << std::endl;
+ // std::vector<FieldVector<double,1>> basisValues_edge1;
+ // std::vector<FieldVector<double,1>> basisValues_edge2;
+ // std::vector<FieldVector<double,1>> basisValues_edge3;
+ // // lagrangeLFE_.localBasis().evaluateFunction({0.5,0.0}, basisValues_edge1);
+ // // lagrangeLFE_.localBasis().evaluateFunction({0.5,0.5}, basisValues_edge2);
+ // // lagrangeLFE_.localBasis().evaluateFunction({0.0,0.5}, basisValues_edge3);
+ // lagrangeLFE_.localBasis().evaluateFunction(geometry.local(edge1Center), basisValues_edge1);
+ // lagrangeLFE_.localBasis().evaluateFunction(geometry.local(edge2Center), basisValues_edge2);
+ // lagrangeLFE_.localBasis().evaluateFunction(geometry.local(edge3Center), basisValues_edge3);
+ // //TEST VERTICES
+ // // lagrangeLFE_.localBasis().evaluateFunction({0.0,0.0}, basisValues_edge1);
+ // // lagrangeLFE_.localBasis().evaluateFunction({1.0,0.0}, basisValues_edge2);
+ // // lagrangeLFE_.localBasis().evaluateFunction({0.0,1.0}, basisValues_edge3);
+
+ // std::cout << "------ print P2-Basis evaluation -----" << std::endl;
+ // for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ // {
+ // std::cout << i << "-th P2-basis function: " << std::endl;
+ // printvector(std::cout, basisValues_edge1[i] , "basisValues_edge1[i] ", "--");
+ // }
+ // for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ // {
+ // std::cout << i << "-th P2-basis function: " << std::endl;
+ // printvector(std::cout, basisValues_edge2[i] , "basisValues_edge2[i] ", "--");
+ // }
+ // for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ // {
+ // std::cout << i << "-th P2-basis function: " << std::endl;
+ // printvector(std::cout, basisValues_edge3[i] , "basisValues_edge3[i] ", "--");
+ // }
+
+ // FieldMatrix<RT, 3, gridDim> discreteGradient_edge1(0);
+ // FieldMatrix<RT, 3, gridDim> discreteGradient_edge2(0);
+ // FieldMatrix<RT, 3, gridDim> discreteGradient_edge3(0);
+
+ // for (int k=0; k<3; k++)
+ // for (int l=0; l<gridDim; l++)
+ // for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ // {
+ // discreteGradient_edge1[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues_edge1[i];
+ // discreteGradient_edge2[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues_edge2[i];
+ // discreteGradient_edge3[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues_edge3[i];
+ // }
+
+ // printmatrix(std::cout, discreteGradient_edge1, "discreteGradient_edge1: ", "--");
+ // printmatrix(std::cout, discreteGradient_edge2, "discreteGradient_edge2: ", "--");
+ // printmatrix(std::cout, discreteGradient_edge3, "discreteGradient_edge3: ", "--");
+ // std::cout << "sin(edge1Center[0]): " << sin(edge1Center[0]) << std::endl;
+ // std::cout << "cos(edge1Center[0]): " << cos(edge1Center[0]) << std::endl;
+ // std::cout << "sin(edge2Center[0]): " << sin(edge2Center[0]) << std::endl;
+ // std::cout << "cos(edge2Center[0]): " << cos(edge2Center[0]) << std::endl;
+ // std::cout << "sin(edge3Center[0]): " << sin(edge3Center[0]) << std::endl;
+ // std::cout << "cos(edge3Center[0]): " << cos(edge3Center[0]) << std::endl;
+ // // std::cout << "sin(geometry.global({0.5,0.0})[0]): " << sin(geometry.global({0.5,0.0})[0]) << std::endl;
+ // // std::cout << "cos(geometry.global({0.5,0.0})[0]): " << cos(geometry.global({0.5,0.0})[0]) << std::endl;
+ // // std::cout << "sin(geometry.global({0.5,0.5})[0]): " << sin(geometry.global({0.5,0.5})[0]) << std::endl;
+ // // std::cout << "cos(geometry.global({0.5,0.5})[0]): " << cos(geometry.global({0.5,0.5})[0]) << std::endl;
+ // // std::cout << "sin(geometry.global({0.0,0.5})[0]): " << sin(geometry.global({0.0,0.5})[0]) << std::endl;
+ // // std::cout << "cos(geometry.global({0.0,0.5})[0]): " << cos(geometry.global({0.0,0.5})[0]) << std::endl;
+ // // TEST VERTICES
+ // // std::cout << "sin(geometry.global({0.0,0.0})[0]): " << sin(geometry.global({0.0,0.0})[0]) << std::endl;
+ // // std::cout << "cos(geometry.global({0.0,0.0})[0]): " << cos(geometry.global({0.0,0.0})[0]) << std::endl;
+ // // std::cout << "sin(geometry.global({1.0,0.0})[0]): " << sin(geometry.global({1.0,0.0})[0]) << std::endl;
+ // // std::cout << "cos(geometry.global({1.0,0.0})[0]): " << cos(geometry.global({1.0,0.0})[0]) << std::endl;
+ // // std::cout << "sin(geometry.global({0.0,1.0})[0]): " << sin(geometry.global({0.0,1.0})[0]) << std::endl;
+ // // std::cout << "cos(geometry.global({0.0,1.0})[0]): " << cos(geometry.global({0.0,1.0})[0]) << std::endl;
+ // exit(0);
+ // if(indexSet.index(element)==1)
+ // exit(0);
+
+
+ // #if 0
+ /**
+ * @brief VERSION - 2 : Mixed version using the binomal formula (a-b)^2 = a^2 - 2ab + b^2
+ * for | II - Beff |^2 and discretize every term individually.
+ */
+ // Gauss-Quadrature:
+ const auto &quadRule = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
+ // Trapezoidal-rule:
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ // Trapezoidal-rule + evaluation on edge centers (only for testing purposes):
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{0.5,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,0.5}, 1.0/6.0}, {{0.5,0.5}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ for (auto&& quadPoint : quadRule)
+ {
+
+ //Get values of the P2-Basis functions on current quadrature point
+ std::vector<FieldVector<double,1>> basisValues;
+ lagrangeLFE_.localBasis().evaluateFunction(quadPoint.position(), basisValues);
+
+
+ // Get Jacobians of the P2-Basis functions on current quadrature point
+ std::vector<FieldMatrix<double, 1, gridDim>> referenceGradients;
+
+ lagrangeLFE_.localBasis().evaluateJacobian(quadPoint.position(), referenceGradients);
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+ const auto integrationElement = geometry.integrationElement(quadPoint.position());
+
+ std::vector<FieldVector<RT, gridDim>> gradients(referenceGradients.size());
+ for (size_t i = 0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+ Tensor3<RT,3,2,2> discreteHessian(0);
+ FieldMatrix<RT, 3, gridDim> discreteGradient(0);
+
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ {
+ discreteGradient[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues[i];
+ discreteHessian[k][l][0] += discreteJacobianCoefficients[i][k][l]*gradients[i][0];
+ discreteHessian[k][l][1] += discreteJacobianCoefficients[i][k][l]*gradients[i][1];
+
+ }
+
+ /**
+ * @brief print the (three) 2x2 slices of the discrete Hessian (check for symmetry)
+ */
+ FieldMatrix<RT, 2, 2> discreteHessian_slice1(0);
+ FieldMatrix<RT, 2, 2> discreteHessian_slice2(0);
+ FieldMatrix<RT, 2, 2> discreteHessian_slice3(0);
+
+ for (int i=0; i<2; i++)
+ for (int l=0; l<gridDim; l++)
+ {
+ discreteHessian_slice1[i][l] = discreteHessian[0][i][l];
+ discreteHessian_slice2[i][l] = discreteHessian[1][i][l];
+ discreteHessian_slice3[i][l] = discreteHessian[2][i][l];
+ }
+
+ // printmatrix(std::cout, discreteHessian_slice1, "discreteHessian_slice1: ", "--");
+ // printmatrix(std::cout, discreteHessian_slice2, "discreteHessian_slice2: ", "--");
+ // printmatrix(std::cout, discreteHessian_slice3, "discreteHessian_slice3: ", "--");
+ // printmatrix(std::cout, sym(discreteHessian_slice1), "SYM discreteHessian_slice1: ", "--");
+ // printmatrix(std::cout, sym(discreteHessian_slice2), "SYM discreteHessian_slice2: ", "--");
+ // printmatrix(std::cout, sym(discreteHessian_slice3), "SYM discreteHessian_slice3: ", "--");
+
+
+ //compute normal vector (cross product of partial derivatives)
+ FieldVector<RT,3> surfaceNormal(0);
+
+ FieldVector<RT,3> partialX = {discreteGradient[0][0], discreteGradient[1][0], discreteGradient[2][0]};
+ FieldVector<RT,3> partialY = {discreteGradient[0][1], discreteGradient[1][1], discreteGradient[2][1]};
+
+ auto normX = partialX.two_norm();
+ auto normY = partialY.two_norm();
+
+ // Normalize
+ // partialX /= normX;
+ // partialY /= normY;
+
+
+ surfaceNormal[0] = partialX[1]*partialY[2] - partialX[2]*partialY[1];
+ surfaceNormal[1] = partialX[2]*partialY[0] - partialX[0]*partialY[2];
+ surfaceNormal[2] = partialX[0]*partialY[1] - partialX[1]*partialY[0];
+ // printvector(std::cout, surfaceNormal, "surfaceNormal", "--");
+ // printvector(std::cout, surfaceNormal.two_norm(), "surfaceNormal.two_norm()", "--");
+
+ // RT surfaceNormalMag = sqrt(surfaceNormal[0]*surfaceNormal[0] + surfaceNormal[1]*surfaceNormal[1] + surfaceNormal[2]*surfaceNormal[2]);
+ // std::cout << "surfaceNormalMag: " << surfaceNormalMag << std::endl;
+
+ // FieldMatrix<RT,2,2> G(0);
+ FieldMatrix<RT,2,2> secondFF(0);
+
+ for(size_t m=0; m<2; m++)
+ for(size_t n=0; n<2; n++)
+ for(size_t k=0; k<3; k++)
+ {
+ secondFF[m][n] += surfaceNormal[k]*discreteHessian[k][n][m];
+ }
+
+ // To get the right sign:
+ secondFF = -1.0*secondFF;
+
+
+ // substract effective prestrain
+ // G = G + Beff;
+
+ //Take symmetric part?
+ // auto symG = 0.5*(G.transposed() + G);
+
+
+
+
+ // printmatrix(std::cout, G , "matrix G: ", "--");
+ // printmatrix(std::cout, symG , "matrix symG: ", "--");
+
+ // std::cout << " discreteHessian.frobenius_norm2(): " << discreteHessian.frobenius_norm2() << std::endl;
+ // std::cout << " G.frobenius_norm2() : " << G.frobenius_norm2() << std::endl;
+ // std::cout << " DIFFERENCE : " << G.frobenius_norm2()- discreteHessian.frobenius_norm2() << std::endl;
+
+ /**
+ * @brief Compute squared Frobenius norm of the discrete Hessian
+ */
+ // auto quadResult = discreteHessian.frobenius_norm2();
+ // std::cout << "quadResult: " << weight*quadResult << std::endl;
+ // auto weight = geometry.volume()/3.0 * element.geometry().integrationElement(geometry.local(geometry.corner(c)));
+
+ // auto weight = geometry.volume()/3.0;
+
+ // std::cout << "geometry.volume()/3.0 : " << geometry.volume()/3.0 << std::endl;
+
+ // we need to use the surface area of the reference element (which is 1/2)
+ // std::cout << "geometry.volume(): " << geometry.volume() << std::endl;
+ // auto weight = (1.0/6.0) * element.geometry().integrationElement(geometry.local(geometry.corner(c)));
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ // std::cout << "weight: " << weight << std::endl;
+
+ // harmonicEnergy += weight * discreteHessian.frobenius_norm2();
+ // harmonicEnergy += weight * symG.frobenius_norm2();
+
+ // auto term1 = 0.5 * weight * discreteHessian.frobenius_norm2();
+ RT term1 = 0.0;
+
+ FieldMatrix<RT,3,3> X = tensorToSymCoefficients(discreteHessian);
+
+ for(int k=0; k<3; k++)
+ {
+ // term1 += Qhom[0][0]*pow((2.0*discreteHessian[k][0][0]),2)+ Qhom[0][1]* ... more efficient (TODO) Qhom is symmetric
+ term1 += Qhom[0][0]*pow(X[k][0],2) + Qhom[0][1]*X[k][1]*X[k][0] + Qhom[0][2]*X[k][2]*X[k][0]
+ + Qhom[1][0]*X[k][0]*X[k][1] + Qhom[1][1]*pow(X[k][1],2) + Qhom[1][2]*X[k][2]*X[k][1]
+ + Qhom[2][0]*X[k][0]*X[k][2] + Qhom[2][1]*X[k][1]*X[k][2] + Qhom[2][2]*pow(X[k][2],2);
+ }
+ // term1 = term1 * 0.5 * weight;
+ term1 = term1 * weight;
+
+
+
+
+ //TODO : Add method "applyQhom"
+
+
+ // Second Term
+
+ // RT term2 = 0.0;
+ // for(size_t i=0; i<2; i++)
+ // for(size_t j=0; j<2; j++)
+ // {
+ // term2 += weight * secondFF[i][j] * Beff[i][j];
+ // }
+
+
+ FieldVector<RT,3> secondFFCoefficients = matrixToSymCoefficients(secondFF);
+ FieldVector<RT,3> BeffCoefficients = matrixToSymCoefficients(Beff);
+
+ RT term2 = 0.0;
+ term2 += Qhom[0][0]*BeffCoefficients[0]*secondFFCoefficients[0] + Qhom[0][1]*BeffCoefficients[1]*secondFFCoefficients[0] + Qhom[0][2]*BeffCoefficients[2]*secondFFCoefficients[0]
+ + Qhom[1][0]*BeffCoefficients[0]*secondFFCoefficients[1] + Qhom[1][1]*BeffCoefficients[1]*secondFFCoefficients[1] + Qhom[1][2]*BeffCoefficients[2]*secondFFCoefficients[1]
+ + Qhom[2][0]*BeffCoefficients[0]*secondFFCoefficients[2] + Qhom[2][1]*BeffCoefficients[1]*secondFFCoefficients[2] + Qhom[2][2]*BeffCoefficients[2]*secondFFCoefficients[2];
+ term2 = 2.0 * term2 * weight;
+
+
+
+
+
+ // auto term3 = 0.5 * weight* Beff.frobenius_norm2();
+
+ RT term3 = 0.0;
+
+ term3 += Qhom[0][0]*BeffCoefficients[0]*BeffCoefficients[0] + Qhom[0][1]*BeffCoefficients[1]*BeffCoefficients[0] + Qhom[0][2]*BeffCoefficients[2]*BeffCoefficients[0]
+ + Qhom[1][0]*BeffCoefficients[0]*BeffCoefficients[1] + Qhom[1][1]*BeffCoefficients[1]*BeffCoefficients[1] + Qhom[1][2]*BeffCoefficients[2]*BeffCoefficients[1]
+ + Qhom[2][0]*BeffCoefficients[0]*BeffCoefficients[2] + Qhom[2][1]*BeffCoefficients[1]*BeffCoefficients[2] + Qhom[2][2]*BeffCoefficients[2]*BeffCoefficients[2];
+
+ // term3 = term3 * 0.5 * weight ;
+ term3 = term3 * weight ;
+
+ // std::cout << "term1: " << term1 << std::endl;
+ // std::cout << "term2: " << term2 << std::endl;
+ // std::cout << "term3: " << term3 << std::endl;
+ // std::cout << "term1 - term2 + term3: " << term1 - term2 + term3 << std::endl;
+
+ // harmonicEnergy += weight * G.frobenius_norm2();
+ harmonicEnergy += term1 - term2 + term3; //eigentlich muss man term2 abziehen bei | II - Z |
+
+ // std::cout << "harmonicEnergy: " << harmonicEnergy << std::endl;
+ }
+ // #endif
+
+
+ #if 0
+ /**
+ * @brief VERSION - 3 : Mixed version using the binomal formula (a-b)^2 = a^2 - 2ab + b^2
+ * for | II - Beff |^2 and discretize every term individually.
+ * *But with |II|^2 instead of |D^2y|^2
+ */
+ // Gauss-Quadrature:
+ const auto &quadRule = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
+ // Trapezoidal-rule:
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ // Trapezoidal-rule + evaluation on edge centers (only for testing purposes):
+ // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{0.5,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,0.5}, 1.0/6.0}, {{0.5,0.5}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
+ for (auto&& quadPoint : quadRule)
+ {
+
+ //Get values of the P2-Basis functions on current quadrature point
+ std::vector<FieldVector<double,1>> basisValues;
+ lagrangeLFE_.localBasis().evaluateFunction(quadPoint.position(), basisValues);
+
+
+ // Get Jacobians of the P2-Basis functions on current quadrature point
+ std::vector<FieldMatrix<double, 1, gridDim>> referenceGradients;
+
+ lagrangeLFE_.localBasis().evaluateJacobian(quadPoint.position(), referenceGradients);
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+ const auto integrationElement = geometry.integrationElement(quadPoint.position());
+
+ std::vector<FieldVector<RT, gridDim>> gradients(referenceGradients.size());
+ for (size_t i = 0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+ Tensor3<RT,3,2,2> discreteHessian(0);
+ FieldMatrix<RT, 3, gridDim> discreteGradient(0);
+
+ for (int k=0; k<3; k++)
+ for (int l=0; l<gridDim; l++)
+ for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
+ {
+ discreteGradient[k][l] += discreteJacobianCoefficients[i][k][l]*basisValues[i];
+ discreteHessian[k][l][0] += discreteJacobianCoefficients[i][k][l]*gradients[i][0];
+ discreteHessian[k][l][1] += discreteJacobianCoefficients[i][k][l]*gradients[i][1];
+
+ }
+
+
+ //compute normal vector (cross product of partial derivatives)
+ FieldVector<RT,3> surfaceNormal(0);
+
+ FieldVector<RT,3> partialX = {discreteGradient[0][0], discreteGradient[1][0], discreteGradient[2][0]};
+ FieldVector<RT,3> partialY = {discreteGradient[0][1], discreteGradient[1][1], discreteGradient[2][1]};
+
+ auto normX = partialX.two_norm();
+ auto normY = partialY.two_norm();
+
+ // Normalize
+ // partialX /= normX;
+ // partialY /= normY;
+
+
+ surfaceNormal[0] = partialX[1]*partialY[2] - partialX[2]*partialY[1];
+ surfaceNormal[1] = partialX[2]*partialY[0] - partialX[0]*partialY[2];
+ surfaceNormal[2] = partialX[0]*partialY[1] - partialX[1]*partialY[0];
+ // printvector(std::cout, surfaceNormal, "surfaceNormal", "--");
+ // printvector(std::cout, surfaceNormal.two_norm(), "surfaceNormal.two_norm()", "--");
+
+ // RT surfaceNormalMag = sqrt(surfaceNormal[0]*surfaceNormal[0] + surfaceNormal[1]*surfaceNormal[1] + surfaceNormal[2]*surfaceNormal[2]);
+ // std::cout << "surfaceNormalMag: " << surfaceNormalMag << std::endl;
+
+
+ FieldMatrix<RT,2,2> secondFF(0);
+
+ for(size_t m=0; m<2; m++)
+ for(size_t n=0; n<2; n++)
+ for(size_t k=0; k<3; k++)
+ {
+ secondFF[m][n] += surfaceNormal[k]*discreteHessian[k][n][m];
+ //Test: transposed version
+ // secondFF[m][n] += discreteHessian[k][m][n]*surfaceNormal[k];
+ // secondFF[0][0] += surfaceNormal[k]*discreteHessian[k][0][0];
+ // secondFF[0][1] += surfaceNormal[k]*discreteHessian[k][1][0];
+ // secondFF[1][0] += surfaceNormal[k]*discreteHessian[k][0][1];
+ // secondFF[1][1] += surfaceNormal[k]*discreteHessian[k][1][1];
+ }
+
+ // substract effective prestrain
+ // G = G + Beff;
+
+ //Take symmetric part?
+ // auto symG = 0.5*(G.transposed() + G);
+
+
+
+
+ // printmatrix(std::cout, G , "matrix G: ", "--");
+ // printmatrix(std::cout, symG , "matrix symG: ", "--");
+
+ // std::cout << " discreteHessian.frobenius_norm2(): " << discreteHessian.frobenius_norm2() << std::endl;
+ // std::cout << " G.frobenius_norm2() : " << G.frobenius_norm2() << std::endl;
+ // std::cout << " DIFFERENCE : " << G.frobenius_norm2()- discreteHessian.frobenius_norm2() << std::endl;
+
+ /**
+ * @brief Compute squared Frobenius norm of the discrete Hessian
+ */
+ // auto quadResult = discreteHessian.frobenius_norm2();
+ // std::cout << "quadResult: " << weight*quadResult << std::endl;
+ // auto weight = geometry.volume()/3.0 * element.geometry().integrationElement(geometry.local(geometry.corner(c)));
+
+ // auto weight = geometry.volume()/3.0;
+
+ // std::cout << "geometry.volume()/3.0 : " << geometry.volume()/3.0 << std::endl;
+
+ // we need to use the surface area of the reference element (which is 1/2)
+ // std::cout << "geometry.volume(): " << geometry.volume() << std::endl;
+ // auto weight = (1.0/6.0) * element.geometry().integrationElement(geometry.local(geometry.corner(c)));
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ // std::cout << "weight: " << weight << std::endl;
+
+ // harmonicEnergy += weight * discreteHessian.frobenius_norm2();
+ // harmonicEnergy += weight * symG.frobenius_norm2();
+
+ // auto term1 = 0.5 * weight * discreteHessian.frobenius_norm2();
+ auto term1 = 0.5 * weight * secondFF.frobenius_norm2();
+
+ RT term2 = 0.0;
+
+ for(size_t i=0; i<2; i++)
+ for(size_t j=0; j<2; j++)
+ {
+ term2 += weight * secondFF[i][j] * Beff[i][j];
+ }
+
+ auto term3 = 0.5 * weight* Beff.frobenius_norm2();
+
+ // std::cout << "term1: " << term1 << std::endl;
+ // std::cout << "term2: " << term2 << std::endl;
+ // std::cout << "term3: " << term3 << std::endl;
+ // std::cout << "term1 - term2 + term3: " << term1 - term2 + term3 << std::endl;
+
+ // harmonicEnergy += weight * G.frobenius_norm2();
+ harmonicEnergy += term1 - term2 + term3; //eigentlich muss man term2 abziehen bei | II - Z |
+
+ // std::cout << "harmonicEnergy: " << harmonicEnergy << std::endl;
+ }
+ #endif
+
+
+ /**
+ * @brief Compute contribution of the force Term
+ *
+ * Integrate the scalar product of the force 'f'
+ * with the local deformation function 'localFunction_'
+ * over the current element.
+ */
+ RT forceTermEnergy = 0;
+
+ if (parameterSet_.get<bool>("assemble_force_term", 1))
+ {
+ for (auto&& quadPoint : quadRule)
+ {
+ auto deformationValue = localFunction_.evaluate(quadPoint.position());
+ auto forceValue = localForce_(quadPoint.position());
+ // printvector(std::cout, deformationValue.globalCoordinates(), "deformationValue", "--");
+ // printvector(std::cout, forceValue, "forceValue", "--");
+ // const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+ auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+ forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
+ // std::cout << "forceTermEnergy:" << forceTermEnergy << std::endl;
+ }
+ }
+
+ return harmonicEnergy - forceTermEnergy;
+ }
+
+ private:
+ LocalDiscreteKirchhoffFunction& localFunction_;
+
+ LocalForce& localForce_;
+
+ P2LocalFiniteElement lagrangeLFE_;
+
+ const Dune::ParameterTree parameterSet_;
+
+ // const PrestrainType& prestrain_;
+
+ // const MatrixType& quadraticForm_;
+ };
+
+} // namespace Dune::GFE
+#endif
\ No newline at end of file
diff --git a/src/macro-problem.cc b/src/macro-problem.cc
index 62526239e0beedc4842380df7b0142f730db7d64..e3102e5334e976c1563d40131af2d655bb989054 100644
--- a/src/macro-problem.cc
+++ b/src/macro-problem.cc
@@ -53,7 +53,8 @@
#include <dune/gfe/energies/discretekirchhoffbendingenergyconforming.hh>
#include <dune/gfe/energies/discretekirchhoffbendingenergyzienkiewicz.hh>
#include <dune/gfe/energies/discretekirchhoffbendingenergyzienkiewiczprojected.hh>
-#include <dune/gfe/energies/discretekirchhoffbendingenergyprestrained.hh>
+
+#include <dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh>
#include <dune/gfe/spaces/stiefelmanifold.hh>
// #include <dune/gfe/spaces/stiefelmatrix.hh>
@@ -309,29 +310,28 @@ int main(int argc, char *argv[])
auto forceGVF = Dune::Functions::makeGridViewFunction(pythonForce, gridView);
auto localForce = localFunction(forceGVF);
-// "def force(x): return [0, 0, 0]"
- /**
- * @brief Read effective prestrain Tensor
- */
- Dune::FieldVector<adouble,3> Beffvec = parameterSet.get<Dune::FieldVector<adouble,3>>("effectivePrestrain", {0.0, 0.0, 0.0});
- Dune::FieldMatrix<adouble,2,2> effectivePrestrain = {{(adouble)Beffvec[0], (adouble)Beffvec[2]}, {(adouble)Beffvec[2],(adouble)Beffvec[1]} };
- printmatrix(std::cout, effectivePrestrain, "effective prestrain (Beff): ", "--");
-
-
- /**
- * @brief Get effective quadratic form Qhom
- *
- * input-vector: [q_1,q_2,q_3,q_12,q_13,q_23]
- * is assembled into a matrix where the off-diagonal entries are divided by 2 (compare with definition in the paper)
- * ( q_1 , 0.5*q_12 , 0.5*q_13 )
- * Q = ( 0.5*q_12 , q_2 , 0.5*q_23 )
- * ( 0.5*q_13 , 0.5*q_23 , q_3 )
- */
- Dune::FieldVector<adouble,6> Qhomvec = parameterSet.get<Dune::FieldVector<adouble,6>>("effectiveQuadraticForm", {1.0, 1.0, 1.0, 0.0, 0.0, 0.0});
- Dune::FieldMatrix<adouble,3,3> Qhom = {{(adouble)Qhomvec[0], (adouble)0.5*Qhomvec[3], (adouble)0.5*Qhomvec[4]},
- {(adouble)0.5*Qhomvec[3], (adouble)Qhomvec[1], (adouble)0.5*Qhomvec[5]},
- {(adouble)0.5*Qhomvec[4], (adouble)0.5*Qhomvec[5], (adouble)Qhomvec[2]}};
- printmatrix(std::cout, Qhom, "effective quadratic form (Qhom): ", "--");
+ // /**
+ // * @brief Read effective prestrain Tensor
+ // */
+ // Dune::FieldVector<adouble,3> Beffvec = parameterSet.get<Dune::FieldVector<adouble,3>>("effectivePrestrain", {0.0, 0.0, 0.0});
+ // Dune::FieldMatrix<adouble,2,2> effectivePrestrain = {{(adouble)Beffvec[0], (adouble)Beffvec[2]}, {(adouble)Beffvec[2],(adouble)Beffvec[1]} };
+ // printmatrix(std::cout, effectivePrestrain, "effective prestrain (Beff): ", "--");
+
+
+ // /**
+ // * @brief Get effective quadratic form Qhom
+ // *
+ // * input-vector: [q_1,q_2,q_3,q_12,q_13,q_23]
+ // * is assembled into a matrix where the off-diagonal entries are divided by 2 (compare with definition in the paper)
+ // * ( q_1 , 0.5*q_12 , 0.5*q_13 )
+ // * Q = ( 0.5*q_12 , q_2 , 0.5*q_23 )
+ // * ( 0.5*q_13 , 0.5*q_23 , q_3 )
+ // */
+ // Dune::FieldVector<adouble,6> Qhomvec = parameterSet.get<Dune::FieldVector<adouble,6>>("effectiveQuadraticForm", {1.0, 1.0, 1.0, 0.0, 0.0, 0.0});
+ // Dune::FieldMatrix<adouble,3,3> Qhom = {{(adouble)Qhomvec[0], (adouble)0.5*Qhomvec[3], (adouble)0.5*Qhomvec[4]},
+ // {(adouble)0.5*Qhomvec[3], (adouble)Qhomvec[1], (adouble)0.5*Qhomvec[5]},
+ // {(adouble)0.5*Qhomvec[4], (adouble)0.5*Qhomvec[5], (adouble)Qhomvec[2]}};
+ // printmatrix(std::cout, Qhom, "effective quadratic form (Qhom): ", "--");
////////////////////////////////////////////////////////////////////////
// Create an assembler for the Discrete Kirchhoff Energy Functional (using ADOL-C)
@@ -349,7 +349,8 @@ int main(int argc, char *argv[])
/**
* @brief Setup energy featuring prestrain
*/
- auto localEnergy_prestrain = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyPrestrained<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce, Qhom, effectivePrestrain, parameterSet);
+ // auto localEnergy_prestrain = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyPrestrained<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce, Qhom, effectivePrestrain, parameterSet);
+ auto localEnergy_prestrain = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyPrestrained<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce, parameterSet, pyModule);
LocalGeodesicFEADOLCStiffness<CoefficientBasis, Coefficient> localGFEADOLCStiffness_conforming(localEnergy_conforming.get());
LocalGeodesicFEADOLCStiffness<CoefficientBasis, Coefficient> localGFEADOLCStiffness_nonconforming(localEnergy_nonconforming.get());