]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add get_function_laplacians() methods in FEValues class, improve performance of get_f...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Dec 2008 23:16:45 +0000 (23:16 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Dec 2008 23:16:45 +0000 (23:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@17954 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
deal.II/doc/doxygen/headers/vector_valued.h
deal.II/doc/news/changes.h
deal.II/examples/step-31/step-31.cc
deal.II/examples/step-32/step-32.cc

index a3b8331285877e726c2614f2170224cfdeca483e..61b36bce72cf963ebb7cb32370663805ccfe7d57 100644 (file)
@@ -281,7 +281,7 @@ namespace FEValuesViews
                                        */
       template <class InputVector>
       void get_function_gradients (const InputVector& fe_function,
-                                  std::vector<gradient_type>& values) const;
+                                  std::vector<gradient_type>& gradients) const;
 
                                       /**
                                        * Return the Hessians of the selected
@@ -302,7 +302,33 @@ namespace FEValuesViews
                                        */
       template <class InputVector>
       void get_function_hessians (const InputVector& fe_function,
-                                 std::vector<hessian_type>& values) const;
+                                 std::vector<hessian_type>& hessians) const;
+
+                                      /**
+                                       * Return the Laplacians of the
+                                       * selected scalar component 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. The Laplacians are
+                                       * the trace of the Hessians.
+                                       *
+                                       * This function is the
+                                       * equivalent of the
+                                       * FEValuesBase::get_function_laplacians
+                                       * function but it only works
+                                       * on the selected scalar
+                                       * component.
+                                       */
+      template <class InputVector>
+      void get_function_laplacians (const InputVector& fe_function,
+                                   std::vector<value_type>& laplacians) const;
       
       
     private:
@@ -431,7 +457,7 @@ namespace FEValuesViews
                                        * finite element, the Hessian is a
                                        * <code>Tensor@<3,dim@></code>.
                                        */
-      typedef Tensor<3,dim>          hessian_type;
+      typedef Tensor<3,spacedim>          hessian_type;
 
                                       /**
                                        * Default constructor. Creates an
@@ -520,7 +546,7 @@ namespace FEValuesViews
       divergence_type
       divergence (const unsigned int shape_function,
                  const unsigned int q_point) const;
-      
+
                                       /**
                                        * Return the Hessian (the tensor of
                                        * rank 2 of all second derivatives) of
@@ -573,7 +599,7 @@ namespace FEValuesViews
                                        */
       template <class InputVector>
       void get_function_gradients (const InputVector& fe_function,
-                                  std::vector<gradient_type>& values) const;
+                                  std::vector<gradient_type>& gradients) const;
 
                                       /**
                                        * Return the symmetrized gradients of
@@ -597,7 +623,7 @@ namespace FEValuesViews
                                        */
       template <class InputVector>
       void get_function_symmetric_gradients (const InputVector& fe_function,
-                                            std::vector<symmetric_gradient_type>& values) const;
+                   std::vector<symmetric_gradient_type>& symmetric_gradients) const;
       
                                       /**
                                        * Return the divergence of the selected
@@ -620,7 +646,7 @@ namespace FEValuesViews
                                        */
       template <class InputVector>
       void get_function_divergences (const InputVector& fe_function,
-                                    std::vector<divergence_type>& values) const;
+                                    std::vector<divergence_type>& divergences) const;
       
                                       /**
                                        * Return the Hessians of the selected
@@ -641,7 +667,33 @@ namespace FEValuesViews
                                        */
       template <class InputVector>
       void get_function_hessians (const InputVector& fe_function,
-                                 std::vector<hessian_type>& values) const;
+                                 std::vector<hessian_type>& hessians) const;
+
+                                      /**
+                                       * Return the Laplacians 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. The Laplacians are
+                                       * the trace of the Hessians.
+                                       *
+                                       * This function is the
+                                       * equivalent of the
+                                       * FEValuesBase::get_function_laplacians
+                                       * function but it only works
+                                       * on the selected vector
+                                       * components.
+                                       */
+      template <class InputVector>
+      void get_function_laplacians (const InputVector& fe_function,
+                                   std::vector<value_type>& laplacians) const;
       
     private:
                                       /**
@@ -972,7 +1024,7 @@ class FEValuesData
                                      */
     std::vector<Point<spacedim> >  normal_vectors;
 
-                                    /** 
+                                    /**
                                      * List of outward vectors normal to the cell  
                                      * surface (line) at the quadrature points
                                      * for the codimension 1 case,
@@ -1848,18 +1900,19 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * BlockVector&lt;T&gt;, or one
                                      * of the sequential PETSc or
                                      * Trilinos vector wrapper
-                                     * classes..It represents a
+                                     * classes. It represents a
                                      * global vector of DoF values
                                      * associated with the DofHandler
                                      * object with which this
                                      * FEValues object was last
                                      * initialized.
                                      *
-                                     * The output are the second derivatives
-                                     * of the function represented by
-                                     * these DoF values, as computed
-                                     * in real space (as opposed to
-                                     * on the unit cell).
+                                     * The output are the second
+                                     * derivatives of the function
+                                     * represented by these DoF
+                                     * values, as computed in real
+                                     * space (as opposed to on the
+                                     * unit cell).
                                      */
     template <class InputVector>
     void
@@ -1919,7 +1972,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     void get_function_hessians (
       const InputVector& fe_function,
       const VectorSlice<const std::vector<unsigned int> >& indices,
-      std::vector<Tensor<2,spacedim> >& result) const;
+      std::vector<Tensor<2,spacedim> >& hessians) const;
 
                                     /**
                                      * Access to the second
@@ -1932,7 +1985,7 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
     void get_function_hessians (
       const InputVector& fe_function,
       const VectorSlice<const std::vector<unsigned int> >& indices,
-      std::vector<std::vector<Tensor<2,spacedim> > >& result,
+      std::vector<std::vector<Tensor<2,spacedim> > >& hessians,
       bool quadrature_points_fastest = false) const;
 
                                     /**
@@ -1953,6 +2006,151 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                  bool = false) const;
     
 
+                                    /**
+                                     * Compute the (scalar) Laplacian
+                                     * of the finite element function
+                                     * characterized by @p
+                                     * fe_function restricted to @p
+                                     * cell at the quadrature
+                                     * points. The Laplacian output
+                                     * vector is equivalent to
+                                     * getting
+                                     * <tt>trace(hessians)</tt>,
+                                     * where <tt>hessian</tt> would
+                                     * be the output of the
+                                     * get_function_hessians()
+                                     * function.
+                                     *
+                                     * The function assumes that the
+                                     * @p laplacians object
+                                     * already has the correct size.
+                                     *
+                                     * This function may only be used if the
+                                     * finite element in use is a scalar one,
+                                     * i.e. has only one vector component. If
+                                     * it is a vector-valued one, then use
+                                     * the other
+                                     * get_function_laplacians()
+                                     * function.
+                                     * 
+                                     * The actual data type of the
+                                     * input vector may be either a
+                                     * Vector&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, or one
+                                     * of the sequential PETSc or
+                                     * Trilinos vector wrapper
+                                     * classes. It represents a
+                                     * global vector of DoF values
+                                     * associated with the DofHandler
+                                     * object with which this
+                                     * FEValues object was last
+                                     * initialized.
+                                     *
+                                     * The output are the traces of
+                                     * the second derivatives
+                                     * (i.e. Laplacians) of the
+                                     * function represented by these
+                                     * DoF values, as computed in
+                                     * real space (as opposed to on
+                                     * the unit cell).
+                                     */
+    template <class InputVector, typename number>
+    void
+    get_function_laplacians (const InputVector& fe_function,
+                            std::vector<number>& laplacians) const;
+
+    
+                                    /**
+                                     * Compute the (scalar) Laplacian
+                                     * of the finite element function
+                                     * characterized by @p
+                                     * fe_function restricted to @p
+                                     * cell at the quadrature
+                                     * points. The Laplacian output
+                                     * vector is equivalent to
+                                     * getting
+                                     * <tt>trace(hessians)</tt>, with
+                                     * <tt>hessian</tt> corresponding
+                                     * to the output of the
+                                     * get_function_hessians()
+                                     * function.
+                                     *
+                                     * The function assumes that the
+                                     * @p laplacians object
+                                     * already has the correct size.
+                                     *
+                                     * This function does the same as
+                                     * the other one with the same
+                                     * name, but applies to
+                                     * vector-valued finite elements.
+                                     *
+                                     * The actual data type of the
+                                     * input vector may be either a
+                                     * Vector&lt;T&gt;,
+                                     * BlockVector&lt;T&gt;, or one
+                                     * of the sequential PETSc or
+                                     * Trilinos vector wrapper
+                                     * classes. It represents a
+                                     * global vector of DoF values
+                                     * associated with the DofHandler
+                                     * object with which this
+                                     * FEValues object was last
+                                     * initialized.
+                                     *
+                                     * The output are the traces of
+                                     * the second derivatives (i.e.
+                                     * Laplacians) of the function
+                                     * represented by these DoF
+                                     * values, as computed in real
+                                     * space (as opposed to on the
+                                     * unit cell).
+                                     */
+    template <class InputVector, typename number>
+    void
+    get_function_laplacians (const InputVector      &fe_function,
+                            std::vector<Vector<number> > &laplacians) const;
+
+                                    /**
+                                     * Access to the second
+                                     * derivatives of a function with
+                                     * more flexibility. see
+                                     * get_function_values() with
+                                     * corresponding arguments.
+                                     */
+    template <class InputVector, typename number>
+    void get_function_laplacians (
+      const InputVector& fe_function,
+      const VectorSlice<const std::vector<unsigned int> >& indices,
+      std::vector<number>& laplacians) const;
+
+                                    /**
+                                     * Access to the second
+                                     * derivatives of a function with
+                                     * more flexibility. see
+                                     * get_function_values() with
+                                     * corresponding arguments.
+                                     */
+    template <class InputVector, typename number>
+    void get_function_laplacians (
+      const InputVector& fe_function,
+      const VectorSlice<const std::vector<unsigned int> >& indices,
+      std::vector<Vector<number> >& laplacians) const;
+
+                                    /**
+                                     * Access to the second
+                                     * derivatives of a function with
+                                     * more flexibility. see
+                                     * get_function_values() with
+                                     * corresponding arguments.
+                                     */
+    template <class InputVector, typename number>
+    void get_function_laplacians (
+      const InputVector& fe_function,
+      const VectorSlice<const std::vector<unsigned int> >& indices,
+      std::vector<std::vector<number> >& laplacians,
+      bool quadrature_points_fastest = false) const;
+
+
                                     //@}
     
                                     /**
@@ -2573,7 +2771,7 @@ class FEValues : public FEValuesBase<dim,spacedim>
                                      * about degrees of
                                      * freedom. These functions are,
                                      * above all, the
-                                     * <tt>get_function_value/gradients/hessians</tt>
+                                     * <tt>get_function_value/gradients/hessians/laplacians</tt>
                                      * functions. If you want to call
                                      * these functions, you have to
                                      * call the @p reinit variants
@@ -2771,7 +2969,7 @@ class FEFaceValuesBase : public FEValuesBase<dim,spacedim>
                                      * function was called.
                                      */
     unsigned int present_face_index;
-    
+
                                     /**
                                      * Store a copy of the quadrature
                                      * formula here.
@@ -3269,122 +3467,6 @@ namespace FEValuesViews
 
 
 
-  template <int dim, int spacedim>
-  template <class InputVector>
-  inline
-  void
-  Scalar<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());
-
-    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)
-       for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-         values[q_point] +=
-           dof_values(shape_function) *
-           fe_values.shape_values(shape_function_data[shape_function]
-                                  .row_index,
-                                  q_point);
-  }
-  
-
-
-
-  template <int dim, int spacedim>  
-  template <class InputVector>
-  inline
-  void
-  Scalar<dim,spacedim>::
-  get_function_gradients (const InputVector &fe_function, 
-                         std::vector<gradient_type> &values) const
-  {
-    typedef FEValuesBase<dim,spacedim> FVB;
-    Assert (fe_values.update_flags & update_gradients,
-           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(), 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)
-       for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-         values[q_point] += 
-           dof_values(shape_function) *
-           fe_values.shape_gradients[shape_function_data[shape_function]
-                                     .row_index][q_point];
-  }
-
-
-
-  template <int dim, int spacedim>  
-  template <class InputVector>
-  inline
-  void
-  Scalar<dim,spacedim>::
-  get_function_hessians (const InputVector &fe_function, 
-                        std::vector<hessian_type> &values) const
-  {
-    typedef FEValuesBase<dim,spacedim> FVB;
-    Assert (fe_values.update_flags & update_hessians,
-           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(), hessian_type());
-
-    for (unsigned int shape_function=0;
-        shape_function<fe_values.fe->dofs_per_cell; ++shape_function)
-      for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-       if (shape_function_data[shape_function].is_nonzero_shape_function_component)
-         values[q_point] += 
-           dof_values(shape_function) *
-           fe_values.shape_hessians[shape_function_data[shape_function]
-                                    .row_index][q_point];
-  }
-  
-
-
   template <int dim, int spacedim>  
   inline
   typename Vector<dim,spacedim>::value_type
@@ -3646,270 +3728,6 @@ namespace FEValuesViews
        return symmetrize(return_value);
       }
   }
-
-
-
-  template <int dim, int spacedim>
-  template <class InputVector>
-  inline
-  void
-  Vector<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());
-
-    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;
-       else if (snc != -1)
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           values[q_point][shape_function_data[shape_function].single_nonzero_component_index]
-             += dof_values(shape_function) * fe_values.shape_values(snc,q_point);
-       else
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           for (unsigned int d=0; d<dim; ++d)    
-             if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
-               values[q_point][d]
-                 += (dof_values(shape_function) *
-                     fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point));
-      }
-  }
-  
-
-
-
-  template <int dim, int spacedim>  
-  template <class InputVector>
-  inline
-  void
-  Vector<dim,spacedim>::
-  get_function_gradients (const InputVector &fe_function, 
-                         std::vector<gradient_type> &values) const
-  {
-    typedef FEValuesBase<dim,spacedim> FVB;
-    Assert (fe_values.update_flags & update_gradients,
-           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(), 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
-         continue;
-       else if (snc != -1)
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           values[q_point][shape_function_data[shape_function].single_nonzero_component_index]
-             += dof_values(shape_function) * fe_values.shape_gradients[snc][q_point];
-       else
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           for (unsigned int d=0; d<dim; ++d)    
-             if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
-               values[q_point][d]
-                 += (dof_values(shape_function) *
-                     fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point]);
-      }
-  }
-
-
-
-  template <int dim, int spacedim>  
-  template <class InputVector>
-  inline
-  void
-  Vector<dim,spacedim>::
-  get_function_symmetric_gradients (const InputVector &fe_function, 
-                                   std::vector<symmetric_gradient_type> &values) const
-  {
-    typedef FEValuesBase<dim,spacedim> FVB;
-    Assert (fe_values.update_flags & update_gradients,
-           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(), 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
-         continue;
-       else if (snc != -1)
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           values[q_point]
-             += dof_values(shape_function) *
-             symmetrize_single_row (shape_function_data[shape_function].single_nonzero_component_index,
-                                    fe_values.shape_gradients[snc][q_point]);
-       else
-         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)    
-               if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
-                 grad[d] = (dof_values(shape_function) *
-                            fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point]);
-             values[q_point] += symmetrize(grad);
-           }
-      }
-  }
-
-
-
-  template <int dim, int spacedim>  
-  template <class InputVector>
-  inline
-  void
-  Vector<dim,spacedim>::
-  get_function_divergences (const InputVector &fe_function, 
-                           std::vector<divergence_type> &values) const
-  {
-    typedef FEValuesBase<dim,spacedim> FVB;
-    Assert (fe_values.update_flags & update_gradients,
-           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(), 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;
-       else if (snc != -1)
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           values[q_point]
-             += dof_values(shape_function) *
-             fe_values.shape_gradients[snc][q_point][shape_function_data[shape_function].single_nonzero_component_index];
-       else
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           {
-             for (unsigned int d=0; d<dim; ++d)    
-               if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
-                 values[q_point] +=
-                   dof_values(shape_function) *
-                   fe_values.shape_gradients[shape_function_data[shape_function].row_index[d]][q_point][d];
-           }
-      }
-  }
-  
-
-
-  template <int dim, int spacedim>  
-  template <class InputVector>
-  inline
-  void
-  Vector<dim,spacedim>::
-  get_function_hessians (const InputVector &fe_function, 
-                        std::vector<hessian_type> &values) const
-  {
-    typedef FEValuesBase<dim,spacedim> FVB;
-    Assert (fe_values.update_flags & update_hessians,
-           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(), hessian_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;
-       else if (snc != -1)
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           values[q_point][shape_function_data[shape_function].single_nonzero_component_index]
-             += dof_values(shape_function) * fe_values.shape_hessians[snc][q_point];
-       else
-         for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
-           for (unsigned int d=0; d<dim; ++d)    
-             if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
-               values[q_point][d]
-                 += (dof_values(shape_function) *
-                     fe_values.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point]);
-      }
-  }
 }
 
 
@@ -4081,7 +3899,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,
index 235a58f586344cd420ea3ed39c6b5557cf4ab8bf..fad84c5413e19be39fa92606c8a73c4a6a3b40d7 100644 (file)
@@ -57,7 +57,7 @@ namespace FEValuesViews
 {
   template <int dim, int spacedim>  
   Scalar<dim,spacedim>::Scalar (const FEValuesBase<dim,spacedim> &fe_values,
-                      const unsigned int       component)
+                               const unsigned int                component)
                  :
                  fe_values (fe_values),
                  component (component),
@@ -222,6 +222,525 @@ namespace FEValuesViews
     Assert (false, ExcInternalError());
     return *this;
   }
+
+
+
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Scalar<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());
+
+    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)
+       {
+         const double value = dof_values(shape_function);
+         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 <class InputVector>
+  void
+  Scalar<dim,spacedim>::
+  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());    
+    Assert (gradients.size() == fe_values.n_quadrature_points,
+           ExcDimensionMismatch(gradients.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 (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)
+       {
+         const double value = dof_values(shape_function);
+         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)
+           gradients[q_point] += value * *shape_gradient_ptr++;
+       }
+  }
+
+
+
+  template <int dim, int spacedim>  
+  template <class InputVector>
+  void
+  Scalar<dim,spacedim>::
+  get_function_hessians (const InputVector &fe_function, 
+                        std::vector<hessian_type> &hessians) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (fe_values.update_flags & update_hessians,
+           typename FVB::ExcAccessToUninitializedField());
+    Assert (hessians.size() == fe_values.n_quadrature_points,
+           ExcDimensionMismatch(hessians.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 (hessians.begin(), hessians.end(), hessian_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)
+       {
+         const double value = dof_values(shape_function);
+         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)
+           hessians[q_point] += value * *shape_hessian_ptr++;
+       }
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Scalar<dim,spacedim>::
+  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()); 
+    Assert (laplacians.size() == fe_values.n_quadrature_points,
+           ExcDimensionMismatch(laplacians.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 (laplacians.begin(), laplacians.end(), value_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)
+       {
+         const double value = dof_values(shape_function);
+         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] += 
+             value * trace(fe_values.shape_hessians[row_index][q_point]);
+       }
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Vector<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());
+
+    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 (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[q_point][comp] += value * *shape_value_ptr++;
+         }
+       else
+         for (unsigned int d=0; d<dim; ++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[q_point][d] += value * *shape_value_ptr++;
+             }
+      }
+  }
+  
+
+
+
+  template <int dim, int spacedim>  
+  template <class InputVector>
+  void
+  Vector<dim,spacedim>::
+  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());    
+    Assert (gradients.size() == fe_values.n_quadrature_points,
+           ExcDimensionMismatch(gradients.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 (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
+         continue;
+
+       const double value = dof_values(shape_function);
+       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)
+             gradients[q_point][comp] += value * *shape_gradient_ptr++;
+         }
+       else
+         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 = 
+                 &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)
+                 gradients[q_point][d] += value * *shape_gradient_ptr++;
+             }
+      }
+  }
+
+
+
+  template <int dim, int spacedim>  
+  template <class InputVector>
+  void
+  Vector<dim,spacedim>::
+  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());    
+    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,
+           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 (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
+         continue;
+
+       const double value = dof_values(shape_function);
+       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)
+             symmetric_gradients[q_point]
+               += value * symmetrize_single_row (comp,*shape_gradient_ptr++);
+         }
+       else
+         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)    
+               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];
+             symmetric_gradients[q_point] += symmetrize(grad);
+           }
+      }
+  }
+
+
+
+  template <int dim, int spacedim>  
+  template <class InputVector>
+  void
+  Vector<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 (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)
+             divergences[q_point] += value * (*shape_gradient_ptr++)[comp];
+         }
+       else
+         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 =
+                 &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)
+                 divergences[q_point] += value * (*shape_gradient_ptr++)[d];
+             }
+      }
+  }
+  
+
+
+  template <int dim, int spacedim>  
+  template <class InputVector>
+  void
+  Vector<dim,spacedim>::
+  get_function_hessians (const InputVector &fe_function, 
+                        std::vector<hessian_type> &hessians) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (fe_values.update_flags & update_hessians,
+           typename FVB::ExcAccessToUninitializedField());
+    Assert (hessians.size() == fe_values.n_quadrature_points,
+           ExcDimensionMismatch(hessians.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 (hessians.begin(), hessians.end(), hessian_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 (snc != -1)
+         {
+           const unsigned int comp = 
+             shape_function_data[shape_function].single_nonzero_component_index;
+           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++;
+         }
+       else
+         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 = 
+                 &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)
+                 hessians[q_point][d] += value * *shape_hessian_ptr++;
+             }
+      }
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Vector<dim,spacedim>::
+  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());    
+    Assert (laplacians.size() == fe_values.n_quadrature_points,
+           ExcDimensionMismatch(laplacians.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 (laplacians.begin(), laplacians.end(), value_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 (snc != -1)
+         {
+           const unsigned int comp =
+             shape_function_data[shape_function].single_nonzero_component_index;
+           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++);
+         }
+       else
+         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 = 
+                 &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)
+                 laplacians[q_point][d] += value * trace(*shape_hessian_ptr++);
+             }
+      }
+  }
 }
 
 
@@ -1222,7 +1741,22 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // deal with scalar finite
                                   // elements, so no need to check
                                   // for non-primitivity of shape
-                                  // functions
+                                  // functions. in order to increase
+                                  // the speed of this function, we
+                                  // directly access the data in the
+                                  // shape_values array, and
+                                  // increment pointers for accessing
+                                  // the data. this saves some lookup
+                                  // time and indexing. moreover, the
+                                  // order of the loops is such that
+                                  // we can access the shape_values
+                                  // data stored contiguously (which
+                                  // is also advantageous because
+                                  // access to dof_values is
+                                  // generally more expensive than
+                                  // access to the std::vector values
+                                  // - so we do the cheaper operation
+                                  // in the innermost loop)
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = dof_values(shape_func);
@@ -1262,7 +1796,23 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // deal with scalar finite
                                   // elements, so no need to check
                                   // for non-primitivity of shape
-                                  // functions
+                                  // functions. in order to increase
+                                  // the speed of this function, we
+                                  // directly access the data in the
+                                  // shape_values array, and
+                                  // increment pointers for accessing
+                                  // the data. this saves some lookup
+                                  // time and indexing. moreover, the
+                                  // order of the loops is such that
+                                  // we can access the shape_values
+                                  // data stored contiguously (which
+                                  // is also advantageous because
+                                  // access to the global vector
+                                  // fe_function is more expensive
+                                  // than access to the small
+                                  // std::vector values - so we do
+                                  // the cheaper operation in the
+                                  // innermost loop)
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = fe_function(indices[shape_func]);
@@ -1318,13 +1868,41 @@ void FEValuesBase<dim,spacedim>::get_function_values (
                                   // shape function is primitive or
                                   // not. if it is, then set its only
                                   // non-zero component, otherwise
-                                  // loop over components
+                                  // loop over components. in order
+                                  // to increase the speed of this
+                                  // function, we directly access the
+                                  // data in the shape_values array,
+                                  // and increment pointers for
+                                  // accessing the data. this saves
+                                  // some lookup time and
+                                  // indexing. moreover, in order of
+                                  // the loops is such that we can
+                                  // access the shape_values data
+                                  // stored contiguously (which is
+                                  // also advantageous because access
+                                  // to the global vector fe_function
+                                  // is more expensive than access to
+                                  // the small std::vector values -
+                                  // so we do the cheaper operation
+                                  // in the innermost loop)
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = dof_values(shape_func);
 
       if (fe->is_primitive(shape_func))
        {
+                                  // compared to the scalar
+                                  // functions, finding the correct
+                                  // index in the shape_value table
+                                  // is more involved, since we have
+                                  // to find the row in shape_values
+                                  // that corresponds to the present
+                                  // shape_func. this is done
+                                  // manually in the same way as in
+                                  // shape_value_component() (that
+                                  // function can't be used because
+                                  // it doesn't return us a pointer
+                                  // to the data).
          const unsigned int 
            row = fe->is_primitive() ? 
                  shape_func : this->shape_function_to_row_table[shape_func];
@@ -1334,6 +1912,8 @@ void FEValuesBase<dim,spacedim>::get_function_values (
          for (unsigned int point=0; point<n_quadrature_points; ++point)
            values[point](comp) += value * *shape_value_ptr++;
        }
+                                  // non-primitive case (vector-valued
+                                  // element)
       else
        for (unsigned int c=0; c<n_components; ++c)
          {
@@ -1522,6 +2102,7 @@ void FEValuesBase<dim,spacedim>::get_function_values (
            const double *shape_value_ptr = &this->shape_values(row, 0);
            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)
                values[comp][point] += value * *shape_value_ptr++;
@@ -1589,11 +2170,29 @@ FEValuesBase<dim,spacedim>::get_function_gradients (
                                   // deal with scalar finite
                                   // elements, so no need to check
                                   // for non-primitivity of shape
-                                  // functions
+                                  // functions. in order to increase
+                                  // the speed of this function, we
+                                  // directly access the data in the
+                                  // shape_gradients array, and
+                                  // increment pointers for accessing
+                                  // the data. this saves some lookup
+                                  // time and indexing. moreover, the
+                                  // order of the loops is such that
+                                  // we can access the
+                                  // shape_gradients data stored
+                                  // contiguously (which is also
+                                  // advantageous because access to
+                                  // the vector dof_values is
+                                  // gerenally more expensive than
+                                  // access to the std::vector
+                                  // gradients - so we do the cheaper
+                                  // operation in the innermost loop)
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = dof_values(shape_func);
-      const Tensor<1,spacedim> *shape_gradient_ptr = &this->shape_gradients[shape_func][0];
+      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++;
     }
@@ -1629,11 +2228,29 @@ void FEValuesBase<dim,spacedim>::get_function_gradients (
                                   // deal with scalar finite
                                   // elements, so no need to check
                                   // for non-primitivity of shape
-                                  // functions
+                                  // functions. in order to increase
+                                  // the speed of this function, we
+                                  // directly access the data in the
+                                  // shape_gradients array, and
+                                  // increment pointers for accessing
+                                  // the data. this saves some lookup
+                                  // time and indexing. moreover, the
+                                  // order of the loops is such that
+                                  // we can access the
+                                  // shape_gradients data stored
+                                  // contiguously (which is also
+                                  // advantageous because access to
+                                  // the global vector fe_function is
+                                  // more expensive than access to
+                                  // the small std::vector gradients
+                                  // - so we do the cheaper operation
+                                  // in the innermost loop)
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = fe_function(indices[shape_func]);
-      const Tensor<1,spacedim> *shape_gradient_ptr = &this->shape_gradients[shape_func][0];
+      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++;
     }
@@ -1865,7 +2482,9 @@ get_function_hessians (const InputVector           &fe_function,
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = dof_values(shape_func);
-      const Tensor<2,spacedim> *shape_hessians_ptr = &this->shape_hessians[shape_func][0];
+      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++;
     }
@@ -1905,7 +2524,9 @@ void FEValuesBase<dim,spacedim>::get_function_hessians (
   for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
     {
       const double value = fe_function(indices[shape_func]);
-      const Tensor<2,spacedim> *shape_hessians_ptr = &this->shape_hessians[shape_func][0];
+      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++;
     }
@@ -2110,6 +2731,387 @@ void FEValuesBase<dim, spacedim>::get_function_hessians (
 
 
 
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+  const InputVector   &fe_function,
+  std::vector<number> &laplacians) const
+{
+  Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
+  Assert (fe->n_components() == 1,
+         ExcDimensionMismatch(fe->n_components(), 1));
+  Assert (laplacians.size() == n_quadrature_points,
+         ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+  Assert (present_cell.get() != 0,
+         ExcMessage ("FEValues object is not reinit'ed to any cell"));
+  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
+  std::fill_n (laplacians.begin(), n_quadrature_points, 0);
+
+                                  // add up contributions of trial
+                                  // functions. note that here we
+                                  // deal with scalar finite
+                                  // elements, so no need to check
+                                  // for non-primitivity of shape
+                                  // functions. note that the
+                                  // laplacian is the trace of the
+                                  // hessian, so we use a pointer to
+                                  // the hessians and get their trace
+  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+    {
+      const double value = dof_values(shape_func);
+      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++);
+    }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+  const InputVector& fe_function,
+  const VectorSlice<const std::vector<unsigned int> >& indices,
+  std::vector<number> &laplacians) const
+{
+  Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
+                                  // This function fills a single
+                                  // component only
+  Assert (fe->n_components() == 1,
+         ExcDimensionMismatch(fe->n_components(), 1));
+                                  // One index for each dof
+  Assert (indices.size() == dofs_per_cell,
+         ExcDimensionMismatch(indices.size(), dofs_per_cell));
+                                  // This vector has one entry for
+                                  // 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
+                                  // elements, so no need to check
+                                  // for non-primitivity of shape
+                                  // functions. note that the
+                                  // laplacian is the trace of the
+                                  // hessian, so we use a pointer to
+                                  // the hessians and get their trace
+  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+    {
+      const double value = fe_function(indices[shape_func]);
+      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++);
+    }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+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.  
+                                  // This vector must correspond to a
+                                  // complete discretization
+//  Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
+//       ExcDimensionMismatch(fe_function.size(),
+//                          present_cell->get_dof_handler().n_dofs()));
+                                  // One entry per quadrature point
+  Assert (present_cell.get() != 0,
+         ExcMessage ("FEValues object is not reinit'ed to any cell"));
+  Assert (laplacians.size() == n_quadrature_points,
+         ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+
+  const unsigned int n_components = fe->n_components();
+                                  // Assert that we can write all
+                                  // components into the result
+                                  // vectors
+  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);
+
+                                  // add up contributions of trial
+                                  // functions. now check whether the
+                                  // shape function is primitive or
+                                  // not. if it is, then set its only
+                                  // non-zero component, otherwise
+                                  // loop over components
+  for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+    {
+      const double value = dof_values(shape_func);
+
+      if (fe->is_primitive(shape_func))
+       {
+         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;
+         for (unsigned int point=0; point<n_quadrature_points; ++point)
+           laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
+       }
+      else
+       for (unsigned int c=0; c<n_components; ++c)
+         {
+           if (fe->get_nonzero_components(shape_func)[c] == false)
+             continue;
+
+           const unsigned int
+             row = (this->shape_function_to_row_table[shape_func]
+                    +
+                    std::count (fe->get_nonzero_components(shape_func).begin(),
+                                fe->get_nonzero_components(shape_func).begin()+c,
+                                true));
+
+           const Tensor<2,spacedim> *shape_hessian_ptr 
+             = &this->shape_hessians[row][0];
+
+           for (unsigned int point=0; point<n_quadrature_points; ++point)
+             laplacians[point](c) += value * trace(*shape_hessian_ptr++);
+         }
+    }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+  const InputVector& fe_function,
+  const VectorSlice<const std::vector<unsigned int> >& indices,
+  std::vector<Vector<number> >& laplacians) const
+{
+                                  // 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
+                                  // function values is generated in
+                                  // each point.
+  Assert (indices.size() % dofs_per_cell == 0,
+         ExcNotMultiple(indices.size(), dofs_per_cell));
+
+                                  // The number of components of the
+                                  // result may be a multiple of the
+                                  // 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));
+
+                                  // 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);
+
+                                  // add up contributions of trial
+                                  // functions. now check whether the
+                                  // shape function is primitive or
+                                  // not. if it is, then set its only
+                                  // non-zero component, otherwise
+                                  // loop over components
+  for (unsigned int mc = 0; mc < component_multiple; ++mc)
+    for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+      {
+       const double value = fe_function(indices[shape_func+mc*dofs_per_cell]);
+
+       if (fe->is_primitive(shape_func))
+         {
+           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 
+                                      + mc * n_components;
+           for (unsigned int point=0; point<n_quadrature_points; ++point)
+             laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
+         }
+       else
+         for (unsigned int c=0; c<n_components; ++c)
+           {
+             if (fe->get_nonzero_components(shape_func)[c] == false)
+               continue;
+
+             const unsigned int
+               row = (this->shape_function_to_row_table[shape_func]
+                      +
+                      std::count (fe->get_nonzero_components(shape_func).begin(),
+                                  fe->get_nonzero_components(shape_func).begin()+c,
+                                  true));
+
+             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++);
+           }
+      }
+}
+
+
+
+template <int dim, int spacedim>
+template <class InputVector, typename number>
+void FEValuesBase<dim,spacedim>::get_function_laplacians (
+  const InputVector& fe_function,
+  const VectorSlice<const std::vector<unsigned int> >& indices,
+  std::vector<std::vector<number> >& 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
+                                  // function values is generated in
+                                  // each point.
+  Assert (indices.size() % dofs_per_cell == 0,
+         ExcNotMultiple(indices.size(), dofs_per_cell));
+
+                                  // The number of components of the
+                                  // result may be a multiple of the
+                                  // 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)
+    {
+      Assert (laplacians.size() == result_components,
+             ExcDimensionMismatch(laplacians.size(), result_components));
+      for (unsigned i=0;i<laplacians.size();++i)
+       Assert (laplacians[i].size() == n_quadrature_points,
+               ExcDimensionMismatch(laplacians[i].size(), n_quadrature_points));
+    }
+  else
+    {
+      Assert(laplacians.size() == n_quadrature_points,
+            ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
+      for (unsigned i=0;i<laplacians.size();++i)
+       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);
+
+                                  // add up contributions of trial
+                                  // functions. now check whether the
+                                  // shape function is primitive or
+                                  // not. if it is, then set its only
+                                  // non-zero component, otherwise
+                                  // loop over components
+  for (unsigned int mc = 0; mc < component_multiple; ++mc)
+    for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+      {
+       const double value = fe_function(indices[shape_func+mc*dofs_per_cell]);
+
+       if (fe->is_primitive(shape_func))
+         {
+           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 
+                                      + 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++);
+           else
+             for (unsigned int point=0; point<n_quadrature_points; ++point)
+               laplacians[point][comp] += value * trace(*shape_hessian_ptr++);
+         }
+       else
+         for (unsigned int c=0; c<n_components; ++c)
+           {
+             if (fe->get_nonzero_components(shape_func)[c] == false)
+               continue;
+
+             const unsigned int
+               row = (this->shape_function_to_row_table[shape_func]
+                      +
+                      std::count (fe->get_nonzero_components(shape_func).begin(),
+                                  fe->get_nonzero_components(shape_func).begin()+c,
+                                  true));
+
+             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++);
+             else
+               for (unsigned int point=0; point<n_quadrature_points; ++point)
+                 laplacians[point][comp] += value * trace(*shape_hessian_ptr++);
+           }
+      }
+}
+
+
+
 template <int dim, int spacedim>
 const std::vector<Point<spacedim> > &
 FEValuesBase<dim,spacedim>::get_cell_normal_vectors () const
index a85fcd5ae1f3edda4a8bf90e1387a35921436803..5009baf9a52b04cbb2a237af3d3f8d2f3925e209 100644 (file)
 //
 //-----------------------------------------------------------------------
 
+// Instantiations of functions in FEValuesViews (explicitly write
+// dealii:: namespace in order to not confuse the compiler with
+// FEValuesViews::Vector).
+for (VEC : SERIAL_VECTORS)
+  {
+    template
+       void FEValuesViews::Scalar<deal_II_dimension>::get_function_values<dealii::VEC>
+       (const dealii::VEC&, std::vector<double>&) const;
+    template
+       void FEValuesViews::Scalar<deal_II_dimension>::get_function_gradients<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension> >&) const;
+    template
+       void FEValuesViews::Scalar<deal_II_dimension>::get_function_hessians<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<2,deal_II_dimension> >&) const;
+    template
+       void FEValuesViews::Scalar<deal_II_dimension>::get_function_laplacians<dealii::VEC>
+       (const dealii::VEC&, std::vector<double>&) const;
+
+    template
+       void FEValuesViews::Vector<deal_II_dimension>::get_function_values<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension> >&) const;
+    template
+       void FEValuesViews::Vector<deal_II_dimension>::get_function_gradients<dealii::VEC>
+       (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;
+    template
+       void FEValuesViews::Vector<deal_II_dimension>::get_function_divergences<dealii::VEC>
+       (const dealii::VEC&, std::vector<double>&) const;
+    template
+       void FEValuesViews::Vector<deal_II_dimension>::get_function_hessians<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<3,deal_II_dimension> >&) const;
+    template
+       void FEValuesViews::Vector<deal_II_dimension>::get_function_laplacians<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension> >&) const;
+
+
+#if deal_II_dimension != 3
+    template
+       void FEValuesViews::Scalar<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_values<dealii::VEC>
+       (const dealii::VEC&, std::vector<value_type>&) const;
+    template
+       void FEValuesViews::Scalar<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_gradients<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension+1> >&) const;
+    template
+       void FEValuesViews::Scalar<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_hessians<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<2,deal_II_dimension+1> >&) const;
+    template
+       void FEValuesViews::Scalar<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_laplacians<dealii::VEC>
+       (const dealii::VEC&, std::vector<double>&) const;
+
+    template
+       void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_values<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension+1> >&) const;
+    template
+       void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_gradients<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<2,deal_II_dimension+1> >&) const;
+    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;
+    template
+       void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_divergences<dealii::VEC>
+       (const dealii::VEC&, std::vector<double>&) const;
+    template
+       void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_hessians<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<3,deal_II_dimension+1> >&) const;
+    template
+       void FEValuesViews::Vector<deal_II_dimension, deal_II_dimension+1>
+       ::get_function_laplacians<dealii::VEC>
+       (const dealii::VEC&, std::vector<Tensor<1,deal_II_dimension+1> >&) const;
+
+#endif
+
+  }
+
+
+
 // Instantiations of functions in FEValuesBase
 
 for (VEC : SERIAL_VECTORS)
-  {  
+  {
     template
        void FEValuesBase<deal_II_dimension>::get_function_values<VEC>
       (const VEC&, std::vector<double>&) const;
@@ -85,6 +172,44 @@ for (VEC : SERIAL_VECTORS)
       (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
        std::vector<std::vector<Tensor<2,deal_II_dimension> > > &, bool) const;
 
+    template
+       void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<double>&) const;
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<float>&) const;
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&, std::vector<double>&) const;
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&, std::vector<float>&) const;
+
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<Vector<double> > &) const;
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<Vector<float> > &) const;
+
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<Vector<double> > &) const;
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<Vector<float> > &) const;
+
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<std::vector<double> > &, bool) const;
+    template
+      void FEValuesBase<deal_II_dimension>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<std::vector<float> > &, bool) const;
+
 
 #if deal_II_dimension != 3
 
@@ -158,6 +283,45 @@ for (VEC : SERIAL_VECTORS)
       (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
        std::vector<std::vector<Tensor<2,deal_II_dimension+1> > > &, bool) const;
 
+    template
+       void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<double>&) const;
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<float>&) const;
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&, std::vector<double>&) const;
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&, std::vector<float>&) const;
+
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<Vector<double> > &) const;
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, std::vector<Vector<float> > &) const;
+
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<Vector<double> > &) const;
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<Vector<float> > &) const;
+
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<std::vector<double> > &, bool) const;
+    template
+      void FEValuesBase<deal_II_dimension,deal_II_dimension+1>::get_function_laplacians<VEC>
+      (const VEC&, const VectorSlice<const std::vector<unsigned int> >&,
+       std::vector<std::vector<float> > &, bool) const;
+
+
 #endif
 
   }
index 08a9800db677da12698bddb6c03e7cf5d1df7034..8c7f8b704519bbb5c8a371ef6769f5607aaed2ec 100644 (file)
@@ -868,6 +868,11 @@ scalar_product (const Tensor<2,dim> &u,
  * FEValuesViews::Vector::get_function_symmetric_gradients and
  * FEValuesViews::Vector::get_function_divergences.
  *
+ * Moreover, there is a shortcut available in case when only the
+ * Laplacians of the solution (which is the trace of the hessians) is
+ * needed, usable for both scalar and vector-valued problems as
+ * FEValuesViews::Scalar::get_function_laplacians and
+ * FEValuesViews::Vector::get_function_laplacians.
  *
  * @ingroup feall feaccess
  */
index 70a46561be45a9517ec598d3ecdddd8f10f6d72b..4020dd77eb3af7b1c46fee89e06af073d801b8bb 100644 (file)
@@ -598,6 +598,16 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li>
+  <p>
+  New: FEValues objects can now directly calculate the Laplacians of
+  solution functions at quadrature points. This is equivalent to the 
+  trace of the Hessians, but faster and more convenient in user code.
+  @ref step_31 "step-31" program has also been updated to use this.
+  <br>
+  (Martin Kronbichler 2008/12/15)
+  </p>
+
   <li>
   <p>
   Improved: The FEValuesViews objects that one gets when writing things
index 2d1b94c356d579b49d586ed5a4092ad8391cfac0..7ef6eb8f98964292c5e9d373d853589f1fd449c6 100644 (file)
@@ -597,8 +597,8 @@ class BoussinesqFlowProblem
                      const std::vector<double>          &old_old_temperature,
                      const std::vector<Tensor<1,dim> >  &old_temperature_grads,
                      const std::vector<Tensor<1,dim> >  &old_old_temperature_grads,
-                     const std::vector<Tensor<2,dim> >  &old_temperature_hessians,
-                     const std::vector<Tensor<2,dim> >  &old_old_temperature_hessians,
+                     const std::vector<double>          &old_temperature_laplacians,
+                     const std::vector<double>          &old_old_temperature_laplacians,
                      const std::vector<Tensor<1,dim> >  &old_velocity_values,
                      const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
                      const std::vector<double>          &gamma_values,
@@ -963,8 +963,8 @@ compute_viscosity (const std::vector<double>          &old_temperature,
                   const std::vector<double>          &old_old_temperature,
                   const std::vector<Tensor<1,dim> >  &old_temperature_grads,
                   const std::vector<Tensor<1,dim> >  &old_old_temperature_grads,
-                  const std::vector<Tensor<2,dim> >  &old_temperature_hessians,
-                  const std::vector<Tensor<2,dim> >  &old_old_temperature_hessians,
+                  const std::vector<double>          &old_temperature_laplacians,
+                  const std::vector<double>          &old_old_temperature_laplacians,
                   const std::vector<Tensor<1,dim> >  &old_velocity_values,
                   const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
                   const std::vector<double>          &gamma_values,
@@ -994,8 +994,8 @@ compute_viscosity (const std::vector<double>          &old_temperature,
                                   old_old_temperature_grads[q]) / 2;
       
       const double kappa_Delta_T = EquationData::kappa
-                                  * (trace(old_temperature_hessians[q]) +
-                                     trace(old_old_temperature_hessians[q])) / 2;
+                                  * (old_temperature_laplacians[q] +
+                                     old_old_temperature_laplacians[q]) / 2;
 
       const double residual
        = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
@@ -1987,8 +1987,8 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
   std::vector<double>         old_old_temperature_values(n_q_points);
   std::vector<Tensor<1,dim> > old_temperature_grads(n_q_points);
   std::vector<Tensor<1,dim> > old_old_temperature_grads(n_q_points);
-  std::vector<Tensor<2,dim> > old_temperature_hessians(n_q_points);
-  std::vector<Tensor<2,dim> > old_old_temperature_hessians(n_q_points);
+  std::vector<double>         old_temperature_laplacians(n_q_points);
+  std::vector<double>         old_old_temperature_laplacians(n_q_points);
 
   EquationData::TemperatureRightHandSide<dim>  temperature_right_hand_side;
   std::vector<double> gamma_values (n_q_points);
@@ -2039,10 +2039,10 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
       temperature_fe_values.get_function_gradients (old_old_temperature_solution,
                                                    old_old_temperature_grads);
       
-      temperature_fe_values.get_function_hessians (old_temperature_solution,
-                                                  old_temperature_hessians);
-      temperature_fe_values.get_function_hessians (old_old_temperature_solution,
-                                                  old_old_temperature_hessians);
+      temperature_fe_values.get_function_laplacians (old_temperature_solution,
+                                                    old_temperature_laplacians);
+      temperature_fe_values.get_function_laplacians (old_old_temperature_solution,
+                                                  old_old_temperature_laplacians);
       
       temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
                                              gamma_values);
@@ -2079,8 +2079,8 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
                             old_old_temperature_values,
                             old_temperature_grads,
                             old_old_temperature_grads,
-                            old_temperature_hessians,
-                            old_old_temperature_hessians,
+                            old_temperature_laplacians,
+                            old_old_temperature_laplacians,
                             old_velocity_values,
                             old_old_velocity_values,
                             gamma_values,
index affab19e45e9c2c0421026a758eccf18116bff76..73ee39434f1e31e96bcc76a2c30285b098e58090 100644 (file)
@@ -284,8 +284,8 @@ class BoussinesqFlowProblem
                      const std::vector<double>          &old_old_temperature,
                      const std::vector<Tensor<1,dim> >  &old_temperature_grads,
                      const std::vector<Tensor<1,dim> >  &old_old_temperature_grads,
-                     const std::vector<Tensor<2,dim> >  &old_temperature_hessians,
-                     const std::vector<Tensor<2,dim> >  &old_old_temperature_hessians,
+                     const std::vector<double>          &old_temperature_laplacians,
+                     const std::vector<double>          &old_old_temperature_laplacians,
                      const std::vector<Tensor<1,dim> >  &present_velocity_values,
                      const std::vector<double>          &gamma_values,
                      const double                        global_u_infty,
@@ -465,7 +465,7 @@ BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
       double min_temperature, max_temperature;
 
       trilinos_communicator.MaxAll(&max_local_temperature, &max_temperature, 1);
-      trilinos_communicator.MaxAll(&min_local_temperature, &min_temperature, 1);
+      trilinos_communicator.MinAll(&min_local_temperature, &min_temperature, 1);
 
       return std::make_pair(min_temperature, max_temperature);
     }
@@ -497,7 +497,7 @@ BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
       double min_temperature, max_temperature;
 
       trilinos_communicator.MaxAll(&max_local_temperature, &max_temperature, 1);
-      trilinos_communicator.MaxAll(&min_local_temperature, &min_temperature, 1);
+      trilinos_communicator.MinAll(&min_local_temperature, &min_temperature, 1);
 
       return std::make_pair(min_temperature, max_temperature);
     }    
@@ -512,8 +512,8 @@ compute_viscosity(const std::vector<double>          &old_temperature,
                  const std::vector<double>          &old_old_temperature,
                  const std::vector<Tensor<1,dim> >  &old_temperature_grads,
                  const std::vector<Tensor<1,dim> >  &old_old_temperature_grads,
-                 const std::vector<Tensor<2,dim> >  &old_temperature_hessians,
-                 const std::vector<Tensor<2,dim> >  &old_old_temperature_hessians,
+                 const std::vector<double>          &old_temperature_laplacians,
+                 const std::vector<double>          &old_old_temperature_laplaceians,
                  const std::vector<Tensor<1,dim> >  &present_velocity_values,
                  const std::vector<double>          &gamma_values,
                  const double                        global_u_infty,
@@ -544,8 +544,8 @@ compute_viscosity(const std::vector<double>          &old_temperature,
                                   old_old_temperature_grads[q]) / 2;
       
       const double kappa_Delta_T = EquationData::kappa
-                                  * (trace(old_temperature_hessians[q]) +
-                                     trace(old_old_temperature_hessians[q])) / 2;
+                                  * (old_temperature_laplacians[q] +
+                                     old_old_temperature_laplacians[q]) / 2;
 
       const double residual
        = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
@@ -1069,8 +1069,8 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
   std::vector<double>         old_old_temperature_values(n_q_points);
   std::vector<Tensor<1,dim> > old_temperature_grads(n_q_points);
   std::vector<Tensor<1,dim> > old_old_temperature_grads(n_q_points);
-  std::vector<Tensor<2,dim> > old_temperature_hessians(n_q_points);
-  std::vector<Tensor<2,dim> > old_old_temperature_hessians(n_q_points);
+  std::vector<double>         old_temperature_laplacians(n_q_points);
+  std::vector<double>         old_old_temperature_laplacians(n_q_points);
 
   
   EquationData::TemperatureRightHandSide<dim>  temperature_right_hand_side;
@@ -1113,9 +1113,9 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
                                                      old_old_temperature_grads);
        
        temperature_fe_values.get_function_hessians (old_temperature_solution,
-                                                    old_temperature_hessians);
+                                                    old_temperature_laplacians);
        temperature_fe_values.get_function_hessians (old_old_temperature_solution,
-                                                    old_old_temperature_hessians);
+                                                    old_old_temperature_laplacians);
        
        temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
                                                gamma_values);
@@ -1128,8 +1128,8 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
                               old_old_temperature_values,
                               old_temperature_grads,
                               old_old_temperature_grads,
-                              old_temperature_hessians,
-                              old_old_temperature_hessians,
+                              old_temperature_laplacians,
+                              old_old_temperature_laplacians,
                               present_velocity_values,
                               gamma_values,
                               global_u_infty,

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.