Skip to content
Snippets Groups Projects
Commit 1f0c95db authored by Lisa Julia Nebel's avatar Lisa Julia Nebel
Browse files

Adjust film-on-substrate.cc and SurfaceCosseratEnergy for the use with a CompositeBasis 

Use GFE::LocalEnergy, GFE::LocalIntegralEnergy, GFE::SumEnergy (and delete GFE::SumCosseratEnergy).
parent 9ff83989
No related branches found
No related tags found
1 merge request!48Migrate from dune-fufem to dune-functions bases
Expand all Hide whitespace changes
Inline Side-by-side
dune/gfe/sumcosseratenergy.hh deleted 100644 → 0
+ 0
− 73
View file @ 9ff83989
#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
#define DUNE_GFE_SUMCOSSERATENERGY_HH
#include <dune/elasticity/assemblers/localfestiffness.hh>
#include <dune/gfe/localenergy.hh>
namespace Dune {
namespace GFE {
template<class Basis, class TargetSpace, class field_type=double>
class SumCosseratEnergy
: public 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;
enum {dim=GridView::dimension};
public:
/** \brief Constructor with elastic energy and cosserat energy
* \param elasticEnergy The elastic energy
* \param cosseratEnergy The cosserat energy
*/
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
SumCosseratEnergy(std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#else
SumCosseratEnergy(std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#endif
std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace>> cosseratEnergy)
: elasticEnergy_(elasticEnergy),
cosseratEnergy_(cosseratEnergy)
{}
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const
{
auto&& element = localView.element();
auto&& localFiniteElement = localView.tree().finiteElement();
std::vector<Dune::FieldVector<field_type,dim>> localConfiguration(localSolution.size());
for (int i = 0; i < localSolution.size(); i++) {
localConfiguration[i] = localSolution[i].r;
}
return elasticEnergy_->energy(element, localFiniteElement, localConfiguration) + cosseratEnergy_->energy(localView, localSolution);
}
private:
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#else
std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#endif
std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace> > cosseratEnergy_;
};
} // namespace GFE
} // namespace Dune
#endif //#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
dune/gfe/surfacecosseratenergy.hh
+ 79
− 21
View file @ 1f0c95db
#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
#define DUNE_GFE_SURFACECOSSERATENERGY_HH
#include <dune/common/indices.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
......@@ -16,14 +17,15 @@
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/vertexnormals.hh>
namespace Dune::GFE {
template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
template<class Basis, class... TargetSpaces>
class SurfaceCosseratEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
: public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>
{
using GridView = typename Basis::GridView;
using DT = typename GridView::ctype ;
using RT = typename TargetSpace::ctype;
using RT = typename Dune::GFE::LocalEnergy<Basis, TargetSpaces...>::RT ;
using Entity = typename GridView::template Codim<0>::Entity ;
using RBM0 = RealTuple<RT,GridView::dimensionworld> ;
using RBM1 = Rotation<RT,GridView::dimensionworld> ;
......@@ -94,13 +96,18 @@ public:
}
RT energy(const typename Basis::LocalView& localView,
const std::vector<RigidBodyMotion<field_type,dimWorld> >& localSolution) const
{
const std::vector<TargetSpaces>&... localSolutions) const
{
static_assert(sizeof...(TargetSpaces) == 2, "SurfaceCosseratEnergy needs exactly two TargetSpace!");
using namespace Dune::Indices;
using TargetSpace0 = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
const std::vector<TargetSpace0>& localSolution0 = std::get<0>(std::forward_as_tuple(localSolutions...));
using TargetSpace1 = typename std::tuple_element<1, std::tuple<TargetSpaces...> >::type;
const std::vector<TargetSpace1>& localSolution1 = std::get<1>(std::forward_as_tuple(localSolutions...));
// 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();
////////////////////////////////////////////////////////////////////////////////////
......@@ -121,8 +128,34 @@ RT energy(const typename Basis::LocalView& localView,
////////////////////////////////////////////////////////////////////////////////////
// Set up the local nonlinear finite element function
////////////////////////////////////////////////////////////////////////////////////
typedef LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
using namespace Dune::TypeTree::Indices;
// The set of shape functions on this element
const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
// to construt a local GFE function, in case they are the shape functions are the same, we can use use one GFE-Function
#if MIXED_SPACE
std::vector<RBM0> localSolutionRBM0(localSolution0.size());
std::vector<RBM1> localSolutionRBM1(localSolution1.size());
for (int i = 0; i < localSolution0.size(); i++)
localSolutionRBM0[i] = localSolution0[i];
for (int i = 0; i< localSolution1.size(); i++)
localSolutionRBM1[i] = localSolution1[i];
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM0> LocalGFEFunctionType0;
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), RBM1> LocalGFEFunctionType1;
LocalGFEFunctionType0 localGeodesicFEFunction0(deformationLocalFiniteElement,localSolutionRBM0);
LocalGFEFunctionType1 localGeodesicFEFunction1(orientationLocalFiniteElement,localSolutionRBM1);
#else
std::vector<RBM> localSolutionRBM(localSolution0.size());
for (int i = 0; i < localSolution0.size(); i++) {
for (int j = 0; j < dimWorld; j++)
localSolutionRBM[i].r[j] = localSolution0[i][j];
localSolutionRBM[i].q = localSolution1[i];
}
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(deformationLocalFiniteElement,localSolutionRBM);
#endif
RT energy = 0;
......@@ -134,8 +167,8 @@ RT energy(const typename Basis::LocalView& localView,
auto id = idSet.subId(it.inside(), it.indexInInside(), 1);
auto boundaryGeometry = geometriesOnShellBoundary_.at(id);
auto quadOrder = (it.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * boundaryDim;
auto quadOrder = (it.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
: deformationLocalFiniteElement.localBasis().order() * boundaryDim;
const auto& quad = Dune::QuadratureRules<DT, boundaryDim>::rule(it.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
......@@ -152,18 +185,42 @@ RT energy(const typename Basis::LocalView& localView,
const DT integrationElement = boundaryGeometry.integrationElement(quad[pt].position());
// The value of the local function
RigidBodyMotion<field_type,dimWorld> value = localGeodesicFEFunction.evaluate(quadPos);
// The derivative of the local function
auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
RBM value;
Dune::FieldMatrix<RT, RBM::embeddedDim, dimWorld> derivative;
Dune::FieldMatrix<RT, RBM::embeddedDim, boundaryDim> derivative2D;
Dune::FieldMatrix<RT,7,2> derivative2D;
for (int i = 0; i < 7; i++) {
for (int j = 0; j < 2; j++) {
derivative2D[i][j] = derivative[i][j];
#if MIXED_SPACE
// The value of the local function
RBM0 value0 = localGeodesicFEFunction0.evaluate(quadPos);
RBM1 value1 = localGeodesicFEFunction1.evaluate(quadPos);
// The derivatives of the local functions
auto derivative0 = localGeodesicFEFunction0.evaluateDerivative(quadPos,value0);
auto derivative1 = localGeodesicFEFunction1.evaluateDerivative(quadPos,value1);
// Put the value and the derivatives together from the separated values
for (int i = 0; i < RBM0::dim; i++)
value.r[i] = value0[i];
value.q = value1;
for (int j = 0; j < dimWorld; j++) {
for (int i = 0; i < RBM0::embeddedDim; i++) {
derivative[i][j] = derivative0[i][j];
if (j < boundaryDim)
derivative2D[i][j] = derivative0[i][j];
}
for (int i = 0; i < RBM1::embeddedDim; i++) {
derivative[RBM0::embeddedDim + i][j] = derivative1[i][j];
if (j < boundaryDim)
derivative2D[RBM0::embeddedDim + i][j] = derivative1[i][j];
}
}
#else
value = localGeodesicFEFunction.evaluate(quadPos);
derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
for (int i = 0; i < RBM::embeddedDim; i++)
for (int j = 0; j < boundaryDim; j++)
derivative2D[i][j] = derivative[i][j];
#endif
//////////////////////////////////////////////////////////
// The rotation and its derivative
// Note: we need it in matrix coordinates
......@@ -343,5 +400,6 @@ private:
double b1_, b2_, b3_;
};
} // namespace Dune::GFE
#endif //#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
src/film-on-substrate.cc
+ 269
− 252
View file @ 1f0c95db
This diff is collapsed.
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