Implement jacobian_apply_volume() analytically

Addresses: #6

src/nonlinearpoissonfem.hh

@@ -29,7 +29,6 @@ namespace PPS {
   template<typename Param>
   class NonlinearPoissonFEM :
     public NumericalJacobianVolume<NonlinearPoissonFEM<Param> >,
-    public NumericalNonlinearJacobianApplyVolume<NonlinearPoissonFEM<Param> >,
     public FullVolumePattern,
     public LocalOperatorDefaultFlags
   {
@@ -164,6 +163,63 @@ namespace PPS {
                    q*phihat[i])*factor;
       }
     }
+
+    //! apply local jacobian of the volume term
+    /**
+     * \param lp The point at which to linearize
+     * \param x  The point to apply the operator to
+     * \param y  Where to add the result to
+     */
+    template<class Element, class LocalBasis, class X, class Y>
+    void jacobian_apply_volume(const Element &elem, const LocalBasis &lb,
+                               const X& lp, const X& x, Y& y) const
+    {
+      // types & dimension
+      const int dim = Element::dimension;
+      using RF = typename LocalBasis::Traits::RangeFieldType;
+
+      // select quadrature rule
+      auto geo = elem.geometry();
+      const int order = incrementorder() + 2*lb.order();
+
+      std::vector<typename LocalBasis::Traits::RangeType> phihat;
+      std::vector<typename LocalBasis::Traits::JacobianType> gradphi(lb.size());
+      std::vector<typename LocalBasis::Traits::JacobianType> gradphihat;
+
+      // loop over quadrature points
+      for (const auto& ip : PPS::quadratureRule(geo, order))
+      {
+        // evaluate basis functions
+        lb.evaluateFunction(ip.position(), phihat);
+
+        // evaluate u
+        auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
+        // w is the function corresponding to coefficients lp
+        auto w = std::inner_product(begin(lp), end(lp), begin(phihat), RF(0));
+
+        // evaluate gradient of shape functions
+        lb.evaluateJacobian(ip.position(), gradphihat);
+
+        // transform gradients of shape functions to real element
+        const auto S = geo.jacobianInverseTransposed(ip.position());
+        for (auto i : Dune::range(lb.size()))
+          S.mv(gradphihat[i][0],gradphi[i][0]);
+
+        // compute gradient of u
+        Dune::FieldVector<RF,dim> gradu(0.0);
+        for (auto i : Dune::range(lb.size()))
+          gradu.axpy(x[i],gradphi[i][0]);
+
+        // integrate (grad u)*grad phi_i + q(u)*phi_i
+        auto factor = ip.weight()*
+          geo.integrationElement(ip.position());
+        auto qprime_u = param().qprime(w)*u;
+        for (auto i : Dune::range(lb.size()))
+          y[i] += (gradu*gradphi[i][0]+
+                   qprime_u*phihat[i])*factor;
+      }
+    }
+
   };
 
 }