From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 11 Dec 2008 19:01:53 +0000 (+0000)
Subject: Modify the FEValuesViews classes so that they compute a bunch of things already at... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f89865812ce7166850c215d71e01f2173a4947f5;p=dealii-svn.git

Modify the FEValuesViews classes so that they compute a bunch of things already at construction time. This should make them much faster when you use them to evaluate stuff. On the other hand, it makes construction much more expensive and so the previous approach of creating them on the fly every time someone wrote fe_values[velocities] doesn't work any more. Rather, create all possible views objects at construction time of the FEValues object and simply return a reference to one of them when a view is requested.

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

diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h
index 5798daaef2..28d39c7b53 100644
--- a/deal.II/deal.II/include/fe/fe_values.h
+++ b/deal.II/deal.II/include/fe/fe_values.h
@@ -184,6 +184,12 @@ namespace FEValuesViews
 					*/
       typedef Tensor<2,spacedim> hessian_type;
 
+				       /**
+					* Default constructor. Creates an
+					* invalid object.
+					*/
+      Scalar ();
+      
 				       /**
 					* Constructor for an object that
 					* represents a single scalar component
@@ -194,6 +200,14 @@ namespace FEValuesViews
       Scalar (const FEValuesBase<dim,spacedim> &fe_values_base,
 	      const unsigned int       component);
 
+				       /**
+					* Copy operator. This is not a
+					* lightweight object so we don't allow
+					* copying and generate an exception if
+					* this function is called.
+					*/
+      Scalar & operator= (const Scalar<dim,spacedim> &);
+      
 				       /**
 					* Return the value of the vector
 					* component selected by this view, for
@@ -241,6 +255,42 @@ namespace FEValuesViews
 					* object.
 					*/
       const unsigned int component;
+
+				       /**
+					* For each shape function, store
+					* whether it is primitive.
+					*/
+      std::vector<bool> is_primitive;
+
+				       /**
+					* For each shape function, 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).
+					*/
+      Table<1,bool> is_nonzero_shape_function_component;
+
+				       /**
+					* For each shape function, 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.
+					*/
+      Table<1,unsigned int> row_index;
   };
   
 
@@ -308,6 +358,12 @@ namespace FEValuesViews
 					*/
       typedef Tensor<3,dim>          hessian_type;
 
+				       /**
+					* Default constructor. Creates an
+					* invalid object.
+					*/
+      Vector ();
+
 				       /**
 					* Constructor for an object that
 					* represents dim components of a
@@ -322,7 +378,15 @@ namespace FEValuesViews
 					*/
       Vector (const FEValuesBase<dim,spacedim> &fe_values_base,
 	      const unsigned int first_vector_component);
-	      
+
+				       /**
+					* Copy operator. This is not a
+					* lightweight object so we don't allow
+					* copying and generate an exception if
+					* this function is called.
+					*/
+      Vector & operator= (const Vector<dim,spacedim> &);
+
 				       /**
 					* Return the value of the vector
 					* components selected by this view,
@@ -407,10 +471,77 @@ namespace FEValuesViews
 					* FEValuesBase object.
 					*/
       const unsigned int first_vector_component;
+
+				       /**
+					* 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).
+					*/
+      Table<2,bool> is_nonzero_shape_function_component;
+
+				       /**
+					* 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.
+					*/
+      Table<2,unsigned int> row_index;
   };
 }
 
 
+namespace internal
+{
+  namespace FEValuesViews
+  {
+				     /**
+				      * A class objects of which store a
+				      * collection of FEValuesViews::Scalar,
+				      * FEValuesViews::Vector, etc object. The
+				      * FEValuesBase class uses it to generate
+				      * all possible Views classes upon
+				      * construction time; we do this at
+				      * construction time since the Views
+				      * classes cache some information and are
+				      * therefore relatively expensive to
+				      * create.
+				      */
+    template <int dim, int spacedim>
+    struct Cache
+    {
+					 /**
+					  * Caches for scalar and
+					  * vector-valued views.
+					  */
+	std::vector<dealii::FEValuesViews::Scalar<dim,spacedim> > scalars;
+	std::vector<dealii::FEValuesViews::Vector<dim,spacedim> > vectors;
+
+					 /**
+					  * Constructor.
+					  */
+	Cache (const FEValuesBase<dim,spacedim> &fe_values);
+    };
+  }
+}
+
+
 
   
 /*!@addtogroup feaccess */
@@ -866,7 +997,7 @@ template <int dim, int spacedim>
 				      * FEValuesViews and in particular
 				      * in the @ref vector_valued module.
 				      */
-    FEValuesViews::Scalar<dim,spacedim>
+    const FEValuesViews::Scalar<dim,spacedim> &
     operator[] (const FEValuesExtractors::Scalar &scalar) const;
 
 				     /**
@@ -879,7 +1010,7 @@ template <int dim, int spacedim>
 				      * the namespace FEValuesViews and in particular
 				      * in the @ref vector_valued module.
 				      */
-    FEValuesViews::Vector<dim,spacedim>
+    const FEValuesViews::Vector<dim,spacedim> &
     operator[] (const FEValuesExtractors::Vector &vector) const;
 
     				     //@}
@@ -2202,7 +2333,7 @@ template <int dim, int spacedim>
                                           */
         const CI cell;
     };
-
+    
 
                                      /**
                                       * Implementation of a derived
@@ -2511,6 +2642,12 @@ template <int dim, int spacedim>
                                       */
     FEValuesBase & operator= (const FEValuesBase &);
 
+				     /**
+				      * A cache for all possible FEValuesViews
+				      * objects.
+				      */
+    internal::FEValuesViews::Cache<dim,spacedim> fe_values_views_cache;
+    
 				     /**
 				      * Make the view classes friends of this
 				      * class, since they access internal
@@ -2532,7 +2669,7 @@ template <int dim, int spacedim>
  * @author Wolfgang Bangerth, 1998, Guido Kanschat, 2001
  */
 template <int dim, int spacedim=dim>
-  class FEValues : public FEValuesBase<dim,spacedim>
+class FEValues : public FEValuesBase<dim,spacedim>
 {
   public:
                                      /**
@@ -3252,20 +3389,6 @@ namespace FEValuesExtractors
 
 namespace FEValuesViews
 {
-  template <int dim, int spacedim>  
-  inline
-  Scalar<dim,spacedim>::Scalar (const FEValuesBase<dim,spacedim> &fe_values,
-		       const unsigned int       component)
-		  :
-		  fe_values (fe_values),
-		  component (component)
-  {
-    Assert (component < fe_values.fe->n_components(),
-	    ExcIndexRange(component, 0, fe_values.fe->n_components()));
-  }
-
-
-
   template <int dim, int spacedim>  
   inline
   typename Scalar<dim,spacedim>::value_type
@@ -3281,37 +3404,15 @@ namespace FEValuesViews
 				     // an adaptation of the
 				     // FEValuesBase::shape_value_component
 				     // function except that here we know the
-				     // component as fixed. see the
-				     // comments there
-				     //
-				     // we can do away with some of the
-				     // assertions since they are already
-				     // taken care of in the constructor of
-				     // this class
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	if (component ==
-	    fe_values.fe->system_to_component_index(shape_function).first)
-	  return fe_values.shape_values(fe_values.shape_function_to_row_table[shape_function],q_point);
-	else
-	  return 0;
-      }
+				     // component as fixed and we have
+				     // pre-computed and cached a bunch of
+				     // information. see the comments there
+    if (is_nonzero_shape_function_component[shape_function])
+      return fe_values.shape_values(row_index[shape_function],q_point);
     else
-      {
-	if (fe_values.fe->get_nonzero_components(shape_function)[component] == false)
-	  return 0;
-
-	const unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin()+
-			     component,
-			     true));
-	return fe_values.shape_values(row, q_point);
-      }
+      return 0;
   }
+  
 
 
 
@@ -3330,36 +3431,13 @@ namespace FEValuesViews
 				     // an adaptation of the
 				     // FEValuesBase::shape_grad_component
 				     // function except that here we know the
-				     // component as fixed. see the
-				     // comments there
-				     //
-				     // we can do away with the assertions
-				     // since they are already taken care of
-				     // in the constructor of this class
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	if (component ==
-	    fe_values.fe->system_to_component_index(shape_function).first)
-	  return fe_values.shape_gradients[fe_values.
-					    shape_function_to_row_table[shape_function]][q_point];
-	else
-	  return gradient_type();
-      }
+				     // component as fixed and we have
+				     // pre-computed and cached a bunch of
+				     // information. see the comments there
+    if (is_nonzero_shape_function_component[shape_function])
+      return fe_values.shape_gradients[row_index[shape_function]][q_point];
     else
-      {
-	if (fe_values.fe->get_nonzero_components(shape_function)[component] == false)
-	  return gradient_type();
-
-	const unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin()+
-			     component,
-			     true));
-	return fe_values.shape_gradients[row][q_point];
-      }
+      return gradient_type();
   }
 
 
@@ -3379,51 +3457,13 @@ namespace FEValuesViews
 				     // an adaptation of the
 				     // FEValuesBase::shape_grad_component
 				     // function except that here we know the
-				     // component as fixed. see the
-				     // comments there
-				     //
-				     // we can do away with the assertions
-				     // since they are already taken care of
-				     // in the constructor of this class
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	if (component ==
-	    fe_values.fe->system_to_component_index(shape_function).first)
-	  return fe_values.shape_hessians[fe_values.
-					   shape_function_to_row_table[shape_function]][q_point];
-	else
-	  return hessian_type();
-      }
+				     // component as fixed and we have
+				     // pre-computed and cached a bunch of
+				     // information. see the comments there
+    if (is_nonzero_shape_function_component[shape_function])
+      return fe_values.shape_hessians[row_index[shape_function]][q_point];
     else
-      {
-	if (fe_values.fe->get_nonzero_components(shape_function)[component] == false)
-	  return hessian_type();
-
-	const unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin()+
-			     component,
-			     true));
-	return fe_values.shape_hessians[row][q_point];
-      }
-  }
-
-
-
-  template <int dim, int spacedim>  
-  inline
-  Vector<dim,spacedim>::Vector (const FEValuesBase<dim,spacedim> &fe_values,
-		       const unsigned int       first_vector_component)
-		  :
-		  fe_values (fe_values),
-		  first_vector_component (first_vector_component)
-  {
-    Assert (first_vector_component+dim-1 < this->fe_values.fe->n_components(),
-	    ExcIndexRange(first_vector_component+dim-1, 0,
-			  this->fe_values.fe->n_components()));
+      return hessian_type();
   }
 
 
@@ -3440,53 +3480,19 @@ namespace FEValuesViews
     Assert (fe_values.update_flags & update_values,
 	    typename FVB::ExcAccessToUninitializedField());    
 
-				     // compared to the scalar case above, we
-				     // can save some work here because we
-				     // know that we are querying a contiguous
-				     // range of components
+				     // same as for the scalar case except
+				     // that we have one more index
+				     //
+				     // for primitive elements we could
+				     // probably do even better than the loop
+				     // below because we then know that only
+				     // for one value of 'd' the
+				     // 'if'-condition is true
     value_type return_value;
-    
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-					 // if this is a primitive shape
-					 // function then at most one element
-					 // of the output vector is
-					 // nonzero. find out if indeed one is
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-	if ((nonzero_component >= first_vector_component)
-	    &&
-	    (nonzero_component < first_vector_component + dim))
-	  return_value[nonzero_component - first_vector_component]
-	    = fe_values.shape_values(fe_values.
-				      shape_function_to_row_table[shape_function],
-				      q_point);
-      }
-    else
-      {
-	unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin() +
-			     first_vector_component,
-			     true));
-	for (unsigned int d=0; d<dim; ++d)
-	  if (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
-	      true)
-	    {
-	      return_value[d] = fe_values.shape_values(row, q_point);
-
-	      if ((d != dim-1)
-		  &&
-		  (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d]
-		   == true))
-		++row;
-	    }
-      }
+    for (unsigned int d=0; d<dim; ++d)    
+      if (is_nonzero_shape_function_component(shape_function,d))
+	return_value[d]
+	  = fe_values.shape_values(row_index(shape_function,d),q_point);
 
     return return_value;
   }
@@ -3508,44 +3514,19 @@ namespace FEValuesViews
     Assert (fe_values.update_flags & update_gradients,
 	    typename FVB::ExcAccessToUninitializedField());    
 
+				     // same as for the scalar case except
+				     // that we have one more index
+				     //
+				     // for primitive elements we could
+				     // probably do even better than the loop
+				     // below because we then know that only
+				     // for one value of 'd' the
+				     // 'if'-condition is true
     gradient_type return_value;
-    
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-	if ((nonzero_component >= first_vector_component)
-	    &&
-	    (nonzero_component < first_vector_component + dim))
-	  return_value[nonzero_component - first_vector_component]
-	    = fe_values.shape_gradients[fe_values.
-					 shape_function_to_row_table[shape_function]][q_point];
-      }
-    else
-      {
-	unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin() +
-			     first_vector_component,
-			     true));
-	for (unsigned int d=0; d<dim; ++d)
-	  if (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
-	      true)
-	    {
-	      return_value[d] = fe_values.shape_gradients[row][q_point];
-
-	      if ((d != dim-1)
-		  &&
-		  (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d]
-		   == true))
-		++row;
-	    }
-      }
+    for (unsigned int d=0; d<dim; ++d)    
+      if (is_nonzero_shape_function_component(shape_function,d))
+	return_value[d]
+	  = fe_values.shape_gradients[row_index(shape_function,d)][q_point];
 
     return return_value;
   }
@@ -3567,48 +3548,21 @@ namespace FEValuesViews
     Assert (fe_values.update_flags & update_gradients,
 	    typename FVB::ExcAccessToUninitializedField());    
 
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-	if ((nonzero_component >= first_vector_component)
-	    &&
-	    (nonzero_component < first_vector_component + dim))
-	  return
-	    fe_values.shape_gradients[fe_values.
-				       shape_function_to_row_table[shape_function]][q_point]
-	    [nonzero_component - first_vector_component];
-	else
-	  return 0;
-      }
-    else
-      {
-	unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin() +
-			     first_vector_component,
-			     true));
-
-	double div = 0;
-	for (unsigned int d=0; d<dim; ++d)
-	  if (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
-	      true)
-	    {
-	      div += fe_values.shape_gradients[row][q_point][d];
-
-	      if ((d != dim-1)
-		  &&
-		  (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d]
-		   == true))
-		++row;
-	    }
-	return div;
-      }
+				     // same as for the scalar case except
+				     // that we have one more index
+				     //
+				     // for primitive elements we could
+				     // probably do even better than the loop
+				     // below because we then know that only
+				     // for one value of 'd' the
+				     // 'if'-condition is true
+    divergence_type return_value = 0;
+    for (unsigned int d=0; d<dim; ++d)    
+      if (is_nonzero_shape_function_component(shape_function,d))
+	return_value
+	  += fe_values.shape_gradients[row_index(shape_function,d)][q_point][d];
+
+    return return_value;
   }
   
 
@@ -3627,285 +3581,31 @@ namespace FEValuesViews
     Assert (fe_values.update_flags & update_hessians,
 	    typename FVB::ExcAccessToUninitializedField());    
 
+				     // same as for the scalar case except
+				     // that we have one more index
+				     //
+				     // for primitive elements we could
+				     // probably do even better than the loop
+				     // below because we then know that only
+				     // for one value of 'd' the
+				     // 'if'-condition is true
     hessian_type return_value;
-    
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-	if ((nonzero_component >= first_vector_component)
-	    &&
-	    (nonzero_component < first_vector_component + dim))
-	  return_value[nonzero_component - first_vector_component]
-	    = fe_values.shape_hessians[fe_values.
-					shape_function_to_row_table[shape_function]][q_point];
-      }
-    else
-      {
-	unsigned int
-	  row = (fe_values.shape_function_to_row_table[shape_function]
-		 +
-		 std::count (fe_values.fe->get_nonzero_components(shape_function).begin(),
-			     fe_values.fe->get_nonzero_components(shape_function).begin() +
-			     first_vector_component,
-			     true));
-	for (unsigned int d=0; d<dim; ++d)
-	  if (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d] ==
-	      true)
-	    {
-	      return_value[d] = fe_values.shape_hessians[row][q_point];
-
-	      if ((d != dim-1)
-		  &&
-		  (fe_values.fe->get_nonzero_components(shape_function)[first_vector_component+d]
-		   == true))
-		++row;
-	    }
-      }
+    for (unsigned int d=0; d<dim; ++d)    
+      if (is_nonzero_shape_function_component(shape_function,d))
+	return_value[d]
+	  = fe_values.shape_hessians[row_index(shape_function,d)][q_point];
 
     return return_value;
   }
 
 
-
-				   // we duplicate the following function for
-				   // each dimension since we need to
-				   // initialize a few arrays of dimension
-				   // dependent size. even though this happens
-				   // in a switch(dim) clause, we get warnings
-				   // from the compiler if we don't separate
-				   // things into different functions
-  template <>
-  inline
-  Vector<1>::symmetric_gradient_type
-  Vector<1>::symmetric_gradient (const unsigned int shape_function,
-				 const unsigned int q_point) const
-  {
-				     // this function works like in the case
-				     // above    
-    Assert (shape_function < fe_values.fe->dofs_per_cell,
-	    ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
-    Assert (fe_values.update_flags & update_gradients,
-	    FEValuesBase<1>::ExcAccessToUninitializedField());    
-
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-	if (nonzero_component == first_vector_component)
-	  {
-					     // first get the one component of
-					     // the nonsymmetrized gradient
-					     // that is not zero
-	    const Tensor<1,1> grad
-	      = fe_values.shape_gradients[fe_values.
-					   shape_function_to_row_table[shape_function]][q_point];
-
-					     // then form a symmetric tensor
-					     // out of it. note that access to
-					     // individual elements of a
-					     // SymmetricTensor object is
-					     // fairly slow. if we implemented
-					     // the following in the naive
-					     // way, we would therefore incur
-					     // a very significant penalty: a
-					     // preliminary version of the
-					     // Stokes tutorial program would
-					     // slow down from 17 to 23
-					     // seconds because of this single
-					     // function!
-					     //
-					     // consequently, we try to be a
-					     // bit smarter by already laying
-					     // out the data in the right
-					     // format and creating a
-					     // symmetric tensor of it
-	    const double array[symmetric_gradient_type::n_independent_components]
-	      = { grad[0] };
-	    return symmetric_gradient_type(array);
-	  }
-	else
-	  return symmetric_gradient_type();
-      }
-    else
-      {
-	return symmetrize (gradient(shape_function, q_point));
-      }
-  }
-
-
-
-  template <>
-  inline
-  Vector<2>::symmetric_gradient_type
-  Vector<2>::symmetric_gradient (const unsigned int shape_function,
-				 const unsigned int q_point) const
-  {
-				     // this function works like in the case
-				     // above    
-    Assert (shape_function < fe_values.fe->dofs_per_cell,
-	    ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
-    Assert (fe_values.update_flags & update_gradients,
-	    FEValuesBase<2>::ExcAccessToUninitializedField());    
-
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-					 // first get the one component of
-					 // the nonsymmetrized gradient
-					 // that is not zero
-	const Tensor<1,2> &grad
-	  = fe_values.shape_gradients[fe_values.
-				       shape_function_to_row_table[shape_function]][q_point];
-
-					 // then form a symmetric tensor
-					 // out of it. note that access to
-					 // individual elements of a
-					 // SymmetricTensor object is
-					 // fairly slow. if we implemented
-					 // the following in the naive
-					 // way, we would therefore incur
-					 // a very significant penalty: a
-					 // preliminary version of the
-					 // Stokes tutorial program would
-					 // slow down from 17 to 23
-					 // seconds because of this single
-					 // function!
-					 //
-					 // consequently, we try to be a
-					 // bit smarter by already laying
-					 // out the data in the right
-					 // format and creating a
-					 // symmetric tensor of it
-	switch (nonzero_component - first_vector_component)
-	  {
-	    case 0:
-	    {
-	      const double array[symmetric_gradient_type::n_independent_components]
-		= { grad[0], 0, grad[1]/2 };
-	      return symmetric_gradient_type(array);
-	    }
-
-	    case 1:
-	    {
-	      const double array[symmetric_gradient_type::n_independent_components]
-		= { 0, grad[1], grad[0]/2 };
-	      return symmetric_gradient_type(array);
-	    }
-
-	    default:
-						   // not a shape
-						   // function that
-						   // shared in the
-						   // components of
-						   // the selected
-						   // vector
-		  return symmetric_gradient_type();
-	  }
-      }
-    else
-      {
-	return symmetrize (gradient(shape_function, q_point));
-      }
-  }
-
-
-
-  template <>
+  template <int dim, int spacedim>
   inline
-  Vector<3>::symmetric_gradient_type
-  Vector<3>::symmetric_gradient (const unsigned int shape_function,
-				 const unsigned int q_point) const
+  typename Vector<dim,spacedim>::symmetric_gradient_type
+  Vector<dim,spacedim>::symmetric_gradient (const unsigned int shape_function,
+					    const unsigned int q_point) const
   {
-				     // this function works like in the case
-				     // above    
-    Assert (shape_function < fe_values.fe->dofs_per_cell,
-	    ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
-    Assert (fe_values.update_flags & update_gradients,
-	    FEValuesBase<3>::ExcAccessToUninitializedField());    
-
-    if (fe_values.fe->is_primitive() ||
-	fe_values.fe->is_primitive(shape_function))
-      {
-	const unsigned int
-	  nonzero_component
-	  = fe_values.fe->system_to_component_index(shape_function).first;
-	
-					 // first get the one component of
-					 // the nonsymmetrized gradient
-					 // that is not zero
-	const Tensor<1,3> &grad
-	  = fe_values.shape_gradients[fe_values.
-				       shape_function_to_row_table[shape_function]][q_point];
-
-					 // then form a symmetric
-					 // tensor out of it. note
-					 // that access to individual
-					 // elements of a
-					 // SymmetricTensor object is
-					 // fairly slow. if we
-					 // implemented the following
-					 // in the naive way, we would
-					 // therefore incur a very
-					 // significant penalty: a
-					 // preliminary version of the
-					 // Stokes tutorial program
-					 // would slow down from 17 to
-					 // 23 seconds because of this
-					 // single function!
-					 //
-					 // consequently, we try to be
-					 // a bit smarter by already
-					 // laying out the data in the
-					 // right format and creating
-					 // a symmetric tensor of it
-	switch (nonzero_component-first_vector_component)
-	  {
-	    case 0:
-	    {
-	      const double array[symmetric_gradient_type::n_independent_components]
-		= { grad[0], 0, 0, grad[1]/2 , grad[2]/2, 0};
-	      return symmetric_gradient_type(array);
-	    }
-
-	    case 1:
-	    {
-	      const double array[symmetric_gradient_type::n_independent_components]
-		= { 0, grad[1], 0, grad[0]/2, 0, grad[2]/2 };
-	      return symmetric_gradient_type(array);
-	    }
-
-	    case 2:
-	    {
-	      const double array[symmetric_gradient_type::n_independent_components]
-		= { 0, 0, grad[2], 0, grad[0]/2, grad[1]/2 };
-	      return symmetric_gradient_type(array);
-	    }
-
-	    default:
-						   // not a shape
-						   // function that
-						   // shared in the
-						   // components of
-						   // the selected
-						   // vector
-		  return symmetric_gradient_type();
-	  }
-      }
-    else
-      {
-	return symmetrize (gradient(shape_function, q_point));
-      }
+    return symmetrize (gradient(shape_function, q_point));
   }
 }
 
@@ -3922,23 +3622,31 @@ namespace FEValuesViews
 
 template <int dim, int spacedim>
 inline
-FEValuesViews::Scalar<dim,spacedim>
+const FEValuesViews::Scalar<dim,spacedim> &
 FEValuesBase<dim,spacedim>::
 operator[] (const FEValuesExtractors::Scalar &scalar) const
 {
-  return FEValuesViews::Scalar<dim,spacedim> (*this, scalar.component);
+  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>  
 inline
-FEValuesViews::Vector<dim,spacedim>
+const FEValuesViews::Vector<dim,spacedim> &
 FEValuesBase<dim,spacedim>::
 operator[] (const FEValuesExtractors::Vector &vector) const
 {
-  return
-    FEValuesViews::Vector<dim,spacedim> (*this, vector.first_vector_component);
+  Assert (vector.first_vector_component <
+	  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];
 }
 
 
diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc
index 808b78ca7f..6bee51751b 100644
--- a/deal.II/deal.II/source/fe/fe_values.cc
+++ b/deal.II/deal.II/source/fe/fe_values.cc
@@ -31,6 +31,217 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+
+namespace
+{
+  template <int dim, int spacedim>
+  inline
+  std::vector<unsigned int>
+  make_shape_function_to_row_table (const FiniteElement<dim,spacedim> &fe)
+  {
+    std::vector<unsigned int> shape_function_to_row_table (fe.dofs_per_cell);
+    unsigned int row = 0;
+    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+      {
+	shape_function_to_row_table[i] = row;
+	row += fe.n_nonzero_components (i);
+      }
+
+    return shape_function_to_row_table;
+  }
+}
+
+
+
+namespace FEValuesViews
+{
+  template <int dim, int spacedim>  
+  Scalar<dim,spacedim>::Scalar (const FEValuesBase<dim,spacedim> &fe_values,
+		       const unsigned int       component)
+		  :
+		  fe_values (fe_values),
+		  component (component),
+		  is_nonzero_shape_function_component (fe_values.fe->dofs_per_cell),
+		  row_index (fe_values.fe->dofs_per_cell)
+  {
+    Assert (component < fe_values.fe->n_components(),
+	    ExcIndexRange(component, 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 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)
+	  is_nonzero_shape_function_component[i]
+	    = (component ==
+	       fe_values.fe->system_to_component_index(i).first);
+	else
+	  is_nonzero_shape_function_component[i]
+	    = (fe_values.fe->get_nonzero_components(i)[component]
+	       == true);
+
+	if (is_nonzero_shape_function_component[i] == true)
+	  {
+	    if (is_primitive == true)
+	      row_index[i] = shape_function_to_row_table[i];
+	    else
+	      row_index[i] = (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
+	  row_index[i] = numbers::invalid_unsigned_int;
+      }
+  }
+
+
+
+  template <int dim, int spacedim>
+  Scalar<dim,spacedim>::Scalar ()
+		  :
+		  fe_values (*static_cast<dealii::FEValuesBase<dim,spacedim>*>(0)),
+		  component (numbers::invalid_unsigned_int)
+  {}
+
+
+  template <int dim, int spacedim>
+  Scalar<dim,spacedim> &
+  Scalar<dim,spacedim>::operator= (const Scalar<dim,spacedim> &)
+  {
+				     // we shouldn't be copying these objects
+    Assert (false, ExcInternalError());
+    return *this;
+  }
+  
+
+
+  template <int dim, int spacedim>  
+  Vector<dim,spacedim>::Vector (const FEValuesBase<dim,spacedim> &fe_values,
+				const unsigned int       first_vector_component)
+		  :
+		  fe_values (fe_values),
+		  first_vector_component (first_vector_component),
+		  is_nonzero_shape_function_component (fe_values.fe->dofs_per_cell,
+						       dim),
+		  row_index (fe_values.fe->dofs_per_cell,
+			     dim)
+  {
+    Assert (first_vector_component+dim-1 < fe_values.fe->n_components(),
+	    ExcIndexRange(first_vector_component+dim-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<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() ||
+				       fe_values.fe->is_primitive(i));
+
+	    if (is_primitive == true)
+	      is_nonzero_shape_function_component[i][d]
+		= (component ==
+		   fe_values.fe->system_to_component_index(i).first);
+	    else
+	      is_nonzero_shape_function_component[i][d]
+		= (fe_values.fe->get_nonzero_components(i)[component]
+		   == true);
+
+	    if (is_nonzero_shape_function_component[i][d] == true)
+	      {
+		if (is_primitive == true)
+		  row_index[i][d] = shape_function_to_row_table[i];
+		else
+		  row_index[i][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
+	      row_index[i][d] = numbers::invalid_unsigned_int;
+	  }
+      }
+  }
+  
+
+  template <int dim, int spacedim>
+  Vector<dim,spacedim>::Vector ()
+		  :
+		  fe_values (*static_cast<dealii::FEValuesBase<dim,spacedim>*>(0)),
+		  first_vector_component (numbers::invalid_unsigned_int)
+  {}
+
+
+
+  template <int dim, int spacedim>
+  Vector<dim,spacedim> &
+  Vector<dim,spacedim>::operator= (const Vector<dim,spacedim> &)
+  {
+				     // we shouldn't be copying these objects
+    Assert (false, ExcInternalError());
+    return *this;
+  }
+}
+
+
+namespace internal
+{
+  namespace FEValuesViews
+  {
+    template <int dim, int spacedim>
+    Cache<dim,spacedim>::Cache (const FEValuesBase<dim,spacedim> &fe_values)
+    {
+      const FiniteElement<dim,spacedim> &fe = fe_values.get_fe();
+
+				       // create the views objects. allocate a
+				       // bunch of default-constructed ones
+				       // then destroy them again and do
+				       // in-place construction of those we
+				       // actually want to use (copying stuff
+				       // is wasteful and we can't do that
+				       // anyway because the class has
+				       // reference members)
+      const unsigned int n_scalars = fe.n_components();
+      scalars.resize (n_scalars);
+      for (unsigned int component=0; component<n_scalars; ++component)
+	{
+	  scalars[component].
+	    dealii::FEValuesViews::Scalar<dim,spacedim>::~Scalar ();
+	  new (&scalars[component])
+	    dealii::FEValuesViews::Scalar<dim,spacedim>(fe_values,
+							component);
+	}
+      
+      const unsigned int n_vectors = (fe.n_components() >= dim ?
+				      fe.n_components()-dim+1 :
+				      0);
+      vectors.resize (n_vectors);
+      for (unsigned int component=0; component<n_vectors; ++component)
+	{
+	  vectors[component].
+	    dealii::FEValuesViews::Vector<dim,spacedim>::~Vector ();
+	  new (&vectors[component])
+	    dealii::FEValuesViews::Vector<dim,spacedim>(fe_values,
+							component);
+	}	
+    }
+  }
+}
+
+
 /* ---------------- FEValuesBase<dim,spacedim>::CellIteratorBase --------- */
 
 template <int dim, int spacedim>
@@ -407,19 +618,16 @@ FEValuesData<dim,spacedim>::initialize (const unsigned int        n_quadrature_p
 				   // from shape function number to
 				   // the rows in the tables denoting
 				   // its first non-zero
-				   // component. with this also count
-				   // the total number of non-zero
-				   // components accumulated over all
-				   // shape functions
-  this->shape_function_to_row_table.resize (fe.dofs_per_cell);
-  unsigned int row = 0;
+				   // component
+  this->shape_function_to_row_table
+    = make_shape_function_to_row_table (fe);
+
+				   // count the total number of non-zero
+				   // components accumulated over all shape
+				   // functions
+  unsigned int n_nonzero_shape_components = 0;
   for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-    {
-      this->shape_function_to_row_table[i] = row;
-      row += fe.n_nonzero_components (i);
-    };
-	 
-  const unsigned int n_nonzero_shape_components = row;
+    n_nonzero_shape_components += fe.n_nonzero_components (i);	 
   Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
 	  ExcInternalError());
 
@@ -483,7 +691,8 @@ FEValuesBase<dim,spacedim>::FEValuesBase (const unsigned int n_q_points,
 		mapping(&mapping),
 		fe(&fe),
 		mapping_data(0),
-		fe_data(0)
+		fe_data(0),
+		fe_values_views_cache (*this)
 {
   this->update_flags = flags;
 }
@@ -1514,11 +1723,11 @@ FEValues<dim,spacedim>::FEValues (const Mapping<dim,spacedim>       &mapping,
 			 const UpdateFlags         update_flags)
 		:
 		FEValuesBase<dim,spacedim> (q.size(),
-				   fe.dofs_per_cell,
-				   update_default,
-				   mapping,
-				   fe),
-  quadrature (q)
+					    fe.dofs_per_cell,
+					    update_default,
+					    mapping,
+					    fe),
+		quadrature (q)
 {
   initialize (update_flags);
 }