From 9ff839890dccb3366e42c8367340c807694d621d Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Tue, 20 Oct 2020 14:29:13 +0200
Subject: [PATCH] Refactor parts of surfacecosseratenergy.hh
---
dune/gfe/surfacecosseratenergy.hh | 137 ++++++++++++------------------
1 file changed, 55 insertions(+), 82 deletions(-)
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index b7b5d122..d6ddb215 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -16,46 +16,20 @@
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/vertexnormals.hh>
-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::LagrangeBasis<GridView,order>,i>
-{
-public:
- static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView,
- std::integral_constant<std::size_t, i> iType)
- -> decltype(localView.tree().finiteElement())
- {
- return localView.tree().finiteElement();
- }
-};
template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
class SurfaceCosseratEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
- // grid types
- typedef typename Basis::GridView GridView;
- typedef typename GridView::ctype ctype;
- typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
- typedef typename GridView::ctype DT;
- typedef typename TargetSpace::ctype RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
- typedef typename TargetSpace::template rebind<adouble>::other ATargetSpace;
-
- enum {dim=GridView::dimension};
- enum {dimworld=GridView::dimensionworld};
+ using GridView = typename Basis::GridView;
+ using DT = typename GridView::ctype ;
+ using RT = typename TargetSpace::ctype;
+ using Entity = typename GridView::template Codim<0>::Entity ;
+ using RBM0 = RealTuple<RT,GridView::dimensionworld> ;
+ using RBM1 = Rotation<RT,GridView::dimensionworld> ;
+ using RBM = RigidBodyMotion<RT,GridView::dimensionworld> ;
+
+ enum {dimWorld=GridView::dimensionworld};
enum {gridDim=GridView::dimension};
static constexpr int boundaryDim = gridDim - 1;
@@ -64,9 +38,9 @@ class SurfaceCosseratEnergy
\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 RigidBodyMotion<field_type,dimworld>& value,
- const Dune::FieldMatrix<field_type,7,boundaryDim>& derivative,
- Tensor3<field_type,3,3,boundaryDim>& DR)
+ static void computeDR(const RBM& value,
+ const Dune::FieldMatrix<RT,RBM::embeddedDim,boundaryDim>& derivative,
+ Tensor3<RT,dimWorld,dimWorld,boundaryDim>& derivative_rotation)
{
// 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
@@ -77,15 +51,15 @@ class SurfaceCosseratEnergy
// corresponding orthogonal matrix, we need to invert the i and j indices
//
// So, DR[i][j][k] contains \partial R_ij / \partial k
- Tensor3<field_type,3 , 3, 4> dd_dq;
+ Tensor3<RT, dimWorld, dimWorld, 4> dd_dq;
value.q.getFirstDerivativesOfDirectors(dd_dq);
- DR = field_type(0);
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
+ derivative_rotation = RT(0);
+ for (int i=0; i<dimWorld; i++)
+ for (int j=0; j<dimWorld; j++)
for (int k=0; k<boundaryDim; k++)
for (int l=0; l<4; l++)
- DR[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
+ derivative_rotation[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
}
@@ -96,11 +70,11 @@ public:
* \param parameters The material parameters
*/
SurfaceCosseratEnergy(const Dune::ParameterTree& parameters,
- const std::vector<UnitVector<double,3> >& vertexNormals,
+ const std::vector<UnitVector<double,dimWorld> >& vertexNormals,
const BoundaryPatch<GridView>* shellBoundary,
- const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dim-1, dim>>& geometriesOnShellBoundary,
- const std::function<double(Dune::FieldVector<double,dim>)> thicknessF,
- const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF)
+ const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dimWorld-1, dimWorld>>& geometriesOnShellBoundary,
+ const std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF,
+ const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF)
: shellBoundary_(shellBoundary),
vertexNormals_(vertexNormals),
geometriesOnShellBoundary_(geometriesOnShellBoundary),
@@ -120,35 +94,34 @@ public:
}
RT energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
+ const std::vector<RigidBodyMotion<field_type,dimWorld> >& localSolution) const
{
// The element geometry
auto element = localView.element();
// The set of shape functions on this element
const auto& localFiniteElement = localView.tree().finiteElement();
+ auto gridView = localView.globalBasis().gridView();
////////////////////////////////////////////////////////////////////////////////////
// Construct a linear (i.e., non-constant!) normal field on each element
////////////////////////////////////////////////////////////////////////////////////
- auto gridView = localView.globalBasis().gridView();
-
- assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
- std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
- for (size_t i=0; i<cornerNormals.size(); i++)
- cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,gridDim)];
-
typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;
//it.type()?
P1FiniteElementCache p1FiniteElementCache;
const auto& p1LocalFiniteElement = p1FiniteElementCache.get(element.type());
+
+ assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
+ std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
+ for (size_t i=0; i<cornerNormals.size(); i++)
+ cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,gridDim)];
Dune::GFE::LocalProjectedFEFunction<gridDim, DT, P1LocalFiniteElement, UnitVector<double,3> > unitNormals(p1LocalFiniteElement, cornerNormals);
////////////////////////////////////////////////////////////////////////////////////
// Set up the local nonlinear finite element function
////////////////////////////////////////////////////////////////////////////////////
- typedef LocalGeodesicFEFunction<gridDim, DT, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
+ typedef LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
RT energy = 0;
@@ -180,7 +153,7 @@ RT energy(const typename Basis::LocalView& localView,
const DT integrationElement = boundaryGeometry.integrationElement(quad[pt].position());
// The value of the local function
- RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ RigidBodyMotion<field_type,dimWorld> value = localGeodesicFEFunction.evaluate(quadPos);
// The derivative of the local function
auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
@@ -196,30 +169,30 @@ RT energy(const typename Basis::LocalView& localView,
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
- Dune::FieldMatrix<field_type,dim,dim> R;
+ Dune::FieldMatrix<RT,dimWorld,dimWorld> R;
value.q.matrix(R);
- auto RT = Dune::GFE::transpose(R);
+ auto rt = Dune::GFE::transpose(R);
- Tensor3<field_type,3,3,boundaryDim> DR;
- computeDR(value, derivative2D, DR);
+ Tensor3<RT,dimWorld,dimWorld,boundaryDim> derivative_rotation;
+ computeDR(value, derivative2D, derivative_rotation);
//////////////////////////////////////////////////////////
// Fundamental forms and curvature
//////////////////////////////////////////////////////////
// First fundamental form
- Dune::FieldMatrix<double,3,3> aCovariant;
+ Dune::FieldMatrix<double,dimWorld,dimWorld> aCovariant;
- // If dimworld==3, then the first two lines of aCovariant are simply the jacobianTransposed
- // of the element. If dimworld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
+ // If dimWorld==3, then the first two lines of aCovariant are simply the jacobianTransposed
+ // of the element. If dimWorld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
const auto jacobianTransposed = boundaryGeometry.jacobianTransposed(quad[pt].position());
// auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
for (int i=0; i<2; i++)
{
- for (int j=0; j<dimworld; j++)
+ for (int j=0; j<dimWorld; j++)
aCovariant[i][j] = jacobianTransposed[i][j];
- for (int j=dimworld; j<3; j++)
+ for (int j=dimWorld; j<3; j++)
aCovariant[i][j] = 0.0;
}
@@ -271,28 +244,28 @@ RT energy(const typename Basis::LocalView& localView,
//////////////////////////////////////////////////////////
// Elastic shell strain
- Dune::FieldMatrix<field_type,3,3> Ee(0);
- Dune::FieldMatrix<field_type,3,3> grad_s_m(0);
+ Dune::FieldMatrix<RT,3,3> Ee(0);
+ Dune::FieldMatrix<RT,3,3> grad_s_m(0);
for (int alpha=0; alpha<boundaryDim; alpha++)
{
- Dune::FieldVector<field_type,3> vec;
+ Dune::FieldVector<RT,3> vec;
for (int i=0; i<3; i++)
vec[i] = derivative[i][alpha];
grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
- Ee = RT * grad_s_m - a;
+ Ee = rt * grad_s_m - a;
// Elastic shell bending-curvature strain
- Dune::FieldMatrix<field_type,3,3> Ke(0);
+ Dune::FieldMatrix<RT,3,3> Ke(0);
for (int alpha=0; alpha<boundaryDim; alpha++)
{
- Dune::FieldMatrix<field_type,3,3> tmp;
+ Dune::FieldMatrix<RT,3,3> tmp;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
- tmp[i][j] = DR[i][j][alpha];
- auto tmp2 = RT * tmp;
- Ke += Dune::GFE::dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
+ tmp[i][j] = derivative_rotation[i][j][alpha];
+ auto tmp2 = rt * tmp;
+ Ke += Dune::GFE::dyadicProduct(SkewMatrix<RT,3>(tmp2).axial(), aContravariant[alpha]);
}
//////////////////////////////////////////////////////////
@@ -318,24 +291,24 @@ RT energy(const typename Basis::LocalView& localView,
return energy;
}
- RT W_m(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
+ RT W_m(const Dune::FieldMatrix<RT,3,3>& S, double mu, double lambda) const
{
return W_mixt(S,S, mu, lambda);
}
- RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T, double mu, double lambda) const
+ RT W_mixt(const Dune::FieldMatrix<RT,3,3>& S, const Dune::FieldMatrix<RT,3,3>& T, double mu, double lambda) const
{
return mu * Dune::GFE::frobeniusProduct(Dune::GFE::sym(S), Dune::GFE::sym(T))
+ mu_c_ * Dune::GFE::frobeniusProduct(Dune::GFE::skew(S), Dune::GFE::skew(T))
+ lambda * mu / (lambda + 2*mu) * Dune::GFE::trace(S) * Dune::GFE::trace(T);
}
- RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S, double mu, double lambda) const
+ RT W_mp(const Dune::FieldMatrix<RT,3,3>& S, double mu, double lambda) const
{
return mu * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda * 0.5 * Dune::GFE::traceSquared(S);
}
- RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S, double mu) const
+ RT W_curv(const Dune::FieldMatrix<RT,3,3>& S, double mu) const
{
return mu * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
+ b2_ * Dune::GFE::skew(S).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(S));
@@ -346,16 +319,16 @@ private:
const BoundaryPatch<GridView>* shellBoundary_;
/** \brief Stress-free geometries of the shell elements*/
- const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dim-1, dim>> geometriesOnShellBoundary_;
+ const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dimWorld-1, dimWorld>> geometriesOnShellBoundary_;
/** \brief The normal vectors at the grid vertices. They are used to compute the reference surface curvature. */
std::vector<UnitVector<double,3> > vertexNormals_;
/** \brief The shell thickness as a function*/
- std::function<double(Dune::FieldVector<double,dim>)> thicknessF_;
+ std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF_;
/** \brief The Lamé-parameters as a function*/
- std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dim>)> lameF_;
+ std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF_;
/** \brief Lame constants */
double mu_, lambda_;
--
GitLab