]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement an entirely new scheme to deal with vector-valued problems. Also provide...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 23 Feb 2008 19:14:34 +0000 (19:14 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 23 Feb 2008 19:14:34 +0000 (19:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@15757 0785d39b-7218-0410-832d-ea1e28bc413d

22 files changed:
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_system.h
deal.II/deal.II/include/fe/fe_values.h
deal.II/doc/doxygen/headers/fe.h
deal.II/doc/doxygen/headers/glossary.h
deal.II/doc/doxygen/headers/vector_valued.h [new file with mode: 0644]
deal.II/doc/news/changes.h
deal.II/examples/step-20/doc/intro.dox
deal.II/examples/step-20/step-20.cc
deal.II/examples/step-21/doc/intro.dox
deal.II/examples/step-21/step-21.cc
deal.II/examples/step-22/doc/intro.dox
deal.II/examples/step-22/step-22.cc
deal.II/examples/step-29/doc/intro.dox
deal.II/examples/step-31/doc/intro.dox
deal.II/examples/step-31/step-31.cc
deal.II/examples/step-8/doc/intro.dox
tests/deal.II/Makefile
tests/deal.II/fe_values_view_01.cc [new file with mode: 0644]
tests/deal.II/fe_values_view_02.cc [new file with mode: 0644]
tests/deal.II/fe_values_view_03.cc [new file with mode: 0644]
tests/deal.II/fe_values_view_04.cc [new file with mode: 0644]

index 4c8006bdf3694bd3442823683bf73a9ed62857c4..6d93d2dab7a66cc029ebed1912e2668704d6563e 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -2684,7 +2684,12 @@ FiniteElement<dim>::is_primitive (const unsigned int i) const
                                   // std::vector<bool> is) is
                                   // probably more expensive than
                                   // just comparing against 1
-  return (n_nonzero_components_table[i] == 1);
+                                  //
+                                  // for good measure, short circuit the test
+                                  // if the entire FE is primitive
+  return ((cached_primitivity == true)
+         ||
+         (n_nonzero_components_table[i] == 1));
 }
 
 
index a8abf92fd789d0d366795ce8c5affa0ed6b3b530..78549296bb3cba7f7d1dea3f4560a309ab0b9e50 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 DEAL_II_NAMESPACE_OPEN
 
 /**
- * This class provides an interface to group several elements together
- * into one. To the outside world, the resulting object looks just
- * like a usual finite element object, which is composed of several
- * other finite elements that are possibly of different type. The
- * result is then a vector-valued finite element.
+ * This class provides an interface to group several elements together into
+ * one. To the outside world, the resulting object looks just like a usual
+ * finite element object, which is composed of several other finite elements
+ * that are possibly of different type. The result is then a vector-valued
+ * finite element. Vector valued elements are discussed in a number of
+ * tutorial programs, for example @ref step_8 "step-8", @ref step_20
+ * "step-20", @ref step_21 "step-21", and in particular in the @ref vector_valued
+ * module.
  *
  * The overall numbering of degrees of freedom is as follows: for each
  * subobject (vertex, line, quad, or hex), the degrees of freedom are
@@ -70,7 +73,7 @@ DEAL_II_NAMESPACE_OPEN
  * multiplicity. The number of blocks of the system is simply the sum
  * of all multiplicities.
  *
- * @ingroup febase fe
+ * @ingroup febase fe vector_valued
  * 
  * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2002, 2003, 2006, Ralf Hartmann 2001.
  */
index f2e9e3421065b90fc5bdea735dd6ef2bed663b32..2f81caf1f0e1853b5ead4c5534691e8209ab7401 100644 (file)
@@ -19,6 +19,7 @@
 #include <base/subscriptor.h>
 #include <base/point.h>
 #include <base/tensor.h>
+#include <base/symmetric_tensor.h>
 #include <base/vector_slice.h>
 #include <base/quadrature.h>
 #include <base/table.h>
 DEAL_II_NAMESPACE_OPEN
 
 template <int dim> class Quadrature;
+template <int dim> class FEValuesBase;
+
+
+
+/**
+ * A namespace in which we declare "extractors", i.e. classes that when used
+ * as subscripts in operator[] expressions on FEValues, FEFaceValues, and
+ * FESubfaceValues objects extract certain components of a vector-valued
+ * element. There are extractors for single scalar components as well as
+ * vector components consisting of <code>dim</code> elements.
+ *
+ * See the description of the @ref vector_valued module for examples how to
+ * use the features of this namespace.
+ *
+ * @ingroup feaccess vector_valued
+ */
+namespace FEValuesExtractors
+{
+                                  /**
+                                   * Extractor for a single scalar component
+                                   * of a vector-valued element. The concept
+                                   * of extractors is defined in the
+                                   * documentation of the namespace
+                                   * FEValuesExtractors and in the @ref
+                                   * vector_valued module.
+                                   *
+                                   * @ingroup feaccess vector_valued
+                                   */
+  struct Scalar
+  {
+                                      /**
+                                       * The selected scalar component of the
+                                       * vector.
+                                       */
+      const unsigned int component;
+
+                                      /**
+                                       * Constructor. Take the selected
+                                       * vector component as argument.
+                                       */
+      Scalar (const unsigned int component);
+  };
+
+  
+                                  /**
+                                   * Extractor for a vector of
+                                   * <code>dim</code> components of a
+                                   * vector-valued element. The value of
+                                   * <code>dim</code> is defined by the
+                                   * FEValues object the extractor is applied
+                                   * to.
+                                   *
+                                   * The concept of
+                                   * extractors is defined in the
+                                   * documentation of the namespace
+                                   * FEValuesExtractors and in the @ref
+                                   * vector_valued module.
+                                   *
+                                   * @ingroup feaccess vector_valued
+                                   */
+  struct Vector
+  {
+                                      /**
+                                       * The first component of the vector view.
+                                       */
+      const unsigned int first_vector_component;
+
+                                      /**
+                                       * Constructor. Take the first
+                                       * component of the selected vector
+                                       * inside the FEValues object as
+                                       * argument.
+                                       */
+      Vector (const unsigned int first_vector_component);
+  };
+}
+
+
+
+/**
+ * A namespace for "views" on a FEValues, FEFaceValues, or FESubfaceValues
+ * object. A view represents only a certain part of the whole: whereas the
+ * FEValues object represents <i>all</i> values, gradients, or second
+ * derivatives of all components of a vector-valued element, views restrict
+ * the attention to only a single component or a subset of components.
+ *
+ * There are classes that present views for single scalar components as well
+ * as vector components consisting of <code>dim</code> elements.
+ *
+ * See the description of the @ref vector_valued module for examples how to
+ * use the features of this namespace.
+ *
+ * @ingroup feaccess vector_valued
+ */
+namespace FEValuesViews
+{
+                                  /**
+                                   * A class representing a view to a single
+                                   * scalar component of a possibly
+                                   * vector-valued finite element. Views are
+                                   * discussed in the @ref vector_valued module.
+                                   *
+                                   * @ingroup feaccess vector_valued
+                                   */
+  template <int dim>
+  class Scalar
+  {
+    public:
+                                      /**
+                                       * A typedef for the data type of
+                                       * values of the view this class
+                                       * represents. Since we deal with a
+                                       * single components, the value type is
+                                       * a scalar double.
+                                       */
+      typedef double        value_type;
+
+                                      /**
+                                       * A typedef for the type of gradients
+                                       * of the view this class
+                                       * represents. Here, for a scalar
+                                       * component of the finite element, the
+                                       * gradient is a
+                                       * <code>Tensor@<1,dim@></code>.
+                                       */
+      typedef Tensor<1,dim> gradient_type;
+
+                                      /**
+                                       * A typedef for the type of second
+                                       * derivatives of the view this class
+                                       * represents. Here, for a scalar
+                                       * component of the finite element, the
+                                       * Hessian is a
+                                       * <code>Tensor@<2,dim@></code>.
+                                       */
+      typedef Tensor<2,dim> hessian_type;
+
+                                      /**
+                                       * Constructor for an object that
+                                       * represents a single scalar component
+                                       * of a FEValuesBase object (or of one
+                                       * of the classes derived from
+                                       * FEValuesBase).
+                                       */
+      Scalar (const FEValuesBase<dim> &fe_values_base,
+             const unsigned int       component);
+
+                                      /**
+                                       * Return the value of the vector
+                                       * component selected by this view, for
+                                       * the shape function and quadrature
+                                       * point selected by the arguments.
+                                       */
+      value_type
+      value (const unsigned int shape_function,
+            const unsigned int q_point) const;
+
+                                      /**
+                                       * Return the gradient (a tensor of
+                                       * rank 1) of the vector component
+                                       * selected by this view, for the shape
+                                       * function and quadrature point
+                                       * selected by the arguments.
+                                       */
+      gradient_type
+      gradient (const unsigned int shape_function,
+               const unsigned int q_point) const;
+
+                                      /**
+                                       * Return the Hessian (the tensor of
+                                       * rank 2 of all second derivatives) of
+                                       * the vector component selected by
+                                       * this view, for the shape function
+                                       * and quadrature point selected by the
+                                       * arguments.
+                                       */
+      hessian_type
+      hessian (const unsigned int shape_function,
+              const unsigned int q_point) const;
+      
+      
+    private:
+                                      /**
+                                       * A pointer to the FEValuesBase object
+                                       * we operator on.
+                                       */
+      const SmartPointer<const FEValuesBase<dim> > fe_values;
+
+                                      /**
+                                       * The single scalar component this
+                                       * view represents of the FEValuesBase
+                                       * object.
+                                       */
+      const unsigned int component;
+  };
+  
 
-//TODO: Add access to mapping values to FEValuesBase
-//TODO: Several FEValuesBase of a system should share Mapping
 
+                                  /**
+                                   * A class representing a view to a set of
+                                   * <code>dim</code> components forming a
+                                   * vector part of a vector-valued finite
+                                   * element. Views are discussed in the
+                                   * @ref vector_valued module.
+                                   *
+                                   * @ingroup feaccess vector_valued
+                                   */
+  template <int dim>
+  class Vector
+  {
+    public:
+                                      /**
+                                       * A typedef for the data type of
+                                       * values of the view this class
+                                       * represents. Since we deal with a set
+                                       * of <code>dim</code> components, the
+                                       * value type is a Tensor<1,dim>.
+                                       */
+      typedef Tensor<1,dim>          value_type;
+
+                                      /**
+                                       * A typedef for the type of gradients
+                                       * of the view this class
+                                       * represents. Here, for a set of
+                                       * <code>dim</code> components of the
+                                       * finite element, the gradient is a
+                                       * <code>Tensor@<2,dim@></code>.
+                                       */
+      typedef Tensor<2,dim>          gradient_type;
+
+                                      /**
+                                       * A typedef for the type of
+                                       * symmetrized gradients of the view
+                                       * this class represents. Here, for a
+                                       * set of <code>dim</code> components
+                                       * of the finite element, the
+                                       * symmetrized gradient is a
+                                       * <code>SymmetricTensor@<2,dim@></code>.
+                                       */
+      typedef SymmetricTensor<2,dim> symmetric_gradient_type;
+
+                                      /**
+                                       * A typedef for the type of the
+                                       * divergence of the view this class
+                                       * represents. Here, for a set of
+                                       * <code>dim</code> components of the
+                                       * finite element, the divergence of
+                                       * course is a scalar.
+                                       */
+      typedef double                 divergence_type;
+      
+                                      /**
+                                       * A typedef for the type of second
+                                       * derivatives of the view this class
+                                       * represents. Here, for a set of
+                                       * <code>dim</code> components of the
+                                       * finite element, the Hessian is a
+                                       * <code>Tensor@<3,dim@></code>.
+                                       */
+      typedef Tensor<3,dim>          hessian_type;
+
+                                      /**
+                                       * Constructor for an object that
+                                       * represents dim components of a
+                                       * FEValuesBase object (or of one of
+                                       * the classes derived from
+                                       * FEValuesBase), representing a
+                                       * vector-valued variable.
+                                       *
+                                       * The second argument denotes the
+                                       * index of the first component of the
+                                       * selected vector.
+                                       */
+      Vector (const FEValuesBase<dim> &fe_values_base,
+             const unsigned int first_vector_component);
+             
+                                      /**
+                                       * Return the value of the vector
+                                       * components selected by this view,
+                                       * for the shape function and
+                                       * quadrature point selected by the
+                                       * arguments. Here, since the view
+                                       * represents a vector-valued part of
+                                       * the FEValues object with
+                                       * <code>dim</code> components, the
+                                       * return type is a tensor of rank 1
+                                       * with <code>dim</code> components.
+                                       */
+      value_type
+      value (const unsigned int shape_function,
+            const unsigned int q_point) const;
+
+                                      /**
+                                       * Return the gradient (a tensor of
+                                       * rank 2) of the vector component
+                                       * selected by this view, for the shape
+                                       * function and quadrature point
+                                       * selected by the arguments.
+                                       */
+      gradient_type
+      gradient (const unsigned int shape_function,
+               const unsigned int q_point) const;
+
+                                      /**
+                                       * Return the symmetric gradient (a
+                                       * symmetric tensor of rank 2) of the
+                                       * vector component selected by this
+                                       * view, for the shape function and
+                                       * quadrature point selected by the
+                                       * arguments.
+                                       *
+                                       * The symmetric gradient is defined as
+                                       * $\frac 12 [(\nabla \phi_i(x_q)) +
+                                       * (\nabla \phi_i(x_q))^T]$, where
+                                       * $\phi_i$ represents the
+                                       * <code>dim</code> components selected
+                                       * from the FEValuesBase object, and
+                                       * $x_q$ is the location of the $q$th
+                                       * quadrature point.
+                                       */
+      symmetric_gradient_type
+      symmetric_gradient (const unsigned int shape_function,
+                         const unsigned int q_point) const;
+
+                                      /**
+                                       * Return the scalar divergence of
+                                       * the vector components selected by
+                                       * this view, for the shape function
+                                       * and quadrature point selected by the
+                                       * arguments.
+                                       */
+      divergence_type
+      divergence (const unsigned int shape_function,
+                 const unsigned int q_point) const;
+      
+                                      /**
+                                       * Return the Hessian (the tensor of
+                                       * rank 2 of all second derivatives) of
+                                       * the vector components selected by
+                                       * this view, for the shape function
+                                       * and quadrature point selected by the
+                                       * arguments.
+                                       */
+      hessian_type
+      hessian (const unsigned int shape_function,
+              const unsigned int q_point) const;
+      
+    private:
+                                      /**
+                                       * A pointer to the FEValuesBase object
+                                       * we operator on.
+                                       */
+      const SmartPointer<const FEValuesBase<dim> > fe_values;
+
+                                      /**
+                                       * The first component of the vector
+                                       * this view represents of the
+                                       * FEValuesBase object.
+                                       */
+      const unsigned int first_vector_component;
+  };
+}
+
+
+
+  
 /*!@addtogroup feaccess */
 /*@{*/
 
+
+
+//TODO: Add access to mapping values to FEValuesBase
+//TODO: Several FEValuesBase of a system should share Mapping
+
 /**
  * Contains all data vectors for FEValues.
  * This class has been extracted from FEValuesBase to be handed
@@ -455,6 +823,39 @@ class FEValuesBase : protected FEValuesData<dim>,
                                      * Destructor.
                                      */
     ~FEValuesBase ();
+
+                                    /// @name Extractors Methods to extract individual components
+                                    //@{
+
+                                    /**
+                                     * Create a view of the current FEValues
+                                     * object that represents a particular
+                                     * scalar component of the possibly
+                                     * vector-valued finite element. The
+                                     * concept of views is explained in the
+                                     * documentation of the namespace
+                                     * FEValuesViews and in particular
+                                     * in the @ref vector_valued module.
+                                     */
+    FEValuesViews::Scalar<dim>
+    operator[] (const FEValuesExtractors::Scalar &scalar) const;
+
+                                    /**
+                                     * Create a view of the current FEValues
+                                     * object that represents a set of
+                                     * <code>dim</code> scalar components
+                                     * (i.e. a vector) of the vector-valued
+                                     * finite element. The concept of views
+                                     * is explained in the documentation of
+                                     * the namespace FEValuesViews and in particular
+                                     * in the @ref vector_valued module.
+                                     */
+    FEValuesViews::Vector<dim>
+    operator[] (const FEValuesExtractors::Vector &vector) const;
+
+                                    //@}
+
+    
                                     /// @name ShapeAccess Access to shape function values
                                     //@{
     
@@ -1793,7 +2194,15 @@ class FEValuesBase : protected FEValuesData<dim>,
                                       * copyable, we make it private,
                                       * and also do not implement it.
                                       */
-    FEValuesBase & operator= (const FEValuesBase &);    
+    FEValuesBase & operator= (const FEValuesBase &);
+
+                                    /**
+                                     * Make the view classes friends of this
+                                     * class, since they access internal
+                                     * data.
+                                     */
+    template <int> friend class FEValuesViews::Scalar;
+    template <int> friend class FEValuesViews::Vector;
 };
 
 
@@ -2487,15 +2896,724 @@ class FESubfaceValues : public FEFaceValuesBase<dim>
 /*@}*/
 
 #ifndef DOXYGEN
+
+
+/*------------------------ Inline functions: namespace FEValuesExtractors --------*/
+
+namespace FEValuesExtractors
+{
+  inline
+  Scalar::Scalar (const unsigned int component)
+                 :
+                 component (component)
+  {}
+
+
+
+  inline
+  Vector::Vector (const unsigned int first_vector_component)
+                 :
+                 first_vector_component (first_vector_component)
+  {}
+}
+
+
+/*------------------------ Inline functions: namespace FEValuesViews --------*/
+
+namespace FEValuesViews
+{
+  template <int dim>
+  inline
+  Scalar<dim>::Scalar (const FEValuesBase<dim> &fe_values,
+                      const unsigned int       component)
+                 :
+                 fe_values (&fe_values),
+                 component (component)
+  {
+    Assert (component < this->fe_values->fe->n_components(),
+           ExcIndexRange(component, 0, this->fe_values->fe->n_components()));
+  }
+
+
+
+  template <int dim>
+  inline
+  typename Scalar<dim>::value_type
+  Scalar<dim>::value (const unsigned int shape_function,
+                     const unsigned int q_point) const
+  {
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_values,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+                                    // an adaptation of the
+                                    // FEValuesBase::shape_value_component
+                                    // function except that here we know the
+                                    // component as fixed. see the
+                                    // comments there
+                                    //
+                                    // we can do away with some of the
+                                    // assertions since they are already
+                                    // taken care of in the constructor of
+                                    // this class
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       if (component ==
+           fe_values->fe->system_to_component_index(shape_function).first)
+         return fe_values->shape_values(fe_values->shape_function_to_row_table[shape_function],q_point);
+       else
+         return 0;
+      }
+    else
+      {
+       if (fe_values->fe->get_nonzero_components(shape_function)[component] == false)
+         return 0;
+
+       const unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin()+
+                            component,
+                            true));
+       return fe_values->shape_values(row, q_point);
+      }
+  }
+
+
+
+  template <int dim>
+  inline
+  typename Scalar<dim>::gradient_type
+  Scalar<dim>::gradient (const unsigned int shape_function,
+                        const unsigned int q_point) const
+  {
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_gradients,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+                                    // an adaptation of the
+                                    // FEValuesBase::shape_grad_component
+                                    // function except that here we know the
+                                    // component as fixed. see the
+                                    // comments there
+                                    //
+                                    // we can do away with the assertions
+                                    // since they are already taken care of
+                                    // in the constructor of this class
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       if (component ==
+           fe_values->fe->system_to_component_index(shape_function).first)
+         return fe_values->shape_gradients[fe_values->
+                                           shape_function_to_row_table[shape_function]][q_point];
+       else
+         return gradient_type();
+      }
+    else
+      {
+       if (fe_values->fe->get_nonzero_components(shape_function)[component] == false)
+         return gradient_type();
+
+       const unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin()+
+                            component,
+                            true));
+       return fe_values->shape_gradients[row][q_point];
+      }
+  }
+
+
+
+  template <int dim>
+  inline
+  typename Scalar<dim>::hessian_type
+  Scalar<dim>::hessian (const unsigned int shape_function,
+                       const unsigned int q_point) const
+  {
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_hessians,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+                                    // an adaptation of the
+                                    // FEValuesBase::shape_grad_component
+                                    // function except that here we know the
+                                    // component as fixed. see the
+                                    // comments there
+                                    //
+                                    // we can do away with the assertions
+                                    // since they are already taken care of
+                                    // in the constructor of this class
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       if (component ==
+           fe_values->fe->system_to_component_index(shape_function).first)
+         return fe_values->shape_hessians[fe_values->
+                                          shape_function_to_row_table[shape_function]][q_point];
+       else
+         return hessian_type();
+      }
+    else
+      {
+       if (fe_values->fe->get_nonzero_components(shape_function)[component] == false)
+         return hessian_type();
+
+       const unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin()+
+                            component,
+                            true));
+       return fe_values->shape_hessians[row][q_point];
+      }
+  }
+
+
+
+  template <int dim>
+  inline
+  Vector<dim>::Vector (const FEValuesBase<dim> &fe_values,
+                      const unsigned int       first_vector_component)
+                 :
+                 fe_values (&fe_values),
+                 first_vector_component (first_vector_component)
+  {
+    Assert (first_vector_component+dim-1 < this->fe_values->fe->n_components(),
+           ExcIndexRange(first_vector_component+dim-1, 0,
+                         this->fe_values->fe->n_components()));
+  }
+
+
+
+  template <int dim>
+  inline
+  typename Vector<dim>::value_type
+  Vector<dim>::value (const unsigned int shape_function,
+                     const unsigned int q_point) const
+  {
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_values,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+                                    // compared to the scalar case above, we
+                                    // can save some work here because we
+                                    // know that we are querying a contiguous
+                                    // range of components
+    value_type return_value;
+    
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+                                        // if this is a primitive shape
+                                        // function then at most one element
+                                        // of the output vector is
+                                        // nonzero. find out if indeed one is
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if ((nonzero_component >= first_vector_component)
+           &&
+           (nonzero_component < first_vector_component + dim))
+         return_value[nonzero_component - first_vector_component]
+           = fe_values->shape_values(fe_values->
+                                     shape_function_to_row_table[shape_function],
+                                     q_point);
+      }
+    else
+      {
+       unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin() +
+                            first_vector_component,
+                            true));
+       for (unsigned int d=0; d<dim; ++d)
+         if (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
+             true)
+           {
+             return_value[d] = fe_values->shape_values(row, q_point);
+
+             if ((d != dim-1)
+                 &&
+                 (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d]
+                  == true))
+               ++row;
+           }
+      }
+
+    return return_value;
+  }
+
+  
+
+  template <int dim>
+  inline
+  typename Vector<dim>::gradient_type
+  Vector<dim>::gradient (const unsigned int shape_function,
+                        const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above    
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_gradients,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+    gradient_type return_value;
+    
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if ((nonzero_component >= first_vector_component)
+           &&
+           (nonzero_component < first_vector_component + dim))
+         return_value[nonzero_component - first_vector_component]
+           = fe_values->shape_gradients[fe_values->
+                                        shape_function_to_row_table[shape_function]][q_point];
+      }
+    else
+      {
+       unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin() +
+                            first_vector_component,
+                            true));
+       for (unsigned int d=0; d<dim; ++d)
+         if (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
+             true)
+           {
+             return_value[d] = fe_values->shape_gradients[row][q_point];
+
+             if ((d != dim-1)
+                 &&
+                 (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d]
+                  == true))
+               ++row;
+           }
+      }
+
+    return return_value;
+  }
+
+
+
+  template <int dim>
+  inline
+  typename Vector<dim>::divergence_type
+  Vector<dim>::divergence (const unsigned int shape_function,
+                          const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above    
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_gradients,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if ((nonzero_component >= first_vector_component)
+           &&
+           (nonzero_component < first_vector_component + dim))
+         return
+           fe_values->shape_gradients[fe_values->
+                                      shape_function_to_row_table[shape_function]][q_point]
+           [nonzero_component - first_vector_component];
+       else
+         return 0;
+      }
+    else
+      {
+       unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin() +
+                            first_vector_component,
+                            true));
+
+       double div = 0;
+       for (unsigned int d=0; d<dim; ++d)
+         if (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
+             true)
+           {
+             div += fe_values->shape_gradients[row][q_point][d];
+
+             if ((d != dim-1)
+                 &&
+                 (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d]
+                  == true))
+               ++row;
+           }
+       return div;
+      }
+  }
+  
+
+
+  template <int dim>
+  inline
+  typename Vector<dim>::hessian_type
+  Vector<dim>::hessian (const unsigned int shape_function,
+                       const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above    
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_hessians,
+           typename FEValuesBase<dim>::ExcAccessToUninitializedField());    
+
+    hessian_type return_value;
+    
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if ((nonzero_component >= first_vector_component)
+           &&
+           (nonzero_component < first_vector_component + dim))
+         return_value[nonzero_component - first_vector_component]
+           = fe_values->shape_hessians[fe_values->
+                                       shape_function_to_row_table[shape_function]][q_point];
+      }
+    else
+      {
+       unsigned int
+         row = (fe_values->shape_function_to_row_table[shape_function]
+                +
+                std::count (fe_values->fe->get_nonzero_components(shape_function).begin(),
+                            fe_values->fe->get_nonzero_components(shape_function).begin() +
+                            first_vector_component,
+                            true));
+       for (unsigned int d=0; d<dim; ++d)
+         if (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
+             true)
+           {
+             return_value[d] = fe_values->shape_hessians[row][q_point];
+
+             if ((d != dim-1)
+                 &&
+                 (fe_values->fe->get_nonzero_components(shape_function)[first_vector_component+d]
+                  == true))
+               ++row;
+           }
+      }
+
+    return return_value;
+  }
+
+
+
+                                  // we duplicate the following function for
+                                  // each dimension since we need to
+                                  // initialize a few arrays of dimension
+                                  // dependent size. even though this happens
+                                  // in a switch(dim) clause, we get warnings
+                                  // from the compiler if we don't separate
+                                  // things into different functions
+  template <>
+  inline
+  Vector<1>::symmetric_gradient_type
+  Vector<1>::symmetric_gradient (const unsigned int shape_function,
+                                const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above    
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_gradients,
+           FEValuesBase<1>::ExcAccessToUninitializedField());    
+
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if (nonzero_component == first_vector_component)
+         {
+                                            // first get the one component of
+                                            // the nonsymmetrized gradient
+                                            // that is not zero
+           const Tensor<1,1> grad
+             = fe_values->shape_gradients[fe_values->
+                                          shape_function_to_row_table[shape_function]][q_point];
+
+                                            // then form a symmetric tensor
+                                            // out of it. note that access to
+                                            // individual elements of a
+                                            // SymmetricTensor object is
+                                            // fairly slow. if we implemented
+                                            // the following in the naive
+                                            // way, we would therefore incur
+                                            // a very significant penalty: a
+                                            // preliminary version of the
+                                            // Stokes tutorial program would
+                                            // slow down from 17 to 23
+                                            // seconds because of this single
+                                            // function!
+                                            //
+                                            // consequently, we try to be a
+                                            // bit smarter by already laying
+                                            // out the data in the right
+                                            // format and creating a
+                                            // symmetric tensor of it
+           const double array[symmetric_gradient_type::n_independent_components]
+             = { grad[0] };
+           return symmetric_gradient_type(array);
+         }
+       else
+         return symmetric_gradient_type();
+      }
+    else
+      {
+       return symmetrize (gradient(shape_function, q_point));
+      }
+  }
+
+
+
+  template <>
+  inline
+  Vector<2>::symmetric_gradient_type
+  Vector<2>::symmetric_gradient (const unsigned int shape_function,
+                                const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above    
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_gradients,
+           FEValuesBase<2>::ExcAccessToUninitializedField());    
+
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if ((nonzero_component >= first_vector_component)
+           &&
+           (nonzero_component < first_vector_component + 2))
+         {
+                                            // first get the one component of
+                                            // the nonsymmetrized gradient
+                                            // that is not zero
+           const Tensor<1,2> grad
+             = fe_values->shape_gradients[fe_values->
+                                          shape_function_to_row_table[shape_function]][q_point];
+
+                                            // then form a symmetric tensor
+                                            // out of it. note that access to
+                                            // individual elements of a
+                                            // SymmetricTensor object is
+                                            // fairly slow. if we implemented
+                                            // the following in the naive
+                                            // way, we would therefore incur
+                                            // a very significant penalty: a
+                                            // preliminary version of the
+                                            // Stokes tutorial program would
+                                            // slow down from 17 to 23
+                                            // seconds because of this single
+                                            // function!
+                                            //
+                                            // consequently, we try to be a
+                                            // bit smarter by already laying
+                                            // out the data in the right
+                                            // format and creating a
+                                            // symmetric tensor of it
+           switch (nonzero_component - first_vector_component)
+             {
+               case 0:
+               {
+                 const double array[symmetric_gradient_type::n_independent_components]
+                   = { grad[0], 0, grad[1]/2 };
+                 return symmetric_gradient_type(array);
+               }
+
+               case 1:
+               {
+                 const double array[symmetric_gradient_type::n_independent_components]
+                   = { 0, grad[1], grad[0]/2 };
+                 return symmetric_gradient_type(array);
+               }
+
+               default:
+                     Assert (false, ExcInternalError());
+                     return symmetric_gradient_type();
+             }
+         }
+       else
+         return symmetric_gradient_type();
+      }
+    else
+      {
+       return symmetrize (gradient(shape_function, q_point));
+      }
+  }
+
+
+
+  template <>
+  inline
+  Vector<3>::symmetric_gradient_type
+  Vector<3>::symmetric_gradient (const unsigned int shape_function,
+                                const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above    
+    Assert (shape_function < fe_values->fe->dofs_per_cell,
+           ExcIndexRange (shape_function, 0, fe_values->fe->dofs_per_cell));
+    Assert (fe_values->update_flags & update_gradients,
+           FEValuesBase<3>::ExcAccessToUninitializedField());    
+
+    if (fe_values->fe->is_primitive() ||
+       fe_values->fe->is_primitive(shape_function))
+      {
+       const unsigned int
+         nonzero_component
+         = fe_values->fe->system_to_component_index(shape_function).first;
+       
+       if ((nonzero_component >= first_vector_component)
+           &&
+           (nonzero_component < first_vector_component + 3))
+         {
+                                            // first get the one component of
+                                            // the nonsymmetrized gradient
+                                            // that is not zero
+           const Tensor<1,3> grad
+             = fe_values->shape_gradients[fe_values->
+                                          shape_function_to_row_table[shape_function]][q_point];
+
+                                            // then form a symmetric tensor
+                                            // out of it. note that access to
+                                            // individual elements of a
+                                            // SymmetricTensor object is
+                                            // fairly slow. if we implemented
+                                            // the following in the naive
+                                            // way, we would therefore incur
+                                            // a very significant penalty: a
+                                            // preliminary version of the
+                                            // Stokes tutorial program would
+                                            // slow down from 17 to 23
+                                            // seconds because of this single
+                                            // function!
+                                            //
+                                            // consequently, we try to be a
+                                            // bit smarter by already laying
+                                            // out the data in the right
+                                            // format and creating a
+                                            // symmetric tensor of it
+           switch (nonzero_component - first_vector_component)
+             {
+               case 0:
+               {
+                 const double array[symmetric_gradient_type::n_independent_components]
+                   = { grad[0], 0, 0, grad[1]/2 , grad[2]/2, 0};
+                 return symmetric_gradient_type(array);
+               }
+
+               case 1:
+               {
+                 const double array[symmetric_gradient_type::n_independent_components]
+                   = { 0, grad[1], 0, grad[0]/2, 0, grad[2]/2 };
+                 return symmetric_gradient_type(array);
+               }
+
+               case 2:
+               {
+                 const double array[symmetric_gradient_type::n_independent_components]
+                   = { 0, 0, grad[2], 0, grad[0]/2, grad[1]/2 };
+                 return symmetric_gradient_type(array);
+               }
+
+               default:
+                     Assert (false, ExcInternalError());
+                     return symmetric_gradient_type();
+             }
+         }
+       else
+         return symmetric_gradient_type();
+      }
+    else
+      {
+       return symmetrize (gradient(shape_function, q_point));
+      }
+  }
+}
+
+
+
+
+
+
 /*------------------------ Inline functions: FEValuesBase ------------------------*/
 
 
+
+
+
+template <int dim>
+inline
+FEValuesViews::Scalar<dim>
+FEValuesBase<dim>::
+operator[] (const FEValuesExtractors::Scalar &scalar) const
+{
+  return FEValuesViews::Scalar<dim> (*this, scalar.component);
+}
+
+
+
+template <int dim>
+inline
+FEValuesViews::Vector<dim>
+FEValuesBase<dim>::
+operator[] (const FEValuesExtractors::Vector &vector) const
+{
+  return
+    FEValuesViews::Vector<dim> (*this, vector.first_vector_component);
+}
+
+
+
 template <int dim>
 inline
 double
 FEValuesBase<dim>::shape_value (const unsigned int i,
                                const unsigned int j) const
 {
+  Assert (i < fe->dofs_per_cell,
+         ExcIndexRange (i, 0, fe->dofs_per_cell));
   Assert (this->update_flags & update_values,
          ExcAccessToUninitializedField());
   Assert (fe->is_primitive (i),
@@ -2527,6 +3645,8 @@ FEValuesBase<dim>::shape_value_component (const unsigned int i,
                                          const unsigned int j,
                                          const unsigned int component) const
 {
+  Assert (i < fe->dofs_per_cell,
+         ExcIndexRange (i, 0, fe->dofs_per_cell));
   Assert (this->update_flags & update_values,
          ExcAccessToUninitializedField());
   Assert (component < fe->n_components(),
@@ -2560,7 +3680,7 @@ FEValuesBase<dim>::shape_value_component (const unsigned int i,
                                       // is non-zero at all within
                                       // this component:
       if (fe->get_nonzero_components(i)[component] == false)
-       return 0.;
+       return 0;
 
                                       // count how many non-zero
                                       // component the shape function
@@ -2577,7 +3697,7 @@ FEValuesBase<dim>::shape_value_component (const unsigned int i,
                           fe->get_nonzero_components(i).begin()+component,
                           true));
       return this->shape_values(row, j);
-    };
+    }
 }
 
 
@@ -2588,6 +3708,8 @@ const Tensor<1,dim> &
 FEValuesBase<dim>::shape_grad (const unsigned int i,
                               const unsigned int j) const
 {
+  Assert (i < fe->dofs_per_cell,
+         ExcIndexRange (i, 0, fe->dofs_per_cell));
   Assert (this->update_flags & update_gradients,
          ExcAccessToUninitializedField());
   Assert (fe->is_primitive (i),
@@ -2623,6 +3745,8 @@ FEValuesBase<dim>::shape_grad_component (const unsigned int i,
                                         const unsigned int j,
                                         const unsigned int component) const
 {
+  Assert (i < fe->dofs_per_cell,
+         ExcIndexRange (i, 0, fe->dofs_per_cell));
   Assert (this->update_flags & update_gradients,
          ExcAccessToUninitializedField());
   Assert (component < fe->n_components(),
@@ -2673,7 +3797,7 @@ FEValuesBase<dim>::shape_grad_component (const unsigned int i,
                           fe->get_nonzero_components(i).begin()+component,
                           true));
       return this->shape_gradients[row][j];
-    };
+    }
 }
 
 
@@ -2682,8 +3806,10 @@ template <int dim>
 inline
 const Tensor<2,dim> &
 FEValuesBase<dim>::shape_hessian (const unsigned int i,
-                                        const unsigned int j) const
+                                 const unsigned int j) const
 {
+  Assert (i < fe->dofs_per_cell,
+         ExcIndexRange (i, 0, fe->dofs_per_cell));
   Assert (this->update_flags & update_hessians,
          ExcAccessToUninitializedField());
   Assert (fe->is_primitive (i),
@@ -2722,13 +3848,16 @@ FEValuesBase<dim>::shape_2nd_derivative (const unsigned int i,
 }
 
 
+
 template <int dim>
 inline
 Tensor<2,dim>
 FEValuesBase<dim>::shape_hessian_component (const unsigned int i,
-                                                  const unsigned int j,
-                                                  const unsigned int component) const
+                                           const unsigned int j,
+                                           const unsigned int component) const
 {
+  Assert (i < fe->dofs_per_cell,
+         ExcIndexRange (i, 0, fe->dofs_per_cell));
   Assert (this->update_flags & update_hessians,
          ExcAccessToUninitializedField());
   Assert (component < fe->n_components(),
@@ -2779,10 +3908,11 @@ FEValuesBase<dim>::shape_hessian_component (const unsigned int i,
                           fe->get_nonzero_components(i).begin()+component,
                           true));
       return this->shape_hessians[row][j];
-    };
+    }
 }
 
 
+
 template <int dim>
 inline
 Tensor<2,dim>
@@ -2794,6 +3924,7 @@ FEValuesBase<dim>::shape_2nd_derivative_component (const unsigned int i,
 }
 
 
+
 template <int dim>
 inline
 const FiniteElement<dim> & 
@@ -2870,19 +4001,24 @@ FEValuesBase<dim>::JxW (const unsigned int i) const
 }
 
 
+
 template <int dim>
 template <class InputVector>
+inline
 void
 FEValuesBase<dim>::get_function_grads (const InputVector           &fe_function,
-                   std::vector<Tensor<1,dim> > &gradients) const
+                                      std::vector<Tensor<1,dim> > &gradients) const
 {
   get_function_gradients(fe_function, gradients);
 }
 
 
+
 template <int dim>
 template <class InputVector>
-void FEValuesBase<dim>::get_function_grads (
+inline
+void
+FEValuesBase<dim>::get_function_grads (
   const InputVector& fe_function,
   const VectorSlice<const std::vector<unsigned int> >& indices,
   std::vector<Tensor<1,dim> > &values) const
@@ -2891,8 +4027,10 @@ void FEValuesBase<dim>::get_function_grads (
 }
 
 
+
 template <int dim>
 template <class InputVector>
+inline
 void
 FEValuesBase<dim>::
 get_function_grads (const InputVector                         &fe_function,
@@ -2902,9 +4040,12 @@ get_function_grads (const InputVector                         &fe_function,
 }
 
 
+
 template <int dim>
 template <class InputVector>
-void FEValuesBase<dim>::get_function_grads (
+inline
+void
+FEValuesBase<dim>::get_function_grads (
   const InputVector& fe_function,
   const VectorSlice<const std::vector<unsigned int> >& indices,
   std::vector<std::vector<Tensor<1,dim> > >& values,
@@ -2917,7 +4058,8 @@ void FEValuesBase<dim>::get_function_grads (
 
 template <int dim>
 template <class InputVector>
-inline void
+inline
+void
 FEValuesBase<dim>::
 get_function_2nd_derivatives (const InputVector           &fe_function,
                              std::vector<Tensor<2,dim> > &hessians) const
@@ -2929,6 +4071,7 @@ get_function_2nd_derivatives (const InputVector           &fe_function,
 
 template <int dim>
 template <class InputVector>
+inline
 void
 FEValuesBase<dim>::
 get_function_2nd_derivatives (const InputVector                         &fe_function,
index 5fbb1c202426bfe2b81a9481adfa4da30db3de1e..7e429b45144f3eb1fb920fae5edb16e74af185fc 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2003, 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -67,7 +67,7 @@
  *
  * The classes in this module are used when one wants to assemble matrices or
  * vectors. They link finite elements, quadrature objects, and mappins: the
- * finite element classes describes a finite element space on a unit cell
+ * finite element classes describe a finite element space on a unit cell
  * (i.e. the unit line segment, square, or cube <tt>[0,1]^d</tt>), the
  * quadrature classes describe where quadrature points are located and what
  * weight they have, and the mapping classes describe how to map a point from
  * functionality as the FEValues class does for cells. Finally, the
  * FESubfaceValues class offers the possibility to ingrate on parts of faces
  * if the neighboring cell is refined and the present cell shares only a part
- * of its face with the neighboring cell.
+ * of its face with the neighboring cell. If vector-valued elements are used,
+ * the FEValues and related classes allow access to all vector components; if
+ * one wants to pick individual components, there are extractor classes that
+ * make this task simpler, as described in the @ref vector_valued module.
  *
  * The last member of this group, the UpdateFlags enumeration, is used as an
  * optimization: instead of letting the FEValues class compute every possible
@@ -89,8 +92,8 @@
  * UpdateFlags enumeration is used to offer symbolic names denoting what you
  * want the FEValues class to compute.
  * 
- * All these classes are used in all tutorial programs from step-3 onward, and
- * are described there in significant detail.
+ * All these classes are used in all @ref Tutorial "tutorial programs" from
+ * step-3 onward, and are described there in significant detail.
  *
  * The actual workings of the FEValues class and friends is
  * complicated because it has to be general yet efficient. The page on
index d925476530c6fc0e94de1e215b448374090ce7e9..d596548693055786a29be5a015b7f6dba07f48ae 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
  * <dd>Blocks were introduced in BlockVector,
  * BlockSparseMatrix and related classes. These are used to reflect the
  * structure of a PDE system in linear algebra, in particular allowing
- * for modular solvers. In DoFHandler, this block structure is
- * prepared by DoFRenumbering::block_wise().
+ * for modular solvers for problems with multiple solution components.
+ * How to implement this is described in more detail in the
+ * @ref vector_valued report and the tutorial programs referenced
+ * therein.
  *
  * Originally, this concept was intermixed with the idea of the vector
  * @ref GlossComponent "component". Since the introduction of
diff --git a/deal.II/doc/doxygen/headers/vector_valued.h b/deal.II/doc/doxygen/headers/vector_valued.h
new file mode 100644 (file)
index 0000000..b16aa27
--- /dev/null
@@ -0,0 +1,740 @@
+//-------------------------------------------------------------------------
+//    $Id: fe.h 15458 2007-11-07 22:59:17Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------------
+
+
+/**
+ * @defgroup vector_valued Handling vector valued problems
+ *
+ *
+ * Vector-valued problems are partial differential equations in which the
+ * solution variable is not a scalar function, but a vector-valued function or
+ * a set of functions. This includes, for example:
+ * <ul>
+ *   <li>The elasticity equation discussed in @ref step_8 "step-8",
+ *       @ref step_17 "step-17", and @ref step_18 "step-18" in which the
+ *       solution is the vector-valued displacement at each point.
+ *   <li>The mixed Laplace equation and its extensions discussed in
+ *       @ref step_20 "step-20", and @ref step_21 "step-21" in which the
+ *       solution is the scalar pressure and the vector-valued velocity
+ *       at each point.
+ *   <li>The Stokes equation and its extensions discussed in
+ *       @ref step_22 "step-22", and @ref step_31 "step-31" in which again
+ *       the solution is the scalar pressure and the vector-valued velocity
+ *       at each point.
+ *   <li>Complex-valued solutions consisting of real and imaginary parts, as
+ *       discussed for example in @ref step_29 "step-29".
+ * </ul>
+ * 
+ * This page gives an overview of how to implement such vector-valued problems
+ * efficiently in deal.II.
+ *
+ * <table class="tutorial" width="50%">
+ * <tr><th><b>%Table of contents</b></th></tr>
+ * <tr><td width="100%" valign="top">
+ * <ol>
+ *  <li> @ref VVPhilosophy "A philosophical view of vector-valued problems"
+ *  <li> @ref VVFEs "Describing finite element spaces"
+ *  <li> @ref VVAssembling "Assembling linear systems"
+ *  <li> @ref VVAlternative "An alternative approach"
+ *  <li> @ref VVBlockSolvers "Block solvers"
+ * </ol> </td> </tr> </table>
+ *
+ *
+ *
+ * @anchor VVPhilosophy
+ * <h3>A philosophical view of vector-valued problems</h3>
+ *
+ * The way one deals systematically with vector-valued problems is not
+ * fundamentally different from scalar problems in that one needs to come up
+ * first with a weak (variational) formulation of the problem that takes into
+ * account all the solution variables. To understand how to do that, let us
+ * consider a simple example, the mixed Laplace equations discussed in
+ * @ref step_20 "step-20":
+@f{eqnarray*}
+  \textbf{u} + \nabla p &=& 0,
+  \\
+  -\textrm{div}\; \textbf{u} &=& f,
+@f}
+ *
+ * Here, we have <code>dim+1</code> solution components: the scalar pressure
+ * $p$ and the vector-valued velocity $\textbf u$ with <code>dim</code> vector
+ * components.
+ *
+ * A systematic way to get a weak or variational form for this and other
+ * vector problems is to first consider it as a problem where the operators
+ * and solution variables are written in vector and matrix form. For the
+ * example, this would read as follows:
+@f{eqnarray*}
+  \left(
+  \begin{array}{cc} \mathbf 1 & \nabla \\ -\nabla^T & 0 \end{array}
+  \right)
+  \left(
+  \begin{array}{c} \mathbf u \\ p \end{array}
+  \right)
+  =
+  \left(
+  \begin{array}{c} \mathbf 0 \\ f \end{array}
+  \right)
+@f}
+ *
+ * This makes it clear that the solution
+@f{eqnarray*}
+  U =
+  \left(
+  \begin{array}{c} \mathbf u \\ p \end{array}
+  \right)
+@f}
+ * indeed has <code>dim+1</code> components. We note that we could change the
+ * ordering of the solution components $\textbf u$ and $p$ inside $U$ if we
+ * also change columns of the matrix operator.
+ *
+ * Next, we need to think about test functions $V$. We want to multiply both
+ * sides of the equation with them, then integrate over $\Omega$. The result
+ * should be a scalar equality. We can achieve this by choosing $V$ also
+ * vector valued as
+@f{eqnarray*}
+  V =
+  \left(
+  \begin{array}{c} \mathbf v \\ q \end{array}
+  \right).
+@f}
+ *
+ * It is convenient to multiply the matrix-vector equation by the test
+ * function from the left, since this way we automatically get the correct
+ * matrix later on (in the linear system, the matrix is also multiplied from
+ * the right with the solution variable, not from the left), whereas if we
+ * multiplied from the right then the matrix so assembled is the transpose of
+ * the one we really want.
+ *
+ * With this in mind, let us multiply by $V$ and integrate to get the
+ * following equation which has to hold for all test functions $V$:
+@f{eqnarray*}
+  \int_\Omega
+  \left(
+  \begin{array}{c} \mathbf v \\ q \end{array}
+  \right)^T
+  \left(
+  \begin{array}{cc} \mathbf 1 & \nabla \\ -\nabla^T & 0 \end{array}
+  \right)
+  \left(
+  \begin{array}{c} \mathbf u \\ p \end{array}
+  \right)
+  \ dx
+  =
+  \int_\Omega
+  \left(
+  \begin{array}{c} \mathbf v \\ q \end{array}
+  \right)^T
+  \left(
+  \begin{array}{c} \mathbf 0 \\ f \end{array}
+  \right)
+  \ dx,
+@f}
+ * or equivalently:
+@f{eqnarray*}
+  (\mathbf v, \mathbf u)
+  +
+  (\mathbf v, \nabla p)
+  -
+  (q, \mathrm{div}\ \mathbf u)
+  =
+  (q,f).
+@f}
+ *
+ * We get the final form by integrating by part the second term:
+@f{eqnarray*}
+  (\mathbf v, \mathbf u)
+  -
+  (\mathrm{div}\ \mathbf v, p)
+  -
+  (q, \mathrm{div}\ \mathbf u)
+  =
+  (q,f) + (\mathbf n\cdot\mathbf v, p)_{\partial\Omega}.
+@f}
+ *
+ * It is this form that we will later use in assembling the discrete weak form
+ * into a matrix and a right hand side vector: the form in which we have
+ * solution and test functions $U,V$ that each consist of a number of vector
+ * components that we can extract.
+ *
+ *
+ * @anchor VVFEs
+ * <h3>Describing finite element spaces</h3>
+ *
+ * Once we have settled on this description, we need to find a way to describe
+ * the vector-valued finite element space from which we draw solution and test
+ * functions. This is where the FESystem class comes in: it composes
+ * vector-valued finite element spaces from simpler one. For example, if we
+ * were to attempt to use $Q_1$ elements for all <code>dim</code> components
+ * of $\mathbf u$ and the one pressure component $p$, we could use the
+ * following object:
+ * @code
+ *   FESystem<dim> finite_element (FE_Q<dim>(1), dim,
+ *                                 FE_Q<dim>(1), 1);
+ * @endcode
+ *
+ * This means that the final finite element will consist of <code>dim</code>
+ * components made up of FE_Q elements of degree 1, and another one also of
+ * degree 1. If course, a simpler (and more efficient) way to achieve the same
+ * is to use the following form instead:
+ * @code
+ *   FESystem<dim> finite_element (FE_Q<dim>(1), dim+1);
+ * @endcode
+ * Another way (also not efficient, but making it clear which components
+ * belong to which element) would be
+ * @code
+ *   FESystem<dim> finite_element (FESystem<dim>(FE_Q<dim>(1), dim), 1,
+ *                                 FE_Q<dim>(1),                     1);
+ * @endcode
+ * meaning that we first create a vector element with <code>dim</code>
+ * components, each consisting of FE_Q elements of order 1. We then couple
+ * one copy of this already vector-valued element to a single scalar
+ * copy of the FE_Q element of order 1 which will describe the pressure
+ * component.
+ *
+ * As it turns out these (equivalent) choices do not lead to a stable scheme
+ * for the mixed Laplace equation. In @ref step_20 "step-20", we therefore use
+ * a Raviart-Thomas element for the velocities. What exactly this means may be
+ * of less concern for now except that the FE_RaviartThomas class describes
+ * elements that already have <code>dim</code> components. For the pressure,
+ * we use piecewise bi-/trilinear elements that may be discontinuous between
+ * cells; this is done using the FE_DGQ class. The combined element will then
+ * be described by
+ * @code
+ *   FESystem<dim> finite_element (FERaviartThomas<dim>(1), 1,
+ *                                 FE_DGQ<dim>(1),          1);
+ * @endcode
+ * i.e. we combine a single copy of the Raviart-Thomas element with a single
+ * copy of the element used for the pressure $p$.
+ *
+ * In this manner, we can combine the whole vector-valued element from its
+ * individual components. On a sidenote it may also be worth mentioning why we
+ * want to also mirror our logical construction of the weak form in our
+ * program: by creating a finite element object as shown above and then
+ * handing it off to a DoFHandler object, we get an enumeration of <i>all</i>
+ * degrees of freedom at once, whether velocity or pressure. Initially, they
+ * are enumerated in a rather random order, but they can be arranged such that
+ * first all velocity components are enumerated and then all pressures by, for
+ * example, calling DoFRenumbering::component_wise . By and large, many things
+ * become much simpler if we can rely on the fact that there is only a single
+ * DoFHandler that describes everything.
+ *
+ *
+ * @anchor VVAssembling
+ * <h3>Assembling linear systems</h3> 
+ *
+ * The next step is to assemble the linear system. How to do this for the
+ * simple case of a scalar problem has been shown in many tutorial programs,
+ * starting with @ref step_3 "step-3". Here we will show how to do it for
+ * vector problems.
+ *
+ * How to do this is possibly best explained by showing an example
+ * illustrating how the local contribution of a cell to the weak form of above
+ * mixed Laplace equations could be assembled. This is essentially how @ref
+ * step_20 "step-20" does it:
+ * @code
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
+
+  ...
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      local_rhs = 0;
+
+      right_hand_side.value_list (fe_values.get_quadrature_points(),
+                                  rhs_values);
+      
+      for (unsigned int q=0; q<n_q_points; ++q) 
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              local_matrix(i,j) += (fe_values[velocities].value (i, q) *
+                                   fe_values[velocities].value (j, q)
+                                    -
+                                   fe_values[velocities].divergence (i, q) *
+                                   fe_values[pressure].value (j, q)
+                                    -
+                                   fe_values[pressure].value (i, q) *
+                                   fe_values[velocities].divergence (j, q)) *
+                                    fe_values.JxW(q);
+
+            local_rhs(i) += - fe_values[pressure].value (i, q)
+                              rhs_values[q] *
+                              fe_values.JxW(q);
+          }
+ * @endcode 
+ *
+ * So here's what is happening:
+ * <ul>
+ *   <li> The first thing we do is to declare "extractors" (see the
+ *        FEValuesExtractors namespace). These are objects that don't
+ *        do much except store which components of a vector-valued finite
+ *        element constitute a single scalar component, or a tensor of
+ *        rank 1 (i.e. what we call a "physical vector", always consisting
+ *        of <code>dim</code> components). Here, we declare
+ *        an object that represents the velocities consisting of
+ *        <code>dim</code> components starting at component zero, and the
+ *        extractor for the pressure, which is a scalar component at
+ *        position <code>dim</code>.
+ *
+ *   <li> We then do our usual loop over all cells, shape functions, and
+ *        quadrature points. In the innermost loop, we compute the local
+ *        contribution of a pair of shape functions to the global matrix
+ *        and right hand side vector. Recall that the cell contributions
+ *        to the bilinear form (i.e. neglecting boundary terms) looked as
+ *        follows, based on shape functions
+ *        $V_i=\left(\begin{array}{c}\mathbf v_i q_i\end{array}\right),
+ *         V_j=\left(\begin{array}{c}\mathbf v_j q_j\end{array}\right)$:
+          @f{eqnarray*}
+            (\mathbf v_i, \mathbf v_j)
+           -
+           (\mathrm{div}\ \mathbf v_i, q_j)
+           -
+           (q_i, \mathrm{div}\ \mathbf v_j)
+          @f}
+ *        whereas the implementation looked like this:
+ *        @code
+              local_matrix(i,j) += (fe_values[velocities].value (i, q) *
+                                   fe_values[velocities].value (j, q)
+                                    -
+                                   fe_values[velocities].divergence (i, q) *
+                                   fe_values[pressure].value (j, q)
+                                    -
+                                   fe_values[pressure].value (i, q) *
+                                   fe_values[velocities].divergence (j, q)) *
+                                    fe_values.JxW(q);
+ *        @endcode
+ *        The similarities are pretty obvious.
+ *
+ *   <li> Essentially, what happens in above code is this: when you do
+ *        <code>fe_values[pressure]</code>, a so-called "view" is created, i.e.
+ *        an object that unlike the full FEValues object represents not all
+ *        components of a finite element, but only the one(s) represented by
+ *        the extractor object <code>pressure</code> or
+ *        <code>velocities</code>.
+ *
+ *   <li> These views can then be asked for information about these individual
+ *        components. For example, when you write
+ *        <code>fe_values[pressure].value(i,q)</code> you get the
+ *        value of the pressure component of the $i$th shape function $V_i$ at
+ *        the $q$th quadrature point. Because the extractor
+ *        <code>pressure</code> represents a scalar component, the results of
+ *        the operator <code>fe_values[pressure].value(i,q)</code> is a scalar
+ *        number. On the other hand, the call
+ *        <code>fe_values[velocities].value(i,q)</code> would produce the
+ *        value of a whole set of <code>dim</code> components, which would
+ *        be of type <code>Tensor@<1,dim@></code>.
+ *
+ *   <li> Other things that can be done with views is to ask for the gradient
+ *        of a particular shape function's components described by an
+ *        extractor. For example, <code>fe_values[pressure].gradient(i,q)</code>
+ *        would represent the gradient of the scalar pressure component, which
+ *        is of type <code>Tensor@<1,dim@></code>, whereas the gradient of the
+ *        velocities components, <code>fe_values[velocities].value(i,q)</code>
+ *        is a <code>Tensor@<2,dim@></code>, i.e. a matrix $G_{ij}$ that consits
+ *        of entries $G_{ij}=\frac{\partial\phi_i}{\partial x_j}$. Finally,
+ *        both scalar and vector views can be asked for the second derivatives
+ *        ("Hessians"); vector views can also be queried for the
+ *        symmetric gradients as well as the divergence (as we do in the
+ *        example), the latter former a tensor of rank 3, the latter a scalar.
+ *        Finally, vector views can be asked for the symmetric gradient,
+ *        defined as $S_{ij}=\frac 12 \left[\frac{\partial\phi_i}{\partial x_j}
+ *        + \frac{\partial\phi_j}{\partial x_i}\right]$.
+ * </ul>
+ * Other examples of using extractors and views are shown in tutorial programs
+ * @ref step_21 "step-21",
+ * @ref step_31 "step-31",
+ * @ref step_22 "step-22" and a few other programs.
+ *
+ *
+ * @anchor VVAlternative
+ * <h3>An alternative approach</h3>
+ *
+ * There are situations in which one can optimize the assembly of a matrix or
+ * right hand side vector a bit, using knowledge of the finite element in
+ * use. Consider, for example, the bilinear form of the elasticity equations
+ * which we are concerned with first in @ref step_8 "step-8":
+ *
+@f[
+  a({\mathbf u}, {\mathbf v}) =
+  \left(
+    \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v}
+  \right)_\Omega
+  +
+  \sum_{i,j}
+  \left(
+    \mu \partial_i u_j, \partial_i v_j
+  \right)_\Omega,
+  +
+  \sum_{i,j}
+  \left(
+    \mu \partial_i u_j, \partial_j v_i
+  \right)_\Omega,
+@f]
+ *
+ * Here, $\mathbf u$ is a vector function with <code>dim</code> components,
+ * $\mathbf v$ the corresponding test function, and $\lambda,\mu$ are material
+ * parameters. Given our discussions above, the obvious way to implement this
+ * bilinear form would be as follows, using an extractor object that interprets
+ * all <code>dim</code> components of the finite element as single vector,
+ * rather than disjoint components:
+ *
+ * @code
+      const FEValuesExtractors::Vector displacements (0);
+
+      ...
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           const Tensor<2,dim> phi_i_grad
+             = fe_values[displacements].gradient (i,q_point);
+           const double phi_i_div
+             = fe_values[displacements].divergence (i,q_point);
+         
+           for (unsigned int j=0; j<dofs_per_cell; ++j) 
+             {
+               const Tensor<2,dim> phi_j_grad
+                 = fe_values[displacements].gradient (j,q_point);
+               const double phi_j_div
+                 = fe_values[displacements].divergence (j,q_point);
+
+               cell_matrix(i,j) 
+                 +=  (phi_i_div * phi_j_div *
+                      lambda_values[q_point]
+                      +
+                      scalar_product(phi_i_grad, phi_j_grad) *
+                      mu_values[q_point]
+                      +
+                      scalar_product(phi_i_grad, transpose(phi_j_grad)) *
+                      mu_values[q_point])
+                     *
+                     fe_values.JxW(q_point);
+             }
+         }
+ * @endcode 
+ *
+ * The scalar product between two tensors used in this bilinear form is
+ * implemented as follows:
+ *
+ * @code
+template <int dim>
+double
+scalar_product (const Tensor<2,dim> &u,
+               const Tensor<2,dim> &v)
+{
+  double tmp = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      tmp += u[i][j] * v[i][j];
+  return tmp;
+}
+ * @endcode 
+ *
+ * Now, this is not the code used in @ref step_8 "step-8". In fact,
+ * if one used the above code over the one implemented in that program,
+ * it would run about 8 per cent slower. It can be improved (bringing
+ * down the penalty to about 4 per cent) by taking a close look at the
+ * bilinear form. In fact, we can transform it as follows:
+@f{eqnarray*} 
+  a({\mathbf u}, {\mathbf v})
+  &=&
+  \left(
+    \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v}
+  \right)_\Omega
+  +
+  \sum_{i,j}
+  \left(
+    \mu \partial_i u_j, \partial_i v_j
+  \right)_\Omega
+  +
+  \sum_{i,j}
+  \left(
+    \mu \partial_i u_j, \partial_j v_i
+  \right)_\Omega
+  \\
+  &=&
+  \left(
+    \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v}
+  \right)_\Omega
+  +
+  2
+  \sum_{i,j}
+  \left(
+    \mu \partial_i u_j, \frac 12[\partial_i v_j + \partial_j v_i]
+  \right)_\Omega  
+  \\
+  &=&
+  \left(
+    \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v}
+  \right)_\Omega
+  +
+  2
+  \sum_{i,j}
+  \left(
+    \mu \frac 12[\partial_i u_j + \partial_j u_i], \frac 12[\partial_i v_j + \partial_j v_i]
+  \right)_\Omega
+  \\
+  &=&
+  \left(
+    \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v}
+  \right)_\Omega
+  +
+  2
+  \sum_{i,j}
+  \left(
+  \mu \varepsilon(\mathbf u), \varepsilon(\mathbf v)
+  \right)_\Omega,
+@f}
+ * where $\varepsilon(\mathbf u) = \frac 12 \left([\nabla\mathbf u] +
+ * [\nabla\mathbf u]^2\right)$ is the symmetrized gradient.
+ * In the second to last step, we used that the scalar product between
+ * an arbitrary tensor $\nabla\mathbf u$ and a symmetric tensor
+ * $\frac 12[\partial_i v_j + \partial_j v_i]$ equals the scalar product
+ * of the symmetric part of the former with the second tensor. Using the
+ * techniques discussed above, the obvious way to implement this goes
+ * like this:
+ *
+ * @code
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           const SymmetricTensor<2,dim> phi_i_symmgrad
+             = fe_values[displacements].symmetric_gradient (i,q_point);
+           const double phi_i_div
+             = fe_values[displacements].divergence (i,q_point);
+         
+           for (unsigned int j=0; j<dofs_per_cell; ++j) 
+             {
+               const SymmetricTensor<2,dim> phi_j_symmgrad
+                 = fe_values[displacements].symmetric_gradient (j,q_point);
+               const double phi_j_div
+                 = fe_values[displacements].divergence (j,q_point);
+
+               cell_matrix(i,j) 
+                 +=  (phi_i_div * phi_j_div *
+                      lambda_values[q_point]
+                      +
+                      2 *
+                      (phi_i_symmgrad * phi_j_symmgrad) *
+                      mu_values[q_point])
+                     *
+                     fe_values.JxW(q_point));
+             }
+         }
+ * @endcode 
+ *
+ * So if, again, this is not the code we use in @ref step_8 "step-8", what do
+ * we do there? The answer rests on the finite element we use. There, we use the
+ * following element:
+ * @code
+ *   FESystem<dim> finite_element (FE_Q<dim>(1), dim);
+ * @endcode
+ * In other words, the finite element we use consists of <code>dim</code> copies
+ * of the same scalar element. This is what we call a @ref GlossPrimitive
+ * "primitive" element: an element that may be vector-valued but where each
+ * shape function has exactly one non-zero component. In other words: if the
+ * $x$-component of a displacement shape function is nonzero, then the $y$-
+ * and $z$-components must be zero and so on. What this means is that also
+ * derived quantities based on shape functions inherit this sparsity property.
+ * For example: the divergence
+ * $\mathrm{div}\ \Phi(x,y,z)=\partial_x\varphi_x(x,y,z) +
+ * \partial_y\varphi_y(x,y,z) + \partial_z\varphi_z(x,y,z)$
+ * of a vector-valued shape function
+ * $\Phi(x,y,z)=(\varphi_x(x,y,z), \varphi_y(x,y,z), \varphi_z(x,y,z))^T$ is,
+ * in the present case, either
+ * $\mathrm{div}\ \Phi(x,y,z)=\partial_x\varphi_x(x,y,z)$,
+ * $\mathrm{div}\ \Phi(x,y,z)=\partial_y\varphi_y(x,y,z)$, or
+ * $\mathrm{div}\ \Phi(x,y,z)=\partial_z\varphi_z(x,y,z)$, because exactly one
+ * of the $\varphi_\ast$ is nonzero. Knowing this means that we can save a
+ * number of computations that, if we were to do them, would only yield
+ * zeros to add up.
+ *
+ * In a similar vein, if only one component of a shape function is nonzero,
+ * then only one row of its gradient $\nabla\Phi$ is nonzero. What this means
+ * for terms like $(\mu \nabla\Phi_i,\nabla\Phi_j)$, where the scalar product
+ * between two tensors is defined as
+ * $(\tau, \gamma)_\Omega=\int_\Omega \sum_{i,j=1}^d \tau_{ij} \gamma_{ij}$,
+ * is that the term is only nonzero if both tensors have their nonzero
+ * entries in the same row, which means that the two shape functions have
+ * to have their single nonzero component in the same location.
+ *
+ * If we use this sort of knowledge, then we can in a first step avoid
+ * computing gradient tensors if we can determine up front that their
+ * scalar product will be nonzero, in a second step avoid
+ * building the entire tensors and only get its nonzero components,
+ * and in a final step simplify the scalar product by only considering
+ * that index $i$ for the one nonzero row, rather than multiplying and
+ * adding up zeros.
+ *
+ * The vehicle for all this is the ability to determine which vector
+ * component is going to be nonzero. This information is provided by the
+ * FiniteElement::system_to_component_index function. What can be done with
+ * it, using the example above, is explained in detail in
+ * @ref step_8 "step-8".
+ *
+ *
+ * @anchor VVBlockSolvers
+ * <h3>Block solvers</h3>
+ *
+ * Using techniques as shown above, it isn't particularly complicated to
+ * assemble the linear system, i.e. matrix and right hand side, for a
+ * vector-valued problem. However, then it also has to be solved. This is more
+ * complicated. Naively, one could just consider the matrix as a whole. For
+ * most problems, this matrix is not going to be definite (except for special
+ * cases like the elasticity equations covered in @ref step_8 "step-8" and
+ * @ref step_17 "step-17"). It will, often, also not be symmetric. This rather
+ * general class of matrices presents problems for iterative solvers: the lack
+ * of structural properties prevents the use of most efficient methods and
+ * preconditioners. While it can be done, the solution process will therefore
+ * most often be slower than necessary.
+ *
+ * The answer to this problem is to make use of the structure of the
+ * problem. For example, for the mixed Laplace equations discussed above, the
+ * operator has the form
+@f{eqnarray*}
+  \left(
+  \begin{array}{cc} \mathbf 1 & \nabla \\ -\nabla^T & 0 \end{array}
+  \right)
+@f}
+ *
+ * It would be nice if this structure could be recovered in the linear system
+ * as well. For example, after discretization, we would like to have a matrix
+ * with the following block structure:
+@f{eqnarray*}
+  \left(
+  \begin{array}{cc} M & B \\ B^T & 0 \end{array}
+  \right),
+@f}
+ * where $M$ represents the mass matrix that results from discretizing the
+ * identity operator $\mathbf 1$ and $B$ the equivalent of the gradient
+ * operator.
+ *
+ * By default, this is not what happens, however. Rather, deal.II assigns
+ * %numbers to degrees of freedom in a rather random manner. Consequently, if
+ * you form a vector out of the values of degrees of freedom will not be
+ * neatly ordered in a vector like
+@f{eqnarray*}
+  \left(
+  \begin{array}{c} U \\ P \end{array}
+  \right).
+@f}
+ * Rather, it will be a permutation of this, with %numbers of degrees of
+ * freedom corresponding to velocities and pressures intermixed. Consequently,
+ * the system matrix will also not have the nice structure mentioned above,
+ * but with the same permutation or rows and columns.
+ *
+ * What is needed is to re-enumerate degrees of freedom so that velocities
+ * come first and pressures last. This can be done using the
+ * DoFRenumbering::component_wise function, as explained in @ref step_20
+ * "step-20", @ref step_21 "step-21", @ref step_22 "step-22", and @ref step_31
+ * "step-31". After this, at least the degrees of freedom are partitioned
+ * properly.
+ *
+ * But then we still have to make use of it, i.e. we have to come up with a
+ * solver that uses the structure. For example, in @ref step_20 "step-20", we
+ * do a block elimination of the linear system
+@f{eqnarray*}
+  \left(
+  \begin{array}{cc} M & B \\ B^T & 0 \end{array}
+  \right)
+  \left(
+  \begin{array}{c} U \\ P \end{array}
+  \right)
+  =
+  \left(
+  \begin{array}{c} F \\ G \end{array}
+  \right).  
+@f}
+ * What this system means, of course, is
+@f{eqnarray*}
+  MU + BP &=& F,\\
+  B^TU  &=& G.
+@f}
+ *
+ * So, if we multiply the first equation by $B^TM^{-1}$ and subtract the
+ * second from the result, we get
+@f{eqnarray*}
+  B^TM^{-1}BP &=& B^TM^{-1}F-G.
+@f}
+ *
+ * This is an equation that now only contains the pressure variables. If we
+ * can solve it, we can in a second step solve for the velocities using
+@f{eqnarray*}
+  MU = F-BP.
+@f}
+ *
+ * This has the advantage that the matrices $B^TM^{-1}B$ and $M$ that we have
+ * to solve with are both symmetric and positive definite, as opposed to the
+ * large whole matrix we had before.
+ *
+ * How a solver like this is implemented is explained in more detail in @ref
+ * step_20 "step-20", @ref step_31 "step-31", and a few other tutorial
+ * programs. What we would like to point out here is that we now need a way to
+ * extract certain parts of a matrix or vector: if we are to multiply, say,
+ * the $U$ part of the solution vector by the $M$ part of the global matrix,
+ * then we need to have a way to access these parts of the whole.
+ *
+ * This is where the BlockVector, BlockSparseMatrix, and similar classes come
+ * in. For all practical purposes, then can be used as regular vectors or
+ * sparse matrices, i.e. they offer element access, provide the usual vector
+ * operations and implement, for example, matrix-vector multiplications. In
+ * other words, assembling matrices and right hand sides works in exactly the
+ * same way as for the non-block versions. That said, internally they store
+ * the elements of vectors and matrices in "blocks"; for example, instead of
+ * using one large array, the BlockVector class stores it as a set of arrays
+ * each of which we call a block. The advantage is that, while the whole thing
+ * can be used as a vector, one can also access an individual block which
+ * then, again, is a vector with all the vector operations.
+ *
+ * To show how to do this, let us consider the second equation $MU=F-BP$ to be
+ * solved above. This can be achieved using the following sequence similar to
+ * what we have in @ref step_20 "step-20":
+ * @code
+    Vector<double> tmp (solution.block(0).size());
+    system_matrix.block(0,1).vmult (tmp, solution.block(1));
+    tmp *= -1;
+    tmp += system_rhs.block(0);
+    
+
+    SolverControl solver_control (solution.block(0).size(),
+                                  1e-8*tmp.l2_norm());
+    SolverCG<> cg (solver_control, vector_memory);
+  
+    cg.solve (system_matrix.block(0,0),
+              solution.block(0),
+             tmp,
+              PreconditionIdentity());        
+ * @endcode
+ *
+ * What's happening here is that we allocate a temporary vector with as many
+ * elements as the first block of the solution vector, i.e. the velocity
+ * component $U$, has. We then set this temporary vector equal to the $(0,1)$
+ * block of the matrix, i.e. $B$, times component 1 of the solution which is
+ * the previously computed pressure $P$. The result is multiplied by $-1$, and
+ * component 0 of the right hand side, $F$ is added to it. The temporary
+ * vector now contains $F-BP$. The rest of the code snippet simply solves a
+ * linear system with $F-BP$ as right hand side and the $(0,0)$ block of the
+ * global matrix, i.e. $M$. Using block vectors and matrices in this way
+ * therefore allows us to quite easily write rather complicated solvers making
+ * use of the block structure of a linear system.
+ *
+ * @ingroup feall feaccess
+ */
index a4674816d9f812884026b84d2fee5637dfafb9f9..c2ea7c8cfef213547c4234599b5aa2fd1e1b7956 100644 (file)
@@ -109,6 +109,15 @@ inconvenience this causes.
   (Moritz Allmaras, 2007/10/31)
   </p>
 
+  <li> <p>New: A significantly simpler way to code the assembly of linear
+  systems for problems with more than one solution variable has been
+  implemented. This is explained in detail in the report on @ref vector_valued
+  and tutorial programs @ref step_20 "step-20" and @ref step_21 "step-21"
+  have been converted as well.
+  <br>
+  (WB, 2008/02/23)
+  </p>
+
   <li> <p>Improved: On Mac OS X, the operating system provides for
   "frameworks", which are essentially collections of shared libraries.
   We now link with the "Accelerate" framework (from Mac OS X 10.4
index 46f0756d032f0863821403b50b7e31f9885a0f39..97bcd8cda161ffe914418e7636ae9dbe2990081e 100644 (file)
@@ -24,6 +24,9 @@ We are going to extend this tutorial program in @ref step_21 "step-21" to
 solve not only the mixed Laplace equation, but add another equation that
 describes the transport of a mixture of two fluids.
 
+The equations covered here fall into the class of vector-valued problems. A
+toplevel overview of this topic can be found in the @ref vector_valued module.
+
 
 <h3>Formulation, weak form, and discrete problem</h3>
 
@@ -202,90 +205,84 @@ This is, at best, tedious, error prone, and not dimension independent. There
 are obvious ways to make things dimension independent, but in the end, the
 code is simply not pretty. What would be much nicer is if we could simply
 extract the ${\mathbf u}$ and $p$ components of a shape function $x_h^i$. In the
-program we do that, by writing functions like this one:
+program we do that in the following way:
 @code
-template <int dim>
-Tensor<1,dim>
-extract_u (const FEValuesBase<dim> &fe_values,
-           const unsigned int i,
-           const unsigned int q)
-{
-  Tensor<1,dim> tmp;
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
 
-  for (unsigned int d=0; d<dim; ++d)
-    tmp[d] = fe_values.shape_value_component (i,q,d);
+  ...
 
-  return tmp;
-}
+  for (unsigned int q=0; q<n_q_points; ++q) 
+    for (unsigned int i=0; i<dofs_per_cell; ++i)
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+        local_matrix(i,j) += (fe_values[velocities].value (i, q) *
+                             k_inverse_values[q] *
+                              fe_values[velocities].value (j, q)
+                              -
+                              fe_values[velocities].divergence (i, q) *
+                              fe_values[pressure].value (j, q)
+                              -
+                              fe_values[pressure].value (i, q) *
+                              fe_values[velocities].divergence (j, q)) *
+                              fe_values.JxW(q);
 @endcode
+This is, in fact, not only the first term of the bilinear form, but the
+whole thing (sans boundary contributions).
 
-What this function does is, given an <code>fe_values</code> object, to extract
+What this piece of code does is, given an <code>fe_values</code> object, to extract
 the values of the first $dim$ components of shape function <code>i</code> at
 quadrature points <code>q</code>, that is the velocity components of that shape
 function. Put differently, if we write shape functions $x_h^i$ as the tuple
 $\{{\mathbf u}_h^i,p_h^i\}$, then the function returns the velocity part of this
 tuple. Note that the velocity is of course a <code>dim</code>-dimensional tensor, and
-that the function returns a corresponding object.
-
-Likewise, we have a function that extracts the pressure component of a shape
-function:
-@code
-template <int dim>
-double extract_p (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  return fe_values.shape_value_component (i,q,dim);
-}
-@endcode
-Finally, the bilinear form contains terms involving the gradients of the
-velocity component of shape functions. To be more precise, we are not really
-interested in the full gradient, but only the divergence of the velocity
-components, i.e. ${\textrm{div}}\ {\mathbf u}_h^i = \sum_{d=0}^{dim-1}
-\frac{\partial}{\partial x_d} ({\mathbf u}_h^i)_d$. Here's a function that returns
-this quantity:
-@code
-template <int dim>
-double
-extract_div_u (const FEValuesBase<dim> &fe_values,
-               const unsigned int i,
-               const unsigned int q)
-{
-  double divergence = 0;
-  for (unsigned int d=0; d<dim; ++d)
-    divergence += fe_values.shape_grad_component (i,q,d)[d];
+that the function returns a corresponding object. Similarly, where we
+subscript with the pressure extractor, we extract the scalar pressure
+component. The whole mechanism is described in more detail in the 
+@ref vector_valued module.
 
-  return divergence;
-}
-@endcode
+In practice, it turns out that we can do a bit better if we evaluate the shape
+functions, their gradients and divergences only once per outermost loop, and
+store the result, as this saves us a few otherwise repeated computations.
+The final result then looks like this, working in every space dimension:
 
-With these three functions, all of which are completely dimension independent
-and will therefore also work in 3d, assembling the local matrix and right hand
-side contributions becomes a charm:
 @code
-for (unsigned int q=0; q<n_q_points; ++q) 
-  for (unsigned int i=0; i<dofs_per_cell; ++i)
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
     {
-      const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q);
-      const double div_phi_i_u    = extract_div_u (fe_values, i, q);
-      const double phi_i_p        = extract_p (fe_values, i, q);
-           
-      for (unsigned int j=0; j<dofs_per_cell; ++j)
-        {
-          const Tensor<1,dim> phi_j_u = extract_u (fe_values, j, q);
-          const double div_phi_j_u    = extract_div_u (fe_values, j, q);
-          const double phi_j_p        = extract_p (fe_values, j, q);
-               
-          local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u
-                                - div_phi_i_u * phi_j_p
-                                - phi_i_p * div_phi_j_u)
-                               * fe_values.JxW(q);
-        }
-
-      local_rhs(i) += -(phi_i_p *
-                        rhs_values[q] *
-                        fe_values.JxW(q));
-    }
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      local_rhs = 0;
+
+      right_hand_side.value_list (fe_values.get_quadrature_points(),
+                                  rhs_values);
+      k_inverse.value_list (fe_values.get_quadrature_points(),
+                            k_inverse_values);
+      
+      for (unsigned int q=0; q<n_q_points; ++q) 
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            const Tensor<1,dim> phi_i_u     = fe_values[velocities].value (i, q);
+            const double        div_phi_i_u = fe_values[velocities].divergence (i, q);
+            const double        phi_i_p     = fe_values[pressure].value (i, q);
+            
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              {
+               const Tensor<1,dim> phi_j_u     = fe_values[velocities].value (j, q);
+               const double        div_phi_j_u = fe_values[velocities].divergence (j, q);
+               const double        phi_j_p     = fe_values[pressure].value (j, q);
+                
+                local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u
+                                      - div_phi_i_u * phi_j_p
+                                      - phi_i_p * div_phi_j_u) *
+                                     fe_values.JxW(q);
+              }
+
+            local_rhs(i) += -phi_i_p *
+                            rhs_values[q] *
+                            fe_values.JxW(q);
+          }
 @endcode
 This very closely resembles the form in which we have originally written down
 the bilinear form and right hand side.
@@ -294,40 +291,29 @@ There is one final term that we have to take care of: the right hand side
 contained the term $(g,{\mathbf v}\cdot {\mathbf n})_{\partial\Omega}$, constituting the
 weak enforcement of pressure boundary conditions. We have already seen in
 @ref step_7 "step-7" how to deal with face integrals: essentially exactly the same as with
-domain integrals, except that we have to use the <code>FEFaceValues</code> class
+domain integrals, except that we have to use the FEFaceValues class
 instead of <code>FEValues</code>. To compute the boundary term we then simply have
-to loop over all boundary faces and integrate there. If you look closely at
-the definitions of the <code>extract_*</code> functions above, you will realize
-that it isn't even necessary to write new functions that extract the velocity
-and pressure components of shape functions from <code>FEFaceValues</code> objects:
-both <code>FEValues</code> and <code>FEFaceValues</code> are derived from a common
-base class, <code>FEValuesBase</code>, and the extraction functions above can
-therefore deal with both in exactly the same way. Assembling the missing
-boundary term then takes on the following form:
+to loop over all boundary faces and integrate there. The mechanism works in
+the same way as above, i.e. the extractor classes also work on FEFaceValues objects:
 @code
-for (unsigned int face_no=0;
-     face_no<GeometryInfo<dim>::faces_per_cell;
-     ++face_no)
-  if (cell->at_boundary(face_no))
-    {
-      fe_face_values.reinit (cell, face_no);
-    
-      pressure_boundary_values
-        .value_list (fe_face_values.get_quadrature_points(),
-                     boundary_values);
-
-      for (unsigned int q=0; q<n_face_q_points; ++q) 
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          {
-            const Tensor<1,dim>
-              phi_i_u = extract_u (fe_face_values, i, q);
-                
-            local_rhs(i) += -(phi_i_u *
-                              fe_face_values.normal_vector(q) *
-                              boundary_values[q] *
-                              fe_face_values.JxW(q));
-        }
-  }
+      for (unsigned int face_no=0;
+          face_no<GeometryInfo<dim>::faces_per_cell;
+          ++face_no)
+       if (cell->at_boundary(face_no))
+         {
+           fe_face_values.reinit (cell, face_no);
+           
+           pressure_boundary_values
+             .value_list (fe_face_values.get_quadrature_points(),
+                          boundary_values);
+
+           for (unsigned int q=0; q<n_face_q_points; ++q) 
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               local_rhs(i) += -(fe_face_values[velocities].value (i, q) *
+                                 fe_face_values.normal_vector(q) *
+                                 boundary_values[q] *
+                                 fe_face_values.JxW(q));
+         }
 @endcode
 
 You will find the exact same code as above in the sources for the present
index 76552ff4fea8c2e35198e92bf040e113e6933bac..0e0df238c148c72c7dc26802cba2f6232441ada4 100644 (file)
@@ -320,54 +320,6 @@ KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
 }
 
 
-                                 // @sect3{extract_u and friends}
-
-                                 // The next three functions are
-                                 // needed for matrix and right hand
-                                 // side assembly. They are described
-                                 // in detail in the introduction to
-                                 // this program, so that we do not
-                                 // need to discuss them here again:
-template <int dim>
-Tensor<1,dim>
-extract_u (const FEValuesBase<dim> &fe_values,
-          const unsigned int i,
-          const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  for (unsigned int d=0; d<dim; ++d)
-    tmp[d] = fe_values.shape_value_component (i,q,d);
-
-  return tmp;
-}
-
-
-
-template <int dim>
-double
-extract_div_u (const FEValuesBase<dim> &fe_values,
-              const unsigned int i,
-              const unsigned int q)
-{
-  double divergence = 0;
-  for (unsigned int d=0; d<dim; ++d)
-    divergence += fe_values.shape_grad_component (i,q,d)[d];
-
-  return divergence;
-}
-
-
-  
-template <int dim>
-double extract_p (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  return fe_values.shape_value_component (i,q,dim);
-}
-
-
 
                                  // @sect3{MixedLaplaceProblem class implementation}
 
@@ -646,6 +598,23 @@ void MixedLaplaceProblem<dim>::assemble_system ()
   std::vector<double> boundary_values (n_face_q_points);
   std::vector<Tensor<2,dim> > k_inverse_values (n_q_points);
 
+                                  // Finally, we need a couple of extractors
+                                  // that we will use to get at the velocity
+                                  // and pressure components of vector-valued
+                                  // shape functions. Their function and use
+                                  // is described in detail in the @ref
+                                  // vector_valued report. Essentially, we
+                                  // will use them as subscripts on the
+                                  // FEValues objects below: the FEValues
+                                  // object describes all vector components
+                                  // of shape functions, while after
+                                  // subscription, it will only refer to the
+                                  // velocities (a set of <code>dim</code>
+                                  // components starting at component zero)
+                                  // or the pressure (a scalar component
+                                  // located at position <code>dim</code>):
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
 
                                    // With all this in place, we can
                                    // go on with the loop over all
@@ -670,15 +639,15 @@ void MixedLaplaceProblem<dim>::assemble_system ()
       for (unsigned int q=0; q<n_q_points; ++q) 
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
-            const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q);
-            const double div_phi_i_u = extract_div_u (fe_values, i, q);
-            const double phi_i_p = extract_p (fe_values, i, q);
+            const Tensor<1,dim> phi_i_u     = fe_values[velocities].value (i, q);
+            const double        div_phi_i_u = fe_values[velocities].divergence (i, q);
+            const double        phi_i_p     = fe_values[pressure].value (i, q);
             
             for (unsigned int j=0; j<dofs_per_cell; ++j)
               {
-                const Tensor<1,dim> phi_j_u = extract_u (fe_values, j, q);
-                const double div_phi_j_u = extract_div_u (fe_values, j, q);
-                const double phi_j_p = extract_p (fe_values, j, q);
+               const Tensor<1,dim> phi_j_u     = fe_values[velocities].value (j, q);
+               const double        div_phi_j_u = fe_values[velocities].divergence (j, q);
+               const double        phi_j_p     = fe_values[pressure].value (j, q);
                 
                 local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * phi_j_u
                                       - div_phi_i_u * phi_j_p
@@ -704,15 +673,10 @@ void MixedLaplaceProblem<dim>::assemble_system ()
 
            for (unsigned int q=0; q<n_face_q_points; ++q) 
              for (unsigned int i=0; i<dofs_per_cell; ++i)
-               {
-                 const Tensor<1,dim>
-                   phi_i_u = extract_u (fe_face_values, i, q);
-
-                 local_rhs(i) += -(phi_i_u *
-                                   fe_face_values.normal_vector(q) *
-                                   boundary_values[q] *
-                                   fe_face_values.JxW(q));
-               }
+               local_rhs(i) += -(fe_face_values[velocities].value (i, q) *
+                                 fe_face_values.normal_vector(q) *
+                                 boundary_values[q] *
+                                 fe_face_values.JxW(q));
          }
 
                                        // The final step in the loop
index 66761658ac34ae20f3a7792a539e9a1af5d68811..d08c107a41168cf5c89376453a355ac13918c1ea 100644 (file)
@@ -10,6 +10,11 @@ equation. This is therefore also the first time-dependent tutorial
 program (besides the somewhat strange time-dependence of @ref step_18
 "step-18").  
 
+The equations covered here are an extension of the material already covered in
+@ref step_20 "step-20". In particular, they fall into the class of
+vector-valued problems. A toplevel overview of this topic can be found in the
+@ref vector_valued module.
+
 
 <h2>The two phase flow problem</h2>
 
index b81f519fcab5a492f4881bd7c2a63ee92ba9efa3..a5d093599bc6e0d85563c9ef5146f250fdec3c7a 100644 (file)
@@ -495,80 +495,6 @@ double f_saturation (const double S,
 
 
 
-                                 // @sect3{extract_u and friends}
-
-                                 // More tools: We need methods to extract the
-                                 // velocity, pressure, and saturation
-                                 // components of finite element shape
-                                 // functions. These functions here are
-                                 // completely analogous to the ones we have
-                                 // already used in step-20:
-template <int dim>
-Tensor<1,dim>
-extract_u (const FEValuesBase<dim> &fe_values,
-           const unsigned int i,
-           const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  for (unsigned int d=0; d<dim; ++d)
-    tmp[d] = fe_values.shape_value_component (i,q,d);
-
-  return tmp;
-}
-
-
-
-template <int dim>
-double
-extract_div_u (const FEValuesBase<dim> &fe_values,
-               const unsigned int i,
-               const unsigned int q)
-{
-  double divergence = 0;
-  for (unsigned int d=0; d<dim; ++d)
-    divergence += fe_values.shape_grad_component (i,q,d)[d];
-
-  return divergence;
-}
-
-
-  
-template <int dim>
-double extract_p (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  return fe_values.shape_value_component (i,q,dim);
-}
-
-
-
-template <int dim>
-double extract_s (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  return fe_values.shape_value_component (i,q,dim+1);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-extract_grad_s (const FEValuesBase<dim> &fe_values,
-                const unsigned int i,
-                const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-  for (unsigned int d=0; d<dim; ++d)
-    tmp[d] = fe_values.shape_grad_component (i,q,dim+1)[d];
-
-  return tmp;
-}
-
-
-
                                  // @sect3{Linear solvers and preconditioners}
 
                                  // The linear solvers we use are also
@@ -855,7 +781,11 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
   std::vector<Vector<double> >      old_solution_values(n_q_points, Vector<double>(dim+2));
   std::vector<std::vector<Tensor<1,dim> > >  old_solution_grads(n_q_points,
                                                                 std::vector<Tensor<1,dim> > (dim+2));
-  
+
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
+  const FEValuesExtractors::Scalar saturation (dim+1);
+
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -908,18 +838,18 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
           {
             const double old_s = old_solution_values[q](dim+1);
 
-            const Tensor<1,dim> phi_i_u      = extract_u (fe_values, i, q);
-            const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
-            const double        phi_i_p      = extract_p (fe_values, i, q);
-            const double        phi_i_s      = extract_s (fe_values, i, q); 
-            const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+            const Tensor<1,dim> phi_i_u      = fe_values[velocities].value (i, q);
+            const double        div_phi_i_u  = fe_values[velocities].divergence (i, q);
+            const double        phi_i_p      = fe_values[pressure].value (i, q);
+            const double        phi_i_s      = fe_values[saturation].value (i, q); 
+            const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient(i, q);
             
             for (unsigned int j=0; j<dofs_per_cell; ++j)
               {
-                const Tensor<1,dim> phi_j_u     = extract_u (fe_values, j, q);
-                const double        div_phi_j_u = extract_div_u (fe_values, j, q);
-                const double        phi_j_p     = extract_p (fe_values, j, q);
-                const double        phi_j_s     = extract_s (fe_values, j, q);
+                const Tensor<1,dim> phi_j_u     = fe_values[velocities].value (j, q);
+                const double        div_phi_j_u = fe_values[velocities].divergence (j, q);
+                const double        phi_j_p     = fe_values[pressure].value (j, q);
+                const double        phi_j_s     = fe_values[saturation].value (j, q);
                 
                 local_matrix(i,j) += (phi_i_u * k_inverse_values[q] *
                                       mobility_inverse(old_s,viscosity) * phi_j_u
@@ -952,7 +882,7 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
               for (unsigned int i=0; i<dofs_per_cell; ++i)
                 {
                   const Tensor<1,dim>
-                    phi_i_u = extract_u (fe_face_values, i, q);
+                    phi_i_u = fe_face_values[velocities].value (i, q);
 
                   local_rhs(i) += -(phi_i_u *
                                     fe_face_values.normal_vector(q) *
@@ -1023,6 +953,8 @@ void TwoPhaseFlowProblem<dim>::assemble_rhs_S ()
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
   SaturationBoundaryValues<dim> saturation_boundary_values;
+
+  const FEValuesExtractors::Scalar saturation (dim+1);
   
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
@@ -1049,8 +981,8 @@ void TwoPhaseFlowProblem<dim>::assemble_rhs_S ()
             for (unsigned int d=0; d<dim; ++d)
               present_u[d] = present_solution_values[q](d);
 
-            const double        phi_i_s      = extract_s(fe_values, i, q);
-            const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+            const double        phi_i_s      = fe_values[saturation].value (i, q);
+            const Tensor<1,dim> grad_phi_i_s = fe_values[saturation].gradient (i, q);
                      
             local_rhs(i) += (time_step *
                              f_saturation(old_s,viscosity) *
@@ -1128,7 +1060,7 @@ void TwoPhaseFlowProblem<dim>::assemble_rhs_S ()
                                               :
                                               neighbor_saturation[q]),
                                              viscosity) *
-                                extract_s(fe_face_values,i,q) *
+                                fe_face_values[saturation].value (i,q) *
                                 fe_face_values.JxW(q);
             }
         }
index 358f177418a1fbf38e931ed9bc1af356402e3686..dbb770c0de8fb62b88126c8ebe13fec14425b6dd 100644 (file)
@@ -92,6 +92,9 @@ mixing over time scales of many million years, a time scale much
 shorter than for the same amount of heat to be distributed by thermal
 conductivity. 
 
+The equations covered here fall into the class of vector-valued problems. A
+toplevel overview of this topic can be found in the @ref vector_valued module.
+
 
 <h2>%Boundary and initial conditions</h2>
 
@@ -264,119 +267,6 @@ temperature field using the just computed velocity field.
 
 
 
-@code
-template <int dim>
-Tensor<1,dim>
-extract_u (const FEValuesBase<dim> &fe_values,
-           const unsigned int i,
-           const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component < dim)
-    tmp[component] = fe_values.shape_value (i,q);
-
-  return tmp;
-}
-
-
-
-template <int dim>
-double
-extract_div_u (const FEValuesBase<dim> &fe_values,
-               const unsigned int i,
-               const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component < dim)
-    return fe_values.shape_grad (i,q)[component];
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-extract_grad_s_u (const FEValuesBase<dim> &fe_values,
-                 const unsigned int i,
-                 const unsigned int q)
-{
-  Tensor<2,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-  
-  if (component < dim)
-    {
-      const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2;
-      
-      for (unsigned int e=0; e<dim; ++e)
-       tmp[component][e] += grad_phi_over_2[e];
-      for (unsigned int d=0; d<dim; ++d)
-       tmp[d][component] += grad_phi_over_2[d];
-    }
-  
-  return tmp;
-}
-
-
-  
-template <int dim>
-double extract_p (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim)
-    return fe_values.shape_value (i,q);
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-double extract_T (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim+1)
-    return fe_values.shape_value (i,q);
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-extract_grad_T (const FEValuesBase<dim> &fe_values,
-                const unsigned int i,
-                const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim+1)
-    tmp = fe_values.shape_grad (i,q);
-
-  return tmp;
-}
-@endcode
-
 @code
   {
     Vector<double> xx(dof_handler.n_dofs());
index fd58cac8eaffc63d43ca65d76f5aec81a081cfc2..c19973e26970b51602202fc43ecd8db001d48274 100644 (file)
@@ -44,7 +44,6 @@
 #include <fe/fe_system.h>
 #include <fe/fe_values.h>
 #include <fe/mapping_q1.h>
-#include <fe/mapping_c1.h>
 
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
@@ -232,122 +231,6 @@ RightHandSide<dim>::vector_value (const Point<dim> &p,
 
 
 
-
-
-template <int dim>
-Tensor<1,dim>
-extract_u (const FEValuesBase<dim> &fe_values,
-           const unsigned int i,
-           const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component < dim)
-    tmp[component] = fe_values.shape_value (i,q);
-
-  return tmp;
-}
-
-
-
-template <int dim>
-double
-extract_div_u (const FEValuesBase<dim> &fe_values,
-               const unsigned int i,
-               const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component < dim)
-    return fe_values.shape_grad (i,q)[component];
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-extract_grad_s_u (const FEValuesBase<dim> &fe_values,
-                 const unsigned int i,
-                 const unsigned int q)
-{
-  Tensor<2,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-  
-  if (component < dim)
-    {
-      const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2;
-      
-      for (unsigned int e=0; e<dim; ++e)
-       tmp[component][e] += grad_phi_over_2[e];
-      for (unsigned int d=0; d<dim; ++d)
-       tmp[d][component] += grad_phi_over_2[d];
-    }
-  
-  return tmp;
-}
-
-
-  
-template <int dim>
-double extract_p (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim)
-    return fe_values.shape_value (i,q);
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-double extract_T (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim+1)
-    return fe_values.shape_value (i,q);
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-extract_grad_T (const FEValuesBase<dim> &fe_values,
-                const unsigned int i,
-                const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim+1)
-    tmp = fe_values.shape_grad (i,q);
-
-  return tmp;
-}
-
-
-
-
 template <class Matrix, class Preconditioner>
 class InverseMatrix : public Subscriptor
 {
@@ -534,19 +417,6 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
 }
 
 
-template <int dim>
-double
-scalar_product (const Tensor<2,dim> &a,
-               const Tensor<2,dim> &b)
-{
-  double tmp = 0;
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=0; j<dim; ++j)
-      tmp += a[i][j] * b[i][j];
-  return tmp;
-}
-
-
 
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_system () 
@@ -592,6 +462,10 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
   const double Raleigh_number = 10;
   
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
+  const FEValuesExtractors::Scalar temperature (dim+1);
+
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -610,24 +484,27 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
 
-             const Tensor<1,dim> phi_i_u      = extract_u (fe_values, i, q);
+             const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
 
              if (rebuild_matrices)
                {
-                 const Tensor<2,dim> phi_i_grads_u= extract_grad_s_u (fe_values, i, q);
-                 const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
-                 const double        phi_i_p      = extract_p (fe_values, i, q);
-                 const double        phi_i_T      = extract_T (fe_values, i, q); 
-                 const Tensor<1,dim> grad_phi_i_T = extract_grad_T(fe_values, i, q);
+                 const SymmetricTensor<2,dim>
+                   phi_i_grads_u = fe_values[velocities].symmetric_gradient (i, q);
+                 const Tensor<1,dim> phi_i_u       = fe_values[velocities].value (i, q);
+                 const double        div_phi_i_u   = fe_values[velocities].divergence (i, q);
+                 const double        phi_i_p       = fe_values[pressure].value (i, q);
+                 const double        phi_i_T       = fe_values[temperature].value (i, q); 
+                 const Tensor<1,dim> grad_phi_i_T  = fe_values[temperature].gradient(i, q);
             
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
                    {
-                     const Tensor<2,dim> phi_j_grads_u     = extract_grad_s_u (fe_values, j, q);
-                     const double        div_phi_j_u = extract_div_u (fe_values, j, q);
-                     const double        phi_j_p     = extract_p (fe_values, j, q);
-                     const double        phi_j_T     = extract_T (fe_values, j, q);
+                     const SymmetricTensor<2,dim>
+                       phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q);
+                     const double        div_phi_j_u = fe_values[velocities].divergence (j, q);
+                     const double        phi_j_p     = fe_values[pressure].value (j, q);
+                     const double        phi_j_T     = fe_values[temperature].value (j, q);
                 
-                     local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_u)
+                     local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u
                                            - div_phi_i_u * phi_j_p
                                            - phi_i_p * div_phi_j_u
                                            + phi_i_p * phi_j_p
@@ -660,7 +537,7 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
               for (unsigned int i=0; i<dofs_per_cell; ++i)
                 {
                   const Tensor<1,dim>
-                    phi_i_u = extract_u (fe_face_values, i, q);
+                    phi_i_u = fe_face_values[velocities].value (i, q);
 
                   local_rhs(i) += -(phi_i_u *
                                     fe_face_values.normal_vector(q) *
@@ -791,6 +668,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
   TemperatureBoundaryValues<dim> temperature_boundary_values;
+  const FEValuesExtractors::Scalar temperature (dim+1);
   
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
@@ -817,8 +695,8 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
             for (unsigned int d=0; d<dim; ++d)
               present_div_u += present_solution_grads[q][d][d];            
            
-            const double        phi_i_T      = extract_T(fe_values, i, q);
-            const Tensor<1,dim> grad_phi_i_T = extract_grad_T(fe_values, i, q);
+            const double        phi_i_T      = fe_values[temperature].value (i, q);
+            const Tensor<1,dim> grad_phi_i_T = fe_values[temperature].gradient (i, q);
 
            const Point<dim> p = fe_values.quadrature_point(q);
            
@@ -901,7 +779,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                   old_solution_values_face[q](dim+1)
                                   :
                                   neighbor_temperature[q]) *
-                                 extract_T(fe_face_values,i,q) *
+                                 fe_face_values[temperature].value (i,q) *
                                  fe_face_values.JxW(q);
              }
          }
@@ -954,7 +832,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                         old_solution_values_face[q](dim+1)
                                         :
                                         neighbor_temperature[q]) *
-                                       extract_T(fe_face_values,i,q) *
+                                       fe_face_values[temperature].value (i,q) *
                                        fe_face_values.JxW(q);
                    }
                }
@@ -1014,7 +892,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                     old_solution_values_face[q](dim+1)
                                     :
                                     neighbor_temperature[q]) *
-                                   extract_T(fe_face_values,i,q) *
+                                   fe_face_values[temperature].value (i,q) *
                                    fe_face_values.JxW(q);
                }
            }      
@@ -1348,7 +1226,7 @@ void BoussinesqFlowProblem<dim>::run ()
        if (timestep_number % 10 == 0)
          refine_mesh ();
     }
-  while (time <= 500);
+  while (time <= 5);
 }
 
     
index 2486a6c5f028de7d45300e996d6eb3a31cca97ae..fedf82c278d2900ed9a9a5591ac0ed6e23e8bc55 100644 (file)
@@ -29,6 +29,10 @@ ParameterHandler class first used in @ref step_19 "step-19", which provides a
 convenient way for reading parameters from a configuration file at runtime
 without the need to recompile the program code.  
 
+The equations covered here fall into the class of vector-valued problems. A
+toplevel overview of this topic can be found in the @ref vector_valued module.
+
+
 <h3>Problem setting</h3>
 
 The original purpose of this program is to simulate the focussing properties
index 34a85e260b27427a83ade3d7f3bd1fb8f0de8e9d..54f260ad8bde9fd2ff82dcea08039de430d623c8 100644 (file)
@@ -29,6 +29,9 @@ To be well-posed, we will have to add boundary conditions to the
 equations. What boundary conditions are readily possible here will
 become clear once we discuss the weak form of the equations.
 
+The equations covered here fall into the class of vector-valued problems. A
+toplevel overview of this topic can be found in the @ref vector_valued module.
+
 
 <h2>Weak form</h2>
 
index b486221bc7608667593127a02f73394860aca38d..a726cb4ed3c11b5cefb9f273e5f563fdd7b9ad21 100644 (file)
@@ -249,8 +249,8 @@ class RightHandSide : public Function<dim>
 
 template <int dim>
 double
-RightHandSide<dim>::value (const Point<dim>  &p,
-                           const unsigned int component) const 
+RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                           const unsigned int /*component*/) const 
 {
   return 0;
 }
@@ -266,120 +266,6 @@ RightHandSide<dim>::vector_value (const Point<dim> &p,
 }
 
 
-
-                        // @sect3{extract_u and friends}
-
-                    // The next four functions are needed for
-                    // the assembly of the system matrix and 
-                    // the right hand side. They are very similar
-                    // to the ones used in step-20, except
-                    // that we are going to use Q(p+1)Qp elements
-                    // instead of divergence-free Raviart-Thomas
-                    // elements, which simplifies this procedure.
-                    // The only function that is new is 
-                    // <code>extract_grad_s_u</code>, which 
-                    // gets the symmetric gradient of u.
-                    // As discussed in the introduction, this
-                    // is a second-rank tensor, formed by
-                    // contributions from the gradient and its
-                    // transpose.
-template <int dim>
-Tensor<1,dim>
-extract_u (const FEValuesBase<dim> &fe_values,
-           const unsigned int i,
-           const unsigned int q)
-{
-  Tensor<1,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component < dim)
-    tmp[component] = fe_values.shape_value (i,q);
-
-  return tmp;
-}
-
-
-
-template <int dim>
-double
-extract_div_u (const FEValuesBase<dim> &fe_values,
-               const unsigned int i,
-               const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component < dim)
-    return fe_values.shape_grad (i,q)[component];
-  else
-    return 0;
-}
-
-
-template <int dim>
-Tensor<2,dim>
-extract_grad_s_u (const FEValuesBase<dim> &fe_values,
-                 const unsigned int i,
-                 const unsigned int q)
-{
-  Tensor<2,dim> tmp;
-
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-  
-  if (component < dim)
-    {
-      const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2;
-      
-      for (unsigned int e=0; e<dim; ++e)
-       tmp[component][e] += grad_phi_over_2[e];
-      for (unsigned int d=0; d<dim; ++d)
-       tmp[d][component] += grad_phi_over_2[d];
-    }
-  
-  return tmp;
-}
-
-
-  
-template <int dim>
-double extract_p (const FEValuesBase<dim> &fe_values,
-                  const unsigned int i,
-                  const unsigned int q)
-{
-  const unsigned int component
-    = fe_values.get_fe().system_to_component_index(i).first;
-
-  if (component == dim)
-    return fe_values.shape_value (i,q);
-  else
-    return 0;
-}
-
-                        // @sect4{Inner product of second-rank tensors}
-                        
-                    // In the assembly process, we will need
-                    // to form inner products of second-rank 
-                    // tensors. The way how to do this was
-                    // discussed in the introduction - just
-                    // take the sum of the product of the
-                    // individual entries.
-template <int dim>
-double
-scalar_product (const Tensor<2,dim> &a,
-                const Tensor<2,dim> &b)
-{
-  double tmp = 0;
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=0; j<dim; ++j)
-      tmp += a[i][j] * b[i][j];
-  return tmp;
-}
-
-
                         // @sect3{Linear solvers and preconditioners}
                         
                     // The linear solvers and preconditioners are 
@@ -746,7 +632,6 @@ void StokesProblem<dim>::setup_dofs ()
                     // hand side, and global
                     // numbers of the degrees of freedom
                     // for the present cell.
-
 template <int dim>
 void StokesProblem<dim>::assemble_system () 
 {
@@ -787,11 +672,18 @@ void StokesProblem<dim>::assemble_system ()
   const PressureBoundaryValues<dim> pressure_boundary_values;
   
   std::vector<double>               boundary_values (n_face_q_points);
-
+  
   const RightHandSide<dim>          right_hand_side;
   std::vector<Vector<double> >      rhs_values (n_q_points,
                                                 Vector<double>(dim+1));
 
+                                  // Next, we need two objects that work as
+                                  // extractors for the FEValues
+                                  // object. Their use is explained in detail
+                                  // in the report on @ref vector_valued :
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
+
                     // This starts the loop over all
                     // cells. With the abbreviations
                     // <code>extract_u</code> etc. 
@@ -810,20 +702,21 @@ void StokesProblem<dim>::assemble_system ()
                                         rhs_values);
       
       for (unsigned int q=0; q<n_q_points; ++q)
-      {
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-        {
-             const Tensor<2,dim> phi_i_grads_u= extract_grad_s_u (fe_values, i, q);
-             const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
-             const double        phi_i_p      = extract_p (fe_values, i, q);
+       {
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           {
+             const SymmetricTensor<2,dim>
+               phi_i_grads_u = fe_values[velocities].symmetric_gradient (i, q);
+             const Tensor<1,dim> phi_i_u       = fe_values[velocities].value (i, q);
+             const double        div_phi_i_u   = fe_values[velocities].divergence (i, q);
+             const double        phi_i_p       = fe_values[pressure].value (i, q);
           
           for (unsigned int j=0; j<dofs_per_cell; ++j)
           {
-            const Tensor<2,dim> phi_j_grads_u= extract_grad_s_u (fe_values, j, q);
-            const double        div_phi_j_u  = extract_div_u (fe_values, j, q);
-            const double        phi_j_p      = extract_p (fe_values, j, q);
-              
-            
+                 const SymmetricTensor<2,dim>
+                   phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q);
+                 const double        div_phi_j_u   = fe_values[velocities].divergence (j, q);
+                 const double        phi_j_p       = fe_values[pressure].value (j, q);
                     // Note the way we write the 
                     // contributions
                     // <code> phi_i_p * phi_j_p </code>,
@@ -837,11 +730,20 @@ void StokesProblem<dim>::assemble_system ()
                     // is only non-zero when all the
                     // other terms vanish (and the other
                     // way around).
-            local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_u)
+                                                  //
+                                                  // Note also that operator*
+                                                  // is overloaded for
+                                                  // symmetric tensors,
+                                                  // yielding the scalar
+                                                  // product between the two
+                                                  // tensors in the first
+                                                  // line:
+            local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u
                       - div_phi_i_u * phi_j_p
                       - phi_i_p * div_phi_j_u
                       + phi_i_p * phi_j_p)
                          * fe_values.JxW(q);     
+
           }
           const unsigned int component_i =
             fe.system_to_component_index(i).first;
@@ -871,15 +773,15 @@ void StokesProblem<dim>::assemble_system ()
 
             for (unsigned int q=0; q<n_face_q_points; ++q)
               for (unsigned int i=0; i<dofs_per_cell; ++i)
-              {
-                const Tensor<1,dim>
-                  phi_i_u = extract_u (fe_face_values, i, q);
-                
-                local_rhs(i) += -(phi_i_u *
-                                  fe_face_values.normal_vector(q) *
-                                  boundary_values[q] *
-                                  fe_face_values.JxW(q));
-              }
+                {
+                  const Tensor<1,dim>
+                    phi_i_u = fe_face_values[velocities].value (i, q);
+
+                  local_rhs(i) += -(phi_i_u *
+                                    fe_face_values.normal_vector(q) *
+                                    boundary_values[q] *
+                                    fe_face_values.JxW(q));
+                }
           }
 
                     // The final step is, as usual,
index c7e916b8dfd78815744a1d3e8b2d4d748d638ab8..7822a7bcc555906e1dc057cc6a35eec58bf109fd 100644 (file)
@@ -4,7 +4,8 @@
 
 In real life, most partial differential equations are really systems
 of equations. Accordingly, the solutions are usually
-vector-valued. The deal.II library supports such problems, and we will show
+vector-valued. The deal.II library supports such problems (see the
+extensive documentation in the @ref vector_valued module), and we will show
 that that is mostly rather simple. The only more complicated problems
 are in assembling matrix and right hand side, but these are easily
 understood as well. 
@@ -56,7 +57,7 @@ and the respective bilinear form is then
   \sum_{i,j}
   \left(
     \mu \partial_i u_j, \partial_i v_j
-  \right)_\Omega,
+  \right)_\Omega
   +
   \sum_{i,j}
   \left(
@@ -74,7 +75,7 @@ or also writing the first term a sum over components:
   \sum_{k,l}
   \left(
     \mu \partial_i u_j, \partial_i v_j
-  \right)_\Omega,
+  \right)_\Omega
   +
   \sum_{i,j}
   \left(
@@ -83,17 +84,27 @@ or also writing the first term a sum over components:
 @f]
 
 
-How do we now assemble the matrix for such an equation? The first thing we
-need is some knowledge about how the shape functions work in the case of
-vector-valued finite elements. Basically, this comes down to the following:
+How do we now assemble the matrix for such an equation? A very long answer
+with a number of different alternatives is given in the documentation of the
+@ref vector_valued module. Historically, the solution shown below was the only
+one available in the early years of the library. It turns out to also be the
+fastest. On the other hand, if a few per cent of compute time do not matter,
+there are simpler and probably more intuitive ways to assemble the linear
+system than the one discussed below but that weren't available until several
+years after this tutorial program was first written; if you are interested in
+them, take a look at the @ref vector_valued module.
+
+Let us go back to the question of how to assemble the linear system. The first
+thing we need is some knowledge about how the shape functions work in the case
+of vector-valued finite elements. Basically, this comes down to the following:
 let $n$ be the number of shape functions for the scalar finite element of
 which we build the vector element (for example, we will use bilinear functions
 for each component of the vector-valued finite element, so the scalar finite
-element is the <code>FE_Q(1)</code> element which we have used in previous examples
-already, and $n=4$ in two space dimensions). Further, let $N$ be the number of
-shape functions for the vector element; in two space dimensions, we need $n$
-shape functions for each component of the vector, so $N=2n$. Then, the $i$th
-shape function of the vector element has the form
+element is the <code>FE_Q(1)</code> element which we have used in previous
+examples already, and $n=4$ in two space dimensions). Further, let $N$ be the
+number of shape functions for the vector element; in two space dimensions, we
+need $n$ shape functions for each component of the vector, so $N=2n$. Then,
+the $i$th shape function of the vector element has the form
 @f[
   \Phi_i({\mathbf x}) = \varphi_{base(i)}({\mathbf x})\ {\mathbf e}_{comp(i)},
 @f]
@@ -311,13 +322,13 @@ Likewise, the contribution of cell $K$ to the right hand side vector is
 This is the form in which we will implement the local stiffness matrix and
 right hand side vectors.
 
-As a final note: in the @ref step_17 "step-17" example program, we will revisit the elastic
-problem laid out here, and will show how to solve it in parallel on a cluster
-of computers. The resulting program will thus be able to solve this problem to
-significantly higher accuracy, and more efficiently if this is
-required. In addition, in @ref step_20 "step-20", @ref step_21 "step-21", and
-@ref step_22 "step-22", we will revisit some 
-vector-valued problems and show a few techniques that may make it
+As a final note: in the @ref step_17 "step-17" example program, we will
+revisit the elastic problem laid out here, and will show how to solve it in
+parallel on a cluster of computers. The resulting program will thus be able to
+solve this problem to significantly higher accuracy, and more efficiently if
+this is required. In addition, in @ref step_20 "step-20", @ref step_21
+"step-21", as well as a few other of the later tutorial programs, we will
+revisit some vector-valued problems and show a few techniques that may make it
 simpler to actually go through all the stuff shown above, with
-<code>FiniteElement::system_to_component_index</code> etc.
+FiniteElement::system_to_component_index etc.
 
index 849142261c4e06ebfa81060193e9258e419bb604..6eea11aaf90226dbe0b4d0040a600591fdd0b893 100644 (file)
@@ -74,6 +74,7 @@ tests_x = block_matrices \
        create_* \
        line_coarsening_3d \
        no_flux_* \
+       fe_values_view_* \
        rt_hessian \
        rt_covariant
 
diff --git a/tests/deal.II/fe_values_view_01.cc b/tests/deal.II/fe_values_view_01.cc
new file mode 100644 (file)
index 0000000..9a6b189
--- /dev/null
@@ -0,0 +1,115 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2007, 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// test the FEValues views and extractor classes. these tests use a primitive
+// finite element and scalar extractors
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <fstream>
+
+
+
+template<int dim>
+void test (const Triangulation<dim>& tr,
+          const FiniteElement<dim>& fe)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  deallog << "FE=" << fe.get_name()
+         << std::endl;
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                          update_values | update_gradients | update_hessians);
+  for (typename DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active(); cell != dof.end(); ++cell)
+    {
+      fe_values.reinit (cell);
+  
+      for (unsigned int c=0; c<fe.n_components(); ++c)
+       {
+         FEValuesExtractors::Scalar single_component (c);
+
+         for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+           for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+             {
+               deallog << "i=" << i << ", q=" << q << std::endl;
+               deallog << "   "
+                       << fe_values[single_component].value (i,q) << ' '
+                       << fe_values[single_component].gradient (i,q) << ' '
+                       << fe_values[single_component].hessian (i,q) << std::endl;
+
+               Assert (fe_values[single_component].value (i,q)
+                       ==
+                       fe_values.shape_value_component (i,q,c),
+                       ExcInternalError());
+
+               Assert (fe_values[single_component].gradient (i,q)
+                       ==
+                       fe_values.shape_grad_component (i,q,c),
+                       ExcInternalError());
+
+               Assert (fe_values[single_component].hessian (i,q)
+                       ==
+                       fe_values.shape_hessian_component (i,q,c),
+                       ExcInternalError());
+             }
+       }
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_sphere()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_ball(tr);
+
+  static const HyperBallBoundary<dim> boundary;
+  tr.set_boundary (0, boundary);
+
+  FESystem<dim> fe (FE_Q<dim>(1), 1,
+                   FE_Q<dim>(2), 2,
+                   FE_DGQ<dim>(3), dim);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("fe_values_view_01/output");
+  deallog << std::setprecision (2);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-12);
+
+  test_hyper_sphere<2>();
+  test_hyper_sphere<3>();
+}
diff --git a/tests/deal.II/fe_values_view_02.cc b/tests/deal.II/fe_values_view_02.cc
new file mode 100644 (file)
index 0000000..23ca64d
--- /dev/null
@@ -0,0 +1,137 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2007, 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// test the FEValues views and extractor classes. these tests use a primitive
+// finite element and vector extractors
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <fstream>
+#include <iostream>
+
+
+template<int dim>
+void test (const Triangulation<dim>& tr,
+          const FiniteElement<dim>& fe)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  deallog << "FE=" << fe.get_name()
+         << std::endl;
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                          update_values | update_gradients | update_hessians);
+  for (typename DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active(); cell != dof.end(); ++cell)
+    {
+      fe_values.reinit (cell);
+  
+      for (unsigned int c=0; c<fe.n_components(); ++c)
+                                        // use a vector extractor if there
+                                        // are sufficiently many components
+                                        // left after the current component
+                                        // 'c'
+       if (c+dim <= fe.n_components())
+         {
+           FEValuesExtractors::Vector vec_components (c);
+
+           for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+             for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+               {
+                 deallog << "i=" << i << ", q=" << q << std::endl;
+                 deallog << "   "
+                         << fe_values[vec_components].value (i,q) << ' '
+                         << fe_values[vec_components].gradient (i,q) << ' '
+                         << fe_values[vec_components].divergence (i,q) << ' '
+                         << fe_values[vec_components].symmetric_gradient (i,q) << ' '
+                         << fe_values[vec_components].hessian (i,q)
+                         << std::endl;
+                 
+                 for (unsigned int d=0; d<dim; ++d)
+                   {
+                     Assert (fe_values[vec_components].value (i,q)[d]
+                             ==
+                             fe_values.shape_value_component (i,q,c+d),
+                             ExcInternalError());
+
+                     Assert (fe_values[vec_components].gradient (i,q)[d]
+                             ==
+                             fe_values.shape_grad_component (i,q,c+d),
+                             ExcInternalError());
+
+                     Assert (fe_values[vec_components].symmetric_gradient (i,q)
+                             ==
+                             (fe_values[vec_components].gradient(i,q) +
+                              transpose(fe_values[vec_components].gradient(i,q)))/2,
+                             ExcInternalError());
+                     
+                     Assert (fe_values[vec_components].hessian (i,q)[d]
+                             ==
+                             fe_values.shape_hessian_component (i,q,c+d),
+                             ExcInternalError());
+                   }
+
+                 Assert (fe_values[vec_components].divergence (i,q)
+                         ==
+                         trace (fe_values[vec_components].gradient (i,q)),
+                         ExcInternalError());
+               }
+         }
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_sphere()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_ball(tr);
+
+  static const HyperBallBoundary<dim> boundary;
+  tr.set_boundary (0, boundary);
+
+  FESystem<dim> fe (FE_Q<dim>(1), 1,
+                   FE_Q<dim>(2), 2,
+                   FE_DGQ<dim>(3), dim);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("fe_values_view_02/output");
+  deallog << std::setprecision (2);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-12);
+
+  test_hyper_sphere<2>();
+  test_hyper_sphere<3>();
+}
diff --git a/tests/deal.II/fe_values_view_03.cc b/tests/deal.II/fe_values_view_03.cc
new file mode 100644 (file)
index 0000000..cbb9c62
--- /dev/null
@@ -0,0 +1,117 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2007, 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// test the FEValues views and extractor classes. these tests use a non-primitive
+// finite element and scalar extractors
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_nedelec.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <fstream>
+
+
+
+template<int dim>
+void test (const Triangulation<dim>& tr,
+          const FiniteElement<dim>& fe)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  deallog << "FE=" << fe.get_name()
+         << std::endl;
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                          update_values | update_gradients | update_hessians);
+  for (typename DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active(); cell != dof.end(); ++cell)
+    {
+      fe_values.reinit (cell);
+  
+      for (unsigned int c=0; c<fe.n_components(); ++c)
+       {
+         FEValuesExtractors::Scalar single_component (c);
+
+         for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+           for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+             {
+               deallog << "i=" << i << ", q=" << q << std::endl;
+               deallog << "   "
+                       << fe_values[single_component].value (i,q) << ' '
+                       << fe_values[single_component].gradient (i,q) << ' '
+                       << fe_values[single_component].hessian (i,q) << std::endl;
+
+               Assert (fe_values[single_component].value (i,q)
+                       ==
+                       fe_values.shape_value_component (i,q,c),
+                       ExcInternalError());
+
+               Assert (fe_values[single_component].gradient (i,q)
+                       ==
+                       fe_values.shape_grad_component (i,q,c),
+                       ExcInternalError());
+
+               Assert (fe_values[single_component].hessian (i,q)
+                       ==
+                       fe_values.shape_hessian_component (i,q,c),
+                       ExcInternalError());
+             }
+       }
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_sphere()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_ball(tr);
+
+  static const HyperBallBoundary<dim> boundary;
+  tr.set_boundary (0, boundary);
+
+  FESystem<dim> fe (FE_Q<dim>(1), 1,
+                   FE_RaviartThomas<dim>(1), 1,
+                   FE_Nedelec<dim>(1), 1);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("fe_values_view_03/output");
+  deallog << std::setprecision (2);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-12);
+
+  test_hyper_sphere<2>();
+  test_hyper_sphere<3>();
+}
diff --git a/tests/deal.II/fe_values_view_04.cc b/tests/deal.II/fe_values_view_04.cc
new file mode 100644 (file)
index 0000000..9368301
--- /dev/null
@@ -0,0 +1,133 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2007, 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// test the FEValues views and extractor classes. these tests use a non-primitive
+// finite element and vector extractors
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_nedelec.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <fstream>
+
+
+
+template<int dim>
+void test (const Triangulation<dim>& tr,
+          const FiniteElement<dim>& fe)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  deallog << "FE=" << fe.get_name()
+         << std::endl;
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                          update_values | update_gradients | update_hessians);
+  for (typename DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active(); cell != dof.end(); ++cell)
+    {
+      fe_values.reinit (cell);
+  
+      for (unsigned int c=0; c<fe.n_components(); ++c)
+                                        // use a vector extractor if there
+                                        // are sufficiently many components
+                                        // left after the current component
+                                        // 'c'
+       if (c+dim <= fe.n_components())
+         {
+           FEValuesExtractors::Vector vec_components (c);
+
+           for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+             for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+               {
+                 deallog << "i=" << i << ", q=" << q << std::endl;
+                 deallog << "   "
+                         << fe_values[vec_components].value (i,q) << ' '
+                         << fe_values[vec_components].gradient (i,q) << ' '
+                         << fe_values[vec_components].symmetric_gradient (i,q)
+                         << fe_values[vec_components].hessian (i,q)
+                         << std::endl;
+
+                 for (unsigned int d=0; d<dim; ++d)
+                   {
+                     Assert (fe_values[vec_components].value (i,q)[d]
+                             ==
+                             fe_values.shape_value_component (i,q,c+d),
+                             ExcInternalError());
+
+                     Assert (fe_values[vec_components].gradient (i,q)[d]
+                             ==
+                             fe_values.shape_grad_component (i,q,c+d),
+                             ExcInternalError());
+
+                     Assert (fe_values[vec_components].symmetric_gradient (i,q)
+                             ==
+                             (fe_values[vec_components].gradient(i,q) +
+                              transpose(fe_values[vec_components].gradient(i,q)))/2,
+                             ExcInternalError());
+                     
+                     Assert (fe_values[vec_components].hessian (i,q)[d]
+                             ==
+                             fe_values.shape_hessian_component (i,q,c+d),
+                             ExcInternalError());
+                   }
+               }
+         }
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_sphere()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_ball(tr);
+
+  static const HyperBallBoundary<dim> boundary;
+  tr.set_boundary (0, boundary);
+
+  FESystem<dim> fe (FE_Q<dim>(1), 1,
+                   FE_RaviartThomas<dim>(1), 1,
+                   FE_Nedelec<dim>(1), 1);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("fe_values_view_04/output");
+  deallog << std::setprecision (2);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-12);
+
+  test_hyper_sphere<2>();
+  test_hyper_sphere<3>();
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.