From f807505b77c859c5d9da22b190348b99f94ad4d6 Mon Sep 17 00:00:00 2001 From: Oliver Sander <oliver.sander@tu-dresden.de> Date: Sun, 27 Dec 2015 22:28:37 +0100 Subject: [PATCH] Hand over a localView, rather than an element plus a localFiniteElement This is how the new dune-functions bases should be used. Also, it is yet another step towards merging CosseratEnergyLocalStiffness and MixedCosseratEnergy. --- dune/gfe/chiralskyrmionenergy.hh | 15 +++++++------- dune/gfe/cosseratenergystiffness.hh | 15 ++++++-------- dune/gfe/geodesicfeassembler.hh | 6 +++--- dune/gfe/harmonicenergystiffness.hh | 20 ++++++++---------- dune/gfe/localgeodesicfeadolcstiffness.hh | 25 ++++++++--------------- dune/gfe/localgeodesicfestiffness.hh | 16 +++++---------- 6 files changed, 39 insertions(+), 58 deletions(-) diff --git a/dune/gfe/chiralskyrmionenergy.hh b/dune/gfe/chiralskyrmionenergy.hh index 943e16ab..7faae8f1 100644 --- a/dune/gfe/chiralskyrmionenergy.hh +++ b/dune/gfe/chiralskyrmionenergy.hh @@ -23,7 +23,6 @@ class ChiralSkyrmionEnergy // various useful types typedef UnitVector<field_type,3> TargetSpace; typedef typename Basis::GridView GridView; - 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; @@ -43,9 +42,8 @@ public: enum { blocksize = TargetSpace::TangentVector::dimension }; /** \brief Assemble the energy for a single element */ - RT energy (const Entity& e, - const LocalFiniteElement& localFiniteElement, - const std::vector<TargetSpace>& localConfiguration) const; + RT energy (const typename Basis::LocalView& localView, + const std::vector<TargetSpace>& localConfiguration) const override; field_type h_; field_type kappa_; @@ -54,15 +52,16 @@ public: template <class Basis, class field_type> typename ChiralSkyrmionEnergy<Basis, field_type>::RT ChiralSkyrmionEnergy<Basis, field_type>:: -energy(const Entity& element, - const LocalFiniteElement& localFiniteElement, +energy(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localConfiguration) const { - assert(element.type() == localFiniteElement.type()); typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; RT energy = 0; - typedef LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> LocalGFEFunctionType; + + const auto element = localView.element(); + const auto& localFiniteElement = localView.tree().finiteElement(); + typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType; LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localConfiguration); int quadOrder = (element.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2 diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh index c85665a9..e51ab6c8 100644 --- a/dune/gfe/cosseratenergystiffness.hh +++ b/dune/gfe/cosseratenergystiffness.hh @@ -26,7 +26,6 @@ class CosseratEnergyLocalStiffness { // grid types typedef typename Basis::GridView GridView; - typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement; typedef typename GridView::ctype DT; typedef RigidBodyMotion<field_type,dim> TargetSpace; typedef typename TargetSpace::ctype RT; @@ -146,9 +145,8 @@ public: } /** \brief Assemble the energy for a single element */ - RT energy (const Entity& e, - const LocalFiniteElement& localFiniteElement, - const std::vector<TargetSpace>& localSolution) const; + RT energy (const typename Basis::LocalView& localView, + const std::vector<TargetSpace>& localSolution) const override; /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in * the first equation of (4.4) in Neff's paper @@ -296,15 +294,14 @@ public: template <class Basis, int dim, class field_type> typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT CosseratEnergyLocalStiffness<Basis,dim,field_type>:: -energy(const Entity& element, - const typename Basis::LocalView::Tree::FiniteElement& localFiniteElement, +energy(const typename Basis::LocalView& localView, const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const { - assert(element.type() == localFiniteElement.type()); - RT energy = 0; - typedef LocalGeodesicFEFunction<gridDim, DT, LocalFiniteElement, TargetSpace> LocalGFEFunctionType; + auto element = localView.element(); + const auto& localFiniteElement = localView.tree().finiteElement(); + typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType; LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution); int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order() diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh index 7a349f9b..ffea522f 100644 --- a/dune/gfe/geodesicfeassembler.hh +++ b/dune/gfe/geodesicfeassembler.hh @@ -149,7 +149,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol, std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct); // setup local matrix and gradient - localStiffness_->assembleGradientAndHessian(*it, localView.tree().finiteElement(), localSolution, localGradient); + localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient); // Add element matrix to global stiffness matrix for(int i=0; i<numOfBaseFct; i++) { @@ -208,7 +208,7 @@ assembleGradient(const std::vector<TargetSpace>& sol, // Assemble local gradient std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs); - localStiffness_->assembleGradient(*it, localView.tree().finiteElement(), localSolution, localGradient); + localStiffness_->assembleGradient(localView, localSolution, localGradient); // Add to global gradient for (size_t i=0; i<nDofs; i++) @@ -249,7 +249,7 @@ computeEnergy(const std::vector<TargetSpace>& sol) const for (size_t i=0; i<nDofs; i++) localSolution[i] = sol[localIndexSet.index(i)[0]]; - energy += localStiffness_->energy(*it, localView.tree().finiteElement(), localSolution); + energy += localStiffness_->energy(localView, localSolution); } diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh index 17097318..ad8f60b6 100644 --- a/dune/gfe/harmonicenergystiffness.hh +++ b/dune/gfe/harmonicenergystiffness.hh @@ -13,7 +13,6 @@ class HarmonicEnergyLocalStiffness { // grid types typedef typename Basis::GridView GridView; - 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; @@ -27,33 +26,32 @@ public: enum { blocksize = TargetSpace::TangentVector::dimension }; /** \brief Assemble the energy for a single element */ - RT energy (const Entity& e, - const LocalFiniteElement& localFiniteElement, - const std::vector<TargetSpace>& localSolution) const; + RT energy (const typename Basis::LocalView& localView, + const std::vector<TargetSpace>& localSolution) const override; }; template <class Basis, class TargetSpace> typename HarmonicEnergyLocalStiffness<Basis, TargetSpace>::RT HarmonicEnergyLocalStiffness<Basis, TargetSpace>:: -energy(const Entity& element, - const LocalFiniteElement& localFiniteElement, +energy(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution) const { - assert(element.type() == localFiniteElement.type()); typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; RT energy = 0; - typedef LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> LocalGFEFunctionType; + + const auto& localFiniteElement = localView.tree().finiteElement(); + typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType; LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution); - int quadOrder = (element.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2 + int quadOrder = (localFiniteElement.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2 : localFiniteElement.localBasis().order() * 2 * gridDim; - + const auto element = localView.element(); const Dune::QuadratureRule<double, gridDim>& quad - = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder); + = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder); for (size_t pt=0; pt<quad.size(); pt++) { diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh index 467db871..5479d276 100644 --- a/dune/gfe/localgeodesicfeadolcstiffness.hh +++ b/dune/gfe/localgeodesicfeadolcstiffness.hh @@ -24,7 +24,6 @@ class LocalGeodesicFEADOLCStiffness { // grid types typedef typename Basis::GridView GridView; - 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; @@ -47,16 +46,14 @@ public: {} /** \brief Compute the energy at the current configuration */ - virtual RT energy (const Entity& e, - const LocalFiniteElement& localFiniteElement, + virtual RT energy (const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution) const; /** \brief Assemble the element gradient of the energy functional This uses the automatic differentiation toolbox ADOL_C. */ - virtual void assembleGradient(const Entity& element, - const LocalFiniteElement& localFiniteElement, + virtual void assembleGradient(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& solution, std::vector<typename TargetSpace::TangentVector>& gradient) const; @@ -64,8 +61,7 @@ public: This uses the automatic differentiation toolbox ADOL_C. */ - virtual void assembleGradientAndHessian(const Entity& e, - const LocalFiniteElement& localFiniteElement, + virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution, std::vector<typename TargetSpace::TangentVector>& localGradient); @@ -77,8 +73,7 @@ public: template <class Basis, class TargetSpace> typename LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::RT LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>:: -energy(const Entity& element, - const LocalFiniteElement& localFiniteElement, +energy(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution) const { double pureEnergy; @@ -107,7 +102,7 @@ energy(const Entity& element, localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble } - energy = localEnergy_->energy(element,localFiniteElement,localASolution); + energy = localEnergy_->energy(localView,localASolution); energy >>= pureEnergy; @@ -126,13 +121,12 @@ energy(const Entity& element, template <class Basis, class TargetSpace> void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>:: -assembleGradient(const Entity& element, - const LocalFiniteElement& localFiniteElement, +assembleGradient(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution, std::vector<typename TargetSpace::TangentVector>& localGradient) const { // Tape energy computation. We may not have to do this every time, but it's comparatively cheap. - energy(element, localFiniteElement, localSolution); + energy(localView, localSolution); // Compute the actual gradient size_t nDofs = localSolution.size(); @@ -174,13 +168,12 @@ assembleGradient(const Entity& element, // /////////////////////////////////////////////////////////// template <class Basis, class TargetSpace> void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>:: -assembleGradientAndHessian(const Entity& element, - const LocalFiniteElement& localFiniteElement, +assembleGradientAndHessian(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution, std::vector<typename TargetSpace::TangentVector>& localGradient) { // Tape energy computation. We may not have to do this every time, but it's comparatively cheap. - energy(element, localFiniteElement, localSolution); + energy(localView, localSolution); ///////////////////////////////////////////////////////////////// // Compute the gradient. It is needed to transform the Hessian diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh index e5cc7938..878c23ff 100644 --- a/dune/gfe/localgeodesicfestiffness.hh +++ b/dune/gfe/localgeodesicfestiffness.hh @@ -10,7 +10,6 @@ class LocalGeodesicFEStiffness { // grid types typedef typename Basis::GridView GridView; - 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; @@ -39,21 +38,18 @@ public: We compute that using a finite difference approximation. */ - virtual void assembleGradientAndHessian(const Entity& e, - const LocalFiniteElement& localFiniteElement, + virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution, std::vector<typename TargetSpace::TangentVector>& localGradient); /** \brief Compute the energy at the current configuration */ - virtual RT energy (const Entity& e, - const LocalFiniteElement& localFiniteElement, + virtual RT energy (const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution) const = 0; /** \brief Assemble the element gradient of the energy functional The default implementation in this class uses a finite difference approximation */ - virtual void assembleGradient(const Entity& element, - const LocalFiniteElement& localFiniteElement, + virtual void assembleGradient(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& solution, std::vector<typename TargetSpace::TangentVector>& gradient) const; @@ -65,8 +61,7 @@ public: template <class Basis, class TargetSpace> void LocalGeodesicFEStiffness<Basis, TargetSpace>:: -assembleGradient(const Entity& element, - const LocalFiniteElement& localFiniteElement, +assembleGradient(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution, std::vector<typename TargetSpace::TangentVector>& localGradient) const { @@ -79,8 +74,7 @@ assembleGradient(const Entity& element, // /////////////////////////////////////////////////////////// template <class Basis, class TargetSpace> void LocalGeodesicFEStiffness<Basis, TargetSpace>:: -assembleGradientAndHessian(const Entity& element, - const LocalFiniteElement& localFiniteElement, +assembleGradientAndHessian(const typename Basis::LocalView& localView, const std::vector<TargetSpace>& localSolution, std::vector<typename TargetSpace::TangentVector>& localGradient) { -- GitLab