diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index ca8e961f93048fa9424d2c8ff449bfc851ae757d..c26c4a03dae062e683e07b714985cd88a0f6b830 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -66,16 +66,35 @@ namespace Dune::GFE {
const auto geometryJacobianInverse = element.geometry().jacobianInverse(quadPos);
- // Function value at the quadrature point
- // This is always needed: Either directly, or to compute the derivative below
- TargetSpace q = localInterpolationRule.evaluate(quadPos);
-
- // The derivative of the finite element function at the quadrature point
- typename LocalInterpolationRule::DerivativeType derivative;
- if (localDensity_->dependsOnDerivative())
- derivative = localInterpolationRule.evaluateDerivative(quadPos, q) * geometryJacobianInverse;
-
- energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,q.globalCoordinates(),derivative);
+ if (localDensity_->dependsOnValue())
+ {
+ if (localDensity_->dependsOnDerivative())
+ {
+ auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos);
+ derivative = derivative * geometryJacobianInverse;
+ energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
+ }
+ else
+ {
+ auto value = localInterpolationRule.evaluate(quadPos);
+ typename LocalInterpolationRule::DerivativeType dummyDerivative;
+ energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
+ }
+ }
+ else
+ {
+ if (localDensity_->dependsOnDerivative())
+ {
+ typename TargetSpace::CoordinateType dummyValue;
+ auto derivative = localInterpolationRule.evaluateDerivative(quadPos);
+ derivative = derivative * geometryJacobianInverse;
+ energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
+ }
+ else
+ {
+ // Density does not depend on anything. That's rather pointless, but not an error.
+ }
+ }
}
}
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 6892e93604dfd22ed4872f30835899739ea1e9d3..cb36d6aa583e3a05e00eab11d88c9ad730938ac5 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -93,6 +93,10 @@ public:
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const;
+ /** \brief Evaluate the value and the derivative of the interpolation function
+ */
+ std::pair<TargetSpace,DerivativeType> evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const;
+
/** \brief Evaluate the derivative of the function value with respect to a coefficient */
void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
@@ -304,6 +308,17 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace
return result;
}
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+std::pair<TargetSpace,typename LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType>
+LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
+evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
+{
+ std::pair<TargetSpace,DerivativeType> result;
+ result.first = evaluate(local);
+ result.second = evaluateDerivative(local,result.first);
+ return result;
+}
+
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index 5eb2c03a0072afa337946018ad62b61d418d6554..5c394feec193d6670e8cf46bcf62bb39256389c2 100644
--- a/dune/gfe/localprojectedfefunction.hh
+++ b/dune/gfe/localprojectedfefunction.hh
@@ -81,10 +81,21 @@ namespace Dune {
/** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
* \param local Local coordinates in the reference element where to evaluate the derivative
* \param q Value of the local gfe function at 'local'. If you provide something wrong here the result will be wrong, too!
+ *
+ * \note This method is only usable in the conforming setting, because it requires the caller
+ * to hand over the interpolation value as a TargetSpace object.
*/
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const;
+ /** \brief Evaluate the value and the derivative of the interpolation function
+ *
+ * \return A std::pair containing the value and the first derivative of the interpolation function.
+ * If the interpolation is conforming then the first member of the pair will be a TargetSpace.
+ * Otherwise it will be a RealTuple.
+ */
+ auto evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const;
+
/** \brief Get the i'th base coefficient. */
TargetSpace coefficient(int i) const
{
@@ -175,6 +186,55 @@ namespace Dune {
return derivativeOfProjection*derivative;
}
+ template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
+ auto
+ LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
+ evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
+ {
+ // Construct the type of the result -- it depends on whether the interpolation
+ // is conforming or not.
+ using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;
+
+ std::pair<Value,DerivativeType> result;
+
+ ///////////////////////////////////////////////////////////
+ // Compute the value of the interpolation function
+ ///////////////////////////////////////////////////////////
+
+ std::vector<Dune::FieldVector<ctype,1> > w;
+ localFiniteElement_.localBasis().evaluateFunction(local,w);
+
+ typename TargetSpace::CoordinateType embeddedInterpolation(0);
+ for (size_t i=0; i<coefficients_.size(); i++)
+ embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());
+
+ if constexpr (conforming)
+ result.first = TargetSpace::projectOnto(embeddedInterpolation);
+ else
+ result.first = (RealTuple<RT, TargetSpace::CoordinateType::dimension>)embeddedInterpolation;
+
+ ///////////////////////////////////////////////////////////
+ // Compute the derivative of the interpolation function
+ ///////////////////////////////////////////////////////////
+
+ // Compute the interpolation in the surrounding space
+ std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
+ localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
+
+ result.second = 0.0;
+ for (size_t i=0; i<embeddedDim; i++)
+ for (size_t j=0; j<dim; j++)
+ for (size_t k=0; k<coefficients_.size(); k++)
+ result.second[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];
+
+ // The derivative of the projection onto the manifold
+ if constexpr(conforming)
+ result.second = TargetSpace::derivativeOfProjection(embeddedInterpolation) * result.second;
+
+ return result;
+ }
+
+
/** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for SO(3)
*
* \tparam dim Dimension of the reference element
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 59a4390e214d06b4b5eab33b91cfb686348533c1..b32a0b3432cf29a957c9b59d9ccbdd84699b8f02 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -39,8 +39,20 @@ foreach(_test ${TESTS})
target_compile_options(${_test} PRIVATE -g)
endforeach()
+dune_add_test(NAME harmonicmaptest-geodesic
+ SOURCES harmonicmaptest.cc
+ COMPILE_DEFINITIONS "GEODESICINTERPOLATION=1; CONFORMING=1"
+ TIMEOUT 1200)
+
+dune_add_test(NAME harmonicmaptest-projected
+ SOURCES harmonicmaptest.cc
+ COMPILE_DEFINITIONS "GEODESICINTERPOLATION=0; CONFORMING=1"
+ TIMEOUT 1200)
+
# Run distributed tests
-dune_add_test(SOURCES harmonicmaptest.cc
+dune_add_test(NAME harmonicmaptest-projected-nonconforming
+ SOURCES harmonicmaptest.cc
+ COMPILE_DEFINITIONS "GEODESICINTERPOLATION=0; CONFORMING=0"
MPI_RANKS 1 4
TIMEOUT 1200
CMAKE_GUARD MPI_FOUND)
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
index 8f16768a156b63f85f6c4eb302c3995733c5a72a..5dcefa4725599969614379eb77facfa756b89145 100644
--- a/test/harmonicmaptest.cc
+++ b/test/harmonicmaptest.cc
@@ -149,12 +149,20 @@ int main (int argc, char *argv[])
//////////////////////////////////////////////////////////////
using ATargetSpace = TargetSpace::rebind<adouble>::other;
- //using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
- using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+
+#if GEODESICINTERPOLATION
+ using InterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+#else
+#if CONFORMING
+ using InterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace,true>;
+#else
+ using InterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace,false>;
+#endif
+#endif
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<GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >(harmonicDensity);
+
+ std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, InterpolationRule, ATargetSpace> >(harmonicDensity);
LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy.get());
@@ -188,7 +196,16 @@ int main (int argc, char *argv[])
x = solver.getSol();
+#if GEODESICINTERPOLATION
+ std::size_t expectedFinalIteration = 8;
+#else
+#if CONFORMING
std::size_t expectedFinalIteration = 10;
+#else
+ std::size_t expectedFinalIteration = 6;
+#endif
+#endif
+
if (solver.getStatistics().finalIteration != expectedFinalIteration)
{
std::cerr << "Trust-region solver did " << solver.getStatistics().finalIteration+1
@@ -196,7 +213,15 @@ int main (int argc, char *argv[])
return 1;
}
+#if GEODESICINTERPOLATION
+ double expectedEnergy = 12.3154833;
+#else
+#if CONFORMING
double expectedEnergy = 12.2927849;
+#else
+ double expectedEnergy = 0.436857464;
+#endif
+#endif
if ( std::abs(solver.getStatistics().finalEnergy - expectedEnergy) > 1e-7)
{
std::cerr << std::setprecision(9);