diff --git a/dune/gfe/localintegralenergy.hh b/dune/gfe/localintegralenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e6d85ef1650633e8fadc113c135c6dc06bebf1f2
--- /dev/null
+++ b/dune/gfe/localintegralenergy.hh
@@ -0,0 +1,105 @@
+#ifndef DUNE_GFE_LOCALINTEGRALENERGY_HH
+#define DUNE_GFE_LOCALINTEGRALENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/localenergy.hh>
+#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::GFE {
+
+ /**
+ \brief Assembles the elastic energy for a single element integrating the localdensity over one element.
+
+ This class works similarly to the class Dune::Elasticity::LocalIntegralEnergy, where Dune::Elasticity::LocalIntegralEnergy extends
+ Dune::Elasticity::LocalEnergy and Dune::GFE::LocalIntegralEnergy extends Dune::GFE::LocalEnergy.
+ */
+template<class Basis, class... TargetSpaces>
+class LocalIntegralEnergy
+ : public Dune::GFE::LocalEnergy<Basis,TargetSpaces...>
+{
+ using LocalView = typename Basis::LocalView;
+ using GridView = typename LocalView::GridView;
+ using DT = typename GridView::Grid::ctype;
+ using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+
+ enum {gridDim=GridView::dimension};
+
+public:
+
+ /** \brief Constructor with a local energy density
+ */
+ LocalIntegralEnergy(const std::shared_ptr<Dune::Elasticity::LocalDensity<gridDim,RT,DT>>& ld)
+ : localDensity_(ld)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ RT energy(const typename Basis::LocalView& localView,
+ const std::vector<TargetSpaces>&... localSolutions) const
+ {
+ static_assert(sizeof...(TargetSpaces) > 0, "LocalIntegralEnergy needs at least one TargetSpace!");
+
+ using TargetSpace = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
+ static_assert( (std::is_same<TargetSpace, RigidBodyMotion<RT,gridDim>>::value)
+ or (std::is_same<TargetSpace, RealTuple<RT,gridDim>>::value), "The first TargetSpace of LocalIntegralEnergy needs to be RigidBodyMotion or RealTuple!" );
+
+ const std::vector<TargetSpace>& localSolution = std::get<0>(std::forward_as_tuple(localSolutions...));
+
+ using namespace Dune::Indices;
+ // powerBasis: grab the finite element of the first child
+ const auto& localFiniteElement = localView.tree().child(_0,0).finiteElement();
+ const auto& element = localView.element();
+
+ RT energy = 0;
+
+ // store gradients of shape functions and base functions
+ std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+ std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+
+ int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+ : localFiniteElement.localBasis().order() * gridDim;
+
+ const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+ for (const auto& qp : quad)
+ {
+ const DT integrationElement = element.geometry().integrationElement(qp.position());
+
+ const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
+
+ // Global position
+ auto x = element.geometry().global(qp.position());
+
+ // Get gradients of shape functions
+ localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
+
+ // compute gradients of base functions
+ for (size_t i=0; i<gradients.size(); ++i)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+ // Deformation gradient
+ FieldMatrix<RT,gridDim,gridDim> deformationGradient(0);
+ for (size_t i=0; i<gradients.size(); i++)
+ for (int j=0; j<gridDim; j++)
+ deformationGradient[j].axpy(localSolution[i][j], gradients[i]);
+ // Integrate the energy density
+ energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+ }
+
+ return energy;
+ }
+
+protected:
+ const std::shared_ptr<Dune::Elasticity::LocalDensity<gridDim,RT,DT>> localDensity_ = nullptr;
+
+};
+
+} // namespace Dune::GFE
+
+#endif //#ifndef DUNE_GFE_LOCALINTEGRALENERGY_HH
diff --git a/dune/gfe/neumannenergy.hh b/dune/gfe/neumannenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dd8f61b1c373c7958ee926244259a782ad5156d8
--- /dev/null
+++ b/dune/gfe/neumannenergy.hh
@@ -0,0 +1,103 @@
+#ifndef DUNE_GFE_NEUMANNENERGY_HH
+#define DUNE_GFE_NEUMANNENERGY_HH
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/elasticity/assemblers/localenergy.hh>
+
+namespace Dune::GFE {
+ /**
+ \brief Assembles the Neumann energy for a single element on the Neumann Boundary using the Neumann Function.
+
+ This class works similarly to the class Dune::Elasticity::NeumannEnergy, where Dune::Elasticity::NeumannEnergy extends
+ Dune::Elasticity::LocalEnergy and Dune::GFE::NeumannEnergy extends Dune::GFE::LocalEnergy.
+ */
+template<class Basis, class... TargetSpaces>
+class NeumannEnergy
+ : public Dune::GFE::LocalEnergy<Basis,TargetSpaces...>
+{
+ using LocalView = typename Basis::LocalView;
+ using GridView = typename LocalView::GridView;
+ using DT = typename GridView::Grid::ctype;
+ using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
+
+ enum {dim=GridView::dimension};
+
+public:
+
+ /** \brief Constructor with a set of material parameters
+ * \param parameters The material parameters
+ */
+ NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView>>& neumannBoundary,
+ std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<DT,dim>)> neumannFunction)
+ : neumannBoundary_(neumannBoundary),
+ neumannFunction_(neumannFunction)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ RT energy(const typename Basis::LocalView& localView,
+ const std::vector<TargetSpaces>&... localSolutions) const
+ {
+ static_assert(sizeof...(TargetSpaces) > 0, "NeumannEnergy needs at least one TargetSpace!");
+
+ using namespace Dune::Indices;
+ using TargetSpace = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
+ const std::vector<TargetSpace>& localSolution = std::get<0>(std::forward_as_tuple(localSolutions...));
+
+ const auto& localFiniteElement = localView.tree().child(_0,0).finiteElement();
+ const auto& element = localView.element();
+
+ RT energy = 0;
+
+ for (auto&& intersection : intersections(neumannBoundary_->gridView(), element)) {
+
+ if (not neumannBoundary_ or not neumannBoundary_->contains(intersection))
+ continue;
+
+ int quadOrder = localFiniteElement.localBasis().order();
+
+ const auto& quad = Dune::QuadratureRules<DT, dim-1>::rule(intersection.type(), quadOrder);
+
+ for (size_t pt=0; pt<quad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<DT,dim>& quadPos = intersection.geometryInInside().global(quad[pt].position());
+
+ const auto integrationElement = intersection.geometry().integrationElement(quad[pt].position());
+
+ // The value of the local function
+ std::vector<Dune::FieldVector<DT,1> > shapeFunctionValues;
+ localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+
+ Dune::FieldVector<RT,dim> value(0);
+ for (size_t i=0; i<localFiniteElement.size(); i++)
+ for (int j=0; j<dim; j++)
+ value[j] += shapeFunctionValues[i] * localSolution[i][j];
+
+ // Value of the Neumann data at the current position
+ auto neumannValue = neumannFunction_( intersection.geometry().global(quad[pt].position()) );
+
+ // Only translational dofs are affected by the Neumann force
+ for (size_t i=0; i<neumannValue.size(); i++)
+ energy += (neumannValue[i] * value[i]) * quad[pt].weight() * integrationElement;
+
+ }
+
+ }
+
+ return energy;
+ }
+
+private:
+ /** \brief The Neumann boundary */
+ const std::shared_ptr<BoundaryPatch<GridView>> neumannBoundary_;
+
+ /** \brief The function implementing the Neumann data */
+ std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<DT,dim>)> neumannFunction_;
+};
+} // namespace Dune::GFE
+
+#endif //#ifndef DUNE_GFE_NEUMANNENERGY_HH
diff --git a/dune/gfe/sumcosseratenergy.hh b/dune/gfe/sumcosseratenergy.hh
deleted file mode 100644
index 51f6241a2090e51977d9fc79f0f426bf2e49911f..0000000000000000000000000000000000000000
--- a/dune/gfe/sumcosseratenergy.hh
+++ /dev/null
@@ -1,73 +0,0 @@
-#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
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index b7b5d1226c2e8279c7c44de3bc62be586cd8a63e..366e233747f0a98c2e986b40fb8c606be0963e4f 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -1,6 +1,7 @@
#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,46 +17,21 @@
#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();
- }
-};
+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...>
{
- // 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 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> ;
+ using RBM = RigidBodyMotion<RT,GridView::dimensionworld> ;
+
+ enum {dimWorld=GridView::dimensionworld};
enum {gridDim=GridView::dimension};
static constexpr int boundaryDim = gridDim - 1;
@@ -64,9 +40,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 +53,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 +72,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,36 +96,66 @@ public:
}
RT energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<field_type,dim> >& 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();
////////////////////////////////////////////////////////////////////////////////////
// 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;
- 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;
@@ -161,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++) {
@@ -179,47 +185,71 @@ 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);
-
- // 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
//////////////////////////////////////////////////////////
- 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 +301,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 +348,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 +376,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_;
@@ -370,5 +400,6 @@ private:
double b1_, b2_, b3_;
};
+} // namespace Dune::GFE
#endif //#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index 9bc4ce139b2090e48002f5bdbb4f4e8ddcb85e0c..fb745548e806d2deb25e33f2e07af4a05c2802c6 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -51,7 +51,7 @@ initialDeformation = "x"
#############################################
# Inner solver, cholmod or multigrid
-solvertype = cholmod
+solvertype = multigrid
# Number of homotopy steps for the Dirichlet boundary conditions
numHomotopySteps = 1
@@ -148,4 +148,4 @@ mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute th
# log: Generalized Rivlin model or polynomial hyperelastic model, using 0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
# square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
-[]
\ No newline at end of file
+[]
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index f0843d81f4e46dca10975d882c9b15ff7e45390a..94c64915ad8ecb0540c32699bf78fd3a569a0f98 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -1,3 +1,4 @@
+#define MIXED_SPACE 0
#include <iostream>
#include <fstream>
@@ -12,59 +13,52 @@
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/bitsetvector.hh>
+#include <dune/common/indices.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/timer.hh>
+#include <dune/common/tuplevector.hh>
#include <dune/common/version.hh>
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-
-#include <dune/grid/io/file/gmshreader.hh>
-#include <dune/grid/io/file/vtk.hh>
-
-#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
#include <dune/elasticity/materials/exphenckydensity.hh>
#include <dune/elasticity/materials/henckydensity.hh>
#include <dune/elasticity/materials/mooneyrivlindensity.hh>
#include <dune/elasticity/materials/neohookedensity.hh>
#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
-#include <dune/elasticity/materials/sumenergy.hh>
-#include <dune/elasticity/materials/localintegralenergy.hh>
-#include <dune/elasticity/materials/neumannenergy.hh>
-#else
-#include <dune/elasticity/materials/exphenckyenergy.hh>
-#include <dune/elasticity/materials/henckyenergy.hh>
-#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
-#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
-#endif
-#include <dune/elasticity/materials/neohookeenergy.hh>
-#include <dune/elasticity/materials/neumannenergy.hh>
-#include <dune/elasticity/materials/sumenergy.hh>
-#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
-#endif
-
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
-#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/dunepython.hh>
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
#include <dune/gfe/cosseratvtkwriter.hh>
-#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/localintegralenergy.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/neumannenergy.hh>
+#include <dune/gfe/surfacecosseratenergy.hh>
+#include <dune/gfe/sumenergy.hh>
+#include <dune/gfe/vertexnormals.hh>
+
+#if MIXED_SPACE
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+#else
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
#include <dune/gfe/riemannianpnsolver.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/vertexnormals.hh>
-#include <dune/gfe/surfacecosseratenergy.hh>
-#include <dune/gfe/sumcosseratenergy.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+#endif
+
+#include <dune/istl/multitypeblockvector.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
@@ -74,7 +68,15 @@
# define WORLD_DIM 3
#endif
const int dim = WORLD_DIM;
-const int order = 2;
+
+const int targetDim = WORLD_DIM;
+
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+#if !MIXED_SPACE
+static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
+#endif
#if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
template<>
@@ -100,8 +102,14 @@ int main (int argc, char *argv[]) try
// initialize MPI, finalize is done automatically on exit
Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
- if (mpiHelper.rank()==0)
- std::cout << "ORDER = " << order << std::endl;
+ if (mpiHelper.rank()==0) {
+ std::cout << "DISPLACEMENTORDER = " << displacementOrder << ", ROTATIONORDER = " << rotationOrder << std::endl;
+#if MIXED_SPACE
+ std::cout << "MIXED_SPACE = 1" << std::endl;
+#else
+ std::cout << "MIXED_SPACE = 0" << std::endl;
+#endif
+ }
// Start Python interpreter
Python::start();
@@ -115,12 +123,6 @@ int main (int argc, char *argv[]) try
<< std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
<< std::endl;
- typedef BlockVector<FieldVector<double,dim> > SolutionType;
- typedef RigidBodyMotion<double,dim> TargetSpace;
- typedef std::vector<TargetSpace> SolutionTypeCosserat;
-
- const int blocksize = TargetSpace::TangentVector::dimension;
-
// parse data file
ParameterTree parameterSet;
@@ -146,15 +148,16 @@ int main (int argc, char *argv[]) try
const bool startFromFile = parameterSet.get<bool>("startFromFile");
const std::string resultPath = parameterSet.get("resultPath", "");
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
+ /////////////////////////////////////////////////////////////
+ // CREATE THE GRID
+ /////////////////////////////////////////////////////////////
typedef UGGrid<dim> GridType;
std::shared_ptr<GridType> grid;
FieldVector<double,dim> lower(0), upper(1);
+
if (parameterSet.get<bool>("structuredGrid")) {
lower = parameterSet.get<FieldVector<double,dim> >("lower");
@@ -206,17 +209,49 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- // FE basis spanning the FE space that we are working in
- typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
- FEBasis feBasis(gridView);
+ /////////////////////////////////////////////////////////////
+ // DATA TYPES
+ /////////////////////////////////////////////////////////////
+
+ using namespace Dune::Indices;
- // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
- typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
- FufemFEBasis fufemFEBasis(feBasis);
+ typedef std::vector<RealTuple<double,dim> > DisplacementVector;
+ typedef std::vector<Rotation<double,dim>> RotationVector;
+ const int dimRotation = Rotation<double,dim>::embeddedDim;
+ typedef TupleVector<DisplacementVector, RotationVector> SolutionType;
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
+ /////////////////////////////////////////////////////////////
+ // FUNCTION SPACE
+ /////////////////////////////////////////////////////////////
+
+ using namespace Dune::Functions::BasisFactory;
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ power<dim>(
+ lagrange<displacementOrder>()
+ ),
+ power<dimRotation>(
+ lagrange<rotationOrder>()
+ )
+ ));
+
+ auto deformationPowerBasis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<displacementOrder>()
+ ));
+ typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
+ DeformationFEBasis deformationFEBasis(gridView);
+ OrientationFEBasis orientationFEBasis(gridView);
+
+ using CompositeBasis = decltype(compositeBasis);
+ using LocalView = typename CompositeBasis::LocalView;
+
+ /////////////////////////////////////////////////////////////
+ // BOUNDARY DATA
+ /////////////////////////////////////////////////////////////
BitSetVector<1> dirichletVertices(gridView.size(dim), false);
BitSetVector<1> neumannVertices(gridView.size(dim), false);
@@ -240,97 +275,64 @@ int main (int argc, char *argv[]) try
auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
+ std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces() << " faces\n";
std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
std::cout << "On rank " << mpiHelper.rank() << ": Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
-
- BitSetVector<1> dirichletNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
- FufemFEBasis fufemFeBasis(feBasis);
- constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
-#else
- constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
-#endif
-
- BitSetVector<1> neumannNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- constructBoundaryDofs(*neumannBoundary,fufemFeBasis,neumannNodes);
-#else
- constructBoundaryDofs(*neumannBoundary,feBasis,neumannNodes);
-#endif
-
- BitSetVector<1> surfaceShellNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- constructBoundaryDofs(surfaceShellBoundary,fufemFeBasis,surfaceShellNodes);
-#else
- constructBoundaryDofs(surfaceShellBoundary,feBasis,surfaceShellNodes);
-#endif
-
- BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
- for (size_t i=0; i<feBasis.size(); i++)
- for (int j=0; j<blocksize; j++) {
- if (dirichletNodes[i][0])
- dirichletDofs[i][j] = true;
+ BitSetVector<1> dirichletNodes(compositeBasis.size({0}),false);
+ BitSetVector<1> surfaceShellNodes(compositeBasis.size({1}),false);
+ constructBoundaryDofs(dirichletBoundary,deformationFEBasis,dirichletNodes);
+ constructBoundaryDofs(surfaceShellBoundary,orientationFEBasis,surfaceShellNodes);
+
+ //Create BitVector matching the tangential space
+ const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+ typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimRotationTangent>> > VectorForBit;
+ typedef Solvers::DefaultBitVector_t<VectorForBit> BitVector;
+
+ BitVector dirichletDofs;
+ dirichletDofs[_0].resize(compositeBasis.size({0}));
+ dirichletDofs[_1].resize(compositeBasis.size({1}));
+ for (size_t i = 0; i < compositeBasis.size({0}); i++){
+ for (int j = 0; j < dim; j++){
+ dirichletDofs[_0][i][j] = dirichletNodes[i][0];
}
- for (size_t i=0; i<feBasis.size(); i++)
- for (int j=3; j<blocksize; j++) {
- if (not surfaceShellNodes[i][0])
- dirichletDofs[i][j] = true;
+ }
+ for (size_t i = 0; i < compositeBasis.size({1}); i++) {
+ for (int j = 0; j < dimRotationTangent; j++){
+ dirichletDofs[_1][i][j] = (not surfaceShellNodes[i][0] or dirichletNodes[i][0]);
}
+ }
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
+ /////////////////////////////////////////////////////////////
+ // INITIAL DATA
+ /////////////////////////////////////////////////////////////
- SolutionTypeCosserat x(feBasis.size());
+ SolutionType x;
+ x[_0].resize(compositeBasis.size({0}));
+ x[_1].resize(compositeBasis.size({1}));
- BlockVector<FieldVector<double,3> > v;
-
//Initial deformation of the underlying substrate
+ BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
- ::Functions::interpolate(fufemFEBasis, v, pythonInitialDeformation);
-
- //Copy over the initial deformation
- for (size_t i=0; i<x.size(); i++) {
- x[i].r = v[i];
- }
-
- ////////////////////////////////////////////////////////
- // Main homotopy loop
- ////////////////////////////////////////////////////////
-
- using namespace Dune::Functions::BasisFactory;
- auto powerBasis = makeBasis(
- gridView,
- power<dim>(
- lagrange<order>(),
- blockedInterleaved()
- ));
+ Dune::Functions::interpolate(deformationPowerBasis, displacement, pythonInitialDeformation);
- BlockVector<FieldVector<double,dim>> identity;
- Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+ BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
+ Dune::Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,dim> x){ return x; });
- BlockVector<FieldVector<double,dim> > displacement(v.size());
- for (int i = 0; i < v.size(); i++) {
- displacement[i] = v[i];
+ for (int i = 0; i < displacement.size(); i++) {
+ x[_0][i] = displacement[i]; //Copy over the initial deformation to the deformation part of x
+ displacement[i] -= identity[i]; //Subtract identity to get the initial displacement as a function
}
- displacement -= identity;
-
- auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
- auto localDisplacementFunction = localFunction(displacementFunction);
-
- FieldVector<double,dim> neumannValues {0,0,0};
-
- if (parameterSet.hasKey("neumannValues"))
- neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
- std::cout << "Neumann values: " << neumannValues << std::endl;
-
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
- vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_0");
+
+ /////////////////////////////////////////////////////////////
+ // INITIAL SURFACE SHELL DATA
+ /////////////////////////////////////////////////////////////
typedef MultiLinearGeometry<double, dim-1, dim> ML;
std::unordered_map<GridType::GlobalIdSet::IdType, ML> geometriesOnShellBoundary;
@@ -339,6 +341,7 @@ int main (int argc, char *argv[]) try
// Read in the grid deformation
if (startFromFile) {
+ // Create a basis of order 1 in order to deform the geometries on the surface shell boundary
const std::string pathToGridDeformationFile = parameterSet.get("pathToGridDeformationFile", "");
// for this, we create a basis of order 1 in order to deform the geometries on the surface shell boundary
auto feBasisOrder1 = makeBasis(
@@ -427,15 +430,30 @@ int main (int argc, char *argv[]) try
}
}
+ /////////////////////////////////////////////////////////////
+ // MAIN HOMOTOPY LOOP
+ /////////////////////////////////////////////////////////////
+ FieldVector<double,dim> neumannValues {0,0,0};
+ if (parameterSet.hasKey("neumannValues"))
+ neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
+ std::cout << "Neumann values: " << neumannValues << std::endl;
+
+ // Vertex Normals for the 3D-Part
+ std::vector<UnitVector<double,dim> > vertexNormals(gridView.size(dim));
+ Dune::FieldVector<double,dim> vertexNormalRaw = {0,0,1};
+ for (int i = 0; i< vertexNormals.size(); i++) {
+ UnitVector vertexNormal(vertexNormalRaw);
+ vertexNormals[i] = vertexNormal;
+ }
+ //Function for the Cosserat material parameters
const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
-
Python::Reference pythonObject = surfaceShellCallable();
PythonFunction<Dune::FieldVector<double, dim>, double> fThickness(pythonObject.get("thickness"));
PythonFunction<Dune::FieldVector<double, dim>, Dune::FieldVector<double, 2>> fLame(pythonObject.get("lame"));
- for (int i=0; i<numHomotopySteps; i++)
+ for (int i = 0; i < numHomotopySteps; i++)
{
double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
@@ -444,17 +462,13 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////////////////////////////////
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
- std::shared_ptr<NeumannFunction> neumannFunction;
- neumannFunction = std::make_shared<NeumannFunction>(neumannValues, homotopyParameter);
-#else
// A constant vector-valued function, for simple Neumann boundary values
std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>> neumannFunctionPtr;
neumannFunctionPtr = std::make_shared<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>>([&](FieldVector<double,dim> ) {
return neumannValues * (-homotopyParameter);
});
-#endif
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
if (mpiHelper.rank() == 0)
{
std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
@@ -462,9 +476,12 @@ int main (int argc, char *argv[]) try
materialParameters.report();
}
- // Assembler using ADOL-C
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2,8)
+
+ // /////////////////////////////////////////////////
+ // Create the energy functional
+ // /////////////////////////////////////////////////
+
std::shared_ptr<Elasticity::LocalDensity<dim,ValueType>> elasticDensity;
if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(materialParameters);
@@ -477,123 +494,96 @@ int main (int argc, char *argv[]) try
if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(materialParameters);
- if(!elasticDensity)
- DUNE_THROW(Exception, "Error: Selected energy not available!");
-
- auto elasticEnergy = std::make_shared<LocalIntegralEnergy<GridView, FEBasis::LocalView::Tree::FiniteElement, ValueType>>(elasticDensity);
- auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(neumannBoundary,neumannFunctionPtr);
- auto elasticAndNeumann = std::make_shared<Dune::SumEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType>>(elasticEnergy, neumannEnergy);
-#else
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
- std::shared_ptr<LocalFEStiffness<GridView,
-#else
- std::shared_ptr<Elasticity::LocalEnergy<GridView,
-#endif
- FEBasis::LocalView::Tree::FiniteElement,
- std::vector<Dune::FieldVector<ValueType, dim>> > > elasticEnergy;
-
- if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
- elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
- if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
- elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-#endif
-
- if (parameterSet.get<std::string>("energy") == "neohooke")
- elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "hencky")
- elasticEnergy = std::make_shared<HenckyEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "exphencky")
- elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-
- if(!elasticEnergy)
+ if(!elasticDensity)
DUNE_THROW(Exception, "Error: Selected energy not available!");
- auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(neumannBoundary.get(),neumannFunction.get());
-
- auto elasticAndNeumann = std::make_shared<SumEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType>>(elasticEnergy, neumannEnergy);
-#endif
-
- using LocalEnergyBase = GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble, dim> >;
-
- std::shared_ptr<LocalEnergyBase> surfaceCosseratEnergy;
- std::vector<UnitVector<double,3> > vertexNormals(gridView.size(3));
- Dune::FieldVector<double,3> vertexNormalRaw = {0,0,1};
- for (int i = 0; i< vertexNormals.size(); i++) {
- UnitVector vertexNormal(vertexNormalRaw);
- vertexNormals[i] = vertexNormal;
- }
-
- surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
-
- std::shared_ptr<LocalEnergyBase> totalEnergy;
- totalEnergy = std::make_shared<GFE::SumCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble>> (elasticAndNeumann, surfaceCosseratEnergy);
-
-
- LocalGeodesicFEADOLCStiffness<FEBasis,
- TargetSpace> localGFEADOLCStiffness(totalEnergy.get());
+ auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(elasticDensity);
+ auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(neumannBoundary,*neumannFunctionPtr);
+ auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> >>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
- GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
+ GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim>> sumEnergy;
+ sumEnergy.addLocalEnergy(neumannEnergy);
+ sumEnergy.addLocalEnergy(elasticEnergy);
+ sumEnergy.addLocalEnergy(surfaceCosseratEnergy);
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > localGFEADOLCStiffness(&sumEnergy);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
////////////////////////////////////////////////////////
// Set Dirichlet values
////////////////////////////////////////////////////////
- Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("dirichletValues"));
+ BitSetVector<RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
+ for(int i = 0; i < dirichletDofs[_0].size(); i++)
+ for (int j = 0; j < RealTuple<double,dim>::TangentVector::dimension; j++)
+ deformationDirichletDofs[i][j] = dirichletDofs[_0][i][j];
+ BitSetVector<Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
+ for(int i = 0; i < dirichletDofs[_1].size(); i++)
+ for (int j = 0; j < Rotation<double,dim>::TangentVector::dimension; j++)
+ orientationDirichletDofs[i][j] = dirichletDofs[_1][i][j];
+ Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("dirichletValues"));
Python::Callable C = dirichletValuesClass.get("DirichletValues");
// Call a constructor.
Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
// Extract object member functions as Dune functions
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
- PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
-
- std::vector<FieldVector<double,3> > ddV;
- ::Functions::interpolate(fufemFEBasis, ddV, deformationDirichletValues, dirichletDofs);
-
- std::vector<FieldMatrix<double,3,3> > dOV;
- ::Functions::interpolate(fufemFEBasis, dOV, orientationDirichletValues, dirichletDofs);
-
- for (size_t j=0; j<x.size(); j++)
- if (dirichletNodes[j][0])
- {
- x[j].r = ddV[j];
- x[j].q.set(dOV[j]);
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,targetDim> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
+ PythonFunction<FieldVector<double,dim>, FieldMatrix<double,targetDim,targetDim> > rotationalDirichletValues(dirichletValuesPythonObject.get("orientation"));
+
+ BlockVector<FieldVector<double,targetDim> > ddV;
+ Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+
+ BlockVector<FieldMatrix<double,targetDim,targetDim> > dOV;
+ Dune::Functions::interpolate(orientationFEBasis, dOV, rotationalDirichletValues, orientationDirichletDofs);
+
+ for (int i = 0; i < compositeBasis.size({0}); i++)
+ if (dirichletDofs[_0][i][0])
+ x[_0][i] = ddV[i];
+ for (int i = 0; i < compositeBasis.size({1}); i++)
+ if (dirichletDofs[_1][i][0])
+ x[_1][i].set(dOV[i]);
+
+#if !MIXED_SPACE
+ //The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
+ //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
+ //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
+ typedef RigidBodyMotion<double, dim> RBM;
+ std::vector<RBM> xRBM(compositeBasis.size({0}));
+ BitSetVector<RBM::TangentVector::dimension> dirichletDofsRBM(compositeBasis.size({0}), false);
+ for (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < dim; j ++) { // Displacement part
+ xRBM[i].r[j] = x[_0][i][j];
+ dirichletDofsRBM[i][j] = dirichletDofs[_0][i][j];
}
- // /////////////////////////////////////////////////
- // Create the solver and solve
- // /////////////////////////////////////////////////
+ xRBM[i].q = x[_1][i]; // Rotation part
+ for (int j = dim; j < RBM::TangentVector::dimension; j ++)
+ dirichletDofsRBM[i][j] = dirichletDofs[_1][i][j-dim];
+ }
+ typedef Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>> GFEAssemblerWrapper;
+ GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+#endif
+
+ ///////////////////////////////////////////////////
+ // Create the chosen solver and solve!
+ ///////////////////////////////////////////////////
if (parameterSet.get<std::string>("solvertype") == "multigrid") {
- RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+#if MIXED_SPACE
+ MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>> solver;
solver.setup(*grid,
- &assembler,
+ &mixedAssembler,
+ deformationFEBasis,
+ orientationFEBasis,
x,
- dirichletDofs,
+ deformationDirichletDofs,
+ orientationDirichletDofs,
tolerance,
maxSolverSteps,
initialTrustRegionRadius,
@@ -607,48 +597,75 @@ int main (int argc, char *argv[]) try
solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
solver.setInitialIterate(x);
solver.solve();
-
x = solver.getSol();
+#else
+ RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+ solver.setup(*grid,
+ &assembler,
+ xRBM,
+ dirichletDofsRBM,
+ tolerance,
+ maxSolverSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setInitialIterate(xRBM);
+ solver.solve();
+ xRBM = solver.getSol();
+ for (int i = 0; i < xRBM.size(); i++) {
+ x[_0][i] = xRBM[i].r;
+ x[_1][i] = xRBM[i].q;
+ }
+#endif
} else { //parameterSet.get<std::string>("solvertype") == "cholmod"
- RiemannianProximalNewtonSolver<FEBasis,TargetSpace> solver;
+
+#if MIXED_SPACE
+ DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+#else
+ RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
- x,
- dirichletDofs,
+ xRBM,
+ dirichletDofsRBM,
tolerance,
maxSolverSteps,
initialRegularization,
instrumented);
- solver.setInitialIterate(x);
+ solver.setInitialIterate(xRBM);
solver.solve();
-
- x = solver.getSol();
+ xRBM = solver.getSol();
+ for (int i = 0; i < xRBM.size(); i++) {
+ x[_0][i] = xRBM[i].r;
+ x[_1][i] = xRBM[i].q;
+ }
+#endif
}
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
std::cout << "Overall calculation took " << overallTimer.elapsed() << " sec." << std::endl;
/////////////////////////////////
// Output result
/////////////////////////////////
- for (size_t i=0; i<x.size(); i++)
- v[i] = x[i].r;
// Compute the displacement
- BlockVector<FieldVector<double,dim> > displacement(v.size());
- for (int i = 0; i < v.size(); i++) {
- displacement[i] = v[i];
+ BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
+ for (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < dim; j++) {
+ displacement[i][j] = x[_0][i][j];
+ }
+ displacement[i] -= identity[i];
}
- displacement -= identity;
-
- auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
- auto localDisplacementFunction = localFunction(displacementFunction);
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
- vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_" + std::to_string(i+1));
}
} catch (Exception& e) {