From 5d742221e0e16b84241aca0ee274cfb9376bb4ad Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Wed, 7 Oct 2020 13:29:34 +0200
Subject: [PATCH] Add localintegralenergy.hh and neumannenergy.hh, they
assemble the energy for a single element and work with a CompositeBasis
These classes work similarly to Dune::Elasticity::LocalIntegralEnergy and Dune::Elasticity::NeumannEnergy,
where the ones in Dune::Elasticity with a power basis (for the displacement function) and
Dune::GFE::LocalIntegralEnergy / Dune::GFE::NeumannEnergy work with a CompositeBasis (for the displacement AND the rotation),
so the indices to access the different parts of the basis are different.
---
dune/gfe/localintegralenergy.hh | 105 ++++++++++++++++++++++++++++++++
dune/gfe/neumannenergy.hh | 103 +++++++++++++++++++++++++++++++
2 files changed, 208 insertions(+)
create mode 100644 dune/gfe/localintegralenergy.hh
create mode 100644 dune/gfe/neumannenergy.hh
diff --git a/dune/gfe/localintegralenergy.hh b/dune/gfe/localintegralenergy.hh
new file mode 100644
index 00000000..e6d85ef1
--- /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 00000000..dd8f61b1
--- /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
--
GitLab