Skip to content
Snippets Groups Projects
Commit f807505b authored by Sander, Oliver's avatar Sander, Oliver
Browse files

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.
parent 44b0720c
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
......@@ -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()
......
......@@ -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);
}
......
......@@ -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++) {
......
......@@ -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
......
......@@ -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)
{
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment