Change Transformator , prepare Morley & Argyris

dune/functions/functionspacebases/argyrisbasis.hh

@@ -22,55 +22,45 @@ namespace Dune
       template <class R>
       class ArgyrisTransformator
       {
+        ArgyrisTransformator()
+        {
+          setupMatrix();
+        }
+
+        template <class Element>
+        void bind(Element const &e)
+        {
+          if (e.geometry().isAffine)
+            fillMatrix(Dune::referenceElement<double, 2>(GeometryTypes::simplex(2)).position(0, 2), e.geometry()); // barycenter, because we need some value.
+        }
+
         template <class Values, class LocalCoordinate, class Geometry>
         static void apply(Values &values, LocalCoordinate const &xi,
                           Geometry const &geometry)
         {
-          Dune::FieldMatrix<R, 21, 21> V;
-          Dune::FieldMatrix<R, 21, 24> E;
-          Dune::FieldMatrix<R, 24, 24> Vc;
-          Dune::FieldMatrix<R, 24, 21> D;
-          // #Todo move this to some static Variable
-          for (std::size_t i = 0; i < 19; ++i)
-            E[i][i] = 1.;
-          E[19][20] = 1;
-          E[20][22] = 1;
-
-          // Construct Vc
-          auto JT = geometry.jacobianTransposed();
-          auto JTinv = geometry.jacobianInverseTransposed();
-
-          // I don't see through without
-          auto const &dxdx = JT[0][0], dxdy = JT[1][0], dydx = JT[0][1],
-                     dydy = JT[1][1];
-
-          Dune::FieldMatrix<R, 3, 3> theta = {
-              {dxdx * dxdx, 2. * dxdx * dydx, dydx * dydx},
-              {dxdy * dxdx, dxdy * dydx + dxdx * dydy, dydx * dydy},
-              {dxdy * dxdy, 2. * dxdy * dydy, dydy * dydy}};
-          for (std::size_t i = 0; i < 3;
-               ++i) // three blocks of vertex evaluations
+
+          if (!geometry.isAffine())
           {
-            Vc[i][i] = 1.;
-            std::size_t shift = 6 * i + 1;
-            for (std::size_t j = 0; j < 2; ++j)
-              for (std::size_t k = 0; k < 2; ++k)
-                Vc[shift + j][shift + k] = JTinv[j][k];
-            shift += 2;
-            for (std::size_t j = 0; j < 3; ++j)
-              for (std::size_t k = 0; k < 3; ++k)
-                Vc[shift + j][shift + k] = theta[j][k];
+            fillMatrix(xi, geometry);
           }
+          Values tmp = values; // needs to be deep copy
+          mat.mv(tmp, values);
+        }
+
+      private:
+        void setupMatrix()
+        {
+          // TODO implement
         }
 
-        template <class Gradients, class LocalCoordinate, class Geometry>
-        static void applyJacobian(Gradients &gradients,
-                                  LocalCoordinate const &xi,
-                                  Geometry const &geometry)
+        template <class LocalCoordinate, class Geometry>
+        void fillMatrix(LocalCoordinate const &xi, Geometry const &geometry)
         {
-          // #Todo implement
+          // TODO implement
         }
+        BCRSMatrix<R> mat_;
 
+      public:
         /** \brief Wrapper around a callable that applies the inverse
          * transformation */
         template <class Function, class LocalCoordinate, class Element>
@@ -83,7 +73,20 @@ namespace Dune
           LocalValuedFunction(const Function &f, const Element &element)
               : f_(f), element_(element) {}
 
-          auto operator()(const LocalCoordinate &xi) const {}
+          auto operator()(const LocalCoordinate &xi) const
+          {
+            // TODO implement
+          }
+
+          friend auto derivative(LocalValuedFunction const &)
+          {
+            // TODO implement
+          }
+
+          auto normalDerivative(std::size_t i)
+          {
+            // TODO implement
+          }
         };
       };
     } // namespace Impl

dune/functions/functionspacebases/lineartransformedlocalfiniteelement.hh

@@ -25,11 +25,15 @@ namespace Dune
      *      - The LinearTransformator classes are expected to implement a method bind(Element)
      *        which sets up the transformation, for example by filling a sparse Matrix with values.
      *      - The LinearTransformator classes are not static.
-     *      - The inner class LocalValuedFunction should work independent of the instance of LinearTransformator,
-     *        as it has all the information needed , i.e. the Element, and should make use of the fact that
-     *        the dune-functions interface mandates the derivative of a LocalFunction to be global-valued.
-     *      - Additionally the inner class LocalValuedFunction shound implement the DifferentiableFunction
-     *        Interface, or be differentiable as often as needed by the Node set.
+     *      - The inner class LocalValuedFunction cannot implement the inverse of the transformation in the current
+     *        implementation, as this would require changes to LinearTransformedInterpolation and also the inverse of
+     *        the transformation (sparse Matrix inversion!). Instead, the LocalValuedFunction is the object, that,
+     *        when evaluated by the LocalInterpolation class, gives the coefficients as if it was evaluated by the
+     *        GlobalValued Interpolation. It can (and should) therefore implement methods used by the, for
+     *        example normalDerivative(size_t) and should make use of the fact that the dune-functions interface
+     *        mandates the derivative of a LocalFunction to be global-valued.
+     *      - Additionally the inner class LocalValuedFunction should be differentiable as often as needed by the Node
+     *        set.
      *
      *    LinearTransformedLocalBasis:
      *      - holds an instance of the Transformator class
@@ -37,7 +41,7 @@ namespace Dune
      *      - higher derivatives are implemented
      *
      *    LinearTransformedLocalInterpolation (defined only as typedef):
-     *      - same as GlobalValuedTransformation
+     *      - same as GlobalValuedInterpolation
      *
      *    LinearTransformedLocalCoefficients (class does not exist):
      *      - same as in the LocalCoefficients of the localvalued finite element

dune/functions/functionspacebases/morleybasis.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_C1ELEMENTS_MORLEYBASIS_HH
+#define DUNE_C1ELEMENTS_MORLEYBASIS_HH
+
+namespace Dune
+{
+  namespace Functions
+  {
+    namespace Impl
+    {
+      /** please doc me */
+      template <class R>
+      class MorleyTransformator
+      {
+        MorleyTransformator()
+        {
+          setupMatrix();
+        }
+
+        template <class Element>
+        void bind(Element const &e)
+        {
+          if (e.geometry().isAffine)
+            fillMatrix(Dune::referenceElement<double, 2>(GeometryTypes::simplex(2)).position(0, 2), e.geometry()); // barycenter, because we need some value.
+        }
+
+        template <class Values, class LocalCoordinate, class Geometry>
+        void apply(Values &values, LocalCoordinate const &xi,
+                   Geometry const &geometry)
+        {
+          if (!geometry.isAffine())
+          {
+            fillMatrix(xi, geometry);
+          }
+          Values tmp = values; // needs to be deep copy
+          mat.mv(tmp, values);
+        }
+
+      private:
+        void setupMatrix()
+        {
+          // TODO implement
+        }
+
+        template <class LocalCoordinate, class Geometry>
+        void fillMatrix(LocalCoordinate const &xi, Geometry const &geometry)
+        {
+          // TODO implement
+        }
+        BCRSMatrix<R> mat_;
+
+      public:
+        /** \brief Wrapper that is evaluated by LocalInterpolation as if it was evaluated by GlobalInterpolation */
+        template <class Function, class LocalCoordinate, class Element>
+        class LocalValuedFunction
+        {
+          const Function &f_;
+          const Element &element_;
+
+        public:
+          LocalValuedFunction(const Function &f, const Element &element)
+              : f_(f), element_(element) {}
+
+          auto operator()(const LocalCoordinate &xi) const
+          {
+            // TODO implement
+          }
+
+          // actually not needed
+          // friend auto derivative(LocalValuedFunction const &)
+          // {
+          //   // TODO implement
+          // }
+
+          auto normalDerivative(std::size_t i)
+          {
+            // TODO implement
+          }
+        };
+      };
+    } // namespace Impl
+
+    template <class GV, class R>
+    class MorleyNode;
+
+    template <class GV, class MI, typename R>
+    class MorleyPreBasis
+    {
+      static const int dim = GV::dimension;
+      static_assert(dim == 2,
+                    "Morley PreBasis only implemented for 2d simplices");
+
+    public:
+      //! The grid view that the FE basis is defined on
+      using GridView = GV;
+
+      //! Type used for indices and size information
+      using size_type = std::size_t;
+
+      //! Template mapping root tree path to type of created tree node
+      using Node = MorleyNode<GridView, R>;
+
+      static constexpr size_type maxMultiIndexSize = 1;
+      static constexpr size_type minMultiIndexSize = 1;
+      static constexpr size_type multiIndexBufferSize = 1;
+
+      //! Constructor for a given grid view object
+      MorleyPreBasis(const GV &gv)
+          : gridView_(gv)
+      {
+      }
+
+      //! Initialize the global indices
+      void initializeIndices()
+      { // TODO check if we need to do something
+      }
+
+      //! Obtain the grid view that the basis is defined on
+      const GridView &gridView() const
+      {
+        return gridView_;
+      }
+
+      //! Update the stored grid view, to be called if the grid has changed
+      void update(const GridView &gv)
+      {
+        gridView_ = gv;
+      }
+
+      /**
+       * \brief Create tree node
+       */
+      Node makeNode() const
+      {
+        return Node();
+      }
+
+      //! Same as size(prefix) with empty prefix
+      size_type size() const
+      {
+        return 6;
+      }
+
+      //! Return number of possible values for next position in multi index
+      template <class SizePrefix>
+      size_type size(const SizePrefix prefix) const
+      {
+        // this basically means this is a leaf node with dimrange 1 right?
+        assert(prefix.size() == 0 || prefix.size() == 1);
+        return (prefix.size() == 0) ? size() : 0;
+      }
+
+      //! Get the total dimension of the space spanned by this basis
+      size_type dimension() const
+      {
+        return size();
+      }
+
+      //! Get the maximal number of DOFs associated to node for any element
+      size_type maxNodeSize() const
+      {
+        return 6;
+      }
+
+      template <typename It>
+      It indices(const Node &node, It it) const
+      {
+        const auto &gridIndexSet = gridView().indexSet();
+        const auto &element = node.element();
+
+        // throw if Element is not simplex
+        if (not(element.type().isSimplex()))
+          DUNE_THROW(Dune::NotImplemented, "Morley Basis only implemented for simplex elements");
+        for (size_type i = 0, end = node.finiteElement().size(); i < end; ++it, ++i)
+        {
+          Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
+          // TODO probably unnecessary
+          if (!gridIndexSet.contains(element))
+            DUNE_THROW(Dune::RangeError, "Element is not in gridIndexSet!");
+
+          *it = {{(size_type)((localKey.codim() == dim) ? gridIndexSet.subIndex(element, localKey.subEntity(), dim) : 3 * gridIndexSet.size(dim) + gridIndexSet.subIndex(element, localKey.subEntity(), 1))}};
+        }
+        return it;
+      }
+
+    protected:
+      GridView gridView_;
+    };
+  };
+
+  template <class GV, class R>
+  class MorleyNode : public LeafBasisNode
+  {
+    using D = typename GV::ctype;
+    using LocalValuedFE = MorleyLocalFiniteElement<D, R>;
+
+  public:
+    using size_type = std::size_t;
+    using Element = typename GV::template Codim<0>::Entity;
+    using FiniteElement = Impl::LinearTransformedLocalFiniteElement<Impl::MorleyTransformator, MorleyLocalFiniteElement<D, R>, Element>;
+
+    MorleyNode()
+        : localValuedFiniteElement_(std::make_shared<LocalValuedFE>()),
+          finiteElement_(std::make_shared<FiniteElement>()),
+          element_(nullptr)
+    {
+      this->setSize(finiteElement_->size());
+    }
+
+    ~MorleyNode()
+    {
+    }
+
+    //! 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)
+    {
+      if (not e.type().isSimplex())
+        DUNE_THROW(Dune::NotImplemented, "MorleyBasis can only be bound to simplex elements");
+      element_ = &e;
+
+      finiteElement_->bind(*localValuedFiniteElement_, *element_);
+    }
+
+  protected:
+    unsigned int order() const { return 3; }
+    std::shared_ptr<LocalValuedFE> localValuedFiniteElement_;
+    std::shared_ptr<FiniteElement> finiteElement_;
+    const Element *element_;
+  };
+
+} // namespace Impl
+} // namespace Functions
+} // namespace Dune
+
+#endif
\ No newline at end of file