]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Andrew McBride: implement extractors for symmetric rank 2 tensors.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 23 Sep 2009 17:48:36 +0000 (17:48 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 23 Sep 2009 17:48:36 +0000 (17:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@19506 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/fe/fe_values.inst.in

index 27fbd4d6254c9620a1870f2b560e04d97ba6f5d2..172e698a69528764129b574087f43a87124e0f8b 100644 (file)
@@ -55,8 +55,10 @@ template <int dim, int spacedim=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.
+ * element. There are extractors for single scalar components,
+ * vector components consisting of <code>dim</code> elements,
+ * and second order symmetric tensors consisting
+ * of <code>(dim*dim + dim)/2<code> components
  *
  * See the description of the @ref vector_valued module for examples how to
  * use the features of this namespace.
@@ -90,7 +92,7 @@ namespace FEValuesExtractors
       Scalar (const unsigned int component);
   };
 
-  
+
                                   /**
                                    * Extractor for a vector of
                                    * <code>dim</code> components of a
@@ -123,10 +125,53 @@ namespace FEValuesExtractors
                                        */
       Vector (const unsigned int first_vector_component);
   };
+
+
+                                  /**
+                                   * Extractor for a symmetric tensor
+                                   * of a rank specified by the
+                                   * template argument. For a second
+                                   * order symmetric tensor, this
+                                   * represents a collection of
+                                   * <code>(dim*dim + dim)/2</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
+                                   *
+                                   * @author Andrew McBride, 2009
+                                   */
+  template <int rank>
+  struct SymmetricTensor
+  {
+                                      /**
+                                       * The first component of the tensor
+                                       * view.
+                                       */
+      unsigned int first_tensor_component;
+
+                                      /**
+                                       * Constructor. Take the first
+                                       * component of the selected tensor
+                                       * inside the FEValues object as
+                                       * argument.
+                                       */
+      SymmetricTensor (const unsigned int first_tensor_component);
+  };
 }
 
 
 
+
 /**
  * A namespace for "views" on a FEValues, FEFaceValues, or FESubfaceValues
  * object. A view represents only a certain part of the whole: whereas the
@@ -134,8 +179,10 @@ namespace FEValuesExtractors
  * 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.
+ * There are classes that present views for single scalar components,
+ * vector components consisting of <code>dim</code> elements, and
+ * symmetric second order tensor components
+ * consisting of <code>(dim*dim + dim)/2</code> elements
  *
  * See the description of the @ref vector_valued module for examples how to
  * use the features of this namespace.
@@ -152,7 +199,7 @@ namespace FEValuesViews
                                    *
                                    * @ingroup feaccess vector_valued
                                    */
-  template <int dim, int spacedim=dim>  
+  template <int dim, int spacedim=dim>
   class Scalar
   {
     public:
@@ -190,7 +237,7 @@ namespace FEValuesViews
                                        * invalid object.
                                        */
       Scalar ();
-      
+
                                       /**
                                        * Constructor for an object that
                                        * represents a single scalar component
@@ -208,7 +255,7 @@ namespace FEValuesViews
                                        * this function is called.
                                        */
       Scalar & operator= (const Scalar<dim,spacedim> &);
-      
+
                                       /**
                                        * Return the value of the vector
                                        * component selected by this view, for
@@ -330,12 +377,12 @@ namespace FEValuesViews
       template <class InputVector>
       void get_function_laplacians (const InputVector& fe_function,
                                    std::vector<value_type>& laplacians) const;
-      
-      
+
+
     private:
                                       /**
                                        * A reference to the FEValuesBase object
-                                       * we operator on.
+                                       * we operate on.
                                        */
       const FEValuesBase<dim,spacedim> &fe_values;
 
@@ -394,7 +441,7 @@ namespace FEValuesViews
                                        */
       std::vector<ShapeFunctionData> shape_function_data;
   };
-  
+
 
 
                                   /**
@@ -406,7 +453,7 @@ namespace FEValuesViews
                                    *
                                    * @ingroup feaccess vector_valued
                                    */
-  template <int dim, int spacedim=dim>  
+  template <int dim, int spacedim=dim>
   class Vector
   {
     public:
@@ -438,7 +485,7 @@ namespace FEValuesViews
                                        * symmetrized gradient is a
                                        * <code>SymmetricTensor@<2,spacedim@></code>.
                                        */
-      typedef SymmetricTensor<2,spacedim> symmetric_gradient_type;
+      typedef dealii::SymmetricTensor<2,spacedim> symmetric_gradient_type;
 
                                       /**
                                        * A typedef for the type of the
@@ -449,7 +496,7 @@ namespace FEValuesViews
                                        * course is a scalar.
                                        */
       typedef double                 divergence_type;
-      
+
                                       /**
                                        * A typedef for the type of second
                                        * derivatives of the view this class
@@ -625,7 +672,7 @@ namespace FEValuesViews
       template <class InputVector>
       void get_function_symmetric_gradients (const InputVector& fe_function,
                    std::vector<symmetric_gradient_type>& symmetric_gradients) const;
-      
+
                                       /**
                                        * Return the divergence of the selected
                                        * vector components of the finite
@@ -648,7 +695,7 @@ namespace FEValuesViews
       template <class InputVector>
       void get_function_divergences (const InputVector& fe_function,
                                     std::vector<divergence_type>& divergences) const;
-      
+
                                       /**
                                        * Return the Hessians of the selected
                                        * vector components of the finite
@@ -695,11 +742,11 @@ namespace FEValuesViews
       template <class InputVector>
       void get_function_laplacians (const InputVector& fe_function,
                                    std::vector<value_type>& laplacians) const;
-      
+
     private:
                                       /**
                                        * A reference to the FEValuesBase object
-                                       * we operator on.
+                                       * we operate on.
                                        */
       const FEValuesBase<dim,spacedim> &fe_values;
 
@@ -717,7 +764,7 @@ namespace FEValuesViews
                                        * much cheaper.
                                        */
       struct ShapeFunctionData
-      {      
+      {
                                           /**
                                            * For each pair (shape
                                            * function,component within
@@ -780,6 +827,251 @@ namespace FEValuesViews
                                        */
       std::vector<ShapeFunctionData> shape_function_data;
   };
+
+
+  template <int rank, int dim, int spacedim = dim>
+  class SymmetricTensor;
+
+                                  /**
+                                   * A class representing a view to a set of
+                                   * <code>(dim*dim + dim)/2</code> components forming a
+                                   * symmetric second-order tensor from a
+                                   * vector-valued finite
+                                   * element. Views are discussed in the
+                                   * @ref vector_valued module.
+                                   *
+                                   * @ingroup feaccess vector_valued
+                                   *
+                                   * @author Andrew McBride, 2009
+                                   */
+  template <int dim, int spacedim>
+  class SymmetricTensor<2,dim,spacedim>
+  {
+    public:
+                                      /**
+                                       * A typedef for the data type of
+                                       * values of the view this class
+                                       * represents. Since we deal with a set
+                                       * of <code>(dim*dim + dim)/2</code> components
+                                       * (i.e. the unique components of a symmetric second-order
+                                       * tensor), the value type is a SymmetricTensor<2,spacedim>.
+                                       */
+      typedef dealii::SymmetricTensor<2, spacedim> value_type;
+
+                                      /**
+                                       * A typedef for the type of the
+                                       * divergence of the view this class
+                                       * represents. Here, for a set of
+                                       * of <code>(dim*dim + dim)/2</code> unique components
+                                       * of the finite element representing a symmetric second-order
+                                       * tensor, the divergence of
+                                       * course is a * <code>Tensor@<1,dim@></code>.
+                                       */
+      typedef Tensor<1, spacedim> divergence_type;
+
+                                      /**
+                                       * Default constructor. Creates an
+                                       * invalid object.
+                                       */
+      SymmetricTensor();
+
+                                      /**
+                                       * Constructor for an object that
+                                       * represents <code>(dim*dim + dim)/2</code>
+                                       * components of a
+                                       * FEValuesBase object (or of one of
+                                       * the classes derived from
+                                       * FEValuesBase), representing the unique
+                                       * components comprising a symmetric second-
+                                       * order tensor valued variable.
+                                       *
+                                       * The second argument denotes the
+                                       * index of the first component of the
+                                       * selected symmetric second order tensor.
+                                       */
+      SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values_base,
+                     const unsigned int first_tensor_component);
+
+                                      /**
+                                       * Copy operator. This is not a
+                                       * lightweight object so we don't allow
+                                       * copying and generate an exception if
+                                       * this function is called.
+                                       */
+      SymmetricTensor & operator=(const SymmetricTensor<2, dim, spacedim> &);
+
+                                      /**
+                                       * 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*dim + dim)/2</code> components
+                                       * (the unique components of a symmetric second-order tensor),
+                                       * the return type is a symmetric tensor of rank 2.
+                                       */
+      value_type
+      value(const unsigned int shape_function,
+           const unsigned int q_point) const;
+
+
+                                      /**
+                                       * Return the vector 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 values of the selected
+                                       * vector components of the finite
+                                       * element function characterized by
+                                       * <tt>fe_function</tt> at the
+                                       * quadrature points of the cell, face
+                                       * or subface selected the last time
+                                       * the <tt>reinit</tt> function of the
+                                       * FEValues object was called, at the
+                                       * quadrature points.
+                                       *
+                                       * This function is the equivalent of
+                                       * the
+                                       * FEValuesBase::get_function_values
+                                       * function but it only works on the
+                                       * selected vector components.
+                                       */
+      template <class InputVector>
+      void get_function_values(const InputVector& fe_function,
+                              std::vector<value_type>& values) const;
+
+                                      /**
+                                       * Return the divergence of the selected
+                                       * vector components of the finite
+                                       * element function characterized by
+                                       * <tt>fe_function</tt> at the
+                                       * quadrature points of the cell, face
+                                       * or subface selected the last time
+                                       * the <tt>reinit</tt> function of the
+                                       * FEValues object was called, at the
+                                       * quadrature points.
+                                       *
+                                       * There is no equivalent function such
+                                       * as
+                                       * FEValuesBase::get_function_gradients
+                                       * in the FEValues classes but the
+                                       * information can be obtained from
+                                       * FEValuesBase::get_function_gradients,
+                                       * of course.
+                                       */
+      template <class InputVector>
+      void get_function_divergences(const InputVector& fe_function,
+                                   std::vector<divergence_type>& divergences) const;
+
+    private:
+                                      /**
+                                       * A reference to the FEValuesBase object
+                                       * we operate on.
+                                       */
+      const FEValuesBase<dim, spacedim> &fe_values;
+
+                                      /**
+                                       * The first component of the vector
+                                       * this view represents of the
+                                       * FEValuesBase object.
+                                       */
+      const unsigned int first_tensor_component;
+
+                                      /**
+                                       * The number of unique
+                                       * components of a symmetric
+                                       * second-order tensor
+                                       */
+      static const unsigned int n_independent_components = (dim*dim + dim)/2;
+
+                                      /**
+                                       * The index of the first-order
+                                       * tensor representation of a
+                                       * symmetric second-order
+                                       * tensor, stored as the
+                                       * components of a tensor
+                                       */
+      std::vector < std::vector<unsigned int> > vector_to_symmetric_tensor_data;
+
+                                      /**
+                                       * A structure where for each shape
+                                       * function we pre-compute a bunch of
+                                       * data that will make later accesses
+                                       * much cheaper.
+                                       */
+      struct ShapeFunctionData
+      {
+                                          /**
+                                           * For each pair (shape
+                                           * function,component within
+                                           * vector), store whether the
+                                           * selected vector component may be
+                                           * nonzero. For primitive shape
+                                           * functions we know for sure
+                                           * whether a certain scalar
+                                           * component of a given shape
+                                           * function is nonzero, whereas for
+                                           * non-primitive shape functions
+                                           * this may not be entirely clear
+                                           * (e.g. for RT elements it depends
+                                           * on the shape of a cell).
+                                           */
+         bool is_nonzero_shape_function_component[n_independent_components];
+
+                                          /**
+                                           * For each pair (shape function,
+                                           * component within vector), store
+                                           * the row index within the
+                                           * shape_values, shape_gradients,
+                                           * and shape_hessians tables (the
+                                           * column index is the quadrature
+                                           * point index). If the shape
+                                           * function is primitive, then we
+                                           * can get this information from
+                                           * the shape_function_to_row_table
+                                           * of the FEValues object;
+                                           * otherwise, we have to work a bit
+                                           * harder to compute this
+                                           * information.
+                                           */
+         unsigned int row_index[n_independent_components];
+
+                                          /**
+                                           * For each shape function say the
+                                           * following: if only a single
+                                           * entry in
+                                           * is_nonzero_shape_function_component
+                                           * for this shape function is
+                                           * nonzero, then store the
+                                           * corresponding value of row_index
+                                           * and
+                                           * single_nonzero_component_index
+                                           * represents the index between 0
+                                           * and (dim^2 + dim)/2 for which it is
+                                           * attained. If multiple components
+                                           * are nonzero, then store -1. If
+                                           * no components are nonzero then
+                                           * store -2.
+                                           */
+         int single_nonzero_component;
+         unsigned int single_nonzero_component_index;
+      };
+
+                                      /**
+                                       * Store the data about shape
+                                       * functions.
+                                       */
+      std::vector<ShapeFunctionData> shape_function_data;
+  };
 }
 
 
@@ -804,10 +1096,13 @@ namespace internal
     {
                                         /**
                                          * Caches for scalar and
-                                         * vector-valued views.
+                                         * vector, and symmetric second-order
+                                          * tensor valued views.
                                          */
        std::vector<dealii::FEValuesViews::Scalar<dim,spacedim> > scalars;
        std::vector<dealii::FEValuesViews::Vector<dim,spacedim> > vectors;
+        std::vector<dealii::FEValuesViews::SymmetricTensor<2,dim,spacedim> >
+               symmetric_second_order_tensors;
 
                                         /**
                                          * Constructor.
@@ -819,7 +1114,7 @@ namespace internal
 
 
 
-  
+
 /*!@addtogroup feaccess */
 /*@{*/
 
@@ -839,7 +1134,7 @@ namespace internal
  *
  * The purpose of this class is discussed
  * on the page on @ref UpdateFlagsEssay.
- * 
+ *
  * @author Guido Kanschat, 2000
  */
 template <int dim, int spacedim=dim>
@@ -979,7 +1274,7 @@ class FEValuesData
     std::vector<Tensor<3,spacedim> > jacobian_grads;
 
                                     /**
-                                     * Array of the inverse Jacobian matrices 
+                                     * Array of the inverse Jacobian matrices
                                      * at the quadrature points.
                                      */
     std::vector<Tensor<2,spacedim> > inverse_jacobians;
@@ -1001,7 +1296,7 @@ class FEValuesData
     std::vector<Point<spacedim> >  normal_vectors;
 
                                     /**
-                                     * List of outward vectors normal to the cell  
+                                     * List of outward vectors normal to the cell
                                      * surface (line) at the quadrature points
                                      * for the codimension 1 case,
                                      * when spacedim=3 (=2).
@@ -1027,14 +1322,14 @@ class FEValuesData
                                      * filled by Mapping.
                                      */
     HessianVector support_jacobians;
-    
+
                                     /**
                                      * Array of the inverse Jacobian of the
                                      * mapping in the support points,
                                      * filled by Mapping.
                                      */
     HessianVector support_inverse_jacobians;
-    
+
                                     /**
                                      * Indicate the first row which a
                                      * given shape function occupies
@@ -1065,7 +1360,7 @@ class FEValuesData
                                      * non-zero components.
                                      */
     std::vector<unsigned int> shape_function_to_row_table;
-    
+
                                      /**
                                      * Original update flags handed
                                      * to the constructor of
@@ -1139,7 +1434,7 @@ class FEValuesData
  *
  *  The functions of this class fall into different cathegories:
  *  <ul>
- *  <li> shape_value(), shape_grad(), etc: return one of the values 
+ *  <li> shape_value(), shape_grad(), etc: return one of the values
  *    of this object at a time. These functions are inlined, so this
  *    is the suggested access to all finite element values. There
  *    should be no loss in performance with an optimizing compiler. If
@@ -1158,13 +1453,13 @@ class FEValuesData
  *    non-zero component, since then the above functions cannot be used,
  *    and you have to walk over all (or only the non-zero) components of
  *    the shape function using this set of functions.
- *   
+ *
  *  <li> get_function_values(), get_function_gradients(), etc.: Compute a
  *    finite element function or its derivative in quadrature points.
  *
  *  <li> reinit: initialize the FEValues object for a certain cell.
  *    This function is not in the present class but only in the derived
- *    classes and has a variable call syntax. 
+ *    classes and has a variable call syntax.
  *    See the docs for the derived classes for more information.
  * </ul>
  *
@@ -1183,7 +1478,7 @@ class FEValuesData
  * give <b>all</b> needed update flags to FEValues.
  *
  * The mechanisms by which this class works is also discussed
- * on the page on @ref UpdateFlagsEssay. 
+ * on the page on @ref UpdateFlagsEssay.
  *
  * @author Wolfgang Bangerth, 1998, 2003, Guido Kanschat, 2001
  */
@@ -1208,7 +1503,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      */
     const unsigned int dofs_per_cell;
 
-    
+
                                     /**
                                      * Constructor. Set up the array
                                      * sizes with <tt>n_q_points</tt>
@@ -1264,12 +1559,26 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     const FEValuesViews::Vector<dim,spacedim> &
     operator[] (const FEValuesExtractors::Vector &vector) const;
 
+                                    /**
+                                     * Create a view of the current FEValues
+                                     * object that represents a set of
+                                     * <code>(dim*dim + dim)/2</code> scalar components
+                                     * (i.e. a symmetric 2nd order tensor)
+                                      * 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.
+                                     */
+    const FEValuesViews::SymmetricTensor<2,dim,spacedim> &
+    operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const;
+
                                     //@}
 
-    
+
                                     /// @name ShapeAccess Access to shape function values
                                     //@{
-    
+
                                     /**
                                      * Value of a shape function at a
                                      * quadrature point on the cell,
@@ -1470,7 +1779,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     shape_hessian_component (const unsigned int function_no,
                             const unsigned int point_no,
                             const unsigned int component) const;
-    
+
                                     /**
                                      * @deprecated Wrapper for shape_hessian_component()
                                      */
@@ -1478,12 +1787,12 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     shape_2nd_derivative_component (const unsigned int function_no,
                                    const unsigned int point_no,
                                    const unsigned int component) const;
-    
+
 
                                     //@}
                                     /// @name FunctionAccess Access to values of global finite element functions
                                     //@{
-    
+
                                     /**
                                      * Returns the values of the
                                      * finite element function
@@ -1512,10 +1821,10 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * it is a vector-valued one, then use
                                      * the other get_function_values()
                                      * function.
-                                     * 
+                                     *
                                      * The function assumes that the
                                      * <tt>values</tt> object already has the
-                                     * correct size. 
+                                     * correct size.
                                      *
                                      * The actual data type of the
                                      * input vector may be either a
@@ -1564,7 +1873,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                     /**
                                      * Generate function values from
                                      * an arbitrary vector.
-                                     * 
+                                     *
                                      * This function offers the
                                      * possibility to extract
                                      * function values in quadrature
@@ -1721,7 +2030,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * it is a vector-valued one, then use
                                      * the other get_function_gradients()
                                      * function.
-                                     * 
+                                     *
                                      * The function assumes that the
                                      * @p gradients object already has the
                                      * right size.
@@ -1799,7 +2108,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     void get_function_gradients (const InputVector               &fe_function,
                                 std::vector<std::vector<Tensor<1,spacedim> > > &gradients) const;
 
-    
+
                                     /**
                                      * @deprecated Use
                                      * get_function_gradients() instead.
@@ -1869,7 +2178,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * the other
                                      * get_function_hessians()
                                      * function.
-                                     * 
+                                     *
                                      * The actual data type of the
                                      * input vector may be either a
                                      * Vector&lt;T&gt;,
@@ -1895,7 +2204,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     get_function_hessians (const InputVector& fe_function,
                           std::vector<Tensor<2,spacedim> >& hessians) const;
 
-    
+
                                     /**
                                      * Compute the tensor of second
                                      * derivatives of the finite
@@ -1980,7 +2289,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     get_function_2nd_derivatives (const InputVector&,
                                  std::vector<std::vector<Tensor<2,spacedim> > >&,
                                  bool = false) const;
-    
+
 
                                     /**
                                      * Compute the (scalar) Laplacian
@@ -2008,7 +2317,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * the other
                                      * get_function_laplacians()
                                      * function.
-                                     * 
+                                     *
                                      * The actual data type of the
                                      * input vector may be either a
                                      * Vector&lt;T&gt;,
@@ -2035,7 +2344,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     get_function_laplacians (const InputVector& fe_function,
                             std::vector<number>& laplacians) const;
 
-    
+
                                     /**
                                      * Compute the (scalar) Laplacian
                                      * of the finite element function
@@ -2128,7 +2437,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
 
 
                                     //@}
-    
+
                                     /**
                                      * Position of the <tt>i</tt>th
                                      * quadrature point in real space.
@@ -2171,7 +2480,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * the values returned by JxW().
                                      */
     const std::vector<double> & get_JxW_values () const;
-    
+
                                     /**
                                      * Return the Jacobian of the
                                      * transformation at the specified
@@ -2212,7 +2521,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
 
                                     /**
                                      * Pointer to the array holding
-                                     * the values returned by 
+                                     * the values returned by
                                      * inverse_jacobian().
                                      */
     const std::vector<Tensor<2,spacedim> > & get_inverse_jacobians () const;
@@ -2235,7 +2544,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * for this object.
                                      */
     UpdateFlags get_update_flags () const;
-    
+
                                     /**
                                      * Return a triangulation
                                      * iterator to the current cell.
@@ -2324,7 +2633,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * @ingroup Exceptions
                                      */
     DeclException0 (ExcFENotPrimitive);
-    
+
   protected:
                                      /**
                                       * Objects of the FEValues
@@ -2385,7 +2694,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                       *
                                       * @author Wolfgang Bangerth, 2003
                                       */
-    class CellIteratorBase 
+    class CellIteratorBase
     {
       public:
                                          /**
@@ -2395,7 +2704,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                           * class.
                                           */
         virtual ~CellIteratorBase ();
-        
+
                                          /**
                                           * Conversion operator to an
                                           * iterator for
@@ -2568,7 +2877,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      */
     template <typename CI> class CellIterator;
     class TriaCellIterator;
-    
+
                                     /**
                                      * Store the cell selected last time the
                                      * reinit() function was called.  This is
@@ -2578,18 +2887,18 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * the extractor classes.
                                      */
     std::auto_ptr<const CellIteratorBase> present_cell;
-    
+
                                     /**
                                      * Storage for the mapping object.
                                      */
     const SmartPointer<const Mapping<dim,spacedim> > mapping;
-    
+
                                     /**
                                      * Store the finite element for later use.
                                      */
     const SmartPointer<const FiniteElement<dim,spacedim> > fe;
 
-    
+
                                     /**
                                      * Internal data of mapping.
                                      */
@@ -2649,7 +2958,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                       * and also do not implement it.
                                       */
     FEValuesBase (const FEValuesBase &);
-    
+
                                      /**
                                       * Copy operator. Since
                                       * objects of this class are not
@@ -2663,7 +2972,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * objects.
                                      */
     internal::FEValuesViews::Cache<dim,spacedim> fe_values_views_cache;
-    
+
                                     /**
                                      * Make the view classes friends of this
                                      * class, since they access internal
@@ -2671,6 +2980,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      */
     template <int, int> friend class FEValuesViews::Scalar;
     template <int, int> friend class FEValuesViews::Vector;
+    template <int, int, int> friend class FEValuesViews::SymmetricTensor;
 };
 
 
@@ -2681,7 +2991,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
  * This function implements the initialization routines for
  * FEValuesBase, if values in quadrature points of a cell are
  * needed. For further documentation see this class.
- * 
+ *
  * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2001
  */
 template <int dim, int spacedim=dim>
@@ -2707,7 +3017,7 @@ class FEValues : public FEValuesBase<dim,spacedim>
                                       * to <code>dim</code>.
                                       */
     static const unsigned int integral_dimension = dim;
-    
+
                                     /**
                                      * Constructor. Gets cell
                                      * independent data from mapping
@@ -2727,7 +3037,7 @@ class FEValues : public FEValuesBase<dim,spacedim>
     FEValues (const FiniteElement<dim,spacedim> &fe,
              const Quadrature<dim>             &quadrature,
              const UpdateFlags                  update_flags);
-    
+
                                     /**
                                      * Reinitialize the gradients,
                                      * Jacobi determinants, etc for
@@ -2807,7 +3117,7 @@ class FEValues : public FEValuesBase<dim,spacedim>
                                      * stored by this object.
                                      */
     const Quadrature<dim> & get_quadrature () const;
-    
+
                                     /**
                                      * Determine an estimate for the
                                      * memory consumption (in bytes)
@@ -2839,7 +3149,7 @@ class FEValues : public FEValuesBase<dim,spacedim>
                                       * FEValues/hpFEValues).
                                       */
     const FEValues<dim,spacedim> & get_present_fe_values () const;
-    
+
   private:
                                     /**
                                      * Store a copy of the quadrature
@@ -2889,7 +3199,7 @@ class FEValues : public FEValuesBase<dim,spacedim>
  *
  *  @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001
  */
-template <int dim, int spacedim=dim>  
+template <int dim, int spacedim=dim>
 class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
 {
   public:
@@ -2933,7 +3243,7 @@ class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
                                      * is normalized to one.
                                      */
     const Point<dim> & normal_vector (const unsigned int i) const;
-    
+
                                     /**
                                      * Boundary form of the
                                      * transformation of the cell at
@@ -2952,7 +3262,7 @@ class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
                                      * <tt>n.ds</tt>.
                                      */
     const Tensor<1,spacedim> & boundary_form (const unsigned int i) const;
-    
+
                                     /**
                                      * Return the list of outward normal
                                      * vectors to the cell at the
@@ -2981,14 +3291,14 @@ class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
                                      * stored by this object.
                                      */
     const Quadrature<dim-1> & get_quadrature () const;
-    
+
                                     /**
                                      * Determine an estimate for the
                                      * memory consumption (in bytes)
                                      * of this object.
                                      */
     unsigned int memory_consumption () const;
-    
+
   protected:
 
                                     /**
@@ -3020,7 +3330,7 @@ class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
  *
  * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001
  */
-template <int dim, int spacedim=dim>  
+template <int dim, int spacedim=dim>
 class FEFaceValues : public FEFaceValuesBase<dim,spacedim>
 {
   public:
@@ -3130,7 +3440,7 @@ class FEFaceValues : public FEFaceValuesBase<dim,spacedim>
                                      */
     void reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                 const unsigned int                                         face_no);
-    
+
                                      /**
                                       * Return a reference to this
                                       * very object.
@@ -3161,7 +3471,7 @@ class FEFaceValues : public FEFaceValuesBase<dim,spacedim>
                                      * Do work common to the two
                                      * constructors.
                                      */
-    void initialize (const UpdateFlags update_flags);    
+    void initialize (const UpdateFlags update_flags);
 
                                      /**
                                       * The reinit() functions do
@@ -3196,7 +3506,7 @@ class FEFaceValues : public FEFaceValuesBase<dim,spacedim>
  *
  * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2000, 2001
  */
-template <int dim, int spacedim=dim>  
+template <int dim, int spacedim=dim>
 class FESubfaceValues : public FEFaceValuesBase<dim,spacedim>
 {
   public:
@@ -3273,7 +3583,7 @@ class FESubfaceValues : public FEFaceValuesBase<dim,spacedim>
     void reinit (const typename hp::DoFHandler<dim,spacedim>::cell_iterator &cell,
                 const unsigned int                    face_no,
                 const unsigned int                    subface_no);
-    
+
                                     /**
                                      * Reinitialize the gradients,
                                      * Jacobi determinants, etc for
@@ -3290,7 +3600,7 @@ class FESubfaceValues : public FEFaceValuesBase<dim,spacedim>
     void reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
                 const unsigned int                    face_no,
                 const unsigned int                    subface_no);
-    
+
                                     /**
                                      * Reinitialize the gradients,
                                      * Jacobi determinants, etc for
@@ -3406,6 +3716,13 @@ namespace FEValuesExtractors
                  :
                  first_vector_component (first_vector_component)
   {}
+
+  template <int rank>
+  inline
+  SymmetricTensor<rank>::SymmetricTensor (const unsigned int first_tensor_component)
+                 :
+                 first_tensor_component (first_tensor_component)
+  {}
 }
 
 
@@ -3413,7 +3730,7 @@ namespace FEValuesExtractors
 
 namespace FEValuesViews
 {
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   inline
   typename Scalar<dim,spacedim>::value_type
   Scalar<dim,spacedim>::value (const unsigned int shape_function,
@@ -3423,7 +3740,7 @@ namespace FEValuesViews
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // an adaptation of the
                                     // FEValuesBase::shape_value_component
@@ -3438,11 +3755,11 @@ namespace FEValuesViews
     else
       return 0;
   }
-  
 
 
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   inline
   typename Scalar<dim,spacedim>::gradient_type
   Scalar<dim,spacedim>::gradient (const unsigned int shape_function,
@@ -3452,7 +3769,7 @@ namespace FEValuesViews
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // an adaptation of the
                                     // FEValuesBase::shape_grad_component
@@ -3469,7 +3786,7 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   inline
   typename Scalar<dim,spacedim>::hessian_type
   Scalar<dim,spacedim>::hessian (const unsigned int shape_function,
@@ -3479,7 +3796,7 @@ namespace FEValuesViews
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // an adaptation of the
                                     // FEValuesBase::shape_grad_component
@@ -3495,7 +3812,7 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   inline
   typename Vector<dim,spacedim>::value_type
   Vector<dim,spacedim>::value (const unsigned int shape_function,
@@ -3505,7 +3822,7 @@ namespace FEValuesViews
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // same as for the scalar case except
                                     // that we have one more index
@@ -3522,28 +3839,28 @@ namespace FEValuesViews
     else
       {
        value_type return_value;
-       for (unsigned int d=0; d<dim; ++d)    
+       for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point);
-       
+
        return return_value;
       }
   }
 
-  
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   inline
   typename Vector<dim,spacedim>::gradient_type
   Vector<dim,spacedim>::gradient (const unsigned int shape_function,
                                  const unsigned int q_point) const
   {
-    typedef FEValuesBase<dim,spacedim> FVB;    
+    typedef FEValuesBase<dim,spacedim> FVB;
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // same as for the scalar case except
                                     // that we have one more index
@@ -3560,7 +3877,7 @@ namespace FEValuesViews
     else
       {
        gradient_type return_value;
-       for (unsigned int d=0; d<dim; ++d)    
+       for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
@@ -3571,7 +3888,7 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   inline
   typename Vector<dim,spacedim>::divergence_type
   Vector<dim,spacedim>::divergence (const unsigned int shape_function,
@@ -3579,11 +3896,11 @@ namespace FEValuesViews
   {
                                     // this function works like in
                                     // the case above
-    typedef FEValuesBase<dim,spacedim> FVB; 
+    typedef FEValuesBase<dim,spacedim> FVB;
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // same as for the scalar case except
                                     // that we have one more index
@@ -3596,7 +3913,7 @@ namespace FEValuesViews
     else
       {
        divergence_type return_value = 0;
-       for (unsigned int d=0; d<dim; ++d)    
+       for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value
              += fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];
@@ -3604,22 +3921,22 @@ namespace FEValuesViews
        return return_value;
       }
   }
-  
 
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   inline
   typename Vector<dim,spacedim>::hessian_type
   Vector<dim,spacedim>::hessian (const unsigned int shape_function,
                                 const unsigned int q_point) const
   {
                                     // this function works like in
-                                    // the case above                          
-    typedef FEValuesBase<dim,spacedim> FVB;    
+                                    // the case above
+    typedef FEValuesBase<dim,spacedim> FVB;
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // same as for the scalar case except
                                     // that we have one more index
@@ -3636,7 +3953,7 @@ namespace FEValuesViews
     else
       {
        hessian_type return_value;
-       for (unsigned int d=0; d<dim; ++d)    
+       for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point];
@@ -3655,20 +3972,20 @@ namespace FEValuesViews
                                      * equal to zero.
                                      */
     inline
-    SymmetricTensor<2,1>
+    dealii::SymmetricTensor<2,1>
     symmetrize_single_row (const unsigned int n,
                           const Tensor<1,1> &t)
     {
       const unsigned int dim = 1;
       Assert (n < dim, ExcIndexRange (n, 0, dim));
-      
+
       const double array[1] = { t[0] };
-      return SymmetricTensor<2,1>(array);
+      return dealii::SymmetricTensor<2,1>(array);
     }
 
 
     inline
-    SymmetricTensor<2,2>
+    dealii::SymmetricTensor<2,2>
     symmetrize_single_row (const unsigned int n,
                           const Tensor<1,2> &t)
     {
@@ -3677,24 +3994,24 @@ namespace FEValuesViews
          case 0:
          {
            const double array[3] = { t[0], 0, t[1]/2 };
-           return SymmetricTensor<2,2>(array);
+           return dealii::SymmetricTensor<2,2>(array);
          }
          case 1:
          {
            const double array[3] = { 0, t[1], t[0]/2 };
-           return SymmetricTensor<2,2>(array);
+           return dealii::SymmetricTensor<2,2>(array);
          }
          default:
          {
            Assert (false, ExcIndexRange (n, 0, 2));
-           return SymmetricTensor<2,2>();
+           return dealii::SymmetricTensor<2,2>();
          }
        }
     }
 
 
     inline
-    SymmetricTensor<2,3>
+    dealii::SymmetricTensor<2,3>
     symmetrize_single_row (const unsigned int n,
                           const Tensor<1,3> &t)
     {
@@ -3703,27 +4020,27 @@ namespace FEValuesViews
          case 0:
          {
            const double array[6] = { t[0], 0, 0, t[1]/2, t[2]/2, 0 };
-           return SymmetricTensor<2,3>(array);
+           return dealii::SymmetricTensor<2,3>(array);
          }
          case 1:
          {
            const double array[6] = { 0, t[1], 0, t[0]/2, 0, t[2]/2 };
-           return SymmetricTensor<2,3>(array);
+           return dealii::SymmetricTensor<2,3>(array);
          }
          case 2:
          {
            const double array[6] = { 0, 0, t[2], 0, t[0]/2, t[1]/2 };
-           return SymmetricTensor<2,3>(array);
+           return dealii::SymmetricTensor<2,3>(array);
          }
          default:
          {
            Assert (false, ExcIndexRange (n, 0, 3));
-           return SymmetricTensor<2,3>();
+           return dealii::SymmetricTensor<2,3>();
          }
        }
     }
   }
-  
+
 
   template <int dim, int spacedim>
   inline
@@ -3731,11 +4048,11 @@ namespace FEValuesViews
   Vector<dim,spacedim>::symmetric_gradient (const unsigned int shape_function,
                                            const unsigned int q_point) const
   {
-    typedef FEValuesBase<dim,spacedim> FVB;    
+    typedef FEValuesBase<dim,spacedim> FVB;
     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 FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
 
                                     // same as for the scalar case except
                                     // that we have one more index
@@ -3748,7 +4065,7 @@ namespace FEValuesViews
     else
       {
        gradient_type return_value;
-       for (unsigned int d=0; d<dim; ++d)    
+       for (unsigned int d=0; d<dim; ++d)
          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
            return_value[d]
              = fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
@@ -3773,13 +4090,13 @@ operator[] (const FEValuesExtractors::Scalar &scalar) const
   Assert (scalar.component < fe_values_views_cache.scalars.size(),
          ExcIndexRange (scalar.component,
                         0, fe_values_views_cache.scalars.size()));
-  
+
   return fe_values_views_cache.scalars[scalar.component];
 }
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const FEValuesViews::Vector<dim,spacedim> &
 FEValuesBase<dim,spacedim>::
@@ -3789,13 +4106,25 @@ operator[] (const FEValuesExtractors::Vector &vector) const
          fe_values_views_cache.vectors.size(),
          ExcIndexRange (vector.first_vector_component,
                         0, fe_values_views_cache.vectors.size()));
-  
+
   return fe_values_views_cache.vectors[vector.first_vector_component];
 }
 
+template <int dim, int spacedim>
+inline
+const FEValuesViews::SymmetricTensor<2,dim,spacedim> &
+FEValuesBase<dim,spacedim>::
+operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const
+{
+  Assert (tensor.first_tensor_component <
+         fe_values_views_cache.symmetric_second_order_tensors.size(),
+         ExcIndexRange (tensor.first_tensor_component,
+                        0, fe_values_views_cache.symmetric_second_order_tensors.size()));
 
+  return fe_values_views_cache.symmetric_second_order_tensors[tensor.first_tensor_component];
+}
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const double &
 FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
@@ -3827,7 +4156,7 @@ FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 double
 FEValuesBase<dim,spacedim>::shape_value_component (const unsigned int i,
@@ -3840,7 +4169,7 @@ FEValuesBase<dim,spacedim>::shape_value_component (const unsigned int i,
          ExcAccessToUninitializedField());
   Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
-                       
+
                                   // if this particular shape
                                   // function is primitive, then we
                                   // can take a short-cut by checking
@@ -3891,7 +4220,7 @@ FEValuesBase<dim,spacedim>::shape_value_component (const unsigned int i,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Tensor<1,spacedim> &
 FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
@@ -3927,7 +4256,7 @@ FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
 
 
 
-template <int dim, int spacedim> 
+template <int dim, int spacedim>
 inline
 Tensor<1,spacedim>
 FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
@@ -3940,7 +4269,7 @@ FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
          ExcAccessToUninitializedField());
   Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
-                       
+
                                   // if this particular shape
                                   // function is primitive, then we
                                   // can take a short-cut by checking
@@ -3991,7 +4320,7 @@ FEValuesBase<dim,spacedim>::shape_grad_component (const unsigned int i,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Tensor<2,spacedim> &
 FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
@@ -4027,7 +4356,7 @@ FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Tensor<2,spacedim> &
 FEValuesBase<dim,spacedim>::shape_2nd_derivative (const unsigned int i,
@@ -4038,7 +4367,7 @@ FEValuesBase<dim,spacedim>::shape_2nd_derivative (const unsigned int i,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 Tensor<2,spacedim>
 FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
@@ -4051,7 +4380,7 @@ FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
          ExcAccessToUninitializedField());
   Assert (component < fe->n_components(),
          ExcIndexRange(component, 0, fe->n_components()));
-                       
+
                                   // if this particular shape
                                   // function is primitive, then we
                                   // can take a short-cut by checking
@@ -4102,7 +4431,7 @@ FEValuesBase<dim,spacedim>::shape_hessian_component (const unsigned int i,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 Tensor<2,spacedim>
 FEValuesBase<dim,spacedim>::shape_2nd_derivative_component (const unsigned int i,
@@ -4114,18 +4443,18 @@ FEValuesBase<dim,spacedim>::shape_2nd_derivative_component (const unsigned int i
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
-const FiniteElement<dim,spacedim> & 
+const FiniteElement<dim,spacedim> &
 FEValuesBase<dim,spacedim>::get_fe () const
 {
   return *fe;
 }
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
-const Mapping<dim,spacedim> & 
+const Mapping<dim,spacedim> &
 FEValuesBase<dim,spacedim>::get_mapping () const
 {
   return *mapping;
@@ -4133,7 +4462,7 @@ FEValuesBase<dim,spacedim>::get_mapping () const
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 UpdateFlags
 FEValuesBase<dim,spacedim>::get_update_flags () const
@@ -4143,7 +4472,7 @@ FEValuesBase<dim,spacedim>::get_update_flags () const
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 const typename Triangulation<dim,spacedim>::cell_iterator
 FEValuesBase<dim,spacedim>::get_cell () const
 {
@@ -4151,7 +4480,7 @@ FEValuesBase<dim,spacedim>::get_cell () const
 }
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const std::vector<Point<spacedim> > &
 FEValuesBase<dim,spacedim>::get_quadrature_points () const
@@ -4162,7 +4491,7 @@ FEValuesBase<dim,spacedim>::get_quadrature_points () const
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const std::vector<double> &
 FEValuesBase<dim,spacedim>::get_JxW_values () const
@@ -4173,7 +4502,7 @@ FEValuesBase<dim,spacedim>::get_JxW_values () const
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const std::vector<Tensor<2,spacedim> >&
 FEValuesBase<dim,spacedim>::get_jacobians () const
@@ -4212,34 +4541,34 @@ FEValuesBase<dim,spacedim>::quadrature_point (const unsigned int i) const
 {
   Assert (this->update_flags & update_quadrature_points, ExcAccessToUninitializedField());
   Assert (i<this->quadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size()));
-  
+
   return this->quadrature_points[i];
 }
 
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 double
 FEValuesBase<dim,spacedim>::JxW (const unsigned int i) const
 {
   Assert (this->update_flags & update_JxW_values, ExcAccessToUninitializedField());
   Assert (i<this->JxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size()));
-  
+
   return this->JxW_values[i];
 }
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Tensor<2,spacedim> &
 FEValuesBase<dim,spacedim>::jacobian (const unsigned int i) const
 {
   Assert (this->update_flags & update_jacobians, ExcAccessToUninitializedField());
   Assert (i<this->jacobians.size(), ExcIndexRange(i, 0, this->jacobians.size()));
-  
+
   return this->jacobians[i];
 }
 
@@ -4252,7 +4581,7 @@ FEValuesBase<dim,spacedim>::jacobian_grad (const unsigned int i) const
 {
   Assert (this->update_flags & update_jacobian_grads, ExcAccessToUninitializedField());
   Assert (i<this->jacobian_grads.size(), ExcIndexRange(i, 0, this->jacobian_grads.size()));
-  
+
   return this->jacobian_grads[i];
 }
 
@@ -4265,7 +4594,7 @@ FEValuesBase<dim,spacedim>::inverse_jacobian (const unsigned int i) const
 {
   Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField());
   Assert (i<this->inverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size()));
-  
+
   return this->inverse_jacobians[i];
 }
 
@@ -4283,7 +4612,7 @@ FEValuesBase<dim,spacedim>::get_function_grads (const InputVector           &fe_
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 template <class InputVector>
 inline
 void
@@ -4297,7 +4626,7 @@ FEValuesBase<dim,spacedim>::get_function_grads (
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 template <class InputVector>
 inline
 void
@@ -4310,7 +4639,7 @@ get_function_grads (const InputVector                         &fe_function,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 template <class InputVector>
 inline
 void
@@ -4325,7 +4654,7 @@ FEValuesBase<dim,spacedim>::get_function_grads (
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 template <class InputVector>
 inline
 void
@@ -4338,7 +4667,7 @@ get_function_2nd_derivatives (const InputVector           &fe_function,
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 template <class InputVector>
 inline
 void
@@ -4352,7 +4681,7 @@ get_function_2nd_derivatives (const InputVector                         &fe_func
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Point<spacedim> &
 FEValuesBase<dim,spacedim>::cell_normal_vector (const unsigned int i) const
@@ -4362,7 +4691,7 @@ FEValuesBase<dim,spacedim>::cell_normal_vector (const unsigned int i) const
          ExcIndexRange(i, 0, this->cell_normal_vectors.size()));
   Assert (this->update_flags & update_cell_normal_vectors,
          typename FVB::ExcAccessToUninitializedField());
-  
+
   return this->cell_normal_vectors[i];
 }
 
@@ -4372,20 +4701,20 @@ FEValuesBase<dim,spacedim>::cell_normal_vector (const unsigned int i) const
 /*------------------------ Inline functions: FEValues ----------------------------*/
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Quadrature<dim> &
-FEValues<dim,spacedim>::get_quadrature () const 
+FEValues<dim,spacedim>::get_quadrature () const
 {
   return quadrature;
 }
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const FEValues<dim,spacedim> &
-FEValues<dim,spacedim>::get_present_fe_values () const 
+FEValues<dim,spacedim>::get_present_fe_values () const
 {
   return *this;
 }
@@ -4394,7 +4723,7 @@ FEValues<dim,spacedim>::get_present_fe_values () const
 /*------------------------ Inline functions: FEFaceValuesBase --------------------*/
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Point<dim> &
 FEFaceValuesBase<dim,spacedim>::normal_vector (const unsigned int i) const
@@ -4404,12 +4733,12 @@ FEFaceValuesBase<dim,spacedim>::normal_vector (const unsigned int i) const
          ExcIndexRange(i, 0, this->normal_vectors.size()));
   Assert (this->update_flags & update_normal_vectors,
          typename FVB::ExcAccessToUninitializedField());
-  
+
   return this->normal_vectors[i];
 }
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 unsigned int
 FEFaceValuesBase<dim,spacedim>::get_face_index () const
@@ -4420,37 +4749,37 @@ FEFaceValuesBase<dim,spacedim>::get_face_index () const
 
 /*------------------------ Inline functions: FE*FaceValues --------------------*/
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Quadrature<dim-1> &
-FEFaceValuesBase<dim,spacedim>::get_quadrature () const 
+FEFaceValuesBase<dim,spacedim>::get_quadrature () const
 {
   return quadrature;
 }
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const FEFaceValues<dim,spacedim> &
-FEFaceValues<dim,spacedim>::get_present_fe_values () const 
+FEFaceValues<dim,spacedim>::get_present_fe_values () const
 {
   return *this;
 }
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const FESubfaceValues<dim,spacedim> &
-FESubfaceValues<dim,spacedim>::get_present_fe_values () const 
+FESubfaceValues<dim,spacedim>::get_present_fe_values () const
 {
   return *this;
 }
 
 
 
-template <int dim, int spacedim>  
+template <int dim, int spacedim>
 inline
 const Tensor<1,spacedim> &
 FEFaceValuesBase<dim,spacedim>::boundary_form (const unsigned int i) const
@@ -4460,7 +4789,7 @@ FEFaceValuesBase<dim,spacedim>::boundary_form (const unsigned int i) const
          ExcIndexRange(i, 0, this->boundary_forms.size()));
   Assert (this->update_flags & update_boundary_forms,
          typename FVB::ExcAccessToUninitializedField());
-  
+
   return this->boundary_forms[i];
 }
 
index 333157f2a36671de89ea43746c734b83884f8c9a..8a3a581771cced4c1c3f3e765e6e74422a30257e 100644 (file)
@@ -55,7 +55,7 @@ namespace
 
 namespace FEValuesViews
 {
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   Scalar<dim,spacedim>::Scalar (const FEValuesBase<dim,spacedim> &fe_values,
                                const unsigned int                component)
                  :
@@ -68,7 +68,7 @@ namespace FEValuesViews
 
     const std::vector<unsigned int> shape_function_to_row_table
       = make_shape_function_to_row_table (*fe_values.fe);
-    
+
     for (unsigned int i=0; i<fe_values.fe->dofs_per_cell; ++i)
       {
        const bool is_primitive = (fe_values.fe->is_primitive() ||
@@ -119,10 +119,10 @@ namespace FEValuesViews
     Assert (false, ExcInternalError());
     return *this;
   }
-  
 
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   Vector<dim,spacedim>::Vector (const FEValuesBase<dim,spacedim> &fe_values,
                                const unsigned int       first_vector_component)
                  :
@@ -136,11 +136,11 @@ namespace FEValuesViews
 
     const std::vector<unsigned int> shape_function_to_row_table
       = make_shape_function_to_row_table (*fe_values.fe);
-    
+
     for (unsigned int d=0; d<dim; ++d)
       {
        const unsigned int component = first_vector_component + d;
-       
+
        for (unsigned int i=0; i<fe_values.fe->dofs_per_cell; ++i)
          {
            const bool is_primitive = (fe_values.fe->is_primitive() ||
@@ -203,7 +203,7 @@ namespace FEValuesViews
          }
       }
   }
-  
+
 
   template <int dim, int spacedim>
   Vector<dim,spacedim>::Vector ()
@@ -224,17 +224,140 @@ namespace FEValuesViews
   }
 
 
+  template <int dim, int spacedim>
+  SymmetricTensor<2, dim, spacedim>::
+  SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values,
+                 const unsigned int first_tensor_component)
+                 :
+                 fe_values(fe_values),
+                 first_tensor_component(first_tensor_component),
+                 shape_function_data(fe_values.fe->dofs_per_cell)
+  {
+
+    Assert(first_tensor_component + n_independent_components - 1 < fe_values.fe->n_components(),
+          ExcIndexRange(first_tensor_component + n_independent_components - 1, 0,
+                        fe_values.fe->n_components()));
+
+    const std::vector<unsigned int> shape_function_to_row_table
+      = make_shape_function_to_row_table(*fe_values.fe);
+
+    for (unsigned int d = 0; d < n_independent_components; ++d) {
+      const unsigned int component = first_tensor_component + d;
+
+      for (unsigned int i = 0; i < fe_values.fe->dofs_per_cell; ++i) {
+       const bool is_primitive = (fe_values.fe->is_primitive() ||
+                                  fe_values.fe->is_primitive(i));
+
+       if (is_primitive == true)
+         shape_function_data[i].is_nonzero_shape_function_component[d]
+           = (component ==
+              fe_values.fe->system_to_component_index(i).first);
+       else
+         shape_function_data[i].is_nonzero_shape_function_component[d]
+           = (fe_values.fe->get_nonzero_components(i)[component]
+              == true);
+
+       if (shape_function_data[i].is_nonzero_shape_function_component[d]
+           == true) {
+         if (is_primitive == true)
+           shape_function_data[i].row_index[d]
+             = shape_function_to_row_table[i];
+         else
+           shape_function_data[i].row_index[d]
+             = (shape_function_to_row_table[i]
+                +
+                std::count(fe_values.fe->get_nonzero_components(i).begin(),
+                            fe_values.fe->get_nonzero_components(i).begin() +
+                            component,
+                            true));
+       } else
+         shape_function_data[i].row_index[d]
+           = numbers::invalid_unsigned_int;
+      }
+    }
+
+    for (unsigned int i = 0; i < fe_values.fe->dofs_per_cell; ++i) {
+      unsigned int n_nonzero_components = 0;
+      for (unsigned int d = 0; d < n_independent_components; ++d)
+       if (shape_function_data[i].is_nonzero_shape_function_component[d]
+           == true)
+         ++n_nonzero_components;
+
+      if (n_nonzero_components == 0)
+       shape_function_data[i].single_nonzero_component = -2;
+      else if (n_nonzero_components > 1)
+       shape_function_data[i].single_nonzero_component = -1;
+      else {
+       for (unsigned int d = 0; d < n_independent_components; ++d)
+         if (shape_function_data[i].is_nonzero_shape_function_component[d]
+             == true) {
+           shape_function_data[i].single_nonzero_component
+             = shape_function_data[i].row_index[d];
+           shape_function_data[i].single_nonzero_component_index
+             = d;
+           break;
+         }
+      }
+    }
+
+    switch(dim) {
+      case(1):
+           vector_to_symmetric_tensor_data[0][0] = 0;
+           break;
+      case(2):
+           vector_to_symmetric_tensor_data[0][0] = 0;
+           vector_to_symmetric_tensor_data[1][1] = 1;
+           vector_to_symmetric_tensor_data[1][2] = 2;
+           vector_to_symmetric_tensor_data[2][1] = 2;
+           break;
+      case(3):
+           vector_to_symmetric_tensor_data[0][0] = 0;
+           vector_to_symmetric_tensor_data[1][1] = 1;
+           vector_to_symmetric_tensor_data[2][2] = 2;
+           vector_to_symmetric_tensor_data[0][1] = 3;
+           vector_to_symmetric_tensor_data[1][0] = 3;
+           vector_to_symmetric_tensor_data[0][2] = 4;
+           vector_to_symmetric_tensor_data[2][0] = 4;
+           vector_to_symmetric_tensor_data[1][2] = 5;
+           vector_to_symmetric_tensor_data[2][1] = 5;
+           break;
+      default:
+           ;
+    }
+
+
+
+  }
+
+  template <int dim, int spacedim>
+  SymmetricTensor<2, dim, spacedim>::SymmetricTensor()
+                 :
+                 fe_values(*static_cast<dealii::FEValuesBase<dim, spacedim>*> (0)),
+                 first_tensor_component(numbers::invalid_unsigned_int)
+  {}
+
+
+  template <int dim, int spacedim>
+  SymmetricTensor<2, dim, spacedim> &
+  SymmetricTensor<2, dim, spacedim>::operator=(const SymmetricTensor<2, dim, spacedim> &)
+  {
+                                    // we shouldn't be copying these objects
+    Assert(false, ExcInternalError());
+    return *this;
+  }
+
+
 
   template <int dim, int spacedim>
   template <class InputVector>
   void
   Scalar<dim,spacedim>::
-  get_function_values (const InputVector &fe_function, 
+  get_function_values (const InputVector &fe_function,
                       std::vector<value_type> &values) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_values,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (values.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -258,25 +381,25 @@ namespace FEValuesViews
          if (value == 0.)
            continue;
 
-         const double * shape_value_ptr = 
+         const double * shape_value_ptr =
            &fe_values.shape_values(shape_function_data[shape_function].row_index, 0);
          for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
            values[q_point] += value * *shape_value_ptr++;
        }
   }
-  
 
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   template <class InputVector>
   void
   Scalar<dim,spacedim>::
-  get_function_gradients (const InputVector &fe_function, 
+  get_function_gradients (const InputVector &fe_function,
                          std::vector<gradient_type> &gradients) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_gradients,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (gradients.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(gradients.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -291,7 +414,7 @@ namespace FEValuesViews
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
 
     std::fill (gradients.begin(), gradients.end(), gradient_type());
-  
+
     for (unsigned int shape_function=0;
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       if (shape_function_data[shape_function].is_nonzero_shape_function_component)
@@ -300,7 +423,7 @@ namespace FEValuesViews
          if (value == 0.)
            continue;
 
-         const Tensor<1,spacedim> * shape_gradient_ptr = 
+         const Tensor<1,spacedim> * shape_gradient_ptr =
            &fe_values.shape_gradients[shape_function_data[shape_function].
                                       row_index][0];
          for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
@@ -310,11 +433,11 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   template <class InputVector>
   void
   Scalar<dim,spacedim>::
-  get_function_hessians (const InputVector &fe_function, 
+  get_function_hessians (const InputVector &fe_function,
                         std::vector<hessian_type> &hessians) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
@@ -327,7 +450,7 @@ namespace FEValuesViews
     Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
            ExcDimensionMismatch(fe_function.size(),
                                 fe_values.present_cell->n_dofs_for_dof_handler()));
-    
+
                                     // get function values of dofs
                                     // on this cell
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
@@ -343,7 +466,7 @@ namespace FEValuesViews
          if (value == 0.)
            continue;
 
-         const Tensor<2,spacedim> * shape_hessian_ptr = 
+         const Tensor<2,spacedim> * shape_hessian_ptr =
            &fe_values.shape_hessians[shape_function_data[shape_function].
                                      row_index][0];
          for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
@@ -357,12 +480,12 @@ namespace FEValuesViews
   template <class InputVector>
   void
   Scalar<dim,spacedim>::
-  get_function_laplacians (const InputVector &fe_function, 
+  get_function_laplacians (const InputVector &fe_function,
                           std::vector<value_type> &laplacians) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_hessians,
-           typename FVB::ExcAccessToUninitializedField()); 
+           typename FVB::ExcAccessToUninitializedField());
     Assert (laplacians.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(laplacians.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -388,7 +511,7 @@ namespace FEValuesViews
 
          const unsigned int row_index = shape_function_data[shape_function].row_index;
          for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           laplacians[q_point] += 
+           laplacians[q_point] +=
              value * trace(fe_values.shape_hessians[row_index][q_point]);
        }
   }
@@ -399,12 +522,12 @@ namespace FEValuesViews
   template <class InputVector>
   void
   Vector<dim,spacedim>::
-  get_function_values (const InputVector &fe_function, 
+  get_function_values (const InputVector &fe_function,
                       std::vector<value_type> &values) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_values,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (values.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -424,7 +547,7 @@ namespace FEValuesViews
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       {
        const int snc = shape_function_data[shape_function].single_nonzero_component;
-       
+
        if (snc == -2)
                                           // shape function is zero for the
                                           // selected components
@@ -453,20 +576,20 @@ namespace FEValuesViews
              }
       }
   }
-  
 
 
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   template <class InputVector>
   void
   Vector<dim,spacedim>::
-  get_function_gradients (const InputVector &fe_function, 
+  get_function_gradients (const InputVector &fe_function,
                          std::vector<gradient_type> &gradients) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_gradients,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (gradients.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(gradients.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -481,12 +604,12 @@ namespace FEValuesViews
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
 
     std::fill (gradients.begin(), gradients.end(), gradient_type());
-  
+
     for (unsigned int shape_function=0;
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       {
        const int snc = shape_function_data[shape_function].single_nonzero_component;
-       
+
        if (snc == -2)
                                           // shape function is zero for the
                                           // selected components
@@ -498,18 +621,18 @@ namespace FEValuesViews
 
        if (snc != -1)
          {
-           const unsigned int comp = 
+           const unsigned int comp =
              shape_function_data[shape_function].single_nonzero_component_index;
-           const Tensor<1,spacedim> * shape_gradient_ptr = 
+           const Tensor<1,spacedim> * shape_gradient_ptr =
              &fe_values.shape_gradients[snc][0];
            for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
              gradients[q_point][comp] += value * *shape_gradient_ptr++;
          }
        else
-         for (unsigned int d=0; d<dim; ++d)    
+         for (unsigned int d=0; d<dim; ++d)
            if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
              {
-               const Tensor<1,spacedim> * shape_gradient_ptr = 
+               const Tensor<1,spacedim> * shape_gradient_ptr =
                  &fe_values.shape_gradients[shape_function_data[shape_function].
                                             row_index[d]][0];
                for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
@@ -520,16 +643,16 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   template <class InputVector>
   void
   Vector<dim,spacedim>::
-  get_function_symmetric_gradients (const InputVector &fe_function, 
+  get_function_symmetric_gradients (const InputVector &fe_function,
                                    std::vector<symmetric_gradient_type> &symmetric_gradients) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_gradients,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (symmetric_gradients.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(symmetric_gradients.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -543,14 +666,14 @@ namespace FEValuesViews
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
 
-    std::fill (symmetric_gradients.begin(), symmetric_gradients.end(), 
+    std::fill (symmetric_gradients.begin(), symmetric_gradients.end(),
               symmetric_gradient_type());
-  
+
     for (unsigned int shape_function=0;
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       {
        const int snc = shape_function_data[shape_function].single_nonzero_component;
-       
+
        if (snc == -2)
                                           // shape function is zero for the
                                           // selected components
@@ -562,9 +685,9 @@ namespace FEValuesViews
 
        if (snc != -1)
          {
-           const unsigned int comp = 
+           const unsigned int comp =
              shape_function_data[shape_function].single_nonzero_component_index;
-           const Tensor<1,spacedim> * shape_gradient_ptr = 
+           const Tensor<1,spacedim> * shape_gradient_ptr =
              &fe_values.shape_gradients[snc][0];
            for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
              symmetric_gradients[q_point]
@@ -574,7 +697,7 @@ namespace FEValuesViews
          for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
            {
              gradient_type grad;
-             for (unsigned int d=0; d<dim; ++d)    
+             for (unsigned int d=0; d<dim; ++d)
                if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
                  grad[d] =  value *
                             fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point];
@@ -585,16 +708,16 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>  
+  template <int dim, int spacedim>
   template <class InputVector>
   void
   Vector<dim,spacedim>::
-  get_function_divergences (const InputVector &fe_function, 
+  get_function_divergences (const InputVector &fe_function,
                            std::vector<divergence_type> &divergences) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_gradients,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (divergences.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(divergences.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -609,12 +732,12 @@ namespace FEValuesViews
     fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
 
     std::fill (divergences.begin(), divergences.end(), divergence_type());
-  
+
     for (unsigned int shape_function=0;
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       {
        const int snc = shape_function_data[shape_function].single_nonzero_component;
-       
+
        if (snc == -2)
                                           // shape function is zero for the
                                           // selected components
@@ -626,7 +749,7 @@ namespace FEValuesViews
 
        if (snc != -1)
          {
-           const unsigned int comp = 
+           const unsigned int comp =
              shape_function_data[shape_function].single_nonzero_component_index;
            const Tensor<1,spacedim> * shape_gradient_ptr =
              &fe_values.shape_gradients[snc][0];
@@ -634,7 +757,7 @@ namespace FEValuesViews
              divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
          }
        else
-         for (unsigned int d=0; d<dim; ++d)    
+         for (unsigned int d=0; d<dim; ++d)
            if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
              {
                const Tensor<1,spacedim> * shape_gradient_ptr =
@@ -645,14 +768,14 @@ namespace FEValuesViews
              }
       }
   }
-  
 
 
-  template <int dim, int spacedim>  
+
+  template <int dim, int spacedim>
   template <class InputVector>
   void
   Vector<dim,spacedim>::
-  get_function_hessians (const InputVector &fe_function, 
+  get_function_hessians (const InputVector &fe_function,
                         std::vector<hessian_type> &hessians) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
@@ -665,7 +788,7 @@ namespace FEValuesViews
     Assert (fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
            ExcDimensionMismatch(fe_function.size(),
                                 fe_values.present_cell->n_dofs_for_dof_handler()));
-    
+
                                     // get function values of dofs
                                     // on this cell
     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
@@ -677,7 +800,7 @@ namespace FEValuesViews
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       {
        const int snc = shape_function_data[shape_function].single_nonzero_component;
-       
+
        if (snc == -2)
                                           // shape function is zero for the
                                           // selected components
@@ -689,9 +812,9 @@ namespace FEValuesViews
 
        if (snc != -1)
          {
-           const unsigned int comp = 
+           const unsigned int comp =
              shape_function_data[shape_function].single_nonzero_component_index;
-           const Tensor<2,spacedim> * shape_hessian_ptr = 
+           const Tensor<2,spacedim> * shape_hessian_ptr =
              &fe_values.shape_hessians[snc][0];
            for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
              hessians[q_point][comp] += value * *shape_hessian_ptr++;
@@ -700,7 +823,7 @@ namespace FEValuesViews
          for (unsigned int d=0; d<dim; ++d)
            if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
              {
-               const Tensor<2,spacedim> * shape_hessian_ptr = 
+               const Tensor<2,spacedim> * shape_hessian_ptr =
                  &fe_values.shape_hessians[shape_function_data[shape_function].
                                            row_index[d]][0];
                for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
@@ -715,12 +838,12 @@ namespace FEValuesViews
   template <class InputVector>
   void
   Vector<dim,spacedim>::
-  get_function_laplacians (const InputVector &fe_function, 
+  get_function_laplacians (const InputVector &fe_function,
                           std::vector<value_type> &laplacians) const
   {
     typedef FEValuesBase<dim,spacedim> FVB;
     Assert (fe_values.update_flags & update_hessians,
-           typename FVB::ExcAccessToUninitializedField());    
+           typename FVB::ExcAccessToUninitializedField());
     Assert (laplacians.size() == fe_values.n_quadrature_points,
            ExcDimensionMismatch(laplacians.size(), fe_values.n_quadrature_points));
     Assert (fe_values.present_cell.get() != 0,
@@ -740,7 +863,7 @@ namespace FEValuesViews
         shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
       {
        const int snc = shape_function_data[shape_function].single_nonzero_component;
-       
+
        if (snc == -2)
                                           // shape function is zero for the
                                           // selected components
@@ -754,7 +877,7 @@ namespace FEValuesViews
          {
            const unsigned int comp =
              shape_function_data[shape_function].single_nonzero_component_index;
-           const Tensor<2,spacedim> * shape_hessian_ptr = 
+           const Tensor<2,spacedim> * shape_hessian_ptr =
              &fe_values.shape_hessians[snc][0];
            for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
              laplacians[q_point][comp] += value * trace(*shape_hessian_ptr++);
@@ -763,7 +886,7 @@ namespace FEValuesViews
          for (unsigned int d=0; d<dim; ++d)
            if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
              {
-               const Tensor<2,spacedim> * shape_hessian_ptr = 
+               const Tensor<2,spacedim> * shape_hessian_ptr =
                  &fe_values.shape_hessians[shape_function_data[shape_function].
                                            row_index[d]][0];
                for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
@@ -771,6 +894,153 @@ namespace FEValuesViews
              }
       }
   }
+
+
+
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  SymmetricTensor<2, dim, spacedim>::
+  get_function_values(const InputVector &fe_function,
+                     std::vector<value_type> &values) const
+  {
+    typedef FEValuesBase<dim, spacedim> FVB;
+    Assert(fe_values.update_flags & update_values,
+          typename FVB::ExcAccessToUninitializedField());
+    Assert(values.size() == fe_values.n_quadrature_points,
+          ExcDimensionMismatch(values.size(), fe_values.n_quadrature_points));
+    Assert(fe_values.present_cell.get() != 0,
+          ExcMessage("FEValues object is not reinit'ed to any cell"));
+    Assert(fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+          ExcDimensionMismatch(fe_function.size(),
+                               fe_values.present_cell->n_dofs_for_dof_handler()));
+
+                                    // get function values of dofs
+                                    // on this cell
+    dealii::Vector<typename InputVector::value_type > dof_values(fe_values.dofs_per_cell);
+    fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+    std::fill(values.begin(), values.end(), value_type());
+
+                                    // the unique components of the
+                                    // second order tensor stored as
+                                    // a vector (i.e. a first-order
+                                    // tensor)
+    typedef Tensor<1, n_independent_components> base_tensor_type;
+
+    std::vector< base_tensor_type > values_in_vector_form(values.size(), base_tensor_type());
+
+    for (unsigned int shape_function = 0;
+        shape_function < fe_values.fe->dofs_per_cell; ++shape_function) {
+      const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+      if (snc == -2)
+                                        // shape function is zero for the
+                                        // selected components
+       continue;
+
+      const double value = dof_values(shape_function);
+      if (value == 0.)
+       continue;
+
+      if (snc != -1) {
+       const unsigned int comp =
+         shape_function_data[shape_function].single_nonzero_component_index;
+       const double * shape_value_ptr = &fe_values.shape_values(snc, 0);
+       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+         values_in_vector_form[q_point][comp] += value * *shape_value_ptr++;
+      } else
+       for (unsigned int d = 0; d < n_independent_components; ++d)
+         if (shape_function_data[shape_function].is_nonzero_shape_function_component[d]) {
+           const double * shape_value_ptr =
+             &fe_values.shape_values(shape_function_data[shape_function].row_index[d], 0);
+           for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+             values_in_vector_form[q_point][d] += value * *shape_value_ptr++;
+         }
+    }
+                                    // copy entries in std::vector to an array as there is no constructor
+                                    // for a second order tensor that take a std::vector
+    double values_array[n_independent_components];
+    for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+      for (unsigned int d = 0; d < n_independent_components; d++)
+       values_array[d] = values_in_vector_form[q_point][d];
+
+      values[q_point] = dealii::SymmetricTensor<2, dim>(values_array);
+    }
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  SymmetricTensor<2, dim, spacedim>::
+  get_function_divergences(const InputVector &fe_function,
+                          std::vector<divergence_type> &divergences) const
+  {
+    typedef FEValuesBase<dim, spacedim> FVB;
+    Assert(fe_values.update_flags & update_gradients,
+          typename FVB::ExcAccessToUninitializedField());
+    Assert(divergences.size() == fe_values.n_quadrature_points,
+          ExcDimensionMismatch(divergences.size(), fe_values.n_quadrature_points));
+    Assert(fe_values.present_cell.get() != 0,
+          ExcMessage("FEValues object is not reinit'ed to any cell"));
+    Assert(fe_function.size() == fe_values.present_cell->n_dofs_for_dof_handler(),
+          ExcDimensionMismatch(fe_function.size(),
+                               fe_values.present_cell->n_dofs_for_dof_handler()));
+
+                                    // get function values of dofs
+                                    // on this cell
+    dealii::Vector<typename InputVector::value_type > dof_values(fe_values.dofs_per_cell);
+    fe_values.present_cell->get_interpolated_dof_values(fe_function, dof_values);
+
+    std::fill(divergences.begin(), divergences.end(), divergence_type());
+
+    for (unsigned int shape_function = 0;
+        shape_function < fe_values.fe->dofs_per_cell; ++shape_function) {
+      const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+      if (snc == -2)
+                                        // shape function is zero for the
+                                        // selected components
+       continue;
+
+      const double value = dof_values(shape_function);
+      if (value == 0.)
+       continue;
+
+      if (snc != -1) {
+       const unsigned int comp =
+         shape_function_data[shape_function].single_nonzero_component_index;
+
+       const Tensor < 1, spacedim> * shape_gradient_ptr =
+         &fe_values.shape_gradients[snc][0];
+       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+            ++q_point, ++shape_gradient_ptr) {
+         for (unsigned int j = 0; j < dim; ++j) {
+           const unsigned int vector_component = vector_to_symmetric_tensor_data[comp][j];
+           divergences[q_point][vector_component] += value * (*shape_gradient_ptr)[j];
+         }
+       }
+      } else
+       for (unsigned int d = 0; d < n_independent_components; ++d)
+         if (shape_function_data[shape_function].is_nonzero_shape_function_component[d]) {
+           const unsigned int comp =
+             shape_function_data[shape_function].single_nonzero_component_index;
+
+           const Tensor < 1, spacedim> * shape_gradient_ptr =
+             &fe_values.shape_gradients[shape_function_data[shape_function].
+                                        row_index[d]][0];
+           for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+                ++q_point, ++shape_gradient_ptr) {
+             for (unsigned int j = 0; j < dim; ++j) {
+               const unsigned int vector_component = vector_to_symmetric_tensor_data[comp][j];
+               divergences[q_point][vector_component] += value * (*shape_gradient_ptr++)[j];
+             }
+           }
+         }
+    }
+  }
 }
 
 
@@ -801,7 +1071,7 @@ namespace internal
            dealii::FEValuesViews::Scalar<dim,spacedim>(fe_values,
                                                        component);
        }
-      
+
       const unsigned int n_vectors = (fe.n_components() >= dim ?
                                      fe.n_components()-dim+1 :
                                      0);
@@ -813,7 +1083,20 @@ namespace internal
          new (&vectors[component])
            dealii::FEValuesViews::Vector<dim,spacedim>(fe_values,
                                                        component);
-       }       
+            }
+
+            const unsigned int n_symmetric_second_order_tensors = (fe.n_components() >= 0.5*(dim*dim + dim) ?
+                    fe.n_components() - 0.5*(dim*dim + dim) + 1 :
+                    0);
+            symmetric_second_order_tensors.resize(n_symmetric_second_order_tensors);
+            for (unsigned int component = 0; component < n_symmetric_second_order_tensors; ++component)
+             {
+               typedef dealii::FEValuesViews::SymmetricTensor<2, dim, spacedim> SymmetricTensorView;
+                symmetric_second_order_tensors[component].SymmetricTensorView::~SymmetricTensorView();
+                new (&symmetric_second_order_tensors[component])
+                 dealii::FEValuesViews::SymmetricTensor<2, dim, spacedim > (fe_values,
+                                                                            component);
+            }
     }
   }
 }
@@ -835,7 +1118,7 @@ FEValuesBase<dim,spacedim>::CellIteratorBase::~CellIteratorBase ()
                                  * interface. See there for a
                                  * description of the use of
                                  * these classes.
-                                 * 
+                                 *
                                  * @author Wolfgang Bangerth, 2003
                                  */
 template <int dim, int spacedim>
@@ -849,7 +1132,7 @@ class FEValuesBase<dim,spacedim>::CellIterator : public FEValuesBase<dim,spacedi
                                      * this class.
                                      */
     CellIterator (const CI &cell);
-        
+
                                     /**
                                      * Conversion operator to an
                                      * iterator for
@@ -865,7 +1148,7 @@ class FEValuesBase<dim,spacedim>::CellIterator : public FEValuesBase<dim,spacedi
                                      */
     virtual
     operator const typename Triangulation<dim,spacedim>::cell_iterator () const;
-        
+
                                     /**
                                      * Return the number of
                                      * degrees of freedom the DoF
@@ -1020,7 +1303,7 @@ class FEValuesBase<dim,spacedim>::CellIterator : public FEValuesBase<dim,spacedi
                                      */
     const CI cell;
 };
-    
+
 
                                 /**
                                  * Implementation of a derived
@@ -1062,7 +1345,7 @@ class FEValuesBase<dim,spacedim>::CellIterator : public FEValuesBase<dim,spacedi
                                  * iterator that allows to
                                  * extract degree of freedom
                                  * information.
-                                 * 
+                                 *
                                  * @author Wolfgang Bangerth, 2003
                                  */
 template <int dim, int spacedim>
@@ -1075,7 +1358,7 @@ class FEValuesBase<dim,spacedim>::TriaCellIterator : public FEValuesBase<dim,spa
                                      * this class.
                                      */
     TriaCellIterator (const typename Triangulation<dim,spacedim>::cell_iterator &cell);
-        
+
                                     /**
                                      * Conversion operator to an
                                      * iterator for
@@ -1263,7 +1546,7 @@ class FEValuesBase<dim,spacedim>::TriaCellIterator : public FEValuesBase<dim,spa
                                      */
     static const char * const message_string;
 };
-    
+
 
 
 
@@ -1644,7 +1927,7 @@ FEValuesData<dim,spacedim>::initialize (const unsigned int        n_quadrature_p
                                   // functions
   unsigned int n_nonzero_shape_components = 0;
   for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-    n_nonzero_shape_components += fe.n_nonzero_components (i);  
+    n_nonzero_shape_components += fe.n_nonzero_components (i);
   Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
          ExcInternalError());
 
@@ -1663,7 +1946,7 @@ FEValuesData<dim,spacedim>::initialize (const unsigned int        n_quadrature_p
   if (flags & update_hessians)
     this->shape_hessians.resize (n_nonzero_shape_components,
                                 std::vector<Tensor<2,spacedim> > (n_quadrature_points));
-  
+
   if (flags & update_quadrature_points)
     this->quadrature_points.resize(n_quadrature_points);
 
@@ -1818,10 +2101,10 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // each quadrature point
   Assert (values.size() == n_quadrature_points,
          ExcDimensionMismatch(values.size(), n_quadrature_points));
-  
+
                                   // initialize with zero
   std::fill_n (values.begin(), n_quadrature_points, 0);
-  
+
                                   // add up contributions of trial
                                   // functions. note that here we
                                   // deal with scalar finite
@@ -1864,7 +2147,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   const InputVector&            fe_function,
   std::vector<Vector<number> >& values) const
 {
-//TODO: Find out how to do this assertion.  
+//TODO: Find out how to do this assertion.
                                   // This vector must correspond to a
                                   // complete discretization
 //  Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
@@ -1883,16 +2166,16 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   for (unsigned i=0;i<values.size();++i)
     Assert (values[i].size() == n_components,
            ExcDimensionMismatch(values[i].size(), n_components));
-  
+
   Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
   Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(),
          ExcDimensionMismatch(fe_function.size(), present_cell->n_dofs_for_dof_handler()));
-    
+
                                   // get function values of dofs
                                   // on this cell
   Vector<typename InputVector::value_type> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  
+
                                   // initialize with zero
   for (unsigned i=0;i<values.size();++i)
     std::fill_n (values[i].begin(), values[i].size(), 0);
@@ -1939,8 +2222,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                           // function can't be used because
                                           // it doesn't return us a pointer
                                           // to the data).
-         const unsigned int 
-           row = fe->is_primitive() ? 
+         const unsigned int
+           row = fe->is_primitive() ?
            shape_func : this->shape_function_to_row_table[shape_func];
 
          const double *shape_value_ptr = &this->shape_values(row, 0);
@@ -1983,9 +2266,9 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // One value per quadrature point
   Assert (n_quadrature_points == values.size(),
          ExcDimensionMismatch(values.size(), n_quadrature_points));
-  
+
   const unsigned int n_components = fe->n_components();
-  
+
                                   // Size of indices must be a
                                   // multiple of dofs_per_cell such
                                   // that an integer number of
@@ -1999,7 +2282,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // number of components of the
                                   // finite element
   const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
-  
+
   for (unsigned i=0;i<values.size();++i)
     Assert (values[i].size() == result_components,
            ExcDimensionMismatch(values[i].size(), result_components));
@@ -2009,9 +2292,9 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // element, we need this number for
                                   // loops filling all components
   const unsigned int component_multiple = result_components / n_components;
-  
+
   Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
-    
+
                                   // initialize with zero
   for (unsigned i=0;i<values.size();++i)
     std::fill_n (values[i].begin(), values[i].size(), 0);
@@ -2031,12 +2314,12 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
        if (fe->is_primitive(shape_func))
          {
-           const unsigned int 
-             row = fe->is_primitive() ? 
+           const unsigned int
+             row = fe->is_primitive() ?
              shape_func : this->shape_function_to_row_table[shape_func];
 
            const double *shape_value_ptr = &this->shape_values(row, 0);
-           const unsigned int comp = fe->system_to_component_index(shape_func).first 
+           const unsigned int comp = fe->system_to_component_index(shape_func).first
                                      + mc * n_components;
            for (unsigned int point=0; point<n_quadrature_points; ++point)
              values[point](comp) += value * *shape_value_ptr++;
@@ -2056,7 +2339,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
              const double *shape_value_ptr = &this->shape_values(row, 0);
              const unsigned int comp = c + mc * n_components;
-           
+
              for (unsigned int point=0; point<n_quadrature_points; ++point)
                values[point](comp) += value * *shape_value_ptr++;
            }
@@ -2074,7 +2357,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
   bool quadrature_points_fastest) const
 {
   const unsigned int n_components = fe->n_components();
-  
+
                                   // Size of indices must be a
                                   // multiple of dofs_per_cell such
                                   // that an integer number of
@@ -2107,15 +2390,15 @@ void FEValuesBase<dim,spacedim>::get_function_values (
        Assert (values[i].size() == result_components,
                ExcDimensionMismatch(values[i].size(), result_components));
     }
-  
+
                                   // If the result has more
                                   // components than the finite
                                   // element, we need this number for
                                   // loops filling all components
   const unsigned int component_multiple = result_components / n_components;
-  
+
   Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
-    
+
                                   // initialize with zero
   for (unsigned i=0;i<values.size();++i)
     std::fill_n (values[i].begin(), values[i].size(), 0);
@@ -2135,12 +2418,12 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
        if (fe->is_primitive(shape_func))
          {
-           const unsigned int 
-             row = fe->is_primitive() ? 
+           const unsigned int
+             row = fe->is_primitive() ?
              shape_func : this->shape_function_to_row_table[shape_func];
 
            const double *shape_value_ptr = &this->shape_values(row, 0);
-           const unsigned int comp = fe->system_to_component_index(shape_func).first 
+           const unsigned int comp = fe->system_to_component_index(shape_func).first
                                      + mc * n_components;
 
            if (quadrature_points_fastest)
@@ -2165,7 +2448,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
 
              const double *shape_value_ptr = &this->shape_values(row, 0);
              const unsigned int comp = c + mc * n_components;
-           
+
              if (quadrature_points_fastest)
                for (unsigned int point=0; point<n_quadrature_points; ++point)
                  values[comp][point] += value * *shape_value_ptr++;
@@ -2233,7 +2516,7 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
       if (value == 0.)
        continue;
 
-      const Tensor<1,spacedim> *shape_gradient_ptr 
+      const Tensor<1,spacedim> *shape_gradient_ptr
        = &this->shape_gradients[shape_func][0];
       for (unsigned int point=0; point<n_quadrature_points; ++point)
        gradients[point] += value * *shape_gradient_ptr++;
@@ -2261,10 +2544,10 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
                                   // each quadrature point
   Assert (gradients.size() == n_quadrature_points,
          ExcDimensionMismatch(gradients.size(), n_quadrature_points));
-  
+
                                   // initialize with zero
   std::fill_n (gradients.begin(), n_quadrature_points, Tensor<1,spacedim>());
-  
+
                                   // add up contributions of trial
                                   // functions. note that here we
                                   // deal with scalar finite
@@ -2293,7 +2576,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
       if (value == 0.)
        continue;
 
-      const Tensor<1,spacedim> *shape_gradient_ptr 
+      const Tensor<1,spacedim> *shape_gradient_ptr
        = &this->shape_gradients[shape_func][0];
       for (unsigned int point=0; point<n_quadrature_points; ++point)
        gradients[point] += value * *shape_gradient_ptr++;
@@ -2348,10 +2631,10 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
 
       if (fe->is_primitive(shape_func))
        {
-         const unsigned int 
-           row = fe->is_primitive() ? 
+         const unsigned int
+           row = fe->is_primitive() ?
            shape_func : this->shape_function_to_row_table[shape_func];
-         const Tensor<1,spacedim> *shape_gradient_ptr 
+         const Tensor<1,spacedim> *shape_gradient_ptr
            = &this->shape_gradients[row][0];
          const unsigned int comp = fe->system_to_component_index(shape_func).first;
          for (unsigned int point=0; point<n_quadrature_points; ++point)
@@ -2370,7 +2653,7 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
                                 fe->get_nonzero_components(shape_func).begin()+c,
                                 true));
 
-           const Tensor<1,spacedim> *shape_gradient_ptr 
+           const Tensor<1,spacedim> *shape_gradient_ptr
              = &this->shape_gradients[row][0];
 
            for (unsigned int point=0; point<n_quadrature_points; ++point)
@@ -2390,7 +2673,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
   bool quadrature_points_fastest) const
 {
   const unsigned int n_components = fe->n_components();
-  
+
                                   // Size of indices must be a
                                   // multiple of dofs_per_cell such
                                   // that an integer number of
@@ -2404,7 +2687,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
                                   // number of components of the
                                   // finite element
   const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
-  
+
                                   // Check if the value argument is
                                   // initialized to the correct sizes
   if (quadrature_points_fastest)
@@ -2429,9 +2712,9 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
                                   // element, we need this number for
                                   // loops filling all components
   const unsigned int component_multiple = result_components / n_components;
-  
+
   Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
-    
+
                                   // initialize with zero
   for (unsigned i=0;i<gradients.size();++i)
     std::fill_n (gradients[i].begin(), gradients[i].size(), Tensor<1,spacedim>());
@@ -2451,10 +2734,10 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
 
        if (fe->is_primitive(shape_func))
          {
-           const unsigned int 
-             row = fe->is_primitive() ? 
+           const unsigned int
+             row = fe->is_primitive() ?
              shape_func : this->shape_function_to_row_table[shape_func];
-           const Tensor<1,spacedim> *shape_gradient_ptr 
+           const Tensor<1,spacedim> *shape_gradient_ptr
              = &this->shape_gradients[row][0];
            const unsigned int comp = fe->system_to_component_index(shape_func).first
                                      + mc * n_components;
@@ -2479,7 +2762,7 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
                                   fe->get_nonzero_components(shape_func).begin()+c,
                                   true));
 
-             const Tensor<1,spacedim> *shape_gradient_ptr 
+             const Tensor<1,spacedim> *shape_gradient_ptr
                = &this->shape_gradients[row][0];
              const unsigned int comp = c + mc * n_components;
 
@@ -2533,7 +2816,7 @@ get_function_hessians (const InputVector           &fe_function,
       if (value == 0.)
        continue;
 
-      const Tensor<2,spacedim> *shape_hessians_ptr 
+      const Tensor<2,spacedim> *shape_hessians_ptr
        = &this->shape_hessians[shape_func][0];
       for (unsigned int point=0; point<n_quadrature_points; ++point)
        hessians[point] += value * *shape_hessians_ptr++;
@@ -2561,10 +2844,10 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
                                   // each quadrature point
   Assert (hessians.size() == n_quadrature_points,
          ExcDimensionMismatch(hessians.size(), n_quadrature_points));
-  
+
                                   // initialize with zero
   std::fill_n (hessians.begin(), n_quadrature_points, Tensor<2,spacedim>());
-  
+
                                   // add up contributions of trial
                                   // functions. note that here we
                                   // deal with scalar finite
@@ -2577,7 +2860,7 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
       if (value == 0.)
        continue;
 
-      const Tensor<2,spacedim> *shape_hessians_ptr 
+      const Tensor<2,spacedim> *shape_hessians_ptr
        = &this->shape_hessians[shape_func][0];
       for (unsigned int point=0; point<n_quadrature_points; ++point)
        hessians[point] += value * *shape_hessians_ptr++;
@@ -2629,11 +2912,11 @@ get_function_hessians (const InputVector                         &fe_function,
 
       if (fe->is_primitive(shape_func))
        {
-         const unsigned int 
-           row = fe->is_primitive() ? 
+         const unsigned int
+           row = fe->is_primitive() ?
            shape_func : this->shape_function_to_row_table[shape_func];
 
-         const Tensor<2,spacedim> *shape_hessian_ptr 
+         const Tensor<2,spacedim> *shape_hessian_ptr
            = &this->shape_hessians[row][0];
          const unsigned int comp = fe->system_to_component_index(shape_func).first;
 
@@ -2657,7 +2940,7 @@ get_function_hessians (const InputVector                         &fe_function,
                                 fe->get_nonzero_components(shape_func).begin()+c,
                                 true));
 
-           const Tensor<2,spacedim> *shape_hessian_ptr 
+           const Tensor<2,spacedim> *shape_hessian_ptr
              = &this->shape_hessians[row][0];
 
            if (quadrature_points_fastest)
@@ -2681,9 +2964,9 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
   bool quadrature_points_fastest) const
 {
   Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField());
-    
+
   const unsigned int n_components = fe->n_components();
-  
+
                                   // Size of indices must be a
                                   // multiple of dofs_per_cell such
                                   // that an integer number of
@@ -2697,7 +2980,7 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
                                   // number of components of the
                                   // finite element
   const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
-  
+
                                   // Check if the value argument is
                                   // initialized to the correct sizes
   if (quadrature_points_fastest)
@@ -2722,7 +3005,7 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
                                   // element, we need this number for
                                   // loops filling all components
   const unsigned int component_multiple = result_components / n_components;
-  
+
                                   // initialize with zero
   for (unsigned i=0;i<hessians.size();++i)
     std::fill_n (hessians[i].begin(), hessians[i].size(), Tensor<2,spacedim>());
@@ -2742,11 +3025,11 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
 
        if (fe->is_primitive(shape_func))
          {
-           const unsigned int 
-             row = fe->is_primitive() ? 
+           const unsigned int
+             row = fe->is_primitive() ?
              shape_func : this->shape_function_to_row_table[shape_func];
 
-           const Tensor<2,spacedim> *shape_hessian_ptr 
+           const Tensor<2,spacedim> *shape_hessian_ptr
              = &this->shape_hessians[row][0];
            const unsigned int comp = fe->system_to_component_index(shape_func).first
                                      + mc * n_components;
@@ -2771,7 +3054,7 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
                                   fe->get_nonzero_components(shape_func).begin()+c,
                                   true));
 
-             const Tensor<2,spacedim> *shape_hessian_ptr 
+             const Tensor<2,spacedim> *shape_hessian_ptr
                = &this->shape_hessians[row][0];
              const unsigned int comp = c + mc * n_components;
 
@@ -2827,7 +3110,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       if (value == 0.)
        continue;
 
-      const Tensor<2,spacedim> *shape_hessian_ptr 
+      const Tensor<2,spacedim> *shape_hessian_ptr
        = &this->shape_hessians[shape_func][0];
       for (unsigned int point=0; point<n_quadrature_points; ++point)
        laplacians[point] += value * trace(*shape_hessian_ptr++);
@@ -2855,10 +3138,10 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                   // each quadrature point
   Assert (laplacians.size() == n_quadrature_points,
          ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
-  
+
                                   // initialize with zero
   std::fill_n (laplacians.begin(), n_quadrature_points, 0);
-  
+
                                   // add up contributions of trial
                                   // functions. note that here we
                                   // deal with scalar finite
@@ -2874,7 +3157,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
       if (value == 0.)
        continue;
 
-      const Tensor<2,spacedim> *shape_hessian_ptr 
+      const Tensor<2,spacedim> *shape_hessian_ptr
        = &this->shape_hessians[shape_func][0];
       for (unsigned int point=0; point<n_quadrature_points; ++point)
        laplacians[point] += value * trace(*shape_hessian_ptr++);
@@ -2889,7 +3172,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   const InputVector&            fe_function,
   std::vector<Vector<number> >& laplacians) const
 {
-//TODO: Find out how to do this assertion.  
+//TODO: Find out how to do this assertion.
                                   // This vector must correspond to a
                                   // complete discretization
 //  Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
@@ -2908,16 +3191,16 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   for (unsigned i=0;i<laplacians.size();++i)
     Assert (laplacians[i].size() == n_components,
            ExcDimensionMismatch(laplacians[i].size(), n_components));
-  
+
   Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
   Assert (fe_function.size() == present_cell->n_dofs_for_dof_handler(),
          ExcDimensionMismatch(fe_function.size(), present_cell->n_dofs_for_dof_handler()));
-    
+
                                   // get function values of dofs
                                   // on this cell
   Vector<typename InputVector::value_type> dof_values (dofs_per_cell);
   present_cell->get_interpolated_dof_values(fe_function, dof_values);
-  
+
                                   // initialize with zero
   for (unsigned i=0;i<laplacians.size();++i)
     std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
@@ -2936,11 +3219,11 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
 
       if (fe->is_primitive(shape_func))
        {
-         const unsigned int 
-           row = fe->is_primitive() ? 
+         const unsigned int
+           row = fe->is_primitive() ?
            shape_func : this->shape_function_to_row_table[shape_func];
 
-         const Tensor<2,spacedim> *shape_hessian_ptr 
+         const Tensor<2,spacedim> *shape_hessian_ptr
            = &this->shape_hessians[row][0];
          const unsigned int comp = fe->system_to_component_index(shape_func).first;
          for (unsigned int point=0; point<n_quadrature_points; ++point)
@@ -2959,7 +3242,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                 fe->get_nonzero_components(shape_func).begin()+c,
                                 true));
 
-           const Tensor<2,spacedim> *shape_hessian_ptr 
+           const Tensor<2,spacedim> *shape_hessian_ptr
              = &this->shape_hessians[row][0];
 
            for (unsigned int point=0; point<n_quadrature_points; ++point)
@@ -2980,9 +3263,9 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                   // One value per quadrature point
   Assert (n_quadrature_points == laplacians.size(),
          ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
-  
+
   const unsigned int n_components = fe->n_components();
-  
+
                                   // Size of indices must be a
                                   // multiple of dofs_per_cell such
                                   // that an integer number of
@@ -2996,7 +3279,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                   // number of components of the
                                   // finite element
   const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
-  
+
   for (unsigned i=0;i<laplacians.size();++i)
     Assert (laplacians[i].size() == result_components,
            ExcDimensionMismatch(laplacians[i].size(), result_components));
@@ -3006,9 +3289,9 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                   // element, we need this number for
                                   // loops filling all components
   const unsigned int component_multiple = result_components / n_components;
-  
+
   Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
-    
+
                                   // initialize with zero
   for (unsigned i=0;i<laplacians.size();++i)
     std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
@@ -3028,13 +3311,13 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
 
        if (fe->is_primitive(shape_func))
          {
-           const unsigned int 
-             row = fe->is_primitive() ? 
+           const unsigned int
+             row = fe->is_primitive() ?
              shape_func : this->shape_function_to_row_table[shape_func];
 
            const Tensor<2,spacedim> *shape_hessian_ptr
              = &this->shape_hessians[row][0];
-           const unsigned int comp = fe->system_to_component_index(shape_func).first 
+           const unsigned int comp = fe->system_to_component_index(shape_func).first
                                      + mc * n_components;
            for (unsigned int point=0; point<n_quadrature_points; ++point)
              laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
@@ -3052,10 +3335,10 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                   fe->get_nonzero_components(shape_func).begin()+c,
                                   true));
 
-             const Tensor<2,spacedim> *shape_hessian_ptr 
+             const Tensor<2,spacedim> *shape_hessian_ptr
                = &this->shape_hessians[row][0];
              const unsigned int comp = c + mc * n_components;
-           
+
              for (unsigned int point=0; point<n_quadrature_points; ++point)
                laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
            }
@@ -3073,7 +3356,7 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
   bool quadrature_points_fastest) const
 {
   const unsigned int n_components = fe->n_components();
-  
+
                                   // Size of indices must be a
                                   // multiple of dofs_per_cell such
                                   // that an integer number of
@@ -3106,15 +3389,15 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
        Assert (laplacians[i].size() == result_components,
                ExcDimensionMismatch(laplacians[i].size(), result_components));
     }
-  
+
                                   // If the result has more
                                   // components than the finite
                                   // element, we need this number for
                                   // loops filling all components
   const unsigned int component_multiple = result_components / n_components;
-  
+
   Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
-    
+
                                   // initialize with zero
   for (unsigned i=0;i<laplacians.size();++i)
     std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
@@ -3134,13 +3417,13 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
 
        if (fe->is_primitive(shape_func))
          {
-           const unsigned int 
-             row = fe->is_primitive() ? 
+           const unsigned int
+             row = fe->is_primitive() ?
              shape_func : this->shape_function_to_row_table[shape_func];
 
-           const Tensor<2,spacedim> *shape_hessian_ptr 
+           const Tensor<2,spacedim> *shape_hessian_ptr
              = &this->shape_hessians[row][0];
-           const unsigned int comp = fe->system_to_component_index(shape_func).first 
+           const unsigned int comp = fe->system_to_component_index(shape_func).first
                                      + mc * n_components;
            if (quadrature_points_fastest)
              for (unsigned int point=0; point<n_quadrature_points; ++point)
@@ -3162,10 +3445,10 @@ void FEValuesBase<dim,spacedim>::get_function_laplacians (
                                   fe->get_nonzero_components(shape_func).begin()+c,
                                   true));
 
-             const Tensor<2,spacedim> *shape_hessian_ptr 
+             const Tensor<2,spacedim> *shape_hessian_ptr
                = &this->shape_hessians[row][0];
              const unsigned int comp = c + mc * n_components;
-           
+
              if (quadrature_points_fastest)
                for (unsigned int point=0; point<n_quadrature_points; ++point)
                  laplacians[comp][point] += value * trace(*shape_hessian_ptr++);
@@ -3222,7 +3505,7 @@ template <int dim, int spacedim>
 UpdateFlags
 FEValuesBase<dim,spacedim>::compute_update_flags (const UpdateFlags update_flags) const
 {
-  
+
                                   // first find out which objects
                                   // need to be recomputed on each
                                   // cell we visit. this we have to
@@ -3243,7 +3526,7 @@ FEValuesBase<dim,spacedim>::compute_update_flags (const UpdateFlags update_flags
 template <int dim, int spacedim>
 inline
 void
-FEValuesBase<dim,spacedim>::check_cell_similarity 
+FEValuesBase<dim,spacedim>::check_cell_similarity
 (const typename Triangulation<dim,spacedim>::cell_iterator &cell)
 {
                                   // case that there has not been any cell
@@ -3260,11 +3543,11 @@ FEValuesBase<dim,spacedim>::check_cell_similarity
     else
       cell_similarity = (cell->is_translation_of
                         (static_cast<const typename Triangulation<dim,spacedim>::cell_iterator>(*this->present_cell))
-                        ? 
-                        CellSimilarity::translation 
-                        : 
+                        ?
+                        CellSimilarity::translation
+                        :
                         CellSimilarity::none);
-  
+
                                   // TODO: here, one could implement other
                                   // checks for similarity, e.g. for
                                   // children of a parallelogram.
@@ -3273,7 +3556,7 @@ FEValuesBase<dim,spacedim>::check_cell_similarity
 
 
 template <int dim, int spacedim>
-CellSimilarity::Similarity 
+CellSimilarity::Similarity
 FEValuesBase<dim,spacedim>::get_cell_similarity () const
 {
   return cell_similarity;
@@ -3346,7 +3629,7 @@ FEValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
 
 
   const UpdateFlags flags = this->compute_update_flags (update_flags);
-  
+
                                   // then get objects into which the
                                   // FE and the Mapping can store
                                   // intermediate data used across
@@ -3381,7 +3664,7 @@ FEValues<dim,spacedim>::reinit (const typename DoFHandler<dim,spacedim>::cell_it
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename DoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -3415,7 +3698,7 @@ FEValues<dim,spacedim>::reinit (const typename hp::DoFHandler<dim,spacedim>::cel
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename hp::DoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -3460,7 +3743,7 @@ FEValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedim>::cell_
                                   // this cell, are the same
 
 //TODO: This was documented out ith the repository. Why?
-  
+
 //  Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
 //       static_cast<const FiniteElementData<dim>&>(cell->get_fe()),
 //       typename FEValuesBase<dim,spacedim>::ExcFEDontMatch());
@@ -3472,7 +3755,7 @@ FEValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedim>::cell_
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename MGDoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -3498,7 +3781,7 @@ void FEValues<dim,spacedim>::reinit (const typename Triangulation<dim,spacedim>:
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::TriaCellIterator (cell));
                                    // this was the part of the work
                                    // that is dependent on the actual
@@ -3554,7 +3837,7 @@ template <int dim, int spacedim>
 FEFaceValuesBase<dim,spacedim>::FEFaceValuesBase (const unsigned int n_q_points,
                                                  const unsigned int dofs_per_cell,
                                                  const UpdateFlags,
-                                                 const Mapping<dim,spacedim> &mapping,      
+                                                 const Mapping<dim,spacedim> &mapping,
                                                  const FiniteElement<dim,spacedim> &fe,
                                                  const Quadrature<dim-1>& quadrature)
                :
@@ -3649,7 +3932,7 @@ void
 FEFaceValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
 {
   const UpdateFlags flags = this->compute_update_flags (update_flags);
-  
+
                                   // then get objects into which the
                                   // FE and the Mapping can store
                                   // intermediate data used across
@@ -3677,16 +3960,16 @@ void FEFaceValues<dim,spacedim>::reinit (const typename DoFHandler<dim,spacedim>
 
   Assert (face_no < GeometryInfo<dim>::faces_per_cell,
          ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-  
+
                                    // set new cell. auto_ptr will take
                                    // care that old object gets
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename DoFHandler<dim,spacedim>::cell_iterator> (cell));
-  
+
                                    // this was the part of the work
                                    // that is dependent on the actual
                                    // data type of the iterator. now
@@ -3713,13 +3996,13 @@ void FEFaceValues<dim,spacedim>::reinit (const typename hp::DoFHandler<dim,space
 
   Assert (face_no < GeometryInfo<dim>::faces_per_cell,
          ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-  
+
                                    // set new cell. auto_ptr will take
                                    // care that old object gets
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename hp::DoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -3769,13 +4052,13 @@ void FEFaceValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedi
 
   Assert (face_no < GeometryInfo<dim>::faces_per_cell,
          ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-  
+
                                    // set new cell. auto_ptr will take
                                    // care that old object gets
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename MGDoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -3800,7 +4083,7 @@ void FEFaceValues<dim,spacedim>::reinit (const typename Triangulation<dim,spaced
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::TriaCellIterator (cell));
 
                                    // this was the part of the work
@@ -3821,7 +4104,7 @@ void FEFaceValues<dim,spacedim>::do_reinit (const unsigned int face_no)
                                   // available)
   const typename Triangulation<dim,spacedim>::cell_iterator cell=*this->present_cell;
   this->present_face_index=cell->face_index(face_no);
-  
+
   this->get_mapping().fill_fe_face_values(*this->present_cell, face_no,
                                          this->quadrature,
                                          *this->mapping_data,
@@ -3892,7 +4175,7 @@ void
 FESubfaceValues<dim,spacedim>::initialize (const UpdateFlags update_flags)
 {
   const UpdateFlags flags = this->compute_update_flags (update_flags);
-  
+
                                   // then get objects into which the
                                   // FE and the Mapping can store
                                   // intermediate data used across
@@ -3947,7 +4230,7 @@ void FESubfaceValues<dim,spacedim>::reinit (const typename DoFHandler<dim,spaced
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename DoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -4001,7 +4284,7 @@ void FESubfaceValues<dim,spacedim>::reinit (const typename hp::DoFHandler<dim,sp
   do_reinit (face_no, subface_no);
 }
 
-  
+
 template <int dim, int spacedim>
 void FESubfaceValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
                                            const unsigned int         face_no,
@@ -4025,7 +4308,7 @@ void FESubfaceValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spac
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::template
      CellIterator<typename MGDoFHandler<dim,spacedim>::cell_iterator> (cell));
 
@@ -4048,13 +4331,13 @@ void FESubfaceValues<dim,spacedim>::reinit (const typename Triangulation<dim,spa
          ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
   Assert (subface_no < cell->face(face_no)->n_children(),
          ExcIndexRange (subface_no, 0, cell->face(face_no)->n_children()));
-  
+
                                    // set new cell. auto_ptr will take
                                    // care that old object gets
                                    // destroyed and also that this
                                    // object gets destroyed in the
                                    // destruction of this class
-  this->present_cell.reset 
+  this->present_cell.reset
     (new typename FEValuesBase<dim,spacedim>::TriaCellIterator (cell));
 
                                    // this was the part of the work
@@ -4139,7 +4422,7 @@ void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
               ExcInternalError());
        this->present_face_index=subface_index;
       }
-  
+
                                   // now ask the mapping and the finite element
                                   // to do the actual work
   this->get_mapping().fill_fe_subface_values(*this->present_cell,
@@ -4150,7 +4433,7 @@ void FESubfaceValues<dim,spacedim>::do_reinit (const unsigned int face_no,
                                             this->JxW_values,
                                             this->boundary_forms,
                                             this->normal_vectors);
-  
+
   this->get_fe().fill_fe_subface_values(this->get_mapping(),
                                        *this->present_cell,
                                         face_no, subface_no,
@@ -4184,8 +4467,8 @@ template class FESubfaceValues<deal_II_dimension>;
 template class FEValuesData<deal_II_dimension,deal_II_dimension+1>;
 template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>;
 template class FEValues<deal_II_dimension,deal_II_dimension+1>;
-template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>:: 
-CellIterator<DoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>; 
+template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
+CellIterator<DoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
 //template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
 //  CellIterator<MGDoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
 
index 5009baf9a52b04cbb2a237af3d3f8d2f3925e209..de06ad83ff2be83ec8aeb6d835b947c40d3cd058 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -37,7 +37,7 @@ for (VEC : SERIAL_VECTORS)
        (const dealii::VEC&, std::vector<Tensor<2,deal_II_dimension> >&) const;
     template
        void FEValuesViews::Vector<deal_II_dimension>::get_function_symmetric_gradients<dealii::VEC>
-       (const dealii::VEC&, std::vector<SymmetricTensor<2,deal_II_dimension> >&) const;
+       (const dealii::VEC&, std::vector<dealii::SymmetricTensor<2,deal_II_dimension> >&) const;
     template
        void FEValuesViews::Vector<deal_II_dimension>::get_function_divergences<dealii::VEC>
        (const dealii::VEC&, std::vector<double>&) const;
@@ -48,6 +48,13 @@ for (VEC : SERIAL_VECTORS)
        void FEValuesViews::Vector<deal_II_dimension>::get_function_laplacians<dealii::VEC>
        (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension> >&) const;
 
+    template
+       void FEValuesViews::SymmetricTensor<2,deal_II_dimension>::get_function_values<dealii::VEC>
+       (const dealii::VEC&, std::vector<dealii::SymmetricTensor<2,deal_II_dimension> >&) const;
+    template
+       void FEValuesViews::SymmetricTensor<2,deal_II_dimension>::get_function_divergences<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension> >&) const;
+
 
 #if deal_II_dimension != 3
     template
@@ -78,7 +85,7 @@ for (VEC : SERIAL_VECTORS)
     template
        void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
        ::get_function_symmetric_gradients<dealii::VEC>
-       (const dealii::VEC&, std::vector<SymmetricTensor<2,deal_II_dimension+1> >&) const;
+       (const dealii::VEC&, std::vector<dealii::SymmetricTensor<2,deal_II_dimension+1> >&) const;
     template
        void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
        ::get_function_divergences<dealii::VEC>
@@ -92,6 +99,15 @@ for (VEC : SERIAL_VECTORS)
        ::get_function_laplacians<dealii::VEC>
        (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension+1> >&) const;
 
+    template
+       void FEValuesViews::SymmetricTensor<2,deal_II_dimension+1>
+        ::get_function_values<dealii::VEC>
+       (const dealii::VEC&, std::vector<dealii::SymmetricTensor<2,deal_II_dimension+1> >&) const;
+    template
+       void FEValuesViews::SymmetricTensor<2,deal_II_dimension+1>
+        ::get_function_divergences<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension+1> >&) const;
+
 #endif
 
   }

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.