Skip to content
Snippets Groups Projects
Commit 39336b40 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

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.
parent f807505b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/cosseratenergystiffness.hh
+ 209
− 2
View file @ 39336b40
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
#include <dune/fufem/boundarypatch.hh> #include <dune/fufem/boundarypatch.hh>
#include "localgeodesicfestiffness.hh" #include "localgeodesicfestiffness.hh"
#include <dune/gfe/mixedlocalgeodesicfestiffness.hh>
#include "localgeodesicfefunction.hh" #include "localgeodesicfefunction.hh"
#include <dune/gfe/rigidbodymotion.hh> #include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/tensor3.hh> #include <dune/gfe/tensor3.hh>
...@@ -19,10 +20,44 @@ ...@@ -19,10 +20,44 @@
//#define QUADRATIC_MEMBRANE_ENERGY //#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> template<class Basis, int dim, class field_type=double>
class CosseratEnergyLocalStiffness 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 // grid types
typedef typename Basis::GridView GridView; typedef typename Basis::GridView GridView;
...@@ -113,6 +148,36 @@ public: // for testing ...@@ -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: public:
/** \brief Constructor with a set of material parameters /** \brief Constructor with a set of material parameters
...@@ -148,6 +213,11 @@ public: ...@@ -148,6 +213,11 @@ public:
RT energy (const typename Basis::LocalView& localView, RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override; 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 /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the first equation of (4.4) in Neff's paper * the first equation of (4.4) in Neff's paper
*/ */
...@@ -300,7 +370,9 @@ energy(const typename Basis::LocalView& localView, ...@@ -300,7 +370,9 @@ energy(const typename Basis::LocalView& localView,
RT energy = 0; RT energy = 0;
auto element = localView.element(); 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; typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution); LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
...@@ -415,5 +487,140 @@ energy(const typename Basis::LocalView& localView, ...@@ -415,5 +487,140 @@ energy(const typename Basis::LocalView& localView,
return energy; 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 #endif //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
dune/gfe/mixedcosseratenergy.hh deleted 100644 → 0
+ 0
− 400
View file @ f807505b
#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
src/mixed-cosserat-continuum.cc
+ 2
− 2
View file @ 39336b40
...@@ -37,7 +37,7 @@ ...@@ -37,7 +37,7 @@
#include <dune/gfe/rigidbodymotion.hh> #include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/mixedlocalgfeadolcstiffness.hh> #include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
#include <dune/gfe/mixedcosseratenergy.hh> #include <dune/gfe/cosseratenergystiffness.hh>
#include <dune/gfe/cosseratvtkwriter.hh> #include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/mixedgfeassembler.hh> #include <dune/gfe/mixedgfeassembler.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh> #include <dune/gfe/mixedriemanniantrsolver.hh>
...@@ -296,7 +296,7 @@ int main (int argc, char *argv[]) try ...@@ -296,7 +296,7 @@ int main (int argc, char *argv[]) try
} }
// Assembler using ADOL-C // Assembler using ADOL-C
MixedCosseratEnergy<decltype(compositeBasis), CosseratEnergyLocalStiffness<decltype(compositeBasis),
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters, 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
&neumannBoundary, &neumannBoundary,
neumannFunction.get()); neumannFunction.get());
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment