diff --git a/dune/gfe/assemblers/CMakeLists.txt b/dune/gfe/assemblers/CMakeLists.txt
index cdfd6ade5b61b5d9181f8c0bc9de747b8fed0dcf..279a37476b2fef992cef984465c698e82edf2374 100644
--- a/dune/gfe/assemblers/CMakeLists.txt
+++ b/dune/gfe/assemblers/CMakeLists.txt
@@ -1,6 +1,5 @@
#install headers
install(FILES
- chiralskyrmionenergy.hh
cosseratenergystiffness.hh
cosseratrodenergy.hh
geodesicfeassembler.hh
diff --git a/dune/gfe/assemblers/chiralskyrmionenergy.hh b/dune/gfe/assemblers/chiralskyrmionenergy.hh
deleted file mode 100644
index d88ef7e3fd1f3a388a893c626d5ffc9726934a10..0000000000000000000000000000000000000000
--- a/dune/gfe/assemblers/chiralskyrmionenergy.hh
+++ /dev/null
@@ -1,134 +0,0 @@
-#ifndef DUNE_GFE_CHIRALSKYRMIONENERGY_HH
-#define DUNE_GFE_CHIRALSKYRMIONENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/gfe/assemblers/localenergy.hh>
-
-namespace Dune {
-
- namespace GFE {
-
- /** \brief Energy of certain chiral Skyrmion
- *
- * 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 LocalInterpolationRule, class field_type>
- class ChiralSkyrmionEnergy
- : public GFE::LocalEnergy<Basis,UnitVector<field_type,3> >
- {
- // various useful types
- typedef UnitVector<field_type,3> TargetSpace;
- typedef typename Basis::GridView GridView;
- typedef typename GridView::ctype DT;
- typedef typename TargetSpace::ctype RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- constexpr static int gridDim = GridView::dimension;
-
- public:
-
- ChiralSkyrmionEnergy(const Dune::ParameterTree& parameters)
- {
- h_ = parameters.template get<double>("h");
- kappa_ = parameters.template get<double>("kappa");
- }
-
- //! Dimension of a tangent space
- constexpr static int blocksize = TargetSpace::TangentVector::dimension;
-
- /** \brief Assemble the energy for a single element */
- RT energy (const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localConfiguration) const override;
-
- RT energy (const typename Basis::LocalView& localView,
- const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
- {
- DUNE_THROW(NotImplemented, "!");
- }
-
- field_type h_;
- field_type kappa_;
- };
-
- 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
- {
- typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
- RT energy = 0;
-
- const auto element = localView.element();
- const auto& localFiniteElement = localView.tree().finiteElement();
- LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
-
- int quadOrder = (element.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
- : localFiniteElement.localBasis().order() * 2 * gridDim;
-
-
-
- const Dune::QuadratureRule<double, gridDim>& quad
- = Dune::QuadratureRules<double, gridDim>::rule(element.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 double integrationElement = element.geometry().integrationElement(quadPos);
-
- const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
- double weight = quad[pt].weight() * integrationElement;
-
- // The value of the function
- auto value = localInterpolationRule.evaluate(quadPos);
-
- // The derivative of the local function defined on the reference element
- typename LocalInterpolationRule::DerivativeType referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos,value);
-
- // The derivative of the function defined on the actual element
- typename LocalInterpolationRule::DerivativeType derivative(0);
-
- for (size_t comp=0; comp<referenceDerivative.N(); comp++)
- jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
-
- //////////////////////////////////////////////////////////////
- // Exchange energy (aka harmonic energy)
- //////////////////////////////////////////////////////////////
- // The Frobenius norm is the correct norm here if the metric of TargetSpace is the identity.
- // (And if the metric of the domain space is the identity, which it always is here.)
- energy += weight * 0.5 * derivative.frobenius_norm2();
-
- //////////////////////////////////////////////////////////////
- // Dzyaloshinskii-Moriya interaction term
- //////////////////////////////////////////////////////////////
-
- // derivative[a][b] contains the partial derivative of m_a in the direction x_b
- Dune::FieldVector<field_type, 3> curl = {derivative[2][1], -derivative[2][0], derivative[1][0]-derivative[0][1]};
-
- FieldVector<field_type, 3> v = value.globalCoordinates();
-
- energy += weight * kappa_ * (v * curl);
-
- //////////////////////////////////////////////////////////////
- // Zeeman interaction term
- //////////////////////////////////////////////////////////////
- v[2] -= 1; // subtract e_3
- energy += weight * 0.5 * h_ * v.two_norm2();
-
- }
-
- return energy;
- }
-
- } // namespace GFE
-
-} // namespace Dune
-#endif
diff --git a/dune/gfe/densities/CMakeLists.txt b/dune/gfe/densities/CMakeLists.txt
index e326ccc61a9c161ae43c3192e51e64f2ac0f1f4c..cea45c0e9d897b47477a9d570a8062154dff5b12 100644
--- a/dune/gfe/densities/CMakeLists.txt
+++ b/dune/gfe/densities/CMakeLists.txt
@@ -1,6 +1,7 @@
#install headers
install(FILES
bulkcosseratdensity.hh
+ chiralskyrmiondensity.hh
cosseratshelldensity.hh
harmonicdensity.hh
localdensity.hh
diff --git a/dune/gfe/densities/chiralskyrmiondensity.hh b/dune/gfe/densities/chiralskyrmiondensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..82e4d8f10989805064c42bc569098cbf7c0501d8
--- /dev/null
+++ b/dune/gfe/densities/chiralskyrmiondensity.hh
@@ -0,0 +1,78 @@
+#ifndef DUNE_GFE_DENSITIES_CHIRALSKYRMIONDENSITY_HH
+#define DUNE_GFE_DENSITIES_CHIRALSKYRMIONDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/gfe/densities/localdensity.hh>
+
+namespace Dune::GFE
+{
+
+ /** \brief Energy density of certain chiral Skyrmion model
+ *
+ * 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 Position, class field_type>
+ class ChiralSkyrmionDensity
+ : public GFE::LocalDensity<Position,UnitVector<field_type,3> >
+ {
+ // various useful types
+ using TargetSpace = UnitVector<field_type,3>;
+ using Derivative = FieldMatrix<field_type, TargetSpace::embeddedDim, Position::size()>;
+
+ public:
+
+ ChiralSkyrmionDensity(const Dune::ParameterTree& parameters)
+ {
+ h_ = parameters.template get<double>("h");
+ kappa_ = parameters.template get<double>("kappa");
+ }
+
+ //! Dimension of a tangent space
+ constexpr static int blocksize = TargetSpace::TangentVector::dimension;
+
+ /** \brief Evaluation with the current position, the deformation function, the deformation gradient, the rotation and the rotation gradient
+ *
+ * \param x The current position
+ * \param value The deformation at the current position
+ * \param derivative The derivative of the deformation at the current position
+ */
+ virtual field_type operator() (const Position& x,
+ const TargetSpace& value,
+ const Derivative& derivative) const override
+ {
+ //////////////////////////////////////////////////////////////
+ // Exchange energy (aka harmonic energy)
+ //////////////////////////////////////////////////////////////
+
+ field_type density = 0.5 * derivative.frobenius_norm2();
+
+ //////////////////////////////////////////////////////////////
+ // Dzyaloshinskii-Moriya interaction term
+ //////////////////////////////////////////////////////////////
+
+ // derivative[a][b] contains the partial derivative of m_a in the direction x_b
+ FieldVector<field_type, 3> curl = {derivative[2][1], -derivative[2][0], derivative[1][0]-derivative[0][1]};
+
+ FieldVector<field_type, 3> v = value.globalCoordinates();
+
+ density += kappa_ * (v * curl);
+
+ //////////////////////////////////////////////////////////////
+ // Zeeman interaction term
+ //////////////////////////////////////////////////////////////
+ v[2] -= 1; // subtract e_3
+ density += 0.5 * h_ * v.two_norm2();
+
+ return density;
+ }
+
+ private:
+ field_type h_;
+ field_type kappa_;
+ };
+
+} // namespace Dune::GFE
+#endif
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 095cc200404c6568a10aa762b6ceb6df667810da..35a50711335bc48011f0de225e6582d1b5f0d341 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -43,8 +43,8 @@
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/assemblers/localintegralenergy.hh>
-#include <dune/gfe/assemblers/chiralskyrmionenergy.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
+#include <dune/gfe/densities/chiralskyrmiondensity.hh>
#include <dune/gfe/densities/harmonicdensity.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
@@ -269,41 +269,37 @@ int main (int argc, char *argv[])
// ////////////////////////////////////////////////////////////
typedef TargetSpace::rebind<adouble>::other ATargetSpace;
- using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
- using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
-
- // Assembler using ADOL-C
- std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
+ // First, the energy density
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 GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace>(harmonicDensity));
- else if (parameterSet["interpolationMethod"] == "projected")
- localEnergy.reset(new GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace>(harmonicDensity));
- else
- DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
+ using LocalCoordinate = GridType::Codim<0>::Entity::Geometry::LocalCoordinate;
+ std::shared_ptr<GFE::LocalDensity<LocalCoordinate,ATargetSpace> > density;
- } else if (energy == "chiral_skyrmion")
+ if (energy == "harmonic")
{
- // // Doesn't work: we are not inside of a template
- // if constexpr (std::is_same<TargetSpace, UnitVector<double,3> >::value)
- // {
- if (parameterSet["interpolationMethod"] == "geodesic")
- localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, GeodesicInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
- else if (parameterSet["interpolationMethod"] == "projected")
- localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, ProjectedInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
- else
- DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
- // } else
- // DUNE_THROW(Exception, "Build program with TargetSpace = UnitVector<3> for the ChiralSkyrmion energy!");
-
+ density = std::make_shared<GFE::HarmonicDensity<LocalCoordinate, ATargetSpace> >();
+ }
+ else if (energy == "chiral_skyrmion")
+ {
+ density = std::make_shared<GFE::ChiralSkyrmionDensity<LocalCoordinate, adouble> >(parameterSet.sub("energyParameters"));
} else
DUNE_THROW(Exception, "Unknown energy type '" << energy << "'");
+ // Next: The local energy, i.e., the integral of the density over one element
+ using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+ using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+
+ std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy;
+
+ if (parameterSet["interpolationMethod"] == "geodesic")
+ localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >(density);
+ else if (parameterSet["interpolationMethod"] == "projected")
+ localEnergy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >(density);
+ else
+ DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
+
+ // Compute local tangent problems by applying ADOL-C directly to the energy on the element
LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy.get());
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, localGFEADOLCStiffness);