From d40fa780570387150264f8bcca631e4ed6e917ea Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sat, 30 Jan 2016 15:24:13 +0100
Subject: [PATCH] Add L2DistanceSquaredEnergy and WeightedSumEnergy
This are needed by the gradientflow code to implement the implicit Euler method
for an L2 gradient flow.
---
dune/gfe/l2distancesquaredenergy.hh | 71 +++++++++++++++++++++++++++++
dune/gfe/weightedsumenergy.hh | 54 ++++++++++++++++++++++
2 files changed, 125 insertions(+)
create mode 100644 dune/gfe/l2distancesquaredenergy.hh
create mode 100644 dune/gfe/weightedsumenergy.hh
diff --git a/dune/gfe/l2distancesquaredenergy.hh b/dune/gfe/l2distancesquaredenergy.hh
new file mode 100644
index 00000000..d08571cd
--- /dev/null
+++ b/dune/gfe/l2distancesquaredenergy.hh
@@ -0,0 +1,71 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_L2_DISTANCE_SQUARED_ENERGY_HH
+#define DUNE_GFE_L2_DISTANCE_SQUARED_ENERGY_HH
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+
+template<class Basis, class TargetSpace>
+class L2DistanceSquaredEnergy
+ : public LocalGeodesicFEStiffness<Basis,TargetSpace>
+{
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+
+ // some other sizes
+ enum {gridDim=GridView::dimension};
+
+public:
+
+ // This is the function that we are computing the L2-distance to
+ std::shared_ptr<VirtualGridViewFunction<GridView,UnitVector<double,3> > > origin_;
+
+ /** \brief Assemble the energy for a single element */
+ RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localSolution) const override
+
+ {
+ RT energy = 0;
+
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
+ LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+
+ // Just guessing an appropriate quadrature order
+ auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
+
+ 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 double integrationElement = element.geometry().integrationElement(quadPos);
+
+ auto weight = quad[pt].weight() * integrationElement;
+
+ // The function value
+ auto value = localGeodesicFEFunction.evaluate(quadPos);
+ UnitVector<double,3> originValue;
+ origin_->evaluateLocal(element,quadPos, originValue);
+
+ // Add the local energy density
+ energy += weight * TargetSpace::distanceSquared(originValue, value);
+
+ }
+
+ return energy;
+ }
+
+};
+
+#endif
+
diff --git a/dune/gfe/weightedsumenergy.hh b/dune/gfe/weightedsumenergy.hh
new file mode 100644
index 00000000..38d5cbef
--- /dev/null
+++ b/dune/gfe/weightedsumenergy.hh
@@ -0,0 +1,54 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_WEIGHTED_SUM_ENERGY_HH
+#define DUNE_GFE_WEIGHTED_SUM_ENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/localgeodesicfestiffness.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+
+template<class Basis, class TargetSpace>
+class WeightedSumEnergy
+ : public LocalGeodesicFEStiffness<Basis,TargetSpace>
+{
+ // grid types
+ typedef typename Basis::GridView GridView;
+ typedef typename GridView::ctype DT;
+ typedef typename TargetSpace::ctype RT;
+
+ // some other sizes
+ enum {gridDim=GridView::dimension};
+
+public:
+
+ std::vector<std::shared_ptr<LocalGeodesicFEStiffness<Basis,TargetSpace> > > addends_;
+
+ std::vector<double> weights_;
+
+ WeightedSumEnergy(std::vector<std::shared_ptr<LocalGeodesicFEStiffness<Basis,TargetSpace> > > addends,
+ std::vector<double> weights)
+ : addends_(addends),
+ weights_(weights)
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ RT energy (const typename Basis::LocalView& localView,
+ const std::vector<TargetSpace>& localConfiguration) const override
+
+ {
+ RT energy = 0;
+
+ assert(weights_.size() == addends_.size());
+
+ for (size_t i=0; i<addends_.size(); i++)
+ energy += weights_[i] * addends_[i]->energy(localView, localConfiguration);
+
+ return energy;
+ }
+
+};
+
+#endif
+
--
GitLab