Also enable computing the H_div seminorm.

This is joint work with Zhen Tao and Arezou Ghesmati.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 9d05c3c28195530af4585140754f44408266ef5f..b955d28707576e8438e8dc80e6b269d4f2532638 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -473,6 +473,13 @@ inconvenience this causes.
   <br>
   (Lei Qiao, 2015/04/19)
   </li>
+  
+  <li> New: The VectorTools::integrate_difference() function can now
+  also compute the $H_\text{div}$ seminorm, using the
+  VectorTools::NormType::Hdiv_seminorm argument.
+  <br>
+  (Zhen Tao, Arezou Ghesmati, Wolfgang Bangerth, 2015/04/17)
+  </li>
 
   <li> New: TrilinosWrappers::BlockSparseMatrix now has member functions
   TrilinosWrappers::BlockSparseMatrix::domain_paritioner() and
@@ -480,7 +487,6 @@ inconvenience this causes.
   vector of the underlying block Epetra_Map.
   <br>
   (Matthias Maier, 2015/04/08)
-
   </li>
 
   <li> Fixed: The class SymmetricTensor<2,dim> is now usable also for dim>3.

diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h
index 77c8958e71c70b8581d3245b5bdf657d57ef8fea..a35cbd384a6f04f82f1d669280d28271e7977b58 100644 (file)
--- a/include/deal.II/numerics/vector_tools.h
+++ b/include/deal.II/numerics/vector_tools.h
@@ -355,6 +355,10 @@ namespace VectorTools
      * #L2_norm of the gradient.
      */
     H1_seminorm,
+    /**
+     * #L2_norm of the divergence of a vector field
+     */
+    Hdiv_seminorm,
     /**
      * The square of this norm is the square of the #L2_norm plus the square
      * of the #H1_seminorm.
@@ -1802,7 +1806,12 @@ namespace VectorTools
    * formula only evaluates the two solutions at these particular points,
    * choosing this quadrature formula may indicate an error far smaller than
    * it actually is.
-   * @param[in] norm The norm $X$ shown above that should be computed.
+   * @param[in] norm The norm $X$ shown above that should be computed. If
+   *   the norm is NormType::Hdiv_seminorm, then the finite element on which this
+   *   function is called needs to have at least dim vector components,
+   *   and the divergence will be computed on the first div components.
+   *   This works, for example, on the finite elements used for the
+   *   mixed Laplace (step-20) and the Stokes equations (step-22).
    * @param[in] weight The additional argument @p weight allows to evaluate
    * weighted norms.  The weight function may be scalar, establishing a
    * spatially variable weight in the domain for all components equally. This

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 919cd6db2ecfca5e85dfba4454a6e3b589f58486..77e7605a48e8d84f648081cdf351f5f2d106aac4 100644 (file)
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -6278,6 +6278,24 @@ namespace VectorTools
           diff = std::sqrt(diff);
           break;
 
+        case Hdiv_seminorm:
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              Assert (n_components >= dim,
+                      ExcMessage ("You can only ask for the Hdiv norm for a finite element "
+                                  "with at least 'dim' components. In that case, this function "
+                                  "will take the divergence of the first 'dim' components."));
+              double sum = 0;
+              // take the trace of the derivatives, square it, multiply it
+              // with the weight function
+              for (unsigned int k=0; k<dim; ++k)
+                sum += (data.psi_grads[q][k][k] * data.psi_grads[q][k][k]) *
+                       data.weight_vectors[q](k);
+              diff += sum * fe_values.JxW(q);
+            }
+          diff = std::sqrt(diff);
+          break;
+
         case W1infty_seminorm:
         case W1infty_norm:
         {
@@ -6336,6 +6354,7 @@ namespace VectorTools
         case L2_norm:
         case H1_seminorm:
         case H1_norm:
+        case Hdiv_seminorm:
           exponent = 2.;
           break;
 
@@ -6352,6 +6371,7 @@ namespace VectorTools
       switch (norm)
         {
         case H1_seminorm:
+        case Hdiv_seminorm:
         case W1p_seminorm:
         case W1infty_seminorm:
           update_flags |= UpdateFlags (update_gradients);