Allow to use finite elements with more than one component.

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

diff --git a/deal.II/deal.II/include/numerics/derivative_approximation.h b/deal.II/deal.II/include/numerics/derivative_approximation.h
index 5efe558ab301298ab181ba316ea088db189d367a..c207ebcab8771ff8deea1eb335ad2a07d8697fb9 100644 (file)
--- a/deal.II/deal.II/include/numerics/derivative_approximation.h
+++ b/deal.II/deal.II/include/numerics/derivative_approximation.h
@@ -173,12 +173,19 @@ class DerivativeApproximation
                                      * receive the cell-wise
                                      * Euclidian norm of the
                                      * approximated gradient.
+                                     *
+                                     * The last parameter denotes the
+                                     * solution component, for which
+                                     * the gradient is to be
+                                     * computed. It defaults to the
+                                     * first component.
                                      */
     template <int dim>
     static void
     approximate_gradient (const DoFHandler<dim> &dof,
                          const Vector<double>  &solution,
-                         Vector<float>         &derivative_norm);
+                         Vector<float>         &derivative_norm,
+                         const unsigned int     component = 0);
 
                                     /**
                                      * This function is the analogue
@@ -195,12 +202,19 @@ class DerivativeApproximation
                                      * derivatives. The spectral norm
                                      * is the matrix norm associated
                                      * to the $l_2$ vector norm.
+                                     *
+                                     * The last parameter denotes the
+                                     * solution component, for which
+                                     * the gradient is to be
+                                     * computed. It defaults to the
+                                     * first component.
                                      */
     template <int dim>
     static void
     approximate_second_derivative (const DoFHandler<dim> &dof,
                                   const Vector<double>  &solution,
-                                  Vector<float>         &derivative_norm);
+                                  Vector<float>         &derivative_norm,
+                                  const unsigned int     component);
     
                                     /**
                                      * Exception
@@ -268,7 +282,8 @@ class DerivativeApproximation
                                          */
        static ProjectedDerivative
        get_projected_derivative (const FEValues<dim>  &fe_values,
-                                 const Vector<double> &solution);
+                                 const Vector<double> &solution,
+                                 const unsigned int    component);
     
                                         /**
                                          * Return the norm of the
@@ -347,7 +362,8 @@ class DerivativeApproximation
                                          */
        static ProjectedDerivative
        get_projected_derivative (const FEValues<dim>  &fe_values,
-                                 const Vector<double> &solution);
+                                 const Vector<double> &solution,
+                                 const unsigned int    component);
        
                                         /**
                                          * Return the norm of the
@@ -404,11 +420,17 @@ class DerivativeApproximation
                                      * administration that is
                                      * independent of the actual
                                      * derivative to be computed.
+                                     *
+                                     * The @p{component} argument
+                                     * denotes which component of the
+                                     * solution vector we are to work
+                                     * on.
                                      */
     template <class DerivativeDescription, int dim>
     static void
     approximate_derivative (const DoFHandler<dim> &dof,
                            const Vector<double>  &solution,
+                           const unsigned int     component,
                            Vector<float>         &derivative_norm);
 
                                     /**
@@ -421,8 +443,9 @@ class DerivativeApproximation
     static void
     approximate (const DoFHandler<dim> &dof,
                 const Vector<double>  &solution,
+                const unsigned int     component,
                 const IndexInterval   &index_interval,
-                            Vector<float>         &derivative_norm);    
+                Vector<float>         &derivative_norm);    
 };
 
 

diff --git a/deal.II/deal.II/source/numerics/derivative_approximation.cc b/deal.II/deal.II/source/numerics/derivative_approximation.cc
index 78433a828b6d217c82f2bfb42f6d0a6e4f356ec2..21185db45136bb261852416897587ff3c744e15c 100644 (file)
--- a/deal.II/deal.II/source/numerics/derivative_approximation.cc
+++ b/deal.II/deal.II/source/numerics/derivative_approximation.cc
@@ -42,11 +42,22 @@ inline
 typename DerivativeApproximation::Gradient<dim>::ProjectedDerivative
 DerivativeApproximation::Gradient<dim>::
 get_projected_derivative (const FEValues<dim>  &fe_values,
-                         const Vector<double> &solution) 
+                         const Vector<double> &solution,
+                         const unsigned int    component)
 {
-  vector<ProjectedDerivative> values (1);
-  fe_values.get_function_values (solution, values);
-  return values[0];
+  if (fe_values.get_fe().n_components() == 1)
+    {
+      vector<ProjectedDerivative> values (1);
+      fe_values.get_function_values (solution, values);
+      return values[0];
+    }
+  else
+    {
+      vector<Vector<double> > values
+       (1, Vector<double>(fe_values.get_fe().n_components()));
+      fe_values.get_function_values (solution, values);
+      return values[0](component);
+    };
 };
 
 
@@ -79,11 +90,22 @@ inline
 typename DerivativeApproximation::SecondDerivative<dim>::ProjectedDerivative
 DerivativeApproximation::SecondDerivative<dim>::
 get_projected_derivative (const FEValues<dim>  &fe_values,
-                         const Vector<double> &solution) 
+                         const Vector<double> &solution,
+                         const unsigned int    component)
 {
-  vector<ProjectedDerivative> values (1);
-  fe_values.get_function_grads (solution, values);
-  return values[0];
+  if (fe_values.get_fe().n_components() == 1)
+    {
+      vector<ProjectedDerivative> values (1);
+      fe_values.get_function_grads (solution, values);
+      return values[0];
+    }
+  else
+    {
+      vector<vector<ProjectedDerivative> > values
+       (1, vector<ProjectedDerivative>(fe_values.get_fe().n_components()));
+      fe_values.get_function_grads (solution, values);
+      return values[0][component];
+    };
 };
 
 
@@ -169,10 +191,12 @@ void
 DerivativeApproximation::
 approximate_gradient (const DoFHandler<dim> &dof_handler,
                      const Vector<double>  &solution,
-                     Vector<float>         &derivative_norm)
+                     Vector<float>         &derivative_norm,
+                     const unsigned int     component)
 {
   approximate_derivative<Gradient<dim>,dim> (dof_handler,
                                             solution,
+                                            component,
                                             derivative_norm);
 };
 
@@ -183,10 +207,12 @@ void
 DerivativeApproximation::
 approximate_second_derivative (const DoFHandler<dim> &dof_handler,
                               const Vector<double>  &solution,
-                              Vector<float>         &derivative_norm)
+                              Vector<float>         &derivative_norm,
+                              const unsigned int     component)
 {
   approximate_derivative<SecondDerivative<dim>,dim> (dof_handler,
                                                     solution,
+                                                    component,
                                                     derivative_norm);
 };
 
@@ -197,6 +223,7 @@ void
 DerivativeApproximation::
 approximate_derivative (const DoFHandler<dim> &dof_handler,
                        const Vector<double>  &solution,
+                       const unsigned int     component,
                        Vector<float>         &derivative_norm)
 {
   Assert (derivative_norm.size() == dof_handler.get_tria().n_active_cells(),
@@ -215,7 +242,7 @@ approximate_derivative (const DoFHandler<dim> &dof_handler,
                    Threads::encapsulate
                    (&DerivativeApproximation::
                     template approximate<DerivativeDescription,dim>)
-                   .collect_args (dof_handler, solution,
+                   .collect_args (dof_handler, solution, component,
                                   index_intervals[i],
                                   derivative_norm));
   thread_manager.wait ();
@@ -227,6 +254,7 @@ template <class DerivativeDescription, int dim>
 void 
 DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
                                      const Vector<double>  &solution,
+                                     const unsigned int     component,
                                      const IndexInterval   &index_interval,
                                      Vector<float>         &derivative_norm)
 {
@@ -278,7 +306,8 @@ DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
       const typename DerivativeDescription::ProjectedDerivative
        this_midpoint_value
        = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
-                                                   solution);
+                                                          solution,
+                                                          component);
                                       // ...and the place where it lives
       const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
 
@@ -364,7 +393,7 @@ DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
          const typename DerivativeDescription::ProjectedDerivative
            neighbor_midpoint_value
            = DerivativeDescription::get_projected_derivative (fe_midpoint_value,
-                                                       solution);
+                                                              solution, component);
          
                                           // ...and the place where it lives
          const Point<dim>
@@ -446,14 +475,16 @@ void
 DerivativeApproximation::
 approximate_gradient (const DoFHandler<deal_II_dimension> &dof_handler,
                      const Vector<double>  &solution,
-                     Vector<float>         &derivative_norm);
+                     Vector<float>         &derivative_norm,
+                     const unsigned int     component);
 
 template
 void 
 DerivativeApproximation::
 approximate_second_derivative (const DoFHandler<deal_II_dimension> &dof_handler,
                               const Vector<double>  &solution,
-                              Vector<float>         &derivative_norm);
+                              Vector<float>         &derivative_norm,
+                              const unsigned int     component);