Q1 local basis free of std::vector

Imported and adapted from dune-localfunctions q1localbasis

src/q1localbasis.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef PACXX_PROJECTSEMINAR_2019_SRC_Q1LOCALBASIS_HH
+#define PACXX_PROJECTSEMINAR_2019_SRC_Q1LOCALBASIS_HH
+
+#include <array>
+#include <cstddef>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/rangeutilities.hh>
+
+namespace PPS {
+  /**
+   * Lagrange shape functions of order 1 on the reference cube.
+   *
+   * Also known as \f$Q^1\f$.
+   *
+   * \tparam D Type to represent the field in the domain.
+   * \tparam R Type to represent the field in the range.
+   * \tparam dim Dimension of the cube
+   */
+  template<class D, class R, int dim>
+  class Q1LocalBasis
+  {
+  public:
+    static constexpr int dimDomain = dim;
+    using DomainField = D;
+    using Domain = Dune::FieldVector<D, dim>;
+    static constexpr int dimRange = 1;
+    using RangeField = R;
+    using Range = Dune::FieldVector<R, 1>;
+    using Jacobian = Dune::FieldMatrix<R,1,dim>;
+
+    //! \brief number of shape functions
+    static constexpr std::size_t size()
+    {
+      return 1<<dim;
+    }
+
+    //! \brief Evaluate all shape functions
+    auto evaluateFunction (const Domain& pos) const
+    {
+      std::array<Range, size()> result;
+
+      for (auto i : Dune::range(size()))
+      {
+        result[i] = 1;
+
+        for(auto j : Dune::range(dim))
+          // if j-th bit of i is set multiply with pos[j], else with 1-pos[j]
+          result[i] *= (i & (1<<j)) ? pos[j] :  1-pos[j];
+      }
+
+      return result;
+    }
+
+    //! \brief Evaluate Jacobian of all shape functions
+    auto evaluateJacobian (const Domain& pos) const
+    {
+      std::array<Jacobian, size()> result;
+
+      // Loop over all shape functions
+      for (auto i : Dune::range(size()))
+        // Loop over all coordinate directions
+        for (auto j : Dune::range(dim))
+        {
+          // Initialize: the overall expression is a product
+          // if j-th bit of i is set to -1, else 1
+          result[i][0][j] = (i & (1<<j)) ? 1 : -1;
+
+          for (auto k : Dune::range(dim))
+            if (j!=k)
+              // if k-th bit of i is set multiply with pos[j], else with 1-pos[j]
+              result[i][0][j] *= (i & (1<<k)) ? pos[k] :  1-pos[k];
+        }
+
+      return result;
+    }
+
+    //! \brief Polynomial order of the shape functions
+    static constexpr unsigned int order()
+    {
+      return 1;
+    }
+  };
+}
+
+#endif // DUNE_Q1_LOCALBASIS_HH