From a425feca08f22718e815529b4156e131aa4f82f1 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 10 May 2016 16:18:32 +0200
Subject: [PATCH] Do not hardwire the local interpolation rule in the energy
implementations
Make them a parameter instead. This allows to switch more easily between
geodesic finite elements and projected finite elements.
---
dune/gfe/chiralskyrmionenergy.hh | 18 ++++++++----------
dune/gfe/harmonicenergystiffness.hh | 16 +++++++---------
src/harmonicmaps.cc | 11 ++++++++---
3 files changed, 23 insertions(+), 22 deletions(-)
diff --git a/dune/gfe/chiralskyrmionenergy.hh b/dune/gfe/chiralskyrmionenergy.hh
index 7faae8f1..bff58bb4 100644
--- a/dune/gfe/chiralskyrmionenergy.hh
+++ b/dune/gfe/chiralskyrmionenergy.hh
@@ -5,7 +5,6 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/localgeodesicfestiffness.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
namespace Dune {
@@ -16,7 +15,7 @@ namespace GFE {
* The energy is discussed in:
* - Christof Melcher, "Chiral skyrmions in the plane", Proc. of the Royal Society, online DOI DOI: 10.1098/rspa.2014.0394
*/
-template<class Basis, class field_type>
+template<class Basis, class LocalInterpolationRule, class field_type>
class ChiralSkyrmionEnergy
: public LocalGeodesicFEStiffness<Basis,UnitVector<field_type,3> >
{
@@ -49,9 +48,9 @@ public:
field_type kappa_;
};
-template <class Basis, class field_type>
-typename ChiralSkyrmionEnergy<Basis, field_type>::RT
-ChiralSkyrmionEnergy<Basis, field_type>::
+template <class Basis, class LocalInterpolationRule, class field_type>
+typename ChiralSkyrmionEnergy<Basis, LocalInterpolationRule, field_type>::RT
+ChiralSkyrmionEnergy<Basis, LocalInterpolationRule, field_type>::
energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localConfiguration) const
{
@@ -61,8 +60,7 @@ energy(const typename Basis::LocalView& localView,
const auto element = localView.element();
const auto& localFiniteElement = localView.tree().finiteElement();
- typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
- LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localConfiguration);
+ LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
int quadOrder = (element.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
: localFiniteElement.localBasis().order() * 2 * gridDim;
@@ -84,13 +82,13 @@ energy(const typename Basis::LocalView& localView,
double weight = quad[pt].weight() * integrationElement;
// The value of the function
- auto value = localGeodesicFEFunction.evaluate(quadPos);
+ auto value = localInterpolationRule.evaluate(quadPos);
// The derivative of the local function defined on the reference element
- typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+ typename LocalInterpolationRule::DerivativeType referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos,value);
// The derivative of the function defined on the actual element
- typename LocalGFEFunctionType::DerivativeType derivative(0);
+ typename LocalInterpolationRule::DerivativeType derivative(0);
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index cc287bc9..b39d882f 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -5,9 +5,8 @@
#include <dune/geometry/quadraturerules.hh>
#include "localgeodesicfestiffness.hh"
-#include "localgeodesicfefunction.hh"
-template<class Basis, class TargetSpace>
+template<class Basis, class LocalInterpolationRule, class TargetSpace>
class HarmonicEnergyLocalStiffness
: public LocalGeodesicFEStiffness<Basis,TargetSpace>
{
@@ -27,17 +26,16 @@ public:
};
-template <class Basis, class TargetSpace>
-typename HarmonicEnergyLocalStiffness<Basis, TargetSpace>::RT
-HarmonicEnergyLocalStiffness<Basis, TargetSpace>::
+template <class Basis, class LocalInterpolationRule, class TargetSpace>
+typename HarmonicEnergyLocalStiffness<Basis, LocalInterpolationRule, TargetSpace>::RT
+HarmonicEnergyLocalStiffness<Basis, LocalInterpolationRule, TargetSpace>::
energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const
{
RT energy = 0;
const auto& localFiniteElement = localView.tree().finiteElement();
- typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
- LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+ LocalInterpolationRule localInterpolationRule(localFiniteElement,localSolution);
int quadOrder = (localFiniteElement.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
: localFiniteElement.localBasis().order() * 2 * gridDim;
@@ -58,10 +56,10 @@ energy(const typename Basis::LocalView& localView,
double weight = quad[pt].weight() * integrationElement;
// The derivative of the local function defined on the reference element
- auto referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
+ auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos);
// The derivative of the function defined on the actual element
- typename LocalGFEFunctionType::DerivativeType derivative(0);
+ typename LocalInterpolationRule::DerivativeType derivative(0);
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index ce16c9a7..eb306a6e 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -35,6 +35,8 @@
#include <dune/gfe/rotation.hh>
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/harmonicenergystiffness.hh>
#include <dune/gfe/chiralskyrmionenergy.hh>
@@ -189,20 +191,23 @@ int main (int argc, char *argv[]) try
// Create an assembler for the Harmonic Energy Functional
// ////////////////////////////////////////////////////////////
- // Assembler using ADOL-C
typedef TargetSpace::rebind<adouble>::other ATargetSpace;
+ typedef LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace> LocalInterpolationRule;
+ //typedef GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace> LocalInterpolationRule;
+
+ // Assembler using ADOL-C
std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > localEnergy;
std::string energy = parameterSet.get<std::string>("energy");
if (energy == "harmonic")
{
- localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, ATargetSpace>);
+ localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, LocalInterpolationRule, ATargetSpace>);
} else if (energy == "chiral_skyrmion")
{
- localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, adouble>(parameterSet.sub("energyParameters")));
+ localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, LocalInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
} else
DUNE_THROW(Exception, "Unknown energy type '" << energy << "'");
--
GitLab