Skip to content
Snippets Groups Projects

WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis

Open WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis
feature/bendingIsometries into master
Open Klaus Böhnlein requested to merge feature/bendingIsometries into master
Viewing commit e526c81d
Show latest version
2 files
+ 435
− 6
Compare changes
  • Side-by-side
  • Inline
Files
2
dune/gfe/energies/discretekirchhoffbendingenergyconforming.hh 0 → 100644
+ 426
− 0
#ifndef DUNE_GFE_ENERGIES_DISCRETEKIRCHHOFFBENDINGENERGYCONFORMING_HH
#define DUNE_GFE_ENERGIES_DISCRETEKIRCHHOFFBENDINGENERGYCONFORMING_HH
#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/localenergy.hh>
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
#include <dune/gfe/polardecomposition.hh>
#include <dune/gfe/rotation.hh>
#include <dune/gfe/localprojectedfefunction.hh>
/**
* Conforming version of 'discretekirchhoffbendingenergy.hh'.
*
* In the future this will be integrated into 'discretekirchhoffbendingenergy.hh'
* with a 'conforming-Flag'.
*
*
*
* \brief Assemble the discrete Kirchhoff bending energy for a single element.
*
* The Kirchhoff bending energy consists of two parts:
*
* 1. The harmonic energy of a discrete Jacobian operator
* (acting as a discrete Hessian operator)
*
* 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 DiscreteKirchhoffBendingEnergyConforming
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
//using Dune::GFE::LocalEnergy<Basis,TargetSpace>::RT;
typedef typename TargetSpace::ctype RT;
typedef typename LocalDiscreteKirchhoffFunction::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:
DiscreteKirchhoffBendingEnergyConforming(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 std::vector<TargetSpace>& localSolution) const
{
RT energy = 0;
/**
* @brief create a projection-based finite element of order 2 with values in SO(3).
* We collect the coefficients along the way and store them in the following vector
*/
std::vector<Rotation<RT,3>> gfeCoefficients(lagrangeLFE_.size());
#if 1
// TODO: More serious order selection
int quadOrder = 3;
#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 &quad = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
const auto element = localView.element();
localFunction_.bind(element);
localForce_.bind(element);
//Update local coefficients
localFunction_.updateLocalCoefficients(localSolution,element);
auto geometry = element.geometry();
auto gridView = localFunction_.gridView();
/**
* @brief Compute energy contribution from the harmonic energy:
*/
RT harmonicEnergy = 0;
/**
* @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;
discreteJacobianCoefficients.resize(lagrangeLFE_.size());
for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
{
// std::cout << lagrangeLFE_.localCoefficients().localKey(i) << std::endl;
// Determine coefficients for basis functions assigned to a vertex
if (lagrangeLFE_.localCoefficients().localKey(i).codim() == 2)
{
size_t nodeIndex = lagrangeLFE_.localCoefficients().localKey(i).subEntity();
// std::cout << "vertex with number:" << nodeIndex << std::endl;
// std::cout << "and coordinates: " << geometry.corner(nodeIndex) << std::endl;
DerivativeType whJacobianValue = localFunction_.evaluateDerivative(geometry.local(geometry.corner(nodeIndex)));
discreteJacobianCoefficients[i] = whJacobianValue;
/**
* save coefficient for GFE:
*
* @brief Convert 3x2 Matrix to Rotation<RT,3>
* TODO: write a routine for this
*/
// extract columns
FieldVector<RT,3> currentDerivative_x(0);
FieldVector<RT,3> currentDerivative_y(0);
for(std::size_t j=0; j<3; j++)
{
currentDerivative_x[j] = whJacobianValue[j][0];
currentDerivative_y[j] = whJacobianValue[j][1];
}
//cross product to get last column of rotation matrix
FieldVector<RT,3> cross = {currentDerivative_x[1]*currentDerivative_y[2] - currentDerivative_x[2]*currentDerivative_y[1],
currentDerivative_x[2]*currentDerivative_y[0] - currentDerivative_x[0]*currentDerivative_y[2],
currentDerivative_x[0]*currentDerivative_y[1] - currentDerivative_x[1]*currentDerivative_y[0]};
// extend to 3x3 matrix:
FieldMatrix<RT,3,3> rotMatrix(0);
for(std::size_t k=0; k<3 ; k++)
{
rotMatrix[0][k] = currentDerivative_x[k];
rotMatrix[1][k] = currentDerivative_y[k];
rotMatrix[2][k] = cross[k];
}
Rotation<RT,3> tmpRot;
tmpRot.set(rotMatrix);
gfeCoefficients[i] = tmpRot;
// std::vector<FieldVector<double, 1>> p2Values;
// lagrangeLFE_.localBasis().evaluateFunction(geometry.local(geometry.corner(nodeIndex)), p2Values);
// std::cout <<"TEST: " << i << "-th basis function value at current node:" << p2Values[i] << std::endl;
}
// Determine coefficients for basis functions assigned to an edge
if (lagrangeLFE_.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 = lagrangeLFE_.localCoefficients().localKey(i).subEntity();
auto edgeGeometry = element.template subEntity<1>(edgeIndex).geometry();
auto edgeEntity = element.template subEntity<1>(edgeIndex);
// std::cout << "Edge with number:" << edgeIndex << std::endl;
// std::cout << "edgeGeometry.center()" << edgeGeometry.center() << std::endl;
// Get the right orientation for tangent and normal vectors on current edge.
const auto &refElement = referenceElement(element);
const auto &indexSet = gridView.indexSet();
// 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);
int edgeOrientation = 1;
if ((localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1))
{
edgeOrientation = -1;
// std::cout << "edge orientation reversed" << std::endl;
}
// std::cout << "Edge vertex 0:" << geometry.corner(localV0) << std::endl;
// std::cout << "Edge vertex 1:" << geometry.corner(localV1) << std::endl;
// for(std::size_t n = 0; n<edgeGeometry.corners(); n++)
// {
// std::cout <<"edge corner number " << n << " : " << edgeGeometry.corner(n) << std::endl;
// }
// Get tangent vector on current edge
FieldVector<RT,2> tmp = (geometry.corner(localV1) - geometry.corner(localV0)) / abs(edgeGeometry.volume());
//convert to FieldMatrix
FieldMatrix<RT,2,1> edgeTangent = {tmp[0], tmp[1]};
// printmatrix(std::cout, edgeTangent , "edgeTangent ", "--");
edgeTangent *= edgeOrientation;
// printmatrix(std::cout, edgeTangent , "edgeTangent with right orientation:", "--");
// rotate edge tangent vector by pi/2 clockwise
FieldMatrix<RT,2,1> edgeNormal = {edgeTangent[1], -1.0*edgeTangent[0]};
// printmatrix(std::cout, edgeNormal , "edgeNormal", "--");
// Evaluate Jacobian of localfunctions_ at edge-vertices and edge-midpoints.
DerivativeType whJacobianValueCenter = localFunction_.evaluateDerivative(geometry.local(edgeGeometry.center()));
DerivativeType whJacobianValueFirstVertex = localFunction_.evaluateDerivative(geometry.local(geometry.corner(localV0)));
DerivativeType whJacobianValueSecondVertex = localFunction_.evaluateDerivative(geometry.local(geometry.corner(localV1)));
// @OS: Is there a elementwise/hadamard matrix multiplication in Dune?
FieldMatrix<RT,3,gridDim> normalComponent(0);
FieldMatrix<RT,3,gridDim> tangentialComponent(0);
auto normalComponentFactor = 0.5 * (whJacobianValueFirstVertex + whJacobianValueSecondVertex) * edgeNormal;
auto tangentialComponentFactor = whJacobianValueCenter * edgeTangent;
// printmatrix(std::cout, normalComponentFactor , "normalComponentFactor ", "--");
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];
}
// printmatrix(std::cout, normalComponent , "normalComponent ", "--");
// printmatrix(std::cout, tangentialComponent , "tangentialComponent ", "--");
discreteJacobianCoefficients[i] = normalComponent + tangentialComponent;
/**
* save coefficient for GFE:
*
* @brief Convert 3x2 Matrix to Rotation<RT,3>
* TODO: write a routine for this
*/
// extract columns
FieldVector<RT,3> currentDerivative_x(0);
FieldVector<RT,3> currentDerivative_y(0);
for(std::size_t j=0; j<3; j++)
{
currentDerivative_x[j] = discreteJacobianCoefficients[i][j][0];
currentDerivative_y[j] = discreteJacobianCoefficients[i][j][1];
}
//cross product to get last column of rotation matrix
FieldVector<RT,3> cross = {currentDerivative_x[1]*currentDerivative_y[2] - currentDerivative_x[2]*currentDerivative_y[1],
currentDerivative_x[2]*currentDerivative_y[0] - currentDerivative_x[0]*currentDerivative_y[2],
currentDerivative_x[0]*currentDerivative_y[1] - currentDerivative_x[1]*currentDerivative_y[0]};
// extend to 3x3 matrix:
FieldMatrix<RT,3,3> rotMatrix(0);
for(std::size_t k=0; k<3 ; k++)
{
rotMatrix[0][k] = currentDerivative_x[k];
rotMatrix[1][k] = currentDerivative_y[k];
rotMatrix[2][k] = cross[k];
}
// The Jacobians on the edge-midpoints are not orthogonal a priori. Here we need to project!
Dune::GFE::PolarDecomposition()(rotMatrix);
Rotation<RT,3> tmpRot;
tmpRot.set(rotMatrix);
gfeCoefficients[i] = tmpRot;
}
}
// std::cout << "------ print Coefficients -----" << std::endl;
// for(std::size_t i=0; i<lagrangeLFE_.size(); i++)
// {
// std::cout << i << "-th Coefficient: " << std::endl;
// printmatrix(std::cout, discreteJacobianCoefficients[i] , "discreteJacobianCoefficients[i] ", "--");
// }
/**
* @brief Now, let's create the geometric finite element
*
*/
GFE::LocalProjectedFEFunction<gridDim, double, P2LocalFiniteElement, Rotation<RT,3>> localGFE(lagrangeLFE_,gfeCoefficients);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
const auto integrationElement = element.geometry().integrationElement(quadPos);
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
auto weight = quad[pt].weight() * integrationElement;
// The derivative of the local function defined on the reference element
auto referenceDerivative = localGFE.evaluateDerivative(quadPos);
// printmatrix(std::cout, referenceDerivative , "referenceDerivative ", "--");
// extract first two columns
FieldMatrix<RT,3,gridDim> derivativeMatrix(0);
for(std::size_t j=0; j<3; j++) // are we eliminating the right row?
{
derivativeMatrix[j][0] = referenceDerivative[j][0];
derivativeMatrix[j][1] = referenceDerivative[j][1];
}
// std::cout << "jacobianInverseTransposed.N() "<< jacobianInverseTransposed.N() << std::endl;
// std::cout << "TargetSpace::embeddedDim"<< TargetSpace::embeddedDim << std::endl;
// std::cout << "jacobianInverseTransposed.M()"<< jacobianInverseTransposed.M() << std::endl;
// Compute the Frobenius norm squared of the derivative. This is the correct term
// if both domain and target space use the metric inherited from an embedding.
for (size_t i=0; i<jacobianInverseTransposed.N(); i++)
for (int j=0; j<3; j++)
{
RT entry = 0;
for (size_t k=0; k<jacobianInverseTransposed.M(); k++)
entry += jacobianInverseTransposed[i][k] * derivativeMatrix[j][k];
harmonicEnergy += weight * entry * entry;
}
}
/**
* @brief Compute the discrete Hessian as a linear combination of
* Jacobians of the P2-basis functions with the coefficients given
* by 'discreteJacobianCoefficients'.
*
* The entries of discrete Hessian are stored in a Fieldmatrix<double,6,2>
* to compute the squared Frobenius norm via the method .frobenius_norm2()
*/
for (auto&& quadPoint : quad)
{
// Get Jacobians of the P2-Basis functions on current quadrature point
std::vector<FieldMatrix<double, 1, gridDim>> referenceGradients;
// std::vector<FieldMatrix<RT, 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]);
// Compute discrete Hessian entries on current quadrature point.
// @OS: is there a better way to do it by using some sort of tensor data structure?
// FieldMatrix<double, 3*gridDim,gridDim> discreteHessian(0);
FieldMatrix<RT, 3*gridDim,gridDim> discreteHessian(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++)
{
discreteHessian[k*gridDim][l] += discreteJacobianCoefficients[i][k][l]*gradients[i][0];
discreteHessian[k*gridDim+1][l] += discreteJacobianCoefficients[i][k][l]*gradients[i][1];
}
auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
/**
* @brief Compute squared Frobenius norm of the discrete Hessian
*/
// auto quadResult = discreteHessian.frobenius_norm2();
// std::cout << "quadResult: " << weight*quadResult << std::endl;
harmonicEnergy += weight * discreteHessian.frobenius_norm2();
// std::cout << "harmonicEnergy: " << harmonicEnergy << std::endl;
}
/**
* @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;
for (auto&& quadPoint : quad)
{
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;
}
// std::cout << "energy output:" << 0.5 * harmonicEnergy - forceTermEnergy << std::endl;
return 0.5 * harmonicEnergy - forceTermEnergy;
}
private:
LocalDiscreteKirchhoffFunction& localFunction_;
LocalForce& localForce_;
P2LocalFiniteElement lagrangeLFE_;
};
} // namespace Dune::GFE
#endif