From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 26 Jun 2014 19:30:26 +0000 (+0000)
Subject: add grad-div residual
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=43f4fc1bbf342df74cefa161477afaf32b2cab75;p=dealii-svn.git

add grad-div residual

git-svn-id: https://svn.dealii.org/trunk@33090 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/integrators/divergence.h b/deal.II/include/deal.II/integrators/divergence.h
index 79438a8fba..4f35e8f3e1 100644
--- a/deal.II/include/deal.II/integrators/divergence.h
+++ b/deal.II/include/deal.II/integrators/divergence.h
@@ -519,6 +519,45 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * The weak form of the grad-div residual
+     * @f[
+     *  \int_Z \nabla\!\cdot\!u \nabla\!\cdot\!v \,dx
+     * @f]
+     *
+     * @author Guido Kanschat
+     * @date 2014
+     */
+    template <int dim, typename number>
+    void grad_div_residual (
+      Vector<number> &result,
+      const FEValuesBase<dim> &fetest,
+      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
+      const double factor = 1.)
+    {
+      const unsigned int n_dofs = fetest.dofs_per_cell;
+
+      AssertDimension(fetest.get_fe().n_components(), dim);
+      AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
+      
+      for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
+        {
+          const double dx = factor * fetest.JxW(k);
+          for (unsigned int i=0; i<n_dofs; ++i)
+	    {
+	      double dv = 0.;
+	      double du = 0.;
+	      for (unsigned int d=0; d<dim; ++d)
+		{
+		  dv += fetest.shape_grad_component(i,k,d)[d];
+		  du += input[d][k][d];
+		}
+	      
+	      result(i) += dx * du * dv;
+	    }
+        }
+    }
+
     /**
      * The jump of the normal component
      * @f[