From 39336b407d708cc7dac2d22c4fa0a0e10ddbee0e Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 28 Dec 2015 21:24:30 +0100
Subject: [PATCH] Merge the MixedCosseratEnergy class into the
CosseratEnergyLocalStiffness class
Because these two classes are for the most part identical.
Currently, this merged version involves some trickery, including some template
specializations and multiple inheritances. I hope to get rid of that, eventually.
---
dune/gfe/cosseratenergystiffness.hh | 211 ++++++++++++++-
dune/gfe/mixedcosseratenergy.hh | 400 ----------------------------
src/mixed-cosserat-continuum.cc | 4 +-
3 files changed, 211 insertions(+), 404 deletions(-)
delete mode 100644 dune/gfe/mixedcosseratenergy.hh
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index e51ab6c8..1bfd818b 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -9,6 +9,7 @@
#include <dune/fufem/boundarypatch.hh>
#include "localgeodesicfestiffness.hh"
+#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
#include "localgeodesicfefunction.hh"
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/tensor3.hh>
@@ -19,10 +20,44 @@
//#define QUADRATIC_MEMBRANE_ENERGY
+/** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
+ *
+ * We instantiate the CosseratEnergyLocalStiffness class with two different kinds of Basis:
+ * A scalar one and a composite one that combines two scalar ones. But code for accessing the
+ * finite elements in the basis tree only work for one kind of basis, not for the other.
+ * To allow both kinds of basis in a single class we need this trickery below.
+ */
+template <class Basis, std::size_t i>
+class LocalFiniteElementFactory
+{
+public:
+ static auto get(const typename Basis::LocalView& localView,
+ std::integral_constant<std::size_t, i> iType)
+ -> decltype(localView.tree().child(iType).finiteElement())
+ {
+ return localView.tree().child(iType).finiteElement();
+ }
+};
+
+/** \brief Specialize for scalar bases, here we cannot call tree().child() */
+template <class GridView, int order, std::size_t i>
+class LocalFiniteElementFactory<Dune::Functions::PQkNodalBasis<GridView,order>,i>
+{
+public:
+ static auto get(const typename Dune::Functions::PQkNodalBasis<GridView,order>::LocalView& localView,
+ std::integral_constant<std::size_t, i> iType)
+ -> decltype(localView.tree().finiteElement())
+ {
+ return localView.tree().finiteElement();
+ }
+};
template<class Basis, int dim, class field_type=double>
class CosseratEnergyLocalStiffness
- : public LocalGeodesicFEStiffness<Basis,RigidBodyMotion<field_type,dim> >
+ : public LocalGeodesicFEStiffness<Basis,RigidBodyMotion<field_type,dim> >,
+ public MixedLocalGeodesicFEStiffness<Basis,
+ RealTuple<field_type,dim>,
+ Rotation<field_type,dim> >
{
// grid types
typedef typename Basis::GridView GridView;
@@ -113,6 +148,36 @@ public: // for testing
}
+ /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
+ \param value Value of the gfe function at a certain point
+ \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
+ \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
+ */
+ static void computeDR(const Rotation<field_type,3>& value,
+ const Dune::FieldMatrix<field_type,4,gridDim>& derivative,
+ Tensor3<field_type,3,3,gridDim>& DR)
+ {
+ // The LocalGFEFunction class gives us the derivatives of the orientation variable,
+ // but as a map into quaternion space. To obtain matrix coordinates we use the
+ // chain rule, which means that we have to multiply the given derivative with
+ // the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
+ // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
+ // However, since the directors of a given unit quaternion are the _columns_ of the
+ // corresponding orthogonal matrix, we need to invert the i and j indices
+ //
+ // So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
+ Tensor3<field_type,3 , 3, 4> dd_dq;
+ value.getFirstDerivativesOfDirectors(dd_dq);
+
+ DR = field_type(0);
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<gridDim; k++)
+ for (int l=0; l<4; l++)
+ DR[i][j][k] += dd_dq[j][i][l] * derivative[l][k];
+
+ }
+
public:
/** \brief Constructor with a set of material parameters
@@ -148,6 +213,11 @@ public:
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override;
+ /** \brief Assemble the energy for a single element */
+ RT energy (const typename Basis::LocalView& localView,
+ const std::vector<RealTuple<field_type,dim> >& localDisplacementConfiguration,
+ const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const override;
+
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the first equation of (4.4) in Neff's paper
*/
@@ -300,7 +370,9 @@ energy(const typename Basis::LocalView& localView,
RT energy = 0;
auto element = localView.element();
- const auto& localFiniteElement = localView.tree().finiteElement();
+
+ using namespace Dune::TypeTree::Indices;
+ const auto& localFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
@@ -415,5 +487,140 @@ energy(const typename Basis::LocalView& localView,
return energy;
}
+template <class Basis, int dim, class field_type>
+typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT
+CosseratEnergyLocalStiffness<Basis,dim,field_type>::
+energy(const typename Basis::LocalView& localView,
+ const std::vector<RealTuple<field_type,dim> >& localDeformationConfiguration,
+ const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const
+{
+ auto element = localView.element();
+
+ RT energy = 0;
+
+ using namespace Dune::TypeTree::Indices;
+ const auto& deformationLocalFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
+ const auto& orientationLocalFiniteElement = LocalFiniteElementFactory<Basis,1>::get(localView,_1);
+
+ typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
+ LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
+
+ typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
+ LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localOrientationConfiguration);
+
+ // \todo Implement smarter quadrature rule selection for more efficiency, i.e., less evaluations of the Rotation GFE function
+ int quadOrder = deformationLocalFiniteElement.localBasis().order() * ((element.type().isSimplex()) ? 1 : gridDim);
+
+ const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+ for (size_t pt=0; pt<quad.size(); pt++)
+ {
+ // Local position of the quadrature point
+ const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+
+ const DT integrationElement = element.geometry().integrationElement(quadPos);
+
+ const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+
+ DT weight = quad[pt].weight() * integrationElement;
+
+ // The value of the local deformation
+ RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
+ Rotation<field_type,dim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
+
+ // The derivative of the local function defined on the reference element
+ typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
+ typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
+
+ // The derivative of the function defined on the actual element
+ typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
+ typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
+
+ for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
+ jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
+
+ for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
+ jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
+
+ /////////////////////////////////////////////////////////
+ // compute U, the Cosserat strain
+ /////////////////////////////////////////////////////////
+ static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
+
+ //
+ Dune::FieldMatrix<field_type,dim,dim> R;
+ orientationValue.matrix(R);
+
+ Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(deformationDerivative,R);
+
+ //////////////////////////////////////////////////////////
+ // Compute the derivative of the rotation
+ // Note: we need it in matrix coordinates
+ //////////////////////////////////////////////////////////
+
+ Tensor3<field_type,3,3,gridDim> DR;
+ computeDR(orientationValue, orientationDerivative, DR);
+
+ // Add the local energy density
+ if (gridDim==2) {
+#ifdef QUADRATIC_MEMBRANE_ENERGY
+ //energy += weight * thickness_ * quadraticMembraneEnergy(U.matrix());
+ energy += weight * thickness_ * longQuadraticMembraneEnergy(U);
+#else
+ energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
+#endif
+ energy += weight * thickness_ * curvatureEnergy(DR);
+ energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
+ } else if (gridDim==3) {
+ energy += weight * quadraticMembraneEnergy(U);
+ energy += weight * curvatureEnergy(DR);
+ } else
+ DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
+
+ }
+
+ //////////////////////////////////////////////////////////////////////////////
+ // Assemble boundary contributions
+ //////////////////////////////////////////////////////////////////////////////
+
+ if (not neumannFunction_)
+ return energy;
+
+ for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
+ {
+ if (not neumannBoundary_ or not neumannBoundary_->contains(it))
+ continue;
+
+ const auto& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
+
+ for (size_t pt=0; pt<quad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
+
+ const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
+
+ // The value of the local function
+ RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
+
+ // Value of the Neumann data at the current position
+ Dune::FieldVector<double,3> neumannValue;
+
+ if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
+ dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
+ else
+ neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+
+ // Only translational dofs are affected by the Neumann force
+ for (size_t i=0; i<neumannValue.size(); i++)
+ energy += thickness_ * (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
+
+ }
+
+ }
+
+ return energy;
+}
+
#endif //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
diff --git a/dune/gfe/mixedcosseratenergy.hh b/dune/gfe/mixedcosseratenergy.hh
deleted file mode 100644
index 050a9269..00000000
--- a/dune/gfe/mixedcosseratenergy.hh
+++ /dev/null
@@ -1,400 +0,0 @@
-#ifndef DUNE_GFE_MIXEDCOSSERATENERGY_HH
-#define DUNE_GFE_MIXEDCOSSERATENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/fufem/functions/virtualgridfunction.hh>
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/tensor3.hh>
-#include <dune/gfe/orthogonalmatrix.hh>
-#include <dune/gfe/cosseratstrain.hh>
-
-#define DONT_USE_CURL
-
-//#define QUADRATIC_MEMBRANE_ENERGY
-
-
-template<class Basis, int dim, class field_type=double>
-class MixedCosseratEnergy
- : public MixedLocalGeodesicFEStiffness<Basis,
- RealTuple<field_type,dim>,
- Rotation<field_type,dim> >
-{
- // grid types
- typedef typename Basis::GridView GridView;
- typedef typename GridView::ctype DT;
- typedef field_type RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- enum {gridDim=GridView::dimension};
-
-
- /** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] + A[j][i]);
- return result;
- }
-
- /** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
- static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- result[i][j] = 0.5 * (A[i][j] - A[j][i]);
- return result;
- }
-
- /** \brief Return the square of the trace of a matrix */
- template <int N>
- static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
- {
- field_type trace = 0;
- for (int i=0; i<N; i++)
- trace += A[i][i];
- return trace*trace;
- }
-
- /** \brief Compute the (row-wise) curl of a matrix R \f$
- \param DR The partial derivatives of the matrix R
- */
- static Dune::FieldMatrix<field_type,dim,dim> curl(const Tensor3<field_type,dim,dim,dim>& DR)
- {
- Dune::FieldMatrix<field_type,dim,dim> result;
-
- for (int i=0; i<dim; i++) {
- result[i][0] = DR[i][2][1] - DR[i][1][2];
- result[i][1] = DR[i][0][2] - DR[i][2][0];
- result[i][2] = DR[i][1][0] - DR[i][0][1];
- }
-
- return result;
- }
-
-public: // for testing
- /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
- \param value Value of the gfe function at a certain point
- \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
- \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
- */
- static void computeDR(const Rotation<field_type,3>& value,
- const Dune::FieldMatrix<field_type,4,gridDim>& derivative,
- Tensor3<field_type,3,3,gridDim>& DR)
- {
- // The LocalGFEFunction class gives us the derivatives of the orientation variable,
- // but as a map into quaternion space. To obtain matrix coordinates we use the
- // chain rule, which means that we have to multiply the given derivative with
- // the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
- // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
- // However, since the directors of a given unit quaternion are the _columns_ of the
- // corresponding orthogonal matrix, we need to invert the i and j indices
- //
- // So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
- Tensor3<field_type,3 , 3, 4> dd_dq;
- value.getFirstDerivativesOfDirectors(dd_dq);
-
- DR = field_type(0);
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- for (int k=0; k<gridDim; k++)
- for (int l=0; l<4; l++)
- DR[i][j][k] += dd_dq[j][i][l] * derivative[l][k];
-
- }
-
-public:
-
- /** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
- */
- MixedCosseratEnergy(const Dune::ParameterTree& parameters,
- const BoundaryPatch<GridView>* neumannBoundary,
- const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction)
- : neumannBoundary_(neumannBoundary),
- neumannFunction_(neumannFunction)
- {
- // The shell thickness
- thickness_ = parameters.template get<double>("thickness");
-
- // Lame constants
- mu_ = parameters.template get<double>("mu");
- lambda_ = parameters.template get<double>("lambda");
-
- // Cosserat couple modulus, preferably 0
- mu_c_ = parameters.template get<double>("mu_c");
-
- // Length scale parameter
- L_c_ = parameters.template get<double>("L_c");
-
- // Curvature exponent
- q_ = parameters.template get<double>("q");
-
- // Shear correction factor
- kappa_ = parameters.template get<double>("kappa");
- }
-
- /** \brief Assemble the energy for a single element */
- RT energy (const typename Basis::LocalView& localView,
- const std::vector<RealTuple<field_type,dim> >& localDisplacementConfiguration,
- const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const;
-
- /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
- * the first equation of (4.4) in Neff's paper
- */
- RT quadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
- {
- Dune::FieldMatrix<field_type,3,3> UMinus1 = U;
- for (int i=0; i<dim; i++)
- UMinus1[i][i] -= 1;
-
- return mu_ * sym(UMinus1).frobenius_norm2()
- + mu_c_ * skew(UMinus1).frobenius_norm2()
- + (mu_*lambda_)/(2*mu_ + lambda_) * traceSquared(sym(UMinus1));
- }
-
- /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
- * the second equation of (4.4) in Neff's paper
- */
- RT longQuadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
- {
- RT result = 0;
-
- // shear-stretch energy
- Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
- for (int i=0; i<dim-1; i++)
- for (int j=0; j<dim-1; j++)
- sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
-
- result += mu_ * sym2x2.frobenius_norm2();
-
- // first order drill energy
- Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
- for (int i=0; i<dim-1; i++)
- for (int j=0; j<dim-1; j++)
- skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
-
- result += mu_c_ * skew2x2.frobenius_norm2();
-
-
- // classical transverse shear energy
- result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
-
- // elongational stretch energy
- result += mu_*lambda_ / (2*mu_ + lambda_) * traceSquared(sym2x2);
-
- return result;
- }
-
- /** \brief Energy for large-deformation problems (private communication by Patrizio Neff)
- */
- RT nonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
- {
- Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
- for (int i=0; i<dim; i++)
- UMinus1[i][i] -= 1;
-
- RT detU = U.determinant();
-
- return mu_ * sym(UMinus1).frobenius_norm2()
- + (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
- }
-
- RT curvatureEnergy(const Tensor3<field_type,3,3,gridDim>& DR) const
- {
-#ifdef DONT_USE_CURL
- return mu_ * std::pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0);
-#else
- return mu_ * std::pow(L_c_ * L_c_ * curl(DR).frobenius_norm2(),q_/2.0);
-#endif
- }
-
- RT bendingEnergy(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
- {
- // left-multiply the derivative of the third director (in DR[][2][]) with R^T
- Dune::FieldMatrix<field_type,3,3> RT_DR3(0);
- for (int i=0; i<3; i++)
- for (int j=0; j<gridDim; j++)
- for (int k=0; k<3; k++)
- RT_DR3[i][j] += R[k][i] * DR[k][2][j];
-
- return mu_ * sym(RT_DR3).frobenius_norm2()
- + mu_c_ * skew(RT_DR3).frobenius_norm2()
- + mu_*lambda_/(2*mu_+lambda_) * traceSquared(RT_DR3);
- }
-
- /** \brief The shell thickness */
- double thickness_;
-
- /** \brief Lame constants */
- double mu_, lambda_;
-
- /** \brief Cosserat couple modulus, preferably 0 */
- double mu_c_;
-
- /** \brief Length scale parameter */
- double L_c_;
-
- /** \brief Curvature exponent */
- double q_;
-
- /** \brief Shear correction factor */
- double kappa_;
-
- /** \brief The Neumann boundary */
- const BoundaryPatch<GridView>* neumannBoundary_;
-
- /** \brief The function implementing the Neumann data */
- const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction_;
-};
-
-template <class Basis, int dim, class field_type>
-typename MixedCosseratEnergy<Basis,dim,field_type>::RT
-MixedCosseratEnergy<Basis,dim,field_type>::
-energy(const typename Basis::LocalView& localView,
- const std::vector<RealTuple<field_type,dim> >& localDeformationConfiguration,
- const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const
-{
- typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
- auto element = localView.element();
-
- RT energy = 0;
-
- using namespace Dune::TypeTree::Indices;
- const auto& deformationLocalFiniteElement = localView.tree().child(_0).finiteElement();
- const auto& orientationLocalFiniteElement = localView.tree().child(_1).finiteElement();
-
- typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
- LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
-
- typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
- LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localOrientationConfiguration);
-
- // \todo Implement smarter quadrature rule selection for more efficiency, i.e., less evaluations of the Rotation GFE function
- int quadOrder = deformationLocalFiniteElement.localBasis().order() * ((element.type().isSimplex()) ? 1 : gridDim);
-
- const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
- const DT integrationElement = element.geometry().integrationElement(quadPos);
-
- const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
- DT weight = quad[pt].weight() * integrationElement;
-
- // The value of the local deformation
- RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
- Rotation<field_type,dim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
-
- // The derivative of the local function defined on the reference element
- typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
- typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
-
- // The derivative of the function defined on the actual element
- typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
- typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
-
- for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
-
- for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
-
- /////////////////////////////////////////////////////////
- // compute U, the Cosserat strain
- /////////////////////////////////////////////////////////
- static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
-
- //
- Dune::FieldMatrix<field_type,dim,dim> R;
- orientationValue.matrix(R);
-
- Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(deformationDerivative,R);
-
- //////////////////////////////////////////////////////////
- // Compute the derivative of the rotation
- // Note: we need it in matrix coordinates
- //////////////////////////////////////////////////////////
-
- Tensor3<field_type,3,3,gridDim> DR;
- computeDR(orientationValue, orientationDerivative, DR);
-
- // Add the local energy density
- if (gridDim==2) {
-#ifdef QUADRATIC_MEMBRANE_ENERGY
- //energy += weight * thickness_ * quadraticMembraneEnergy(U.matrix());
- energy += weight * thickness_ * longQuadraticMembraneEnergy(U);
-#else
- energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
-#endif
- energy += weight * thickness_ * curvatureEnergy(DR);
- energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
- } else if (gridDim==3) {
- energy += weight * quadraticMembraneEnergy(U);
- energy += weight * curvatureEnergy(DR);
- } else
- DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
-
- }
-
- //////////////////////////////////////////////////////////////////////////////
- // Assemble boundary contributions
- //////////////////////////////////////////////////////////////////////////////
-
- if (not neumannFunction_)
- return energy;
-
- for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
- {
- if (not neumannBoundary_ or not neumannBoundary_->contains(it))
- continue;
-
- const Dune::QuadratureRule<DT, gridDim-1>& quad
- = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
-
- const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
-
- // The value of the local function
- RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
-
- // Value of the Neumann data at the current position
- Dune::FieldVector<double,3> neumannValue;
-
- if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
- dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
- else
- neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
-
- // Only translational dofs are affected by the Neumann force
- for (size_t i=0; i<neumannValue.size(); i++)
- energy += thickness_ * (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
-
- }
-
- }
-
- return energy;
-}
-
-#endif //#ifndef DUNE_GFE_MIXEDCOSSERATENERGY_HH
-
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
index 57f6211b..d7189d7f 100644
--- a/src/mixed-cosserat-continuum.cc
+++ b/src/mixed-cosserat-continuum.cc
@@ -37,7 +37,7 @@
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
-#include <dune/gfe/mixedcosseratenergy.hh>
+#include <dune/gfe/cosseratenergystiffness.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/mixedgfeassembler.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
@@ -296,7 +296,7 @@ int main (int argc, char *argv[]) try
}
// Assembler using ADOL-C
- MixedCosseratEnergy<decltype(compositeBasis),
+ CosseratEnergyLocalStiffness<decltype(compositeBasis),
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
&neumannBoundary,
neumannFunction.get());
--
GitLab