Argyris LocalInterpolation & bugfixes

dune/localfunctions/argyris.hh

@@ -16,9 +16,9 @@
 #include <dune/localfunctions/common/localinterpolation.hh>
 #include <dune/localfunctions/common/localkey.hh>
 
+#include <dune/localfunctions/morley.hh>
 #include <dune/localfunctions/polynomialbasiscoefficients.hh>
 #include <dune/localfunctions/utility/polynomialbasis.hh>
-
 namespace Dune
 {
   namespace Impl
@@ -33,6 +33,7 @@ namespace Dune
        */
 
     template <class D, class R>
+    class ArgyrisLocalBasis:
     class ArgyrisLocalBasis
         : public PolynomialBasisWithMatrix<
               StandardEvaluator<VirtualMonomialBasis<2, D>>, SparseCoeffMatrix<D, 1>, D, R>
@@ -68,8 +69,13 @@ namespace Dune
 
       ArgyrisLocalBasis(std::bitset<numberOfEdges> edgeOrientation) : ArgyrisLocalBasis()
       {
-        for (std::size_t i = 0; i < edgeOrientation_.size(); i++)
-          edgeOrientation_[i] *= edgeOrientation[i] ? -1.0 : 1.0;
+        for (std::size_t i = 0; i < edgeOrientation_.size(); ++i)
+          edgeOrientation_[i] = edgeOrientation[i] ? -1.0 : 1.0;
+      }
+
+      ArgyrisLocalBasis(ArgyrisLocalBasis const &other) : ArgyrisLocalBasis()
+      {
+        edgeOrientation_ = other.edgeOrientation_;
       }
 
       static unsigned int size() { return coeffSize; }
@@ -82,6 +88,7 @@ namespace Dune
       void evaluateFunction(const typename Traits::DomainType &in,
                             std::vector<typename Traits::RangeType> &out) const
       {
+        out.resize(size());
         Base::evaluateFunction(in, out);
         for (std::size_t i = 0; i < numberOfEdges; i++)
           out[18 + i] *= edgeOrientation_[i];
@@ -95,6 +102,8 @@ namespace Dune
       void evaluateJacobian(const typename Traits::DomainType &in,
                             std::vector<typename Traits::JacobianType> &out) const
       {
+        out.resize(size());
+
         Base::evaluateJacobian(in, out);
         for (std::size_t i = 0; i < numberOfEdges; i++)
           out[18 + i] *= edgeOrientation_[i];
@@ -110,6 +119,8 @@ namespace Dune
                    const typename Traits::DomainType &in,
                    std::vector<typename Traits::RangeType> &out) const
       {
+        out.resize(size());
+
         Base::partial(order, in, out);
         for (std::size_t i = 0; i < numberOfEdges; i++)
           out[18 + i] *= edgeOrientation_[i];
@@ -160,36 +171,20 @@ namespace Dune
      *  \tparam LocalBasis The corresponding set of shape functions
      */
     template <class LocalBasis>
-    class ArgyrisLocalInterpolation
+    class ArgyrisLocalInterpolation: MorleyLocalInterpolation<LocalBasis>
     {
 
       using D = typename LocalBasis::Traits::DomainFieldType;
       using LocalCoordinate = typename LocalBasis::Traits::DomainType;
       static constexpr unsigned int numberOfEdges = 3;
+      using Base = MorleyLocalInterpolation<LocalBasis>;
+      using Base::matrixToVector;
+      using Base::normalDerivative;
+      using Base::s_;
+      using Base::v_;
 
     public:
-      ArgyrisLocalInterpolation(std::bitset<numberOfEdges> s = 0)
-      {
-        auto refElement = Dune::referenceElement<double, 2>(GeometryTypes::simplex(2));
-
-        for (std::size_t i = 0; i < numberOfEdges; i++)
-          m_[i] = refElement.position(i, 1);
-
-        for (std::size_t i = 0; i < numberOfEdges; i++)
-        {
-          auto vertexIterator = refElement.subEntities(i, 1, 2).begin();
-          v_[i] = refElement.position(*vertexIterator, 2);
-          auto v0 = *vertexIterator;
-          auto v1 = *(++vertexIterator);
-          // By default, edges point from the vertex with the smaller index
-          // to the vertex with the larger index.
-          if (v0 > v1)
-            std::swap(v0, v1);
-          edge_[i] = refElement.position(v1, 2) - refElement.position(v0, 2);
-          edge_[i] *= (s[i]) ? -1.0 : 1.0;
-          normal_[i] = {-edge[i][1], edge[i][0]}; // rotation by 90 degree
-        }
-      }
+      ArgyrisLocalInterpolation(std::bitset<numberOfEdges> s = 0) : Base(s) {}
 
       /** \brief Evaluate a given function at the Lagrange nodes
        *
@@ -205,8 +200,12 @@ namespace Dune
 
         auto &&f = Impl::makeFunctionWithCallOperator<typename LocalBasis::Traits::DomainType>(ff);
 
-        auto df = derivative(ff);
-        auto hessf = derivative(derivative(ff));
+        // auto df = derivative(ff);
+        // auto hessf = derivative(derivative(ff));
+
+        auto df = makeDerivative<LocalCoordinate>(ff);
+        auto hessf = makeDerivative<LocalCoordinate>(df);
+
         out.resize(LocalBasis::size());
 
         // Iterate over vertices, 6 dofs per vertex
@@ -214,33 +213,22 @@ namespace Dune
         {
           auto const &derivativeValue = df(v_[i]);
           auto const &hessianValue = hessf(v_[i]);
+
           out[i * 6 + 0] = f(v_[i]);
-          out[i * 6 + 1] = derivativeValue[0];
-          out[i * 6 + 2] = derivativeValue[1];
-          out[i * 6 + 3] = hessianValue[0][0];
-          out[i * 6 + 4] = hessianValue[1][0];
-          out[i * 6 + 5] = hessianValue[1][1];
+          out[i * 6 + 1] = matrixToVector(derivativeValue)[0];
+          out[i * 6 + 2] = matrixToVector(derivativeValue)[1];
+          out[i * 6 + 3] = matrixToVector(hessianValue[0])[0];
+          out[i * 6 + 4] = matrixToVector(hessianValue[0])[1];
+          out[i * 6 + 5] = matrixToVector(hessianValue[1])[1];
         }
 
         // iterate over edges, one dof per edge
-        for (unsigned int i = 0; i < 3; ++i)
+        for (unsigned int i = 0; i < numberOfEdges; ++i)
         {
-          auto const &edgeCenter = m_[i];
-          out[shiftVertices + i] = df(edgeCenter).dot(normal_[i]);
+          out[18 + i] = normalDerivative(ff, i) * (s_[i] ? -1 : 1);
         }
         return;
       }
-
-    private:
-      // Vertices of the reference simplex
-      std::array<typename LocalBasis::Traits::DomainType, numberOfEdges> v_;
-      // Edge midpoints of the reference simplex
-      std::array<typename LocalBasis::Traits::DomainType, numberOfEdges> m_;
-      // Edges of the reference simplex
-      std::array<typename LocalBasis::Traits::DomainType, numberOfEdges>
-          edge_; // TODO evtl remove member
-      // normals of the reference simplex
-      std::array<typename LocalBasis::Traits::DomainType, numberOfEdges> normal_;
     };
   } // namespace Impl
   /**
@@ -259,10 +247,7 @@ namespace Dune
     ArgyrisLocalFiniteElement(std::bitset<3> s = 0)
         : basis_(s), interpolation_(s) {}
 
-    const typename Traits::LocalBasisType &localBasis() const
-    {
-      return localBasis_;
-    }
+    const typename Traits::LocalBasisType &localBasis() const { return basis_; }
 
     /** \brief Returns the assignment of the degrees of freedom to the element subentities
      */
@@ -292,7 +277,7 @@ namespace Dune
     }
 
   private:
-    typename Traits::LocalBasisType localBasis_;
+    typename Traits::LocalBasisType basis_;
     typename Traits::LocalCoefficientsType coefficients_;
     typename Traits::LocalInterpolationType interpolation_;
   };