]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Modify the FEValuesViews classes so that they compute a bunch of things already at...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 11 Dec 2008 19:01:53 +0000 (19:01 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 11 Dec 2008 19:01:53 +0000 (19:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@17913 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 5798daaef20d46270f7c2902c567779d0ecd8ab2..28d39c7b53957d522332c450daa4d072170ea1ad 100644 (file)
@@ -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];
 }
 
 
index 808b78ca7f8bea187ce0a2d5b3f3e33cdc0e9438..6bee51751b3a498003fc56e8ae185646eee6eddf 100644 (file)
 
 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);
 }

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.