diff --git a/dune/gfe/assemblers/CMakeLists.txt b/dune/gfe/assemblers/CMakeLists.txt
index 2a6d415715361df50741cecb6e5a2df1fd250545..cdfd6ade5b61b5d9181f8c0bc9de747b8fed0dcf 100644
--- a/dune/gfe/assemblers/CMakeLists.txt
+++ b/dune/gfe/assemblers/CMakeLists.txt
@@ -5,7 +5,6 @@ install(FILES
cosseratrodenergy.hh
geodesicfeassembler.hh
geodesicfeassemblerwrapper.hh
- harmonicenergy.hh
l2distancesquaredenergy.hh
localenergy.hh
localfirstordermodel.hh
diff --git a/dune/gfe/assemblers/harmonicenergy.hh b/dune/gfe/assemblers/harmonicenergy.hh
deleted file mode 100644
index 90c47a93ca7627f286ace775b9639bae253ef40e..0000000000000000000000000000000000000000
--- a/dune/gfe/assemblers/harmonicenergy.hh
+++ /dev/null
@@ -1,84 +0,0 @@
-#ifndef DUNE_GFE_HARMONICENERGY_HH
-#define DUNE_GFE_HARMONICENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/densities/harmonicdensity.hh>
-
-template<class Basis, class LocalInterpolationRule, class TargetSpace>
-class HarmonicEnergy
- : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
-{
- // grid types
- typedef typename Basis::GridView GridView;
- typedef typename GridView::ctype DT;
- typedef typename TargetSpace::ctype RT;
-
- // some other sizes
- constexpr static int gridDim = GridView::dimension;
-
-public:
-
- /** \brief Assemble the energy for a single element */
- RT energy (const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution) const override;
-
- virtual RT energy (const typename Basis::LocalView& localView,
- const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
- {
- DUNE_THROW(Dune::NotImplemented, "!");
- }
-
-};
-
-template <class Basis, class LocalInterpolationRule, class TargetSpace>
-typename HarmonicEnergy<Basis, LocalInterpolationRule, TargetSpace>::RT
-HarmonicEnergy<Basis, LocalInterpolationRule, TargetSpace>::
-energy(const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution) const
-{
- Dune::GFE::HarmonicDensity<Dune::FieldVector<DT,gridDim>,TargetSpace> density;
-
- RT energy = 0;
-
- const auto& localFiniteElement = localView.tree().finiteElement();
- LocalInterpolationRule localInterpolationRule(localFiniteElement,localSolution);
-
- int quadOrder = (localFiniteElement.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
- : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
-
- const auto element = localView.element();
-
- const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
- const auto integrationElement = element.geometry().integrationElement(quadPos);
-
- const auto jacobianInverse = element.geometry().jacobianInverse(quadPos);
-
- auto weight = quad[pt].weight() * integrationElement;
-
- // The derivative of the local function defined on the reference element
-
- // Function value at the point where we are evaluating the derivative
- // (not needed by the energy, but needed to compute the derivative)
- TargetSpace q = localInterpolationRule.evaluate(quadPos);
-
- // Compute the derivative
- auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos, q);
- auto derivative = referenceDerivative * jacobianInverse;
-
- energy += weight * density(quadPos,q,derivative);
- }
-
- return energy;
-}
-
-#endif
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index b1bee967ac5ccf0256b835a2d03d26526d52bb0a..4317850f67f38eaf5282af2e3bec7ab010a7b9ca 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -51,86 +51,125 @@ namespace Dune::GFE {
/** \brief Assemble the energy for a single element */
RT energy(const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolutions) const
+ const std::vector<TargetSpace>& localConfiguration) const
{
- DUNE_THROW(NotImplemented, "!");
+ RT energy = 0;
+
+ if constexpr (not Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
+ {
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
+
+ int quadOrder = (localFiniteElement.type().isSimplex())
+ ? (localFiniteElement.localBasis().order()-1) * 2
+ : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
+
+ const auto element = localView.element();
+
+ const auto& quad = QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
+
+ for (auto&& qp : quad)
+ {
+ // Local position of the quadrature point
+ const auto& quadPos = qp.position();
+
+ const auto integrationElement = element.geometry().integrationElement(quadPos);
+
+ const auto jacobianInverse = element.geometry().jacobianInverse(quadPos);
+
+ // Function value at the quadrature point
+ TargetSpace q = localInterpolationRule.evaluate(quadPos);
+
+ // The derivative of the finite element function at the quadrature point
+ auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos, q);
+ auto derivative = referenceDerivative * jacobianInverse;
+
+ energy += qp.weight() * integrationElement * (*localDensityGFE_)(quadPos,q,derivative);
+ }
+ }
+
+ return energy;
}
RT energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
- // TODO: Cosserat materials are hard-wired here for historical reasons.
- static_assert(TargetSpace::size() == 2, "LocalGeodesicIntegralEnergy needs two TargetSpaces!");
- using TargetSpaceDeformation = typename std::tuple_element<0, TargetSpace>::type;
- using TargetSpaceRotation = typename std::tuple_element<1, TargetSpace>::type;
+ RT energy = 0;
- const auto& element = localView.element();
+ if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
+ {
+ // TODO: Cosserat materials are hard-wired here for historical reasons.
+ static_assert(TargetSpace::size() == 2, "LocalGeodesicIntegralEnergy needs two TargetSpaces!");
+ using TargetSpaceDeformation = typename std::tuple_element<0, TargetSpace>::type;
+ using TargetSpaceRotation = typename std::tuple_element<1, TargetSpace>::type;
- using namespace Indices;
- const std::vector<TargetSpaceDeformation>& localDeformationConfiguration = coefficients[_0];
- const std::vector<TargetSpaceRotation>& localOrientationConfiguration = coefficients[_1];
+ const auto& element = localView.element();
- // composite Basis: grab the finite element of the first child
- const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
- const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
+ using namespace Indices;
+ const std::vector<TargetSpaceDeformation>& localDeformationConfiguration = coefficients[_0];
+ const std::vector<TargetSpaceRotation>& localOrientationConfiguration = coefficients[_1];
- using LocalDeformationGFEFunctionType = typename std::tuple_element<0, LocalInterpolationRule>::type;
- using LocalOrientationGFEFunctionType = typename std::tuple_element<1, LocalInterpolationRule>::type;
+ // composite Basis: grab the finite element of the first child
+ const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
+ const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
- LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
- LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localOrientationConfiguration);
+ using LocalDeformationGFEFunctionType = typename std::tuple_element<0, LocalInterpolationRule>::type;
+ using LocalOrientationGFEFunctionType = typename std::tuple_element<1, LocalInterpolationRule>::type;
- RT energy = 0;
+ LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
+ LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localOrientationConfiguration);
- int quadOrder = (element.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
+
+ int quadOrder = (element.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
: deformationLocalFiniteElement.localBasis().order() * gridDim;
- const auto& quad = QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+ const auto& quad = QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
- for (size_t pt=0; pt<quad.size(); pt++)
- {
- // Local position of the quadrature point
- const FieldVector<DT,gridDim>& quadPos = quad[pt].position();
+ for (size_t pt=0; pt<quad.size(); pt++)
+ {
+ // Local position of the quadrature point
+ const FieldVector<DT,gridDim>& quadPos = quad[pt].position();
- auto x = element.geometry().global(quadPos);
+ auto x = element.geometry().global(quadPos);
- const DT integrationElement = element.geometry().integrationElement(quadPos);
+ const DT integrationElement = element.geometry().integrationElement(quadPos);
- const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+ const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
- DT weightWithintegrationElement = quad[pt].weight() * integrationElement;
+ DT weightWithintegrationElement = quad[pt].weight() * integrationElement;
- // The value of the local deformation
- RealTuple<RT,gridDim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
+ // The value of the local deformation
+ RealTuple<RT,gridDim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
- // The derivative of the deformation defined on the reference element
- typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
+ // The derivative of the deformation defined on the reference element
+ typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
- // The derivative of the deformation defined on the actual element
- typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
- for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
+ // The derivative of the deformation defined on the actual element
+ typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
+ for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
+ jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
- // Integrate the energy density
- if (localDensityElasticity_)
- energy += weightWithintegrationElement * (*localDensityElasticity_)(x, deformationDerivative);
- else if (localDensityGFE_) {
- // The value of the local rotation
- Rotation<RT,gridDim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
+ // Integrate the energy density
+ if (localDensityElasticity_)
+ energy += weightWithintegrationElement * (*localDensityElasticity_)(x, deformationDerivative);
+ else if (localDensityGFE_) {
+ // The value of the local rotation
+ Rotation<RT,gridDim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
- // The derivative of the rotation defined on the reference element
- typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
+ // The derivative of the rotation defined on the reference element
+ typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
- // The derivative of the rotation defined on the actual element
- typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
- for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
+ // The derivative of the rotation defined on the actual element
+ typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
+ for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
+ jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
- energy += weightWithintegrationElement * (*localDensityGFE_)(x,
- deformationValue,
- deformationDerivative,
- orientationValue,
- orientationDerivative);
+ energy += weightWithintegrationElement * (*localDensityGFE_)(x,
+ deformationValue,
+ deformationDerivative,
+ orientationValue,
+ orientationDerivative);
+ }
}
}
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 005bfed4a8a2246e8ac54f018a2d57a2df3d2828..c55aed10ff833454277f082a93aa3d27f7a417f3 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -39,9 +39,10 @@
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/globalgfefunction.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
-#include <dune/gfe/assemblers/harmonicenergy.hh>
+#include <dune/gfe/assemblers/localintegralenergy.hh>
#include <dune/gfe/assemblers/l2distancesquaredenergy.hh>
#include <dune/gfe/assemblers/weightedsumenergy.hh>
+#include <dune/gfe/densities/harmonicdensity.hh>
#include <dune/gfe/spaces/unitvector.hh>
// grid dimension
@@ -208,7 +209,9 @@ int main (int argc, char *argv[]) try
std::vector<std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > > addends(2);
using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
- addends[0] = std::make_shared<HarmonicEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >();
+
+ auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity::Geometry::LocalCoordinate, ATargetSpace> >();
+ addends[0] = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(harmonicDensity);
addends[1] = l2DistanceSquaredEnergy;
std::vector<double> weights = {1.0, 1.0/(2*timeStepSize)};
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index f6e7fff752109ffcd57675da06fbe4770e5c6870..095cc200404c6568a10aa762b6ceb6df667810da 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -42,9 +42,10 @@
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
-#include <dune/gfe/assemblers/harmonicenergy.hh>
+#include <dune/gfe/assemblers/localintegralenergy.hh>
#include <dune/gfe/assemblers/chiralskyrmionenergy.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
+#include <dune/gfe/densities/harmonicdensity.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/spaces/realtuple.hh>
@@ -277,10 +278,12 @@ int main (int argc, char *argv[])
std::string energy = parameterSet.get<std::string>("energy");
if (energy == "harmonic")
{
+ auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity::Geometry::LocalCoordinate, ATargetSpace> >();
+
if (parameterSet["interpolationMethod"] == "geodesic")
- localEnergy.reset(new HarmonicEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace>);
+ localEnergy.reset(new GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace>(harmonicDensity));
else if (parameterSet["interpolationMethod"] == "projected")
- localEnergy.reset(new HarmonicEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace>);
+ localEnergy.reset(new GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace>(harmonicDensity));
else
DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 0789635fe108f492e94993c243034f5f49eea67c..9141ecd85da90b0e57eaa023a7b69f5d9144b179 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -5,7 +5,6 @@ set(TESTS
cosseratrodenergytest
cosseratrodtest
embeddedglobalgfefunctiontest
- harmonicenergytest
localgeodesicfefunctiontest
localgfetestfunctiontest
localprojectedfefunctiontest
diff --git a/test/harmonicenergytest.cc b/test/harmonicenergytest.cc
deleted file mode 100644
index bdfa68dfa1ff24b7fe1b73313999ed135d823644..0000000000000000000000000000000000000000
--- a/test/harmonicenergytest.cc
+++ /dev/null
@@ -1,130 +0,0 @@
-#include "config.h"
-
-#include <dune/grid/uggrid.hh>
-
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
-
-#include <dune/gfe/assemblers/harmonicenergy.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/spaces/unitvector.hh>
-
-#include "multiindex.hh"
-#include "valuefactory.hh"
-
-
-
-using namespace Dune;
-
-/** \brief Computes the diameter of a set */
-template <class TargetSpace>
-double diameter(const std::vector<TargetSpace>& v)
-{
- double d = 0;
- for (size_t i=0; i<v.size(); i++)
- for (size_t j=0; j<v.size(); j++)
- d = std::max(d, TargetSpace::distance(v[i],v[j]));
- return d;
-}
-
-template <class Basis, class TargetSpace>
-void testEnergy(const Basis& basis, const std::vector<TargetSpace>& coefficients)
-{
- using GridView = typename Basis::GridView;
- GridView gridView = basis.gridView();
-
- using GeodesicInterpolationRule = LocalGeodesicFEFunction<Basis::GridView::dimension, double, typename Basis::LocalView::Tree::FiniteElement, TargetSpace>;
-
- HarmonicEnergy<Basis,GeodesicInterpolationRule,TargetSpace> assembler;
- std::vector<TargetSpace> rotatedCoefficients(coefficients.size());
-
- auto localView = basis.localView();
- localView.bind(*gridView.template begin<0>());
-
- for (int i=0; i<10; i++) {
-
- Rotation<double,3> rotation(FieldVector<double,3>(1), double(i));
-
- FieldMatrix<double,3,3> matrix;
- rotation.matrix(matrix);
- for (size_t j=0; j<coefficients.size(); j++) {
- FieldVector<double,3> tmp;
- matrix.mv(coefficients[j].globalCoordinates(), tmp);
- rotatedCoefficients[j] = tmp;
- }
-
- std::cout << "energy: " << assembler.energy(localView,
- rotatedCoefficients) << std::endl;
-
- }
-
-}
-
-
-template <int domainDim, class TargetSpace>
-void test()
-{
- // ////////////////////////////////////////////////////////
- // Make a test grid consisting of a single simplex
- // ////////////////////////////////////////////////////////
-
- typedef UGGrid<domainDim> GridType;
-
- GridFactory<GridType> factory;
-
- FieldVector<double,domainDim> pos(0);
- factory.insertVertex(pos);
-
- for (int i=0; i<domainDim; i++) {
- pos = 0;
- pos[i] = 1;
- factory.insertVertex(pos);
- }
-
- std::vector<unsigned int> v(domainDim+1);
- for (int i=0; i<domainDim+1; i++)
- v[i] = i;
- factory.insertElement(GeometryTypes::simplex(domainDim), v);
-
- auto grid = std::unique_ptr<const GridType>(factory.createGrid());
- using GridView = typename GridType::LeafGridView;
- auto gridView = grid->leafGridView();
-
- // A function space basis
- using Basis = Functions::LagrangeBasis<GridView,1>;
- Basis basis(gridView);
-
- // //////////////////////////////////////////////////////////
- // Test whether the energy is invariant under isometries
- // //////////////////////////////////////////////////////////
-
- std::vector<TargetSpace> testPoints;
- ValueFactory<TargetSpace>::get(testPoints);
-
- // Set up elements of S^2
- std::vector<TargetSpace> coefficients(domainDim+1);
-
- ::MultiIndex index(domainDim+1, testPoints.size());
- int numIndices = index.cycle();
-
- for (int i=0; i<numIndices; i++, ++index) {
-
- for (int j=0; j<domainDim+1; j++)
- coefficients[j] = testPoints[index[j]];
-
- // This may be overly restrictive, but see the TODO in targetspacetrsolver.cc
- if (diameter(coefficients) > TargetSpace::convexityRadius)
- continue;
-
- testEnergy<Basis>(basis, coefficients);
-
- }
-
-}
-
-
-int main(int argc, char** argv)
-{
- MPIHelper::instance(argc, argv);
-
- test<2,UnitVector<double,3> >();
-}
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
index 6dd18b70ff65a383171eccc6688b25b73089a188..3550b42664173380e98bb79f58daeb2143569cf3 100644
--- a/test/harmonicmaptest.cc
+++ b/test/harmonicmaptest.cc
@@ -25,8 +25,9 @@
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
-#include <dune/gfe/assemblers/harmonicenergy.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
+#include <dune/gfe/assemblers/localintegralenergy.hh>
+#include <dune/gfe/densities/harmonicdensity.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/spaces/unitvector.hh>
@@ -151,8 +152,9 @@ int main (int argc, char *argv[])
//using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+ auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<GridType::Codim<0>::Entity::Geometry::LocalCoordinate, ATargetSpace> >();
//std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy = std::make_shared<HarmonicEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >();
- std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy = std::make_shared<HarmonicEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >();
+ std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >(harmonicDensity);
LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy.get());