new function compute_mean_value in VectorTools

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

diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h
index 3e010209b49fedc569a9a6c04665818f247f4246..3f51c7e1fb41d9b24d60385a6cc7dd87ee799fb5 100644 (file)
--- a/deal.II/deal.II/include/numerics/vectors.h
+++ b/deal.II/deal.II/include/numerics/vectors.h
@@ -504,13 +504,16 @@ class VectorTools
                                     /**
                                      * Mean-value filter for Stokes.
                                      * The pressure in Stokes'
-                                     * equations is determined up to a
-                                     * constant only. This function
-                                     * allows to subtract the mean
-                                     * value of the pressure. It is
-                                     * usually called in a
-                                     * preconditioner and generates
-                                     * updates with mean value zero.
+                                     * equations with only Dirichlet
+                                     * boundaries for the velocities
+                                     * is determined up to a constant
+                                     * only. This function allows to
+                                     * subtract the mean value of the
+                                     * pressure. It is usually called
+                                     * in a preconditioner and
+                                     * generates updates with mean
+                                     * value zero. The mean value is
+                                     * understood in the l1-sense.
                                      *
                                      * Apart from the vector @p{v} to
                                      * operate on, this function
@@ -522,8 +525,39 @@ class VectorTools
                                      */
     static void subtract_mean_value(Vector<double>     &v,
                                    const vector<bool> &p_select);
-
-
+    
+                                    /**
+                                     * Compute the mean value of one
+                                     * component of the solution.
+                                     *
+                                     * This function integrates the
+                                     * chosen component over the
+                                     * whole domain and returns the
+                                     * result.
+                                     *
+                                     * Subtracting this mean value
+                                     * from the node vector does not
+                                     * generally yield the desired
+                                     * result of a finite element
+                                     * function with mean value
+                                     * zero. In fact, it only works
+                                     * for Lagrangian
+                                     * elements. Therefore, it is
+                                     * necessary to compute the mean
+                                     * value and subtract it in the
+                                     * evaluation routine.
+                                     *
+                                     * So far, this is needed only in
+                                     * the error evaluation for
+                                     * Stokes with complete Dirichlet
+                                     * boundaries for the velocities.
+                                     */
+    template <int dim>
+    static double compute_mean_value (const DoFHandler<dim> &dof,
+                                     const Quadrature<dim> &quadrature,
+                                     Vector<double>        &v,
+                                     const unsigned int component);
+    
                                     /**
                                      * Exception
                                      */

diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc
index 38f8712d7db7eaa865d228f7253674aadae0f21d..5c439e143b113d6173d493f3551de5dc71e7e8df 100644 (file)
--- a/deal.II/deal.II/source/numerics/vectors.cc
+++ b/deal.II/deal.II/source/numerics/vectors.cc
@@ -1110,6 +1110,45 @@ VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
 };
 
 
+
+template <int dim>
+double
+VectorTools::compute_mean_value (const DoFHandler<dim> &dof,
+                                const Quadrature<dim> &quadrature,
+                                Vector<double>        &v,
+                                const unsigned int component)
+{
+  Assert (component < dof.get_fe().n_components(),
+         ExcIndexRange(component, 0, dof.get_fe().n_components()));
+  
+  FEValues<dim> fe(dof.get_fe(), quadrature,
+                  UpdateFlags(update_JxW_values
+                              | update_values));
+
+  DoFHandler<dim>::active_cell_iterator c;
+  vector<Vector<double> > values(quadrature.n_quadrature_points,
+                                Vector<double> (dof.get_fe().n_components()));
+  
+  double mean = 0.;
+  double area = 0.;
+                                  // Compute mean value
+  for (c = dof.begin_active();
+       c != dof.end();
+       ++c)
+    {
+      fe.reinit (c);
+      fe.get_function_values(v, values);
+      for (unsigned int k=0; k< quadrature.n_quadrature_points; ++k)
+       {
+         mean += fe.JxW(k) * values[k](component);
+         area += fe.JxW(k);
+       }
+    }
+  return (mean/area);
+}
+
+
+
 // explicit instantiations
 template
 void VectorTools::integrate_difference (const DoFHandler<deal_II_dimension> &,
@@ -1141,6 +1180,14 @@ void VectorTools::interpolate (const DoFHandler<deal_II_dimension> &,
                               const Function<deal_II_dimension>   &,
                               Vector<double>                      &);
 
+template
+double VectorTools::compute_mean_value (const DoFHandler<deal_II_dimension> &dof,
+                                       const Quadrature<deal_II_dimension> &quadrature,
+                                       Vector<double>        &v,
+                                       const unsigned int component);
+
+
+
 // the following two functions are not derived from a template in 1d
 // and thus need no explicit instantiation
 #if deal_II_dimension > 1