Implement jacobian_apply_boundary() analytically

src/nitschenonlinearpoissonfem.hh

@@ -30,9 +30,7 @@ namespace PPS {
   template<typename Param>
   class NitscheNonlinearPoissonFEM :
     public NonlinearPoissonFEM<Param>,
-    public NumericalJacobianBoundary<NitscheNonlinearPoissonFEM<Param> >,
-    public NumericalNonlinearJacobianApplyBoundary
-             <NitscheNonlinearPoissonFEM<Param> >
+    public NumericalJacobianBoundary<NitscheNonlinearPoissonFEM<Param> >
   {
     using Real = typename Param::value_type;
     Real stab_;
@@ -118,6 +116,76 @@ namespace PPS {
              +stab_*(u-g)*phihat[i])*factor;
       }
     }
+
+    //! volume integral depending on test and ansatz functions
+    template<typename Intersection, class LocalBasis, typename X, typename Y>
+    void jacobian_apply_boundary(const Intersection& isect,
+                                 const LocalBasis &lb,
+                                 const X& lp, const X& x, Y& y) const
+    {
+      // evaluate boundary condition type
+      auto localgeo = isect.geometryInInside();
+      auto facecenterlocal =
+        referenceElement(localgeo).position(0,0);
+      bool isdirichlet=this->param().b(isect, facecenterlocal);
+
+      // skip rest if we are _not_ on Dirichlet boundary
+      if (!isdirichlet) return;
+
+      // types & dimension
+      const int dim = Intersection::Entity::dimension;
+      using RF = typename LocalBasis::Traits::RangeFieldType;
+
+      // select quadrature rule
+      auto geo = isect.geometry();
+      const int order = this->incrementorder() + 2*lb.order();
+
+      // geometry of inside cell
+      auto geo_inside = isect.inside().geometry();
+
+      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))
+      {
+        // quadrature point in local coordinates of element
+        auto local = localgeo.global(ip.position());
+
+        // evaluate basis functions
+        lb.evaluateFunction(local, phihat);
+
+        // evaluate u
+        auto u = std::inner_product(begin(x), end(x), begin(phihat), RF(0));
+
+        // evaluate gradient of shape functions
+        lb.evaluateJacobian(local, gradphihat);
+
+        // transform gradients of shape functions to real element
+        const auto S = geo_inside.jacobianInverseTransposed(local);
+        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]);
+
+        // get unit outer normal vector
+        auto n = isect.unitOuterNormal(ip.position());
+
+        // integrate
+        // -(grad u)*n * phi_i - u * (grad phi_i)*n + stab*u*phi_i
+        auto factor = ip.weight()*
+          geo.integrationElement(ip.position());
+        for (auto i : Dune::range(lb.size()))
+          y[i] +=
+            (-(gradu*n)*phihat[i]
+             -u*(gradphi[i][0]*n)
+             +stab_*u*phihat[i])*factor;
+      }
+    }
   };
 
 }