From afd79096e3dff81c904c02151270d7794fd67659 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 8 Feb 2016 14:01:10 +0100
Subject: [PATCH] Allow the possibility to have a periodic 1d function space if
the grid is 1d
---
dune/gfe/periodic1dpq1nodalbasis.hh | 256 ++++++++++++++++++++++++++++
src/gradient-flow.cc | 2 +
2 files changed, 258 insertions(+)
create mode 100644 dune/gfe/periodic1dpq1nodalbasis.hh
diff --git a/dune/gfe/periodic1dpq1nodalbasis.hh b/dune/gfe/periodic1dpq1nodalbasis.hh
new file mode 100644
index 00000000..df31bc5a
--- /dev/null
+++ b/dune/gfe/periodic1dpq1nodalbasis.hh
@@ -0,0 +1,256 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_PERIODIC_1D_PQ1NODALBASIS_HH
+#define DUNE_GFE_PERIODIC_1D_PQ1NODALBASIS_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/typetree/leafnode.hh>
+
+#include <dune/functions/functionspacebases/nodes.hh>
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/flatmultiindex.hh>
+
+
+namespace Dune {
+namespace Functions {
+
+// *****************************************************************************
+// This is the reusable part of the basis. It contains
+//
+// PQ1NodeFactory
+// PQ1NodeIndexSet
+// PQ1Node
+//
+// The factory allows to create the others and is the owner of possible shared
+// state. These three components do _not_ depend on the global basis or index
+// set and can be used without a global basis.
+// *****************************************************************************
+
+template<typename GV, typename ST, typename TP>
+class Periodic1DPQ1Node;
+
+template<typename GV, class MI, class TP, class ST>
+class Periodic1DPQ1NodeIndexSet;
+
+template<typename GV, class MI, class ST>
+class Periodic1DPQ1NodeFactory;
+
+template<typename GV, class MI, class ST>
+class Periodic1DPQ1NodeFactory
+{
+ static const int dim = GV::dimension;
+
+public:
+
+ /** \brief The grid view that the FE space is defined on */
+ using GridView = GV;
+ using size_type = ST;
+
+ template<class TP>
+ using Node = Periodic1DPQ1Node<GV, size_type, TP>;
+
+ template<class TP>
+ using IndexSet = Periodic1DPQ1NodeIndexSet<GV, MI, TP, ST>;
+
+ /** \brief Type used for global numbering of the basis vectors */
+ using MultiIndex = MI;
+
+ using SizePrefix = Dune::ReservedVector<size_type, 2>;
+
+ /** \brief Constructor for a given grid view object */
+ Periodic1DPQ1NodeFactory(const GridView& gv) :
+ gridView_(gv)
+ {}
+
+ void initializeIndices()
+ {}
+
+ /** \brief Obtain the grid view that the basis is defined on
+ */
+ const GridView& gridView() const
+ {
+ return gridView_;
+ }
+
+ template<class TP>
+ Node<TP> node(const TP& tp) const
+ {
+ return Node<TP>{tp};
+ }
+
+ template<class TP>
+ IndexSet<TP> indexSet() const
+ {
+ return IndexSet<TP>{*this};
+ }
+
+ size_type size() const
+ {
+ return gridView_.size(dim)-1;
+ }
+
+ //! Return number possible values for next position in multi index
+ size_type size(const SizePrefix prefix) const
+ {
+ if (prefix.size() == 0)
+ return size();
+ if (prefix.size() == 1)
+ return 0;
+ DUNE_THROW(RangeError, "Method size() can only be called for prefixes of length up to one");
+ }
+
+ /** \todo This method has been added to the interface without prior discussion. */
+ size_type dimension() const
+ {
+ return size()-1;
+ }
+
+ size_type maxNodeSize() const
+ {
+ return StaticPower<2,GV::dimension>::power;
+ }
+
+//protected:
+ const GridView gridView_;
+};
+
+
+
+template<typename GV, typename ST, typename TP>
+class Periodic1DPQ1Node :
+ public LeafBasisNode<ST, TP>
+{
+ static const int dim = GV::dimension;
+ static const int maxSize = StaticPower<2,GV::dimension>::power;
+
+ using Base = LeafBasisNode<ST,TP>;
+ using FiniteElementCache = typename Dune::PQkLocalFiniteElementCache<typename GV::ctype, double, dim, 1>;
+
+public:
+
+ using size_type = ST;
+ using TreePath = TP;
+ using Element = typename GV::template Codim<0>::Entity;
+ using FiniteElement = typename FiniteElementCache::FiniteElementType;
+
+ Periodic1DPQ1Node(const TreePath& treePath) :
+ Base(treePath),
+ finiteElement_(nullptr),
+ element_(nullptr)
+ {}
+
+ //! Return current element, throw if unbound
+ const Element& element() const
+ {
+ return *element_;
+ }
+
+ /** \brief Return the LocalFiniteElement for the element we are bound to
+ *
+ * The LocalFiniteElement implements the corresponding interfaces of the dune-localfunctions module
+ */
+ const FiniteElement& finiteElement() const
+ {
+ return *finiteElement_;
+ }
+
+ //! Bind to element.
+ void bind(const Element& e)
+ {
+ element_ = &e;
+ finiteElement_ = &(cache_.get(element_->type()));
+ this->setSize(finiteElement_->size());
+ }
+
+protected:
+
+ FiniteElementCache cache_;
+ const FiniteElement* finiteElement_;
+ const Element* element_;
+};
+
+
+
+template<typename GV, class MI, class TP, class ST>
+class Periodic1DPQ1NodeIndexSet
+{
+ enum {dim = GV::dimension};
+
+public:
+
+ using size_type = ST;
+
+ /** \brief Type used for global numbering of the basis vectors */
+ using MultiIndex = MI;
+
+ using NodeFactory = Periodic1DPQ1NodeFactory<GV, MI, ST>;
+
+ using Node = typename NodeFactory::template Node<TP>;
+
+ Periodic1DPQ1NodeIndexSet(const NodeFactory& nodeFactory) :
+ nodeFactory_(&nodeFactory)
+ {}
+
+ /** \brief Bind the view to a grid element
+ *
+ * Having to bind the view to an element before being able to actually access any of its data members
+ * offers to centralize some expensive setup code in the 'bind' method, which can save a lot of run-time.
+ */
+ void bind(const Node& node)
+ {
+ node_ = &node;
+ }
+
+ /** \brief Unbind the view
+ */
+ void unbind()
+ {
+ node_ = nullptr;
+ }
+
+ /** \brief Size of subtree rooted in this node (element-local)
+ */
+ size_type size() const
+ {
+ return node_->finiteElement().size();
+ }
+
+ //! Maps from subtree index set [0..size-1] to a globally unique multi index in global basis
+ MultiIndex index(size_type i) const
+ {
+ Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
+ const auto& gridIndexSet = nodeFactory_->gridView().indexSet();
+ const auto& element = node_->element();
+
+ //return {{ gridIndexSet.subIndex(element,localKey.subEntity(),dim) }};
+
+ MultiIndex idx{{gridIndexSet.subIndex(element,localKey.subEntity(),dim) }};
+
+ // make periodic
+ if (idx == gridIndexSet.size(dim)-1)
+ idx = {{0}};
+
+ return idx;
+ }
+
+protected:
+ const NodeFactory* nodeFactory_;
+
+ const Node* node_;
+};
+
+/** \brief Nodal basis of a scalar first-order Lagrangian finite element space
+ *
+ * \tparam GV The GridView that the space is defined on
+ * \tparam ST The type used for local indices; global indices are FlatMultiIndex<ST>
+ */
+template<typename GV, class ST = std::size_t>
+using Periodic1DPQ1NodalBasis = DefaultGlobalBasis<Periodic1DPQ1NodeFactory<GV, FlatMultiIndex<ST>, ST> >;
+
+} // end namespace Functions
+} // end namespace Dune
+
+#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_PQ1NODALBASIS_HH
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 57de84e0..957881a7 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -44,6 +44,7 @@
#include <dune/gfe/harmonicenergystiffness.hh>
#include <dune/gfe/l2distancesquaredenergy.hh>
#include <dune/gfe/weightedsumenergy.hh>
+#include <dune/gfe/periodic1dpq1nodalbasis.hh>
// grid dimension
const int dim = 1;
@@ -142,6 +143,7 @@ int main (int argc, char *argv[]) try
//////////////////////////////////////////////////////////////////////////////////
typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView, order> FEBasis;
+ //typedef Dune::Functions::Periodic1DPQ1NodalBasis<typename GridType::LeafGridView> FEBasis;
FEBasis feBasis(grid->leafGridView());
///////////////////////////////////////////
--
GitLab