]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add one more access function. Reindent comments.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 14 Dec 2012 12:13:31 +0000 (12:13 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 14 Dec 2012 12:13:31 +0000 (12:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@27823 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/matrix_free/fe_evaluation.h

index 573dd3ea3234f686584d36b8208f874a1e79ace7..8164f470486e0ae22bc447230d603d329035a59c 100644 (file)
@@ -96,40 +96,29 @@ public:
    */
   //@{
   /**
-   * Initializes the operation pointer to the
-   * current cell. Unlike the FEValues::reinit
-   * function, where the information related to
-   * a particular cell is generated in the
-   * reinit call, this function is very cheap
-   * since all data is pre-computed in @p
-   * matrix_free, and only a few indices
-   * have to be set appropriately.
+   * Initializes the operation pointer to the current cell. Unlike the
+   * FEValues::reinit function, where the information related to a particular
+   * cell is generated in the reinit call, this function is very cheap since
+   * all data is pre-computed in @p matrix_free, and only a few indices have
+   * to be set appropriately.
    */
   void reinit (const unsigned int cell);
 
   /**
-   * For the transformation information stored
-   * in MappingInfo, this function returns the
-   * index which belongs to the current cell as
-   * specified in @p reinit. Note that
-   * MappingInfo has different fields for
-   * Cartesian cells, cells with affine mapping
-   * and with general mappings, so in order to
-   * access the correct data, this interface
-   * must be used together with get_cell_type.
+   * For the transformation information stored in MappingInfo, this function
+   * returns the index which belongs to the current cell as specified in @p
+   * reinit. Note that MappingInfo has different fields for Cartesian cells,
+   * cells with affine mapping and with general mappings, so in order to
+   * access the correct data, this interface must be used together with
+   * get_cell_type.
    */
   unsigned int get_cell_data_number() const;
 
   /**
-   * Returns the type of the cell the @p
-   * reinit function has been called
-   * for. Valid values are @p cartesian
-   * for Cartesian cells (which allows
-   * for considerable data compression),
-   * @p affine for cells with affine
-   * mappings, and @p general for
-   * general cells without any
-   * compressed storage applied.
+   * Returns the type of the cell the @p reinit function has been called
+   * for. Valid values are @p cartesian for Cartesian cells (which allows for
+   * considerable data compression), @p affine for cells with affine mappings,
+   * and @p general for general cells without any compressed storage applied.
    */
   internal::MatrixFreeFunctions::CellType get_cell_type() const;
 
@@ -140,186 +129,136 @@ public:
    */
   //@{
   /**
-   * For the vector @p src, read out the values
-   * on the degrees of freedom of the current
-   * cell, and store them internally. Similar
-   * functionality as the function
-   * DoFAccessor::get_interpolated_dof_values
-   * when no constraints are present, but it
-   * also includes constraints from hanging
-   * nodes, so one can see it as a similar
-   * function to
-   * ConstraintMatrix::read_dof_values as
-   * well. Note that if vectorization is
-   * enabled, the DoF values for several cells
-   * are set.
+   * For the vector @p src, read out the values on the degrees of freedom of
+   * the current cell, and store them internally. Similar functionality as the
+   * function DoFAccessor::get_interpolated_dof_values when no constraints are
+   * present, but it also includes constraints from hanging nodes, so one can
+   * see it as a similar function to ConstraintMatrix::read_dof_values as
+   * well. Note that if vectorization is enabled, the DoF values for several
+   * cells are set.
    *
-   * If some constraints on the vector are
-   * inhomogeneous, use the function
-   * read_dof_values_plain instead and provide
-   * the vector with useful data also in
-   * constrained positions by calling
-   * ConstraintMatrix::distribute. When
-   * accessing vector entries during the
-   * solution of linear systems, the temporary
-   * solution should always have homogeneous
-   * constraints and this method is the correct
-   * one.
+   * If some constraints on the vector are inhomogeneous, use the function
+   * read_dof_values_plain instead and provide the vector with useful data
+   * also in constrained positions by calling
+   * ConstraintMatrix::distribute. When accessing vector entries during the
+   * solution of linear systems, the temporary solution should always have
+   * homogeneous constraints and this method is the correct one.
    */
   template <typename VectorType>
   void read_dof_values (const VectorType &src);
 
   /**
-   * For a collection of several vector @p src,
-   * read out the values on the degrees of
-   * freedom of the current cell for @p
-   * n_components (template argument), starting
-   * at @p first_index, and store them
-   * internally. Similar functionality as the
-   * function ConstraintMatrix::read_dof_values.
-   * Note that if vectorization is enabled, the
-   * DoF values for several cells are set.
+   * For a collection of several vector @p src, read out the values on the
+   * degrees of freedom of the current cell for @p n_components (template
+   * argument), starting at @p first_index, and store them internally. Similar
+   * functionality as the function ConstraintMatrix::read_dof_values.  Note
+   * that if vectorization is enabled, the DoF values for several cells are
+   * set.
    */
   template <typename VectorType>
   void read_dof_values (const std::vector<VectorType> &src,
                         const unsigned int             first_index=0);
 
   /**
-   * Reads data from several vectors. Same as
-   * other function with std::vector, but
-   * accepts a vector of pointers to vectors.
+   * Reads data from several vectors. Same as other function with std::vector,
+   * but accepts a vector of pointers to vectors.
    */
   template <typename VectorType>
   void read_dof_values (const std::vector<VectorType *> &src,
                         const unsigned int              first_index=0);
 
   /**
-   * For the vector @p src, read out the values
-   * on the degrees of freedom of the current
-   * cell, and store them internally. Similar
-   * functionality as the function
-   * DoFAccessor::get_interpolated_dof_values. As
-   * opposed to the read_dof_values function,
-   * this function reads out the plain entries
-   * from vectors, without taking stored
-   * constraints into account. This way of
-   * access is appropriate when the constraints
-   * have been distributed on the vector by a
-   * call to ConstraintMatrix::distribute
-   * previously. This function is also necessary
-   * when inhomogeneous constraints are to be
-   * used, as MatrixFree can only handle
-   * homogeneous constraints. Note that if
-   * vectorization is enabled, the DoF values
-   * for several cells are set.
+   * For the vector @p src, read out the values on the degrees of freedom of
+   * the current cell, and store them internally. Similar functionality as the
+   * function DoFAccessor::get_interpolated_dof_values. As opposed to the
+   * read_dof_values function, this function reads out the plain entries from
+   * vectors, without taking stored constraints into account. This way of
+   * access is appropriate when the constraints have been distributed on the
+   * vector by a call to ConstraintMatrix::distribute previously. This
+   * function is also necessary when inhomogeneous constraints are to be used,
+   * as MatrixFree can only handle homogeneous constraints. Note that if
+   * vectorization is enabled, the DoF values for several cells are set.
    */
   template <typename VectorType>
   void read_dof_values_plain (const VectorType &src);
 
   /**
-   * For a collection of several vector @p src,
-   * read out the values on the degrees of
-   * freedom of the current cell for @p
-   * n_components (template argument), starting
-   * at @p first_index, and store them
-   * internally. Similar functionality as the
-   * function DoFAccessor::read_dof_values.
-   * Note that if vectorization is enabled, the
-   * DoF values for several cells are set.
+   * For a collection of several vector @p src, read out the values on the
+   * degrees of freedom of the current cell for @p n_components (template
+   * argument), starting at @p first_index, and store them internally. Similar
+   * functionality as the function DoFAccessor::read_dof_values.  Note that if
+   * vectorization is enabled, the DoF values for several cells are set.
    */
   template <typename VectorType>
   void read_dof_values_plain (const std::vector<VectorType> &src,
                               const unsigned int             first_index=0);
 
   /**
-   * Reads data from several vectors without
-   * resolving constraints. Same as other
-   * function with std::vector, but accepts a
-   * vector of pointers to vectors.
+   * Reads data from several vectors without resolving constraints. Same as
+   * other function with std::vector, but accepts a vector of pointers to
+   * vectors.
    */
   template <typename VectorType>
   void read_dof_values_plain (const std::vector<VectorType *> &src,
                               const unsigned int              first_index=0);
 
   /**
-   * Takes the values stored internally on dof
-   * values of the current cell and sums them
-   * into the vector @p dst. The function also
-   * applies constraints during the write
-   * operation. The functionality is hence
-   * similar to the function
-   * ConstraintMatrix::distribute_local_to_global.
-   * Note that if vectorization is enabled, the
-   * DoF values for several cells are used.
+   * Takes the values stored internally on dof values of the current cell and
+   * sums them into the vector @p dst. The function also applies constraints
+   * during the write operation. The functionality is hence similar to the
+   * function ConstraintMatrix::distribute_local_to_global.  Note that if
+   * vectorization is enabled, the DoF values for several cells are used.
    */
   template<typename VectorType>
   void distribute_local_to_global (VectorType &dst) const;
 
   /**
-   * Takes the values stored internally on dof
-   * values of the current cell for a
-   * vector-valued problem consisting of @p
-   * n_components (template argument) and sums
-   * them into the collection of vectors vector
-   * @p dst, starting at index @p
-   * first_index. The function also applies
-   * constraints during the write operation. The
-   * functionality is hence similar to the
-   * function
-   * ConstraintMatrix::distribute_local_to_global.
-   * Note that if vectorization is enabled, the
-   * DoF values for several cells are used.
+   * Takes the values stored internally on dof values of the current cell for
+   * a vector-valued problem consisting of @p n_components (template argument)
+   * and sums them into the collection of vectors vector @p dst, starting at
+   * index @p first_index. The function also applies constraints during the
+   * write operation. The functionality is hence similar to the function
+   * ConstraintMatrix::distribute_local_to_global.  Note that if vectorization
+   * is enabled, the DoF values for several cells are used.
    */
   template<typename VectorType>
   void distribute_local_to_global (std::vector<VectorType> &dst,
                                    const unsigned int       first_index=0) const;
 
   /**
-   * Writes data to several vectors. Same as
-   * other function with std::vector, but
-   * accepts a vector of pointers to vectors.
+   * Writes data to several vectors. Same as other function with std::vector,
+   * but accepts a vector of pointers to vectors.
    */
   template<typename VectorType>
   void distribute_local_to_global (std::vector<VectorType *> &dst,
                                    const unsigned int       first_index=0) const;
 
   /**
-   * Takes the values stored internally on dof
-   * values of the current cell and sums them
-   * into the vector @p dst. The function also
-   * applies constraints during the write
-   * operation. The functionality is hence
-   * similar to the function
-   * ConstraintMatrix::distribute_local_to_global.
-   * Note that if vectorization is enabled, the
-   * DoF values for several cells are used.
+   * Takes the values stored internally on dof values of the current cell and
+   * sums them into the vector @p dst. The function also applies constraints
+   * during the write operation. The functionality is hence similar to the
+   * function ConstraintMatrix::distribute_local_to_global.  Note that if
+   * vectorization is enabled, the DoF values for several cells are used.
    */
   template<typename VectorType>
   void set_dof_values (VectorType &dst) const;
 
   /**
-   * Takes the values stored internally on dof
-   * values of the current cell for a
-   * vector-valued problem consisting of @p
-   * n_components (template argument) and sums
-   * them into the collection of vectors vector
-   * @p dst, starting at index @p
-   * first_index. The function also applies
-   * constraints during the write operation. The
-   * functionality is hence similar to the
-   * function
-   * ConstraintMatrix::distribute_local_to_global.
-   * Note that if vectorization is enabled, the
-   * DoF values for several cells are used.
+   * Takes the values stored internally on dof values of the current cell for
+   * a vector-valued problem consisting of @p n_components (template argument)
+   * and sums them into the collection of vectors vector @p dst, starting at
+   * index @p first_index. The function also applies constraints during the
+   * write operation. The functionality is hence similar to the function
+   * ConstraintMatrix::distribute_local_to_global.  Note that if vectorization
+   * is enabled, the DoF values for several cells are used.
    */
   template<typename VectorType>
   void set_dof_values (std::vector<VectorType> &dst,
                        const unsigned int       first_index=0) const;
 
   /**
-   * Writes data to several vectors. Same as
-   * other function with std::vector, but
-   * accepts a vector of pointers to vectors.
+   * Writes data to several vectors. Same as other function with std::vector,
+   * but accepts a vector of pointers to vectors.
    */
   template<typename VectorType>
   void set_dof_values (std::vector<VectorType *> &dst,
@@ -332,185 +271,128 @@ public:
    */
   //@{
   /**
-   * Returns the value stored for the local
-   * degree of freedom with index @p dof. If the
-   * object is vector-valued, a vector-valued
-   * return argument is given. Note that when
-   * vectorization is enabled, values from
-   * several cells are grouped together. If @p
-   * set_dof_values was called last, the value
-   * corresponds to the one set there. If @p
-   * integrate was called last, it instead
-   * corresponds to the value of the integrated
-   * function with the test function of the
-   * given index.
+   * Returns the value stored for the local degree of freedom with index @p
+   * dof. If the object is vector-valued, a vector-valued return argument is
+   * given. Note that when vectorization is enabled, values from several cells
+   * are grouped together. If @p set_dof_values was called last, the value
+   * corresponds to the one set there. If @p integrate was called last, it
+   * instead corresponds to the value of the integrated function with the test
+   * function of the given index.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   value_type get_dof_value (const unsigned int dof) const;
 
   /**
-   * Write a value to the field containing the
-   * degrees of freedom with component @p
-   * dof. Writes to the same field as is
-   * accessed through @p
-   * get_dof_value. Therefore, the original data
-   * that was read from a vector is overwritten
-   * as soon as a value is submitted.
+   * Write a value to the field containing the degrees of freedom with
+   * component @p dof. Writes to the same field as is accessed through @p
+   * get_dof_value. Therefore, the original data that was read from a vector
+   * is overwritten as soon as a value is submitted.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   void submit_dof_value (const value_type   val_in,
                          const unsigned int dof);
 
   /**
-   * Returns the value of a finite
-   * element function at quadrature
-   * point number @p q_point after a
-   * call to @p evaluate(true,...), or
-   * the value that has been stored
-   * there with a call to @p
-   * submit_value. If the object is
-   * vector-valued, a vector-valued
-   * return argument is given. Note that
-   * when vectorization is enabled,
-   * values from several cells are
-   * grouped together.
+   * Returns the value of a finite element function at quadrature point number
+   * @p q_point after a call to @p evaluate(true,...), or the value that has
+   * been stored there with a call to @p submit_value. If the object is
+   * vector-valued, a vector-valued return argument is given. Note that when
+   * vectorization is enabled, values from several cells are grouped together.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   value_type get_value (const unsigned int q_point) const;
 
   /**
-   * Write a value to the field containing the
-   * values on quadrature points with component
-   * @p q_point. Access to the same field as
-   * through @p get_value. If applied before the
-   * function @p integrate(true,...) is
-   * called, this specifies the value which is
-   * tested by all basis function on the current
-   * cell and integrated over.
+   * Write a value to the field containing the values on quadrature points
+   * with component @p q_point. Access to the same field as through @p
+   * get_value. If applied before the function @p integrate(true,...) is
+   * called, this specifies the value which is tested by all basis function on
+   * the current cell and integrated over.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   void submit_value (const value_type   val_in,
                      const unsigned int q_point);
 
   /**
-   * Returns the gradient of a finite element
-   * function at quadrature point number @p
-   * q_point after a call to @p
-   * evaluate(...,true,...), or the value
-   * that has been stored there with a call to
-   * @p submit_gradient.
+   * Returns the gradient of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true,...), or the value
+   * that has been stored there with a call to @p submit_gradient.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   gradient_type get_gradient (const unsigned int q_point) const;
 
   /**
-   * Write a contribution that is tested
-   * by the gradient to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. Access to the same field
-   * as through @p get_gradient. If
-   * applied before the function @p
-   * integrate(...,true) is called, this
-   * specifies what is tested by all
-   * basis function gradients on the
-   * current cell and integrated over.
+   * Write a contribution that is tested by the gradient to the field
+   * containing the values on quadrature points with component @p
+   * q_point. Access to the same field as through @p get_gradient. If applied
+   * before the function @p integrate(...,true) is called, this specifies what
+   * is tested by all basis function gradients on the current cell and
+   * integrated over.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   void submit_gradient(const gradient_type grad_in,
                        const unsigned int  q_point);
 
   /**
-   * Returns the Hessian of a finite element
-   * function at quadrature point number @p
-   * q_point after a call to @p
-   * evaluate(...,true). If only the
-   * diagonal or even the trace of the Hessian,
-   * the Laplacian, is needed, use the other
-   * functions below.
+   * Returns the Hessian of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true). If only the
+   * diagonal or even the trace of the Hessian, the Laplacian, is needed, use
+   * the other functions below.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >
   get_hessian (const unsigned int q_point) const;
 
   /**
-   * Returns the diagonal of the Hessian of a
-   * finite element function at quadrature point
-   * number @p q_point after a call to @p
-   * evaluate(...,true).
+   * Returns the diagonal of the Hessian of a finite element function at
+   * quadrature point number @p q_point after a call to @p evaluate(...,true).
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   gradient_type get_hessian_diagonal (const unsigned int q_point) const;
 
   /**
-   * Returns the Laplacian (i.e., the trace of
-   * the Hessian) of a finite element function
-   * at quadrature point number @p q_point after
-   * a call to @p evaluate(...,true). Compared
-   * to the case when computing the full
-   * Hessian, some operations can be saved when
-   * only the Laplacian is requested.
+   * Returns the Laplacian (i.e., the trace of the Hessian) of a finite
+   * element function at quadrature point number @p q_point after a call to @p
+   * evaluate(...,true). Compared to the case when computing the full Hessian,
+   * some operations can be saved when only the Laplacian is requested.
    *
-   * Note that the derived class
-   * FEEvaluationAccess overloads this operation
-   * with specializations for the scalar case
-   * (n_components == 1) and for the
+   * Note that the derived class FEEvaluationAccess overloads this operation
+   * with specializations for the scalar case (n_components == 1) and for the
    * vector-valued case (n_components == dim).
    */
   value_type get_laplacian (const unsigned int q_point) const;
 
   /**
-   * Takes values on quadrature points,
-   * multiplies by the Jacobian determinant and
-   * quadrature weights (JxW) and sums the
-   * values for all quadrature points on the
-   * cell. The result is a scalar, representing
-   * the integral over the function over the
-   * cell. If a vector-element is used, the
-   * resulting components are still
-   * separated. Moreover, if vectorization is
-   * enabled, the integral values of several
-   * cells are represented together.
+   * Takes values on quadrature points, multiplies by the Jacobian determinant
+   * and quadrature weights (JxW) and sums the values for all quadrature
+   * points on the cell. The result is a scalar, representing the integral
+   * over the function over the cell. If a vector-element is used, the
+   * resulting components are still separated. Moreover, if vectorization is
+   * enabled, the integral values of several cells are represented together.
    */
   value_type integrate_value () const;
 
@@ -521,122 +403,100 @@ public:
    */
   //@{
   /**
-   * Returns a read-only pointer to the first
-   * field of function values on quadrature
-   * points. First come the function values on
-   * all quadrature points for the first
-   * component, then all values for the second
-   * component, and so on. This is related to
-   * the internal data structures used in this
-   * class. The raw data after a call to @p
-   * evaluate only contains unit cell
-   * operations, so possible transformations,
-   * quadrature weights etc. must be applied
-   * manually. In general, it is safer to use
-   * the get_value() function instead, which
-   * does all the transformation internally.
+   * Returns a read-only pointer to the first field of the dof values. This is
+   * the data field the read_dof_values() functions write into. First come the
+   * the dof values for the first component, then all values for the second
+   * component, and so on. This is related to the internal data structures
+   * used in this class. In general, it is safer to use the get_dof_value()
+   * function instead.
+   */
+  const VectorizedArray<Number> *begin_dof_values () const;
+
+  /**
+   * Returns a read and write pointer to the first field of the dof
+   * values. This is the data field the read_dof_values() functions write
+   * into. First come the the dof values for the first component, then all
+   * values for the second component, and so on. This is related to the
+   * internal data structures used in this class. In general, it is safer to
+   * use the get_dof_value() function instead.
+   */
+  VectorizedArray<Number> *begin_dof_values ();
+
+  /**
+   * Returns a read-only pointer to the first field of function values on
+   * quadrature points. First come the function values on all quadrature
+   * points for the first component, then all values for the second component,
+   * and so on. This is related to the internal data structures used in this
+   * class. The raw data after a call to @p evaluate only contains unit cell
+   * operations, so possible transformations, quadrature weights etc. must be
+   * applied manually. In general, it is safer to use the get_value() function
+   * instead, which does all the transformation internally.
    */
   const VectorizedArray<Number> *begin_values () const;
 
   /**
-   * Returns a read and write pointer to the
-   * first field of function values on
-   * quadrature points. First come the function
-   * values on all quadrature points for the
-   * first component, then all values for the
-   * second component, and so on. This is
-   * related to the internal data structures
-   * used in this class. The raw data after a
-   * call to @p evaluate only contains unit
-   * cell operations, so possible
-   * transformations, quadrature weights
-   * etc. must be applied manually. In general,
-   * it is safer to use the get_value() function
-   * instead, which does all the transformation
-   * internally.
+   * Returns a read and write pointer to the first field of function values on
+   * quadrature points. First come the function values on all quadrature
+   * points for the first component, then all values for the second component,
+   * and so on. This is related to the internal data structures used in this
+   * class. The raw data after a call to @p evaluate only contains unit cell
+   * operations, so possible transformations, quadrature weights etc. must be
+   * applied manually. In general, it is safer to use the get_value() function
+   * instead, which does all the transformation internally.
    */
   VectorizedArray<Number> *begin_values ();
 
   /**
-   * Returns a read-only pointer to the first
-   * field of function gradients on quadrature
-   * points. First comes the x-component of the
-   * gradient for the first component on all
-   * quadrature points, then the y-component,
-   * and so on. Next comes the x-component of
-   * the second component, and so on. This is
-   * related to the internal data structures
-   * used in this class. The raw data after a
-   * call to @p evaluate only contains unit
-   * cell operations, so possible
-   * transformations, quadrature weights
-   * etc. must be applied manually. In general,
-   * it is safer to use the get_gradient() function
-   * instead, which does all the transformation
-   * internally.
+   * Returns a read-only pointer to the first field of function gradients on
+   * quadrature points. First comes the x-component of the gradient for the
+   * first component on all quadrature points, then the y-component, and so
+   * on. Next comes the x-component of the second component, and so on. This
+   * is related to the internal data structures used in this class. The raw
+   * data after a call to @p evaluate only contains unit cell operations, so
+   * possible transformations, quadrature weights etc. must be applied
+   * manually. In general, it is safer to use the get_gradient() function
+   * instead, which does all the transformation internally.
    */
   const VectorizedArray<Number> *begin_gradients () const;
 
   /**
-   * Returns a read and write pointer to the
-   * first field of function gradients on
-   * quadrature points. First comes the
-   * x-component of the gradient for the first
-   * component on all quadrature points, then
-   * the y-component, and so on. Next comes the
-   * x-component of the second component, and so
-   * on. This is related to the internal data
-   * structures used in this class. The raw data
-   * after a call to @p evaluate only
-   * contains unit cell operations, so possible
-   * transformations, quadrature weights
-   * etc. must be applied manually. In general,
-   * it is safer to use the get_gradient()
-   * function instead, which does all the
-   * transformation internally.
+   * Returns a read and write pointer to the first field of function gradients
+   * on quadrature points. First comes the x-component of the gradient for the
+   * first component on all quadrature points, then the y-component, and so
+   * on. Next comes the x-component of the second component, and so on. This
+   * is related to the internal data structures used in this class. The raw
+   * data after a call to @p evaluate only contains unit cell operations, so
+   * possible transformations, quadrature weights etc. must be applied
+   * manually. In general, it is safer to use the get_gradient() function
+   * instead, which does all the transformation internally.
    */
   VectorizedArray<Number> *begin_gradients ();
 
   /**
-   * Returns a read-only pointer to the first
-   * field of function hessians on quadrature
-   * points. First comes the xx-component of the
-   * hessian for the first component on all
-   * quadrature points, then the yy-component,
-   * zz-component in (3D), then the
-   * xy-component, and so on. Next comes the
-   * xx-component of the second component, and
-   * so on. This is related to the internal data
-   * structures used in this class. The raw data
-   * after a call to @p evaluate only
-   * contains unit cell operations, so possible
-   * transformations, quadrature weights
-   * etc. must be applied manually. In general,
-   * it is safer to use the get_laplacian() or
-   * get_hessian() functions instead, which does
-   * all the transformation internally.
+   * Returns a read-only pointer to the first field of function hessians on
+   * quadrature points. First comes the xx-component of the hessian for the
+   * first component on all quadrature points, then the yy-component,
+   * zz-component in (3D), then the xy-component, and so on. Next comes the
+   * xx-component of the second component, and so on. This is related to the
+   * internal data structures used in this class. The raw data after a call to
+   * @p evaluate only contains unit cell operations, so possible
+   * transformations, quadrature weights etc. must be applied manually. In
+   * general, it is safer to use the get_laplacian() or get_hessian()
+   * functions instead, which does all the transformation internally.
    */
   const VectorizedArray<Number> *begin_hessians () const;
 
   /**
-   * Returns a read and write pointer to the
-   * first field of function hessians on
-   * quadrature points. First comes the
-   * xx-component of the hessian for the first
-   * component on all quadrature points, then
-   * the yy-component, zz-component in (3D),
-   * then the xy-component, and so on. Next
-   * comes the xx-component of the second
-   * component, and so on. This is related to
-   * the internal data structures used in this
-   * class. The raw data after a call to @p
-   * evaluate only contains unit cell
-   * operations, so possible transformations,
-   * quadrature weights etc. must be applied
-   * manually. In general, it is safer to use
-   * the get_laplacian() or get_hessian()
-   * functions instead, which does all the
-   * transformation internally.
+   * Returns a read and write pointer to the first field of function hessians
+   * on quadrature points. First comes the xx-component of the hessian for the
+   * first component on all quadrature points, then the yy-component,
+   * zz-component in (3D), then the xy-component, and so on. Next comes the
+   * xx-component of the second component, and so on. This is related to the
+   * internal data structures used in this class. The raw data after a call to
+   * @p evaluate only contains unit cell operations, so possible
+   * transformations, quadrature weights etc. must be applied manually. In
+   * general, it is safer to use the get_laplacian() or get_hessian()
+   * functions instead, which does all the transformation internally.
    */
   VectorizedArray<Number> *begin_hessians ();
 
@@ -645,122 +505,94 @@ public:
 protected:
 
   /**
-   * Constructor. Made protected to prevent
-   * users from directly using this class. Takes
-   * all data stored in MatrixFree. If applied
-   * to problems with more than one finite
-   * element or more than one quadrature formula
-   * selected during construction of @p
-   * matrix_free, @p fe_no and @p quad_no allow
-   * to select the appropriate components.
+   * Constructor. Made protected to prevent users from directly using this
+   * class. Takes all data stored in MatrixFree. If applied to problems with
+   * more than one finite element or more than one quadrature formula selected
+   * during construction of @p matrix_free, @p fe_no and @p quad_no allow to
+   * select the appropriate components.
    */
   FEEvaluationBase (const MatrixFree<dim,Number> &matrix_free,
                     const unsigned int            fe_no   = 0,
                     const unsigned int            quad_no = 0);
 
   /**
-   * A unified function to read from and write
-   * into vectors based on the given template
-   * operation. It can perform the operation for
-   * @p read_dof_values, @p
-   * distribute_local_to_global, and @p
-   * set_dof_values. It performs the operation
-   * for several vectors at a time.
+   * A unified function to read from and write into vectors based on the given
+   * template operation. It can perform the operation for @p read_dof_values,
+   * @p distribute_local_to_global, and @p set_dof_values. It performs the
+   * operation for several vectors at a time.
    */
   template<typename VectorType, typename VectorOperation>
   void read_write_operation (const VectorOperation &operation,
                              VectorType            *vectors[]) const;
 
   /**
-   * For a collection of several vector @p src,
-   * read out the values on the degrees of
-   * freedom of the current cell for @p
-   * n_components (template argument), and store
-   * them internally. Similar functionality as
-   * the function
-   * DoFAccessor::read_dof_values. Note
-   * that if vectorization is enabled, the DoF
-   * values for several cells are set.
+   * For a collection of several vector @p src, read out the values on the
+   * degrees of freedom of the current cell for @p n_components (template
+   * argument), and store them internally. Similar functionality as the
+   * function DoFAccessor::read_dof_values. Note that if vectorization is
+   * enabled, the DoF values for several cells are set.
    */
   template<typename VectorType>
   void read_dof_values_plain (const VectorType *src_data[]);
 
   /**
-   * Internal data fields that store the
-   * values. Since all array lengths are known
-   * at compile time and since they are rarely
-   * more than a few kilobytes, allocate them on
-   * the stack. This makes it possible to
-   * cheaply set up a FEEvaluation object and
-   * write thread-safe programs by letting each
-   * thread own a private object of this type.
+   * Internal data fields that store the values. Since all array lengths are
+   * known at compile time and since they are rarely more than a few
+   * kilobytes, allocate them on the stack. This makes it possible to cheaply
+   * set up a FEEvaluation object and write thread-safe programs by letting
+   * each thread own a private object of this type.
    *
-   * This field stores the values for local
-   * degrees of freedom (e.g. after reading out
-   * from a vector but before applying unit cell
-   * transformations or before distributing them
-   * into a result vector). The methods
-   * get_dof_value() and submit_dof_value()
-   * read from or write to this field.
+   * This field stores the values for local degrees of freedom (e.g. after
+   * reading out from a vector but before applying unit cell transformations
+   * or before distributing them into a result vector). The methods
+   * get_dof_value() and submit_dof_value() read from or write to this field.
    */
   VectorizedArray<Number> values_dofs[n_components][dofs_per_cell>0?dofs_per_cell:1];
 
   /**
-   * This field stores the values of the finite
-   * element function on quadrature points after
-   * applying unit cell transformations or
-   * before integrating. The methods get_value()
-   * and submit_value() access this field.
+   * This field stores the values of the finite element function on quadrature
+   * points after applying unit cell transformations or before
+   * integrating. The methods get_value() and submit_value() access this
+   * field.
    */
   VectorizedArray<Number> values_quad[n_components][n_q_points>0?n_q_points:1];
 
   /**
-   * This field stores the gradients of the
-   * finite element function on quadrature
-   * points after applying unit cell
-   * transformations or before integrating. The
-   * methods get_gradient() and
-   * submit_gradient() (as well as some
-   * specializations like
-   * get_symmetric_gradient() or
-   * get_divergence()) access this field.
+   * This field stores the gradients of the finite element function on
+   * quadrature points after applying unit cell transformations or before
+   * integrating. The methods get_gradient() and submit_gradient() (as well as
+   * some specializations like get_symmetric_gradient() or get_divergence())
+   * access this field.
    */
   VectorizedArray<Number> gradients_quad[n_components][dim][n_q_points>0?n_q_points:1];
 
   /**
-   * This field stores the Hessians of the
-   * finite element function on quadrature
-   * points after applying unit cell
-   * transformations. The methods get_hessian(),
-   * get_laplacian(), get_hessian_diagonal()
-   * access this field.
+   * This field stores the Hessians of the finite element function on
+   * quadrature points after applying unit cell transformations. The methods
+   * get_hessian(), get_laplacian(), get_hessian_diagonal() access this field.
    */
   VectorizedArray<Number> hessians_quad[n_components][(dim*(dim+1))/2][n_q_points>0?n_q_points:1];
 
   /**
-   * Stores the number of the quadrature formula
-   * of the present cell.
+   * Stores the number of the quadrature formula of the present cell.
    */
   const unsigned int quad_no;
 
   /**
-   * Stores the number of components in the
-   * finite element as detected in the
-   * MatrixFree storage class for comparison
-   * with the template argument.
+   * Stores the number of components in the finite element as detected in the
+   * MatrixFree storage class for comparison with the template argument.
    */
   const unsigned int n_fe_components;
 
   /**
-   * Stores the active fe index for this class
-   * for efficient indexing in the hp case.
+   * Stores the active fe index for this class for efficient indexing in the
+   * hp case.
    */
   const unsigned int active_fe_index;
 
   /**
-   * Stores the active quadrature index for this
-   * class for efficient indexing in the hp
-   * case.
+   * Stores the active quadrature index for this class for efficient indexing
+   * in the hp case.
    */
   const unsigned int active_quad_index;
 
@@ -770,183 +602,146 @@ protected:
   const MatrixFree<dim,Number>   &matrix_info;
 
   /**
-   * Stores a reference to the underlying DoF
-   * indices and constraint description for the
-   * component specified at construction. Also
-   * contained in matrix_info, but it simplifies
-   * code if we store a reference to it.
+   * Stores a reference to the underlying DoF indices and constraint
+   * description for the component specified at construction. Also contained
+   * in matrix_info, but it simplifies code if we store a reference to it.
    */
   const internal::MatrixFreeFunctions::DoFInfo      &dof_info;
 
   /**
-   * Stores a reference to the underlying
-   * transformation data from unit to real cells
-   * for the given quadrature formula specified
-   * at construction.  Also contained in
-   * matrix_info, but it simplifies code if we
-   * store a reference to it.
+   * Stores a reference to the underlying transformation data from unit to
+   * real cells for the given quadrature formula specified at construction.
+   * Also contained in matrix_info, but it simplifies code if we store a
+   * reference to it.
    */
   const internal::MatrixFreeFunctions::MappingInfo<dim,Number> &mapping_info;
 
   /**
-   * Stores a reference to the unit cell data,
-   * i.e., values, gradients and Hessians in 1D
-   * at the quadrature points that constitute
-   * the tensor product. Also contained in
-   * matrix_info, but it simplifies code if we
+   * Stores a reference to the unit cell data, i.e., values, gradients and
+   * Hessians in 1D at the quadrature points that constitute the tensor
+   * product. Also contained in matrix_info, but it simplifies code if we
    * store a reference to it.
    */
   const internal::MatrixFreeFunctions::ShapeInfo<Number> &data;
 
   /**
-   * A pointer to the Cartesian Jacobian
-   * information of the present cell. Only set
-   * to a useful value if on a Cartesian cell,
-   * otherwise zero.
+   * A pointer to the Cartesian Jacobian information of the present cell. Only
+   * set to a useful value if on a Cartesian cell, otherwise zero.
    */
   const Tensor<1,dim,VectorizedArray<Number> > *cartesian_data;
 
   /**
-   * A pointer to the Jacobian information of
-   * the present cell. Only set to a useful
-   * value if on a non-Cartesian cell.
+   * A pointer to the Jacobian information of the present cell. Only set to a
+   * useful value if on a non-Cartesian cell.
    */
   const Tensor<2,dim,VectorizedArray<Number> > *jacobian;
 
   /**
-   * A pointer to the Jacobian determinant of
-   * the present cell. If on a Cartesian cell or
-   * on a cell with constant Jacobian, this is
-   * just the Jacobian determinant, otherwise
-   * the Jacobian determinant times the
+   * A pointer to the Jacobian determinant of the present cell. If on a
+   * Cartesian cell or on a cell with constant Jacobian, this is just the
+   * Jacobian determinant, otherwise the Jacobian determinant times the
    * quadrature weight.
    */
   const VectorizedArray<Number> *J_value;
 
   /**
-   * A pointer to the quadrature weights of the
-   * underlying quadrature formula.
+   * A pointer to the quadrature weights of the underlying quadrature formula.
    */
   const VectorizedArray<Number> *quadrature_weights;
 
   /**
-   * A pointer to the quadrature points on the
-   * present cell.
+   * A pointer to the quadrature points on the present cell.
    */
   const Point<dim,VectorizedArray<Number> > *quadrature_points;
 
   /**
-   * A pointer to the diagonal part of the
-   * Jacobian gradient on the present
-   * cell. Only set to a useful value if on a
-   * general cell with non-constant Jacobian.
+   * A pointer to the diagonal part of the Jacobian gradient on the present
+   * cell. Only set to a useful value if on a general cell with non-constant
+   * Jacobian.
    */
   const Tensor<2,dim,VectorizedArray<Number> > *jacobian_grad;
 
   /**
-   * A pointer to the upper diagonal part of the
-   * Jacobian gradient on the present cell. Only
-   * set to a useful value if on a general cell
-   * with non-constant Jacobian.
+   * A pointer to the upper diagonal part of the Jacobian gradient on the
+   * present cell. Only set to a useful value if on a general cell with
+   * non-constant Jacobian.
    */
   const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > * jacobian_grad_upper;
 
   /**
-   * After a call to reinit(), stores the number
-   * of the cell we are currently working with.
+   * After a call to reinit(), stores the number of the cell we are currently
+   * working with.
    */
   unsigned int cell;
 
   /**
-   * Stores the type of the cell we are
-   * currently working with after a call to
-   * reinit(). Valid values are @p cartesian, @p
-   * affine and @p general, which have different
-   * implications on how the Jacobian
-   * transformations are stored internally in
-   * MappingInfo.
+   * Stores the type of the cell we are currently working with after a call to
+   * reinit(). Valid values are @p cartesian, @p affine and @p general, which
+   * have different implications on how the Jacobian transformations are
+   * stored internally in MappingInfo.
    */
   internal::MatrixFreeFunctions::CellType cell_type;
 
   /**
-   * The stride to access the correct data in
-   * MappingInfo.
+   * The stride to access the correct data in MappingInfo.
    */
   unsigned int cell_data_number;
 
   /**
-   * If the present cell chunk for vectorization
-   * is not completely filled up with data, this
-   * field stores how many physical cells are
-   * underlying. Is between 1 and
-   * VectorizedArray<Number>::n_array_elements-1
-   * (inclusive).
+   * If the present cell chunk for vectorization is not completely filled up
+   * with data, this field stores how many physical cells are underlying. Is
+   * between 1 and VectorizedArray<Number>::n_array_elements-1 (inclusive).
    */
   unsigned int n_irreg_components_filled;
 
   /**
-   * Stores whether the present cell chunk used
-   * in vectorization is not completely filled
-   * up with physical cells. E.g. if
-   * vectorization dictates that four cells
-   * should be worked with but only three
-   * physical cells are left, this flag will be
-   * set to true, otherwise to false. Mainly
-   * used for internal checking when reading
-   * from vectors or writing to vectors.
+   * Stores whether the present cell chunk used in vectorization is not
+   * completely filled up with physical cells. E.g. if vectorization dictates
+   * that four cells should be worked with but only three physical cells are
+   * left, this flag will be set to true, otherwise to false. Mainly used for
+   * internal checking when reading from vectors or writing to vectors.
    */
   bool at_irregular_cell;
 
   /**
-   * Debug information to track whether dof
-   * values have been initialized before
-   * accessed. Used to control exceptions when
-   * uninitialized data is used.
+   * Debug information to track whether dof values have been initialized
+   * before accessed. Used to control exceptions when uninitialized data is
+   * used.
    */
   bool dof_values_initialized;
 
   /**
-   * Debug information to track whether values
-   * on quadrature points have been initialized
-   * before accessed. Used to control exceptions
-   * when uninitialized data is used.
+   * Debug information to track whether values on quadrature points have been
+   * initialized before accessed. Used to control exceptions when
+   * uninitialized data is used.
    */
   bool values_quad_initialized;
 
   /**
-   * Debug information to track whether
-   * gradients on quadrature points have been
-   * initialized before accessed. Used to
-   * control exceptions when uninitialized data
-   * is used.
+   * Debug information to track whether gradients on quadrature points have
+   * been initialized before accessed. Used to control exceptions when
+   * uninitialized data is used.
    */
   bool gradients_quad_initialized;
 
   /**
-   * Debug information to track whether
-   * Hessians on quadrature points have been
-   * initialized before accessed. Used to
-   * control exceptions when uninitialized data
-   * is used.
+   * Debug information to track whether Hessians on quadrature points have
+   * been initialized before accessed. Used to control exceptions when
+   * uninitialized data is used.
    */
   bool hessians_quad_initialized;
 
   /**
-   * Debug information to track whether values
-   * on quadrature points have been submitted
-   * for integration before the integration is
-   * actually stared. Used to control exceptions
-   * when uninitialized data is used.
+   * Debug information to track whether values on quadrature points have been
+   * submitted for integration before the integration is actually stared. Used
+   * to control exceptions when uninitialized data is used.
    */
   bool values_quad_submitted;
 
   /**
-   * Debug information to track whether
-   * gradients on quadrature points have been
-   * submitted for integration before the
-   * integration is actually stared. Used to
-   * control exceptions when uninitialized data
-   * is used.
+   * Debug information to track whether gradients on quadrature points have
+   * been submitted for integration before the integration is actually
+   * stared. Used to control exceptions when uninitialized data is used.
    */
   bool gradients_quad_submitted;
 };
@@ -979,14 +774,11 @@ public:
 
 protected:
   /**
-   * Constructor. Made protected to prevent
-   * initialization in user code. Takes all data
-   * stored in MatrixFree. If applied to
-   * problems with more than one finite element
-   * or more than one quadrature formula
-   * selected during construction of @p
-   * matrix_free, @p fe_no and @p quad_no allow
-   * to select the appropriate components.
+   * Constructor. Made protected to prevent initialization in user code. Takes
+   * all data stored in MatrixFree. If applied to problems with more than one
+   * finite element or more than one quadrature formula selected during
+   * construction of @p matrix_free, @p fe_no and @p quad_no allow to select
+   * the appropriate components.
    */
   FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
                       const unsigned int            fe_no   = 0,
@@ -1018,141 +810,98 @@ public:
   typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,1,Number> BaseClass;
 
   /**
-   * Returns the value stored for the local
-   * degree of freedom with index @p dof. If the
-   * object is vector-valued, a vector-valued
-   * return argument is given. Note that when
-   * vectorization is enabled, values from
-   * several cells are grouped together. If @p
-   * set_dof_values was called last, the value
-   * corresponds to the one set there. If @p
-   * integrate was called last, it instead
-   * corresponds to the value of the integrated
-   * function with the test function of the
-   * given index.
+   * Returns the value stored for the local degree of freedom with index @p
+   * dof. If the object is vector-valued, a vector-valued return argument is
+   * given. Note that when vectorization is enabled, values from several cells
+   * are grouped together. If @p set_dof_values was called last, the value
+   * corresponds to the one set there. If @p integrate was called last, it
+   * instead corresponds to the value of the integrated function with the test
+   * function of the given index.
    */
   value_type get_dof_value (const unsigned int dof) const;
 
   /**
-   * Write a value to the field containing the
-   * degrees of freedom with component @p
-   * dof. Access to the same field as through @p
-   * get_dof_value.
+   * Write a value to the field containing the degrees of freedom with
+   * component @p dof. Access to the same field as through @p get_dof_value.
    */
   void submit_dof_value (const value_type   val_in,
                          const unsigned int dof);
 
   /**
-   * Returns the value of a finite element
-   * function at quadrature point number @p
-   * q_point after a call to @p
-   * evaluate(true,...), or the value that
-   * has been stored there with a call to @p
-   * submit_value. If the object is
-   * vector-valued, a vector-valued return
-   * argument is given. Note that when
-   * vectorization is enabled, values from
-   * several cells are grouped together.
+   * Returns the value of a finite element function at quadrature point number
+   * @p q_point after a call to @p evaluate(true,...), or the value that has
+   * been stored there with a call to @p submit_value. If the object is
+   * vector-valued, a vector-valued return argument is given. Note that when
+   * vectorization is enabled, values from several cells are grouped together.
    */
   value_type get_value (const unsigned int q_point) const;
 
   /**
-   * Write a value to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. Access to the same field
-   * as through @p get_value. If applied
-   * before the function @p
-   * integrate(true,...) is called, this
-   * specifies the value which is tested
-   * by all basis function on the
-   * current cell and integrated over.
+   * Write a value to the field containing the values on quadrature points
+   * with component @p q_point. Access to the same field as through @p
+   * get_value. If applied before the function @p integrate(true,...) is
+   * called, this specifies the value which is tested by all basis function on
+   * the current cell and integrated over.
    */
   void submit_value (const value_type   val_in,
                      const unsigned int q_point);
 
   /**
-   * Returns the gradient of a finite
-   * element function at quadrature
-   * point number @p q_point after a
-   * call to @p evaluate(...,true,...),
-   * or the value that has been stored
-   * there with a call to @p
-   * submit_gradient.
+   * Returns the gradient of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true,...), or the value
+   * that has been stored there with a call to @p submit_gradient.
    */
   gradient_type get_gradient (const unsigned int q_point) const;
 
   /**
-   * Write a contribution that is tested
-   * by the gradient to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. Access to the same field
-   * as through @p get_gradient. If
-   * applied before the function @p
-   * integrate(...,true) is called, this
-   * specifies what is tested by all
-   * basis function gradients on the
-   * current cell and integrated over.
+   * Write a contribution that is tested by the gradient to the field
+   * containing the values on quadrature points with component @p
+   * q_point. Access to the same field as through @p get_gradient. If applied
+   * before the function @p integrate(...,true) is called, this specifies what
+   * is tested by all basis function gradients on the current cell and
+   * integrated over.
    */
   void submit_gradient(const gradient_type grad_in,
                        const unsigned int  q_point);
 
   /**
-   * Returns the Hessian of a finite
-   * element function at quadrature
-   * point number @p q_point after a
-   * call to @p evaluate(...,true). If
-   * only the diagonal part of the
-   * Hessian or its trace, the
-   * Laplacian, are needed, use the
-   * respective functions below.
+   * Returns the Hessian of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true). If only the
+   * diagonal part of the Hessian or its trace, the Laplacian, are needed, use
+   * the respective functions below.
    */
   Tensor<2,dim,VectorizedArray<Number> >
   get_hessian (unsigned int q_point) const;
 
   /**
-   * Returns the diagonal of the Hessian
-   * of a finite element function at
-   * quadrature point number @p q_point
-   * after a call to @p
-   * evaluate(...,true).
+   * Returns the diagonal of the Hessian of a finite element function at
+   * quadrature point number @p q_point after a call to @p evaluate(...,true).
    */
   gradient_type get_hessian_diagonal (const unsigned int q_point) const;
 
   /**
-   * Returns the Laplacian of a finite
-   * element function at quadrature
-   * point number @p q_point after a
-   * call to @p evaluate(...,true).
+   * Returns the Laplacian of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true).
    */
   value_type get_laplacian (const unsigned int q_point) const;
 
   /**
-   * Takes values on quadrature points,
-   * multiplies by the Jacobian determinant and
-   * quadrature weights (JxW) and sums the
-   * values for all quadrature points on the
-   * cell. The result is a scalar, representing
-   * the integral over the function over the
-   * cell. If a vector-element is used, the
-   * resulting components are still
-   * separated. Moreover, if vectorization is
-   * enabled, the integral values of several
-   * cells are represented together.
+   * Takes values on quadrature points, multiplies by the Jacobian determinant
+   * and quadrature weights (JxW) and sums the values for all quadrature
+   * points on the cell. The result is a scalar, representing the integral
+   * over the function over the cell. If a vector-element is used, the
+   * resulting components are still separated. Moreover, if vectorization is
+   * enabled, the integral values of several cells are represented together.
    */
   value_type integrate_value () const;
 
 protected:
   /**
-   * Constructor. Made protected to avoid
-   * initialization in user code. Takes all data
-   * stored in MatrixFree. If applied to
-   * problems with more than one finite element
-   * or more than one quadrature formula
-   * selected during construction of @p
-   * matrix_free, @p fe_no and @p quad_no allow
-   * to select the appropriate components.
+   * Constructor. Made protected to avoid initialization in user code. Takes
+   * all data stored in MatrixFree. If applied to problems with more than one
+   * finite element or more than one quadrature formula selected during
+   * construction of @p matrix_free, @p fe_no and @p quad_no allow to select
+   * the appropriate components.
    */
   FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
                       const unsigned int            fe_no   = 0,
@@ -1185,150 +934,106 @@ public:
   typedef FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,dim,Number> BaseClass;
 
   /**
-   * Returns the gradient of a finite element
-   * function at quadrature point number @p
-   * q_point after a call to @p
-   * evaluate(...,true,...).
+   * Returns the gradient of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true,...).
    */
   gradient_type get_gradient (const unsigned int q_point) const;
 
   /**
-   * Returns the divergence of a vector-valued
-   * finite element at quadrature point number
-   * @p q_point after a call to @p
-   * evaluate(...,true,...).
+   * Returns the divergence of a vector-valued finite element at quadrature
+   * point number @p q_point after a call to @p evaluate(...,true,...).
    */
   VectorizedArray<Number> get_divergence (const unsigned int q_point) const;
 
   /**
-   * Returns the symmetric gradient of a
-   * vector-valued finite element at
-   * quadrature point number @p q_point
-   * after a call to @p
-   * evaluate(...,true,...). It
-   * corresponds to <tt>0.5
+   * Returns the symmetric gradient of a vector-valued finite element at
+   * quadrature point number @p q_point after a call to @p
+   * evaluate(...,true,...). It corresponds to <tt>0.5
    * (grad+grad<sup>T</sup>)</tt>.
    */
   SymmetricTensor<2,dim,VectorizedArray<Number> >
   get_symmetric_gradient (const unsigned int q_point) const;
 
   /**
-   * Returns the curl of the vector field,
-   * $nabla \times v$ after a call to @p
+   * Returns the curl of the vector field, $nabla \times v$ after a call to @p
    * evaluate(...,true,...).
    */
   Tensor<1,dim==2?1:dim,VectorizedArray<Number> >
   get_curl (const unsigned int q_point) const;
 
   /**
-   * Returns the Hessian of a finite
-   * element function at quadrature
-   * point number @p q_point after a
-   * call to @p evaluate(...,true). If
-   * only the diagonal of the Hessian or
-   * its trace, the Laplacian, is
-   * needed, use the respective
-   * functions.
+   * Returns the Hessian of a finite element function at quadrature point
+   * number @p q_point after a call to @p evaluate(...,true). If only the
+   * diagonal of the Hessian or its trace, the Laplacian, is needed, use the
+   * respective functions.
    */
   Tensor<3,dim,VectorizedArray<Number> >
   get_hessian (const unsigned int q_point) const;
 
   /**
-   * Returns the diagonal of the Hessian
-   * of a finite element function at
-   * quadrature point number @p q_point
-   * after a call to @p
-   * evaluate(...,true).
+   * Returns the diagonal of the Hessian of a finite element function at
+   * quadrature point number @p q_point after a call to @p evaluate(...,true).
    */
   gradient_type get_hessian_diagonal (const unsigned int q_point) const;
 
   /**
-   * Write a contribution that is tested
-   * by the gradient to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. Access to the same field
-   * as through @p get_gradient. If
-   * applied before the function @p
-   * integrate(...,true) is called, this
-   * specifies what is tested by all
-   * basis function gradients on the
-   * current cell and integrated over.
+   * Write a contribution that is tested by the gradient to the field
+   * containing the values on quadrature points with component @p
+   * q_point. Access to the same field as through @p get_gradient. If applied
+   * before the function @p integrate(...,true) is called, this specifies what
+   * is tested by all basis function gradients on the current cell and
+   * integrated over.
    */
   void submit_gradient(const gradient_type grad_in,
                        const unsigned int  q_point);
 
   /**
-   * Write a contribution that is tested
-   * by the gradient to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. This function is an
-   * alternative to the other
-   * submit_gradient function when using
-   * a system of fixed number of
-   * equations which happens to coincide
-   * with the dimension for some
-   * dimensions, but not all. To allow
-   * for dimension-independent
-   * programming, this function can be
-   * used instead.
+   * Write a contribution that is tested by the gradient to the field
+   * containing the values on quadrature points with component @p
+   * q_point. This function is an alternative to the other submit_gradient
+   * function when using a system of fixed number of equations which happens
+   * to coincide with the dimension for some dimensions, but not all. To allow
+   * for dimension-independent programming, this function can be used instead.
    */
   void submit_gradient(const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > > grad_in,
                        const unsigned int q_point);
 
   /**
-   * Write a constribution that is
-   * tested by the divergence to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. Access to the same field
-   * as through @p get_gradient. If
-   * applied before the function @p
-   * integrate(...,true) is called, this
-   * specifies what is tested by all
-   * basis function gradients on the
-   * current cell and integrated over.
+   * Write a constribution that is tested by the divergence to the field
+   * containing the values on quadrature points with component @p
+   * q_point. Access to the same field as through @p get_gradient. If applied
+   * before the function @p integrate(...,true) is called, this specifies what
+   * is tested by all basis function gradients on the current cell and
+   * integrated over.
    */
   void submit_divergence (const VectorizedArray<Number> div_in,
                           const unsigned int q_point);
 
   /**
-   * Write a contribution that is tested
-   * by the gradient to the field
-   * containing the values on quadrature
-   * points with component @p
-   * q_point. Access to the same field
-   * as through @p get_gradient. If
-   * applied before the function @p
-   * integrate(...,true) is called, this
-   * specifies the gradient which is
-   * tested by all basis function
-   * gradients on the current cell and
-   * integrated over.
+   * Write a contribution that is tested by the gradient to the field
+   * containing the values on quadrature points with component @p
+   * q_point. Access to the same field as through @p get_gradient. If applied
+   * before the function @p integrate(...,true) is called, this specifies the
+   * gradient which is tested by all basis function gradients on the current
+   * cell and integrated over.
    */
   void submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray<Number> > grad_in,
                                  const unsigned int      q_point);
 
   /**
-   * Write the components of a curl containing
-   * the values on quadrature point @p
-   * q_point. Access to the same data field as
-   * through @p get_gradient.
+   * Write the components of a curl containing the values on quadrature point
+   * @p q_point. Access to the same data field as through @p get_gradient.
    */
   void submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl_in,
                     const unsigned int q_point);
 
 protected:
   /**
-   * Constructor. Made protected to avoid
-   * initialization in user code. Takes all data
-   * stored in MatrixFree. If applied to
-   * problems with more than one finite element
-   * or more than one quadrature formula
-   * selected during construction of @p
-   * matrix_free, @p fe_no and @p quad_no allow
-   * to select the appropriate components.
+   * Constructor. Made protected to avoid initialization in user code. Takes
+   * all data stored in MatrixFree. If applied to problems with more than one
+   * finite element or more than one quadrature formula selected during
+   * construction of @p matrix_free, @p fe_no and @p quad_no allow to select
+   * the appropriate components.
    */
   FEEvaluationAccess (const MatrixFree<dim,Number> &matrix_free,
                       const unsigned int            fe_no   = 0,
@@ -1398,49 +1103,38 @@ public:
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
   /**
-   * Constructor. Takes all data stored in
-   * MatrixFree. If applied to problems with
-   * more than one finite element or more than
-   * one quadrature formula selected during
-   * construction of @p matrix_free, @p
-   * fe_no and @p quad_no allow to select the
-   * appropriate components.
+   * Constructor. Takes all data stored in MatrixFree. If applied to problems
+   * with more than one finite element or more than one quadrature formula
+   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
+   * allow to select the appropriate components.
    */
   FEEvaluationGeneral (const MatrixFree<dim,Number> &matrix_free,
                        const unsigned int            fe_no   = 0,
                        const unsigned int            quad_no = 0);
 
   /**
-   * Evaluates the function values, the
-   * gradients, and the Laplacians of the FE
-   * function given at the DoF values in the
-   * input vector at the quadrature points.  The
-   * function arguments specify which parts
-   * shall actually be computed. Needs to be
-   * called before the functions @p get_value(),
-   * @p get_gradient() or @p get_laplacian
-   * return useful information.
+   * Evaluates the function values, the gradients, and the Laplacians of the
+   * FE function given at the DoF values in the input vector at the quadrature
+   * points.  The function arguments specify which parts shall actually be
+   * computed. Needs to be called before the functions @p get_value(), @p
+   * get_gradient() or @p get_laplacian return useful information.
    */
   void evaluate (const bool evaluate_val,
                  const bool evaluate_grad,
                  const bool evaluate_hess = false);
 
   /**
-   * This function takes the values and/or
-   * gradients that are stored on quadrature
-   * points, tests them by all the basis
-   * functions/gradients on the cell and
-   * performs the cell integration. The two
-   * function arguments @p integrate_val and @p
-   * integrate_grad are used to enable/disable
-   * some of values or gradients.
+   * This function takes the values and/or gradients that are stored on
+   * quadrature points, tests them by all the basis functions/gradients on the
+   * cell and performs the cell integration. The two function arguments @p
+   * integrate_val and @p integrate_grad are used to enable/disable some of
+   * values or gradients.
    */
   void integrate (const bool integrate_val,
                   const bool integrate_grad);
 
   /**
-   * Returns the q-th quadrature point stored in
-   * MappingInfo.
+   * Returns the q-th quadrature point stored in MappingInfo.
    */
   Point<dim,VectorizedArray<Number> >
   quadrature_point (const unsigned int q_point) const;
@@ -1448,51 +1142,37 @@ public:
 protected:
 
   /**
-   * Internal function that applies the function
-   * values of the tensor product in a given
-   * coordinate direction (first template
-   * argument), from polynomials to values on
-   * quadrature points (second flag set to true)
-   * or in an integration loop from values on
-   * quadrature points to values tested by
-   * different test function (second flag set to
-   * false), and if the result is to be added to
-   * previous content in the data fields or
-   * not.
+   * Internal function that applies the function values of the tensor product
+   * in a given coordinate direction (first template argument), from
+   * polynomials to values on quadrature points (second flag set to true) or
+   * in an integration loop from values on quadrature points to values tested
+   * by different test function (second flag set to false), and if the result
+   * is to be added to previous content in the data fields or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_values (const VectorizedArray<Number> in [],
                      VectorizedArray<Number> out []);
 
   /**
-   * Internal function that applies the gradient
-   * operation of the tensor product in a given
-   * coordinate direction (first template
-   * argument), from polynomials to values on
-   * quadrature points (second flag set to true)
-   * or in an integration loop from values on
-   * quadrature points to values tested by
-   * different test function (second flag set to
-   * false), and if the result is to be added to
-   * previous content in the data fields or
-   * not.
+   * Internal function that applies the gradient operation of the tensor
+   * product in a given coordinate direction (first template argument), from
+   * polynomials to values on quadrature points (second flag set to true) or
+   * in an integration loop from values on quadrature points to values tested
+   * by different test function (second flag set to false), and if the result
+   * is to be added to previous content in the data fields or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_gradients (const VectorizedArray<Number> in [],
                         VectorizedArray<Number> out []);
 
   /**
-   * Internal function that applies the second
-   * derivative operation (Hessian) of the
-   * tensor product in a given coordinate
-   * direction (first template argument), from
-   * polynomials to values on quadrature points
-   * (second flag set to true) or in an
-   * integration loop from values on quadrature
-   * points to values tested by different test
-   * function (second flag set to false), and if
-   * the result is to be added to previous
-   * content in the data fields or not.
+   * Internal function that applies the second derivative operation (Hessian)
+   * of the tensor product in a given coordinate direction (first template
+   * argument), from polynomials to values on quadrature points (second flag
+   * set to true) or in an integration loop from values on quadrature points
+   * to values tested by different test function (second flag set to false),
+   * and if the result is to be added to previous content in the data fields
+   * or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_hessians (const VectorizedArray<Number> in [],
@@ -1566,44 +1246,33 @@ public:
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
   /**
-   * Constructor. Takes all data stored in
-   * MatrixFree. If applied to problems with
-   * more than one finite element or more than
-   * one quadrature formula selected during
-   * construction of @p matrix_free, @p
-   * fe_no and @p quad_no allow to select the
-   * appropriate components.
+   * Constructor. Takes all data stored in MatrixFree. If applied to problems
+   * with more than one finite element or more than one quadrature formula
+   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
+   * allow to select the appropriate components.
    */
   FEEvaluation (const MatrixFree<dim,Number> &matrix_free,
                 const unsigned int            fe_no   = 0,
                 const unsigned int            quad_no = 0);
 
   /**
-   * Evaluates the function values, the
-   * gradients, and the Laplacians of the FE
-   * function given at the DoF values in the
-   * input vector at the quadrature points on
-   * the unit cell.  The function arguments
-   * specify which parts shall actually be
-   * computed. Needs to be called before the
-   * functions @p get_value(), @p get_gradient()
-   * or @p get_laplacian give useful information
-   * (unless these values have been set
-   * manually).
+   * Evaluates the function values, the gradients, and the Laplacians of the
+   * FE function given at the DoF values in the input vector at the quadrature
+   * points on the unit cell.  The function arguments specify which parts
+   * shall actually be computed. Needs to be called before the functions @p
+   * get_value(), @p get_gradient() or @p get_laplacian give useful
+   * information (unless these values have been set manually).
    */
   void evaluate (const bool evaluate_val,
                  const bool evaluate_grad,
                  const bool evaluate_hess = false);
 
   /**
-   * This function takes the values and/or
-   * gradients that are stored on quadrature
-   * points, tests them by all the basis
-   * functions/gradients on the cell and
-   * performs the cell integration. The two
-   * function arguments @p integrate_val and @p
-   * integrate_grad are used to enable/disable
-   * some of values or gradients.
+   * This function takes the values and/or gradients that are stored on
+   * quadrature points, tests them by all the basis functions/gradients on the
+   * cell and performs the cell integration. The two function arguments @p
+   * integrate_val and @p integrate_grad are used to enable/disable some of
+   * values or gradients.
    */
   void integrate (const bool integrate_val,
                   const bool integrate_grad);
@@ -1611,51 +1280,37 @@ public:
 protected:
 
   /**
-   * Internal function that applies the function
-   * values of the tensor product in a given
-   * coordinate direction (first template
-   * argument), from polynomials to values on
-   * quadrature points (second flag set to true)
-   * or in an integration loop from values on
-   * quadrature points to values tested by
-   * different test function (second flag set to
-   * false), and if the result is to be added to
-   * previous content in the data fields or
-   * not.
+   * Internal function that applies the function values of the tensor product
+   * in a given coordinate direction (first template argument), from
+   * polynomials to values on quadrature points (second flag set to true) or
+   * in an integration loop from values on quadrature points to values tested
+   * by different test function (second flag set to false), and if the result
+   * is to be added to previous content in the data fields or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_values (const VectorizedArray<Number> in [],
                      VectorizedArray<Number> out []);
 
   /**
-   * Internal function that applies the gradient
-   * operation of the tensor product in a given
-   * coordinate direction (first template
-   * argument), from polynomials to values on
-   * quadrature points (second flag set to true)
-   * or in an integration loop from values on
-   * quadrature points to values tested by
-   * different test function (second flag set to
-   * false), and if the result is to be added to
-   * previous content in the data fields or
-   * not.
+   * Internal function that applies the gradient operation of the tensor
+   * product in a given coordinate direction (first template argument), from
+   * polynomials to values on quadrature points (second flag set to true) or
+   * in an integration loop from values on quadrature points to values tested
+   * by different test function (second flag set to false), and if the result
+   * is to be added to previous content in the data fields or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_gradients (const VectorizedArray<Number> in [],
                         VectorizedArray<Number> out []);
 
   /**
-   * Internal function that applies the second
-   * derivative operation (Hessian) of the
-   * tensor product in a given coordinate
-   * direction (first template argument), from
-   * polynomials to values on quadrature points
-   * (second flag set to true) or in an
-   * integration loop from values on quadrature
-   * points to values tested by different test
-   * function (second flag set to false), and if
-   * the result is to be added to previous
-   * content in the data fields or not.
+   * Internal function that applies the second derivative operation (Hessian)
+   * of the tensor product in a given coordinate direction (first template
+   * argument), from polynomials to values on quadrature points (second flag
+   * set to true) or in an integration loop from values on quadrature points
+   * to values tested by different test function (second flag set to false),
+   * and if the result is to be added to previous content in the data fields
+   * or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_hessians (const VectorizedArray<Number> in [],
@@ -1724,60 +1379,45 @@ public:
   static const unsigned int n_q_points    = BaseClass::n_q_points;
 
   /**
-   * Constructor. Takes all data stored in
-   * MatrixFree. If applied to problems with
-   * more than one finite element or more than
-   * one quadrature formula selected during
-   * construction of @p matrix_free, @p
-   * fe_no and @p quad_no allow to select the
-   * appropriate components.
+   * Constructor. Takes all data stored in MatrixFree. If applied to problems
+   * with more than one finite element or more than one quadrature formula
+   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
+   * allow to select the appropriate components.
    */
   FEEvaluationGL (const MatrixFree<dim,Number> &matrix_free,
                   const unsigned int          fe_no   = 0,
                   const unsigned int          quad_no = 0);
 
   /**
-   * Evaluates the function values, the
-   * gradients, and the Hessians of the FE
-   * function given at the DoF values in the
-   * input vector at the quadrature points of
-   * the unit cell. The function arguments
-   * specify which parts shall actually be
-   * computed. Needs to be called before the
-   * functions @p get_value(), @p get_gradient()
-   * or @p get_laplacian give useful information
-   * (unless these values have been set
-   * manually).
+   * Evaluates the function values, the gradients, and the Hessians of the FE
+   * function given at the DoF values in the input vector at the quadrature
+   * points of the unit cell. The function arguments specify which parts shall
+   * actually be computed. Needs to be called before the functions @p
+   * get_value(), @p get_gradient() or @p get_laplacian give useful
+   * information (unless these values have been set manually).
    */
   void evaluate (const bool evaluate_val,
                  const bool evaluate_grad,
                  const bool evaluate_lapl = false);
 
   /**
-   * This function takes the values and/or
-   * gradients that are stored on quadrature
-   * points, tests them by all the basis
-   * functions/gradients on the cell and
-   * performs the cell integration. The two
-   * function arguments @p integrate_val and @p
-   * integrate_grad are used to enable/disable
-   * some of values or gradients.
+   * This function takes the values and/or gradients that are stored on
+   * quadrature points, tests them by all the basis functions/gradients on the
+   * cell and performs the cell integration. The two function arguments @p
+   * integrate_val and @p integrate_grad are used to enable/disable some of
+   * values or gradients.
    */
   void integrate (const bool integrate_val,
                   const bool integrate_grad);
 
 protected:
   /**
-   * Internal function that applies the gradient
-   * operation of the tensor product in a given
-   * coordinate direction (first template
-   * argument), from polynomials to values on
-   * quadrature points (second flag set to true)
-   * or in an integration loop from values on
-   * quadrature points to values tested by
-   * different test function (second flag set to
-   * false), and if the result is to be added to
-   * some previous results or not.
+   * Internal function that applies the gradient operation of the tensor
+   * product in a given coordinate direction (first template argument), from
+   * polynomials to values on quadrature points (second flag set to true) or
+   * in an integration loop from values on quadrature points to values tested
+   * by different test function (second flag set to false), and if the result
+   * is to be added to some previous results or not.
    */
   template <int direction, bool dof_to_quad, bool add>
   void apply_gradients (const VectorizedArray<Number> in [],
@@ -1842,9 +1482,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                       "as the number of underlying vector components."));
 
 
-  // do not check for correct dimensions of data
-  // fields here, should be done in derived
-  // classes
+  // do not check for correct dimensions of data fields here, should be done
+  // in derived classes
 }
 
 
@@ -1955,8 +1594,7 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 
 namespace internal
 {
-  // write access to generic vectors that have
-  // operator ().
+  // write access to generic vectors that have operator ().
   template <typename VectorType>
   inline
   typename VectorType::value_type &
@@ -1968,8 +1606,7 @@ namespace internal
 
 
 
-  // read access to generic vectors that have
-  // operator ().
+  // read access to generic vectors that have operator ().
   template <typename VectorType>
   inline
   typename VectorType::value_type
@@ -1981,9 +1618,8 @@ namespace internal
 
 
 
-  // write access to distributed MPI vectors
-  // that have a local_element(uint) method to access data in
-  // local index space, which is what we use in
+  // write access to distributed MPI vectors that have a local_element(uint)
+  // method to access data in local index space, which is what we use in
   // DoFInfo and hence in read_dof_values etc.
   template <typename Number>
   inline
@@ -1996,9 +1632,8 @@ namespace internal
 
 
 
-  // read access to distributed MPI vectors that
-  // have a local_element(uint) method to access data in local
-  // index space, which is what we use in
+  // read access to distributed MPI vectors that have a local_element(uint)
+  // method to access data in local index space, which is what we use in
   // DoFInfo and hence in read_dof_values etc.
   template <typename Number>
   inline
@@ -2011,10 +1646,8 @@ namespace internal
 
 
 
-  // this is to make sure that the parallel
-  // partitioning in the
-  // parallel::distributed::Vector is really the
-  // same as stored in MatrixFree
+  // this is to make sure that the parallel partitioning in the
+  // parallel::distributed::Vector is really the same as stored in MatrixFree
   template <typename VectorType>
   inline
   void check_vector_compatibility (const VectorType                             &vec,
@@ -2036,8 +1669,7 @@ namespace internal
                        "compatible vector."));
   }
 
-  // A class to use the same code to read from
-  // and write to vector
+  // A class to use the same code to read from and write to vector
   template <typename Number>
   struct VectorReader
   {
@@ -2076,8 +1708,7 @@ namespace internal
     }
   };
 
-  // A class to use the same code to read from
-  // and write to vector
+  // A class to use the same code to read from and write to vector
   template <typename Number>
   struct VectorDistributorLocalToGlobal
   {
@@ -2115,8 +1746,7 @@ namespace internal
   };
 
 
-  // A class to use the same code to read from
-  // and write to vector
+  // A class to use the same code to read from and write to vector
   template <typename Number>
   struct VectorSetter
   {
@@ -2165,25 +1795,18 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 ::read_write_operation (const VectorOperation &operation,
                         VectorType            *src[]) const
 {
-  // This functions processes all the functions
-  // read_dof_values,
-  // distribute_local_to_global, and
-  // set_dof_values with the same code. The
-  // distinction between these three cases is
-  // made by the input VectorOperation that
-  // either reads values from a vector and puts
-  // the data into the local data field or write
-  // local data into the vector. Certain
-  // operations are no-ops for the given use
-  // case.
+  // This functions processes all the functions read_dof_values,
+  // distribute_local_to_global, and set_dof_values with the same code. The
+  // distinction between these three cases is made by the input
+  // VectorOperation that either reads values from a vector and puts the data
+  // into the local data field or write local data into the vector. Certain
+  // operations are no-ops for the given use case.
 
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
 
-  // loop over all local dofs. ind_local holds
-  // local number on cell, index iterates over
-  // the elements of index_local_to_global and
-  // dof_indices points to the global indices
-  // stored in index_local_to_global
+  // loop over all local dofs. ind_local holds local number on cell, index
+  // iterates over the elements of index_local_to_global and dof_indices
+  // points to the global indices stored in index_local_to_global
   const unsigned int *dof_indices = dof_info.begin_indices(cell);
   const std::pair<unsigned short,unsigned short> *indicators =
     dof_info.begin_indicators(cell);
@@ -2191,9 +1814,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
     dof_info.end_indicators(cell);
   unsigned int ind_local = 0;
 
-  // scalar case (or case when all components
-  // have the same degrees of freedom and sit on
-  // a different vector each)
+  // scalar case (or case when all components have the same degrees of freedom
+  // and sit on a different vector each)
   if (n_fe_components == 1)
     {
       const unsigned int n_local_dofs =
@@ -2205,12 +1827,11 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
         local_data[comp] =
           const_cast<Number *>(&values_dofs[comp][0][0]);
 
-      // standard case where there are sufficiently
-      // many cells to fill all vectors
+      // standard case where there are sufficiently many cells to fill all
+      // vectors
       if (at_irregular_cell == false)
         {
-          // check whether there is any constraint on
-          // the current cell
+          // check whether there is any constraint on the current cell
           if (indicators != indicators_end)
             {
               for ( ; indicators != indicators_end; ++indicators)
@@ -2224,9 +1845,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                   ind_local += indicators->first;
                   dof_indices   += indicators->first;
 
-                  // constrained case: build the local value as
-                  // a linear combination of the global value
-                  // according to constraints
+                  // constrained case: build the local value as a linear
+                  // combination of the global value according to constraints
                   Number value [n_components];
                   for (unsigned int comp=0; comp<n_components; ++comp)
                     operation.pre_constraints (local_data[comp][ind_local],
@@ -2248,8 +1868,7 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                   ind_local++;
                 }
 
-              // get the dof values past the last
-              // constraint
+              // get the dof values past the last constraint
               for (; ind_local < n_local_dofs; ++dof_indices, ++ind_local)
                 {
                   for (unsigned int comp=0; comp<n_components; ++comp)
@@ -2259,8 +1878,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             }
           else
             {
-              // no constraint at all: loop bounds are
-              // known, compiler can unroll without checks
+              // no constraint at all: loop bounds are known, compiler can
+              // unroll without checks
               AssertDimension (dof_info.end_indices(cell)-dof_indices,
                                static_cast<int>(n_local_dofs));
               for (unsigned int j=0; j<n_local_dofs; ++j)
@@ -2270,10 +1889,9 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             }
         }
 
-      // non-standard case: need to fill in zeros
-      // for those components that are not present
-      // (a bit more expensive), but there is not
-      // more than one such cell
+      // non-standard case: need to fill in zeros for those components that
+      // are not present (a bit more expensive), but there is not more than
+      // one such cell
       else
         {
           Assert (n_irreg_components_filled > 0, ExcInternalError());
@@ -2281,15 +1899,13 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             {
               for (unsigned int j=0; j<indicators->first; ++j)
                 {
-                  // non-constrained case: copy the data from
-                  // the global vector, src, to the local one,
-                  // local_src.
+                  // non-constrained case: copy the data from the global
+                  // vector, src, to the local one, local_src.
                   for (unsigned int comp=0; comp<n_components; ++comp)
                     operation.process_dof (dof_indices[j], *src[comp],
                                            local_data[comp][ind_local]);
 
-                  // here we jump over all the components that
-                  // are artificial
+                  // here we jump over all the components that are artificial
                   ++ind_local;
                   while (ind_local % VectorizedArray<Number>::n_array_elements
                          >= n_irreg_components_filled)
@@ -2301,9 +1917,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                 }
               dof_indices += indicators->first;
 
-              // constrained case: build the local value as
-              // a linear combination of the global value
-              // according to constraint
+              // constrained case: build the local value as a linear
+              // combination of the global value according to constraint
               Number value [n_components];
               for (unsigned int comp=0; comp<n_components; ++comp)
                 operation.pre_constraints (local_data[comp][ind_local],
@@ -2336,9 +1951,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
               Assert (dof_indices != dof_info.end_indices(cell),
                       ExcInternalError());
 
-              // non-constrained case: copy the data from
-              // the global vector, src, to the local one,
-              // local_dst.
+              // non-constrained case: copy the data from the global vector,
+              // src, to the local one, local_dst.
               for (unsigned int comp=0; comp<n_components; ++comp)
                 operation.process_dof (*dof_indices, *src[comp],
                                        local_data[comp][ind_local]);
@@ -2354,13 +1968,10 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
         }
     }
   else
-    // case with vector-valued finite elements
-    // where all components are included in one
-    // single vector. Assumption: first come all
-    // entries to the first component, then all
-    // entries to the second one, and so on. This
-    // is ensured by the way MatrixFree reads
-    // out the indices.
+    // case with vector-valued finite elements where all components are
+    // included in one single vector. Assumption: first come all entries to
+    // the first component, then all entries to the second one, and so
+    // on. This is ensured by the way MatrixFree reads out the indices.
     {
       internal::check_vector_compatibility (*src[0], dof_info);
       Assert (n_fe_components == n_components_, ExcNotImplemented());
@@ -2370,8 +1981,7 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
         const_cast<Number *>(&values_dofs[0][0][0]);
       if (at_irregular_cell == false)
         {
-          // check whether there is any constraint on
-          // the current cell
+          // check whether there is any constraint on the current cell
           if (indicators != indicators_end)
             {
               for ( ; indicators != indicators_end; ++indicators)
@@ -2383,9 +1993,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                   ind_local += indicators->first;
                   dof_indices   += indicators->first;
 
-                  // constrained case: build the local value as
-                  // a linear combination of the global value
-                  // according to constraints
+                  // constrained case: build the local value as a linear
+                  // combination of the global value according to constraints
                   Number value;
                   operation.pre_constraints (local_data[ind_local], value);
 
@@ -2402,8 +2011,7 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                   ind_local++;
                 }
 
-              // get the dof values past the last
-              // constraint
+              // get the dof values past the last constraint
               for (; ind_local<n_local_dofs; ++dof_indices, ++ind_local)
                 operation.process_dof (*dof_indices, *src[0],
                                        local_data[ind_local]);
@@ -2412,8 +2020,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             }
           else
             {
-              // no constraint at all: loop bounds are
-              // known, compiler can unroll without checks
+              // no constraint at all: loop bounds are known, compiler can
+              // unroll without checks
               AssertDimension (dof_info.end_indices(cell)-dof_indices,
                                static_cast<int>(n_local_dofs));
               for (unsigned int j=0; j<n_local_dofs; ++j)
@@ -2422,10 +2030,9 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             }
         }
 
-      // non-standard case: need to fill in zeros
-      // for those components that are not present
-      // (a bit more expensive), but there is not
-      // more than one such cell
+      // non-standard case: need to fill in zeros for those components that
+      // are not present (a bit more expensive), but there is not more than
+      // one such cell
       else
         {
           Assert (n_irreg_components_filled > 0, ExcInternalError());
@@ -2433,14 +2040,12 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             {
               for (unsigned int j=0; j<indicators->first; ++j)
                 {
-                  // non-constrained case: copy the data from
-                  // the global vector, src, to the local one,
-                  // local_src.
+                  // non-constrained case: copy the data from the global
+                  // vector, src, to the local one, local_src.
                   operation.process_dof (dof_indices[j], *src[0],
                                          local_data[ind_local]);
 
-                  // here we jump over all the components that
-                  // are artificial
+                  // here we jump over all the components that are artificial
                   ++ind_local;
                   while (ind_local % VectorizedArray<Number>::n_array_elements
                          >= n_irreg_components_filled)
@@ -2451,9 +2056,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                 }
               dof_indices += indicators->first;
 
-              // constrained case: build the local value as
-              // a linear combination of the global value
-              // according to constraint
+              // constrained case: build the local value as a linear
+              // combination of the global value according to constraint
               Number value;
               operation.pre_constraints (local_data[ind_local], value);
 
@@ -2480,9 +2084,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
               Assert (dof_indices != dof_info.end_indices(cell),
                       ExcInternalError());
 
-              // non-constrained case: copy the data from
-              // the global vector, src, to the local one,
-              // local_dst.
+              // non-constrained case: copy the data from the global vector,
+              // src, to the local one, local_dst.
               operation.process_dof (*dof_indices, *src[0],
                                      local_data[ind_local]);
               ++ind_local;
@@ -2509,9 +2112,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 {
   AssertDimension (n_components_, n_fe_components);
 
-  // only need one component, but to silent
-  // compiler warnings, use n_components copies
-  // here (but these will not be used)
+  // only need one component, but to silent compiler warnings, use
+  // n_components copies here (but these will not be used)
   VectorType *src_data[n_components];
   for (unsigned int d=0; d<n_components; ++d)
     src_data[d] = const_cast<VectorType *>(&src);
@@ -2593,9 +2195,9 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 ::read_dof_values_plain (const VectorType &src)
 {
   AssertDimension (n_components_, n_fe_components);
-  // only need one component, but to avoid
-  // compiler warnings, use n_components copies
-  // here (but these will not be used)
+
+  // only need one component, but to avoid compiler warnings, use n_components
+  // copies here (but these will not be used)
   const VectorType *src_data[n_components];
   for (unsigned int d=0; d<n_components; ++d)
     src_data[d] = &src;
@@ -2660,9 +2262,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
   Assert (dof_values_initialized==true,
           internal::ExcAccessToUninitializedField());
 
-  // only need one component, but to avoid
-  // compiler warnings, use n_components copies
-  // here (but these will not be used)
+  // only need one component, but to avoid compiler warnings, use n_components
+  // copies here (but these will not be used)
   VectorType *dst_data [n_components];
   for (unsigned int d=0; d<n_components; ++d)
     dst_data[d] = &dst;
@@ -2739,9 +2340,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
   Assert (dof_values_initialized==true,
           internal::ExcAccessToUninitializedField());
 
-  // only need one component, but to avoid
-  // compiler warnings, use n_components copies
-  // here (but these will not be used)
+  // only need one component, but to avoid compiler warnings, use n_components
+  // copies here (but these will not be used)
   VectorType *dst_data [n_components];
   for (unsigned int d=0; d<n_components; ++d)
     dst_data[d] = &dst;
@@ -2816,23 +2416,18 @@ void
 FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 ::read_dof_values_plain (const VectorType *src[])
 {
-  // this is different from the other three
-  // operations because we do not use
-  // constraints here, so this is a separate
-  // function.
+  // this is different from the other three operations because we do not use
+  // constraints here, so this is a separate function.
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
 
-  // loop over all local dofs. ind_local holds
-  // local number on cell, index iterates over
-  // the elements of index_local_to_global and
-  // dof_indices points to the global indices
-  // stored in index_local_to_global
+  // loop over all local dofs. ind_local holds local number on cell, index
+  // iterates over the elements of index_local_to_global and dof_indices
+  // points to the global indices stored in index_local_to_global
   const unsigned int *dof_indices = dof_info.begin_indices_plain(cell);
 
-  // scalar case (or case when all components
-  // have the same degrees of freedom and sit on
-  // a different vector each)
+  // scalar case (or case when all components have the same degrees of freedom
+  // and sit on a different vector each)
   if (n_fe_components == 1)
     {
       const unsigned int n_local_dofs =
@@ -2843,8 +2438,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
       for (unsigned int comp=0; comp<n_components; ++comp)
         local_src_number[comp] = &values_dofs[comp][0][0];
 
-      // standard case where there are sufficiently
-      // many cells to fill all vectors
+      // standard case where there are sufficiently many cells to fill all
+      // vectors
       if (at_irregular_cell == false)
         {
           for (unsigned int j=0; j<n_local_dofs; ++j)
@@ -2853,19 +2448,17 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                 internal::vector_access (*src[comp], dof_indices[j]);
         }
 
-      // non-standard case: need to fill in zeros
-      // for those components that are not present
-      // (a bit more expensive), but there is not
-      // more than one such cell
+      // non-standard case: need to fill in zeros for those components that
+      // are not present (a bit more expensive), but there is not more than
+      // one such cell
       else
         {
           Assert (n_irreg_components_filled > 0, ExcInternalError());
           for (unsigned int ind_local=0; ind_local<n_local_dofs;
                ++dof_indices)
             {
-              // non-constrained case: copy the data from
-              // the global vector, src, to the local one,
-              // local_dst.
+              // non-constrained case: copy the data from the global vector,
+              // src, to the local one, local_dst.
               for (unsigned int comp=0; comp<n_components; ++comp)
                 local_src_number[comp][ind_local] =
                   internal::vector_access (*src[comp], *dof_indices);
@@ -2880,13 +2473,10 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
         }
     }
   else
-    // case with vector-valued finite elements
-    // where all components are included in one
-    // single vector. Assumption: first come all
-    // entries to the first component, then all
-    // entries to the second one, and so on. This
-    // is ensured by the way MatrixFree reads
-    // out the indices.
+    // case with vector-valued finite elements where all components are
+    // included in one single vector. Assumption: first come all entries to
+    // the first component, then all entries to the second one, and so
+    // on. This is ensured by the way MatrixFree reads out the indices.
     {
       internal::check_vector_compatibility (*src[0], dof_info);
       Assert (n_fe_components == n_components_, ExcNotImplemented());
@@ -2900,18 +2490,16 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
               internal::vector_access (*src[0], dof_indices[j]);
         }
 
-      // non-standard case: need to fill in zeros
-      // for those components that are not present
-      // (a bit more expensive), but there is not
-      // more than one such cell
+      // non-standard case: need to fill in zeros for those components that
+      // are not present (a bit more expensive), but there is not more than
+      // one such cell
       else
         {
           Assert (n_irreg_components_filled > 0, ExcInternalError());
           for (unsigned int ind_local=0; ind_local<n_local_dofs; ++dof_indices)
             {
-              // non-constrained case: copy the data from
-              // the global vector, src, to the local one,
-              // local_dst.
+              // non-constrained case: copy the data from the global vector,
+              // src, to the local one, local_dst.
               local_src_number[ind_local] =
                 internal::vector_access (*src[0], *dof_indices);
               ++ind_local;
@@ -2934,6 +2522,34 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 
 /*------------------------------ access to data fields ----------------------*/
 
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+         int n_components, typename Number>
+inline
+const VectorizedArray<Number> *
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
+begin_dof_values () const
+{
+  return &values_dofs[0][0];
+}
+
+
+
+template <int dim, int dofs_per_cell_, int n_q_points_,
+         int n_components, typename Number>
+inline
+VectorizedArray<Number> *
+FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
+begin_dof_values ()
+{
+#ifdef DEBUG
+  dof_values_initialized = true;
+#endif
+  return &values_dofs[0][0];
+}
+
+
+
 template <int dim, int dofs_per_cell_, int n_q_points_,
          int n_components, typename Number>
 inline
@@ -3097,9 +2713,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 
 namespace internal
 {
-  // compute tmp = hess_unit(u) * J^T. do this
-  // manually because we do not store the lower
-  // diagonal because of symmetry
+  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
+  // store the lower diagonal because of symmetry
   template <int dim, int n_q_points, typename Number>
   inline
   void
@@ -3198,15 +2813,13 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
             & jac_grad_UT = jacobian_grad_upper[q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
-          // compute laplacian before the gradient
-          // because it needs to access unscaled
-          // gradient data
+          // compute laplacian before the gradient because it needs to access
+          // unscaled gradient data
           VectorizedArray<Number> tmp[dim][dim];
           internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
                                             q_point, tmp);
 
-          // compute first part of hessian,
-          // J * tmp = J * hess_unit(u) * J^T
+          // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
             for (unsigned int e=d; e<dim; ++e)
               {
@@ -3234,22 +2847,19 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
               hessian_out[comp][e][d] = hessian_out[comp][d][e];
         }
     }
-  // cell with general Jacobian, but constant
-  // within the cell
+  // cell with general Jacobian, but constant within the cell
   else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
     {
       const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
-          // compute laplacian before the gradient
-          // because it needs to access unscaled
-          // gradient data
+          // compute laplacian before the gradient because it needs to access
+          // unscaled gradient data
           VectorizedArray<Number> tmp[dim][dim];
           internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
                                             q_point, tmp);
 
-          // compute first part of hessian,
-          // J * tmp = J * hess_unit(u) * J^T
+          // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
             for (unsigned int e=d; e<dim; ++e)
               {
@@ -3258,9 +2868,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                   hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
               }
 
-          // no J' * grad(u) part here because the
-          // Jacobian is constant throughout the cell
-          // and hence, its derivative is zero
+          // no J' * grad(u) part here because the Jacobian is constant
+          // throughout the cell and hence, its derivative is zero
 
           // take symmetric part
           for (unsigned int d=0; d<dim; ++d)
@@ -3304,15 +2913,14 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
       const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
-          // compute laplacian before the gradient
-          // because it needs to access unscaled
-          // gradient data
+          // compute laplacian before the gradient because it needs to access
+          // unscaled gradient data
           VectorizedArray<Number> tmp[dim][dim];
           internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
                                             q_point, tmp);
 
-          // compute only the trace part of hessian,
-          // J * tmp = J * hess_unit(u) * J^T
+          // compute only the trace part of hessian, J * tmp = J *
+          // hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
             {
               hessian_out[comp][d] = jac[d][0] * tmp[0][d];
@@ -3326,22 +2934,20 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
                                        this->gradients_quad[comp][e][q_point]);
         }
     }
-  // cell with general Jacobian, but constant
-  // within the cell
+  // cell with general Jacobian, but constant within the cell
   else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
     {
       const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
       for (unsigned int comp=0; comp<n_components; comp++)
         {
-          // compute laplacian before the gradient
-          // because it needs to access unscaled
-          // gradient data
+          // compute laplacian before the gradient because it needs to access
+          // unscaled gradient data
           VectorizedArray<Number> tmp[dim][dim];
           internal::hessian_unit_times_jac (jac, this->hessians_quad[comp],
                                             q_point, tmp);
 
-          // compute only the trace part of hessian,
-          // J * tmp = J * hess_unit(u) * J^T
+          // compute only the trace part of hessian, J * tmp = J *
+          // hess_unit(u) * J^T
           for (unsigned int d=0; d<dim; ++d)
             {
               hessian_out[comp][d] = jac[d][0] * tmp[0][d];
@@ -3558,9 +3164,8 @@ Tensor<1,dim,VectorizedArray<Number> >
 FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number>
 ::get_gradient (const unsigned int q_point) const
 {
-  // could use the base class gradient, but that
-  // involves too many inefficient
-  // initialization
+  // could use the base class gradient, but that involves too many inefficient
+  // initialization operations on tensors
 
   Assert (this->gradients_quad_initialized==true,
           internal::ExcAccessToUninitializedField());
@@ -3793,8 +3398,7 @@ SymmetricTensor<2,dim,VectorizedArray<Number> >
 FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
 ::get_symmetric_gradient (const unsigned int q_point) const
 {
-  // copy from generic function into
-  // dim-specialization function
+  // copy from generic function into dim-specialization function
   const Tensor<2,dim,VectorizedArray<Number> > grad = get_gradient(q_point);
   VectorizedArray<Number> symmetrized [(dim*dim+dim)/2];
   VectorizedArray<Number> half = make_vectorized_array (0.5);
@@ -3830,8 +3434,7 @@ Tensor<1,dim==2?1:dim,VectorizedArray<Number> >
 FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
 ::get_curl (const unsigned int q_point) const
 {
-  // copy from generic function into
-  // dim-specialization function
+  // copy from generic function into dim-specialization function
   const Tensor<2,dim,VectorizedArray<Number> > grad = get_gradient(q_point);
   Tensor<1,dim==2?1:dim,VectorizedArray<Number> > curl (false);
   switch (dim)
@@ -3963,11 +3566,9 @@ FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
                             sym_grad,
                             const unsigned int q_point)
 {
-  // could have used base class operator, but
-  // that involves some overhead which is
-  // inefficient. it is nice to have the
-  // symmetric tensor because that saves some
-  // operations
+  // could have used base class operator, but that involves some overhead
+  // which is inefficient. it is nice to have the symmetric tensor because
+  // that saves some operations
 #ifdef DEBUG
   Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   AssertIndexRange (q_point, n_q_points);
@@ -4069,8 +3670,8 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
   BaseClass (data_in, fe_no, quad_no_in)
 {
 #ifdef DEBUG
-  // print error message when the dimensions do
-  // not match. Propose a possible fix
+  // print error message when the dimensions do not match. Propose a possible
+  // fix
   if (dofs_per_cell != this->data.dofs_per_cell ||
       n_q_points != this->data.n_q_points)
     {
@@ -4083,8 +3684,8 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
       message += Utilities::int_to_string(fe_no) + ", ";
       message += Utilities::int_to_string(quad_no_in) + ")\n";
 
-      // check whether some other vector component
-      // has the correct number of points
+      // check whether some other vector component has the correct number of
+      // points
       unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
                    proposed_quad_comp = numbers::invalid_unsigned_int;
       for (unsigned int no=0; no<this->matrix_info.n_components(); ++no)
@@ -4120,9 +3721,8 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
             correct_pos += "  \n";
           message += "                                                   " + correct_pos;
         }
-      // ok, did not find the numbers specified by
-      // the template arguments in the given
-      // list. Suggest correct template arguments
+      // ok, did not find the numbers specified by the template arguments in
+      // the given list. Suggest correct template arguments
       const unsigned int proposed_fe_degree = static_cast<unsigned int>(std::pow(1.001*this->data.dofs_per_cell,1./dim))-1;
       const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(std::pow(1.001*this->data.n_q_points,1./dim));
       message += "Wrong template arguments:\n";
@@ -4159,10 +3759,9 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
 
 namespace internal
 {
-  // evaluates the given shape data in 1d-3d
-  // using the tensor product form. does not use
-  // the tensor product form and corresponds to
-  // a usual matrix-matrix product
+  // evaluates the given shape data in 1d-3d using the tensor product
+  // form. does not use the tensor product form and corresponds to a usual
+  // matrix-matrix product
   template <int dim, int fe_degree, int n_q_points_1d, typename Number,
            int direction, bool dof_to_quad, bool add>
   inline
@@ -4205,10 +3804,9 @@ namespace internal
                   out[stride*col]        += res0;
               }
 
-            // increment: in regular case, just go to the
-            // next point in x-direction. If we are at the
-            // end of one chunk in x-dir, need to jump
-            // over to the next layer in z-direction
+            // increment: in regular case, just go to the next point in
+            // x-direction. If we are at the end of one chunk in x-dir, need
+            // to jump over to the next layer in z-direction
             switch (direction)
               {
               case 0:
@@ -4234,14 +3832,10 @@ namespace internal
 
 
 
-  // This performs the evaluation of function
-  // values, gradients and Hessians for
-  // tensor-product finite elements. The
-  // operation is used for both
-  // FEEvaluationGeneral and FEEvaluation, which
-  // provide different functions apply_values,
-  // apply_gradients in the individual
-  // coordinate directions
+  // This performs the evaluation of function values, gradients and Hessians
+  // for tensor-product finite elements. The operation is used for both
+  // FEEvaluationGeneral and FEEvaluation, which provide different functions
+  // apply_values, apply_gradients in the individual coordinate directions
   template <typename FEEval>
   inline
   void
@@ -4596,10 +4190,9 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
           ExcNotInitialized());
   AssertIndexRange (q, n_q_points);
 
-  // Cartesian mesh: not all quadrature points
-  // are stored, only the diagonal. Hence, need
-  // to find the tensor product index and
-  // retrieve the value from that
+  // Cartesian mesh: not all quadrature points are stored, only the
+  // diagonal. Hence, need to find the tensor product index and retrieve the
+  // value from that
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       Point<dim,VectorizedArray<Number> > point (false);
@@ -4621,8 +4214,7 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
           return point;
         }
     }
-  // all other cases: just return the respective
-  // data as it is fully stored
+  // all other cases: just return the respective data as it is fully stored
   else
     return this->quadrature_points[q];
 }
@@ -4707,9 +4299,8 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                                                 -j-1][0]) < zero_tol,
               ExcMessage(error_message));
 
-  // shape values should be zero at for all
-  // basis functions except for one where they
-  // are one in the middle
+  // shape values should be zero at for all basis functions except for one
+  // where they are one in the middle
   if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
     {
       for (int i=0; i<static_cast<int>(n_dofs_1d/2); ++i)
@@ -4721,8 +4312,8 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
               ExcMessage(error_message));
     }
 
-  // skew-symmetry for gradient, zero of middle
-  // basis function in middle quadrature point
+  // skew-symmetry for gradient, zero of middle basis function in middle
+  // quadrature point
   for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
     for (unsigned int j=0; j<n_q_points_1d; ++j)
       Assert (std::fabs(this->data.shape_gradients[i*n_q_points_1d+j][0] +
@@ -4794,15 +4385,11 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
   const int stride    = Utilities::fixed_int_power<nn,direction>::value;
 
-  // This loop specializes the general
-  // application of tensor-product based
-  // elements for "symmetric" finite elements,
-  // i.e., when the shape functions are
-  // symmetric about 0.5 and the quadrature
-  // points are, too. In that case, the 1D shape
-  // values read (sorted lexicographically, rows
-  // run over 1D dofs, columns over quadrature
-  // points):
+  // This loop specializes the general application of tensor-product based
+  // elements for "symmetric" finite elements, i.e., when the shape functions
+  // are symmetric about 0.5 and the quadrature points are, too. In that case,
+  // the 1D shape values read (sorted lexicographically, rows run over 1D
+  // dofs, columns over quadrature points):
   // Q2 --> [ 0.687  0 -0.087 ]
   //        [ 0.4    1  0.4   ]
   //        [-0.087  0  0.687 ]
@@ -4816,10 +4403,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   //        [ 0.176  0.039  0  1.059  0.608 ]
   //        [-0.032 -0.007  0  0.022  0.658 ]
   //
-  // In these matrices, we want to use avoid
-  // computations involving zeros and ones and
-  // in addition use the symmetry in entries to
-  // reduce the number of read operations.
+  // In these matrices, we want to use avoid computations involving zeros and
+  // ones and in addition use the symmetry in entries to reduce the number of
+  // read operations.
   const VectorizedArray<Number> *shape_values = this->data.shape_values.begin();
   for (int i2=0; i2<n_blocks2; ++i2)
     {
@@ -4960,10 +4546,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                 out[stride*n_cols] += res0;
             }
 
-          // increment: in regular case, just go to the
-          // next point in x-direction. If we are at the
-          // end of one chunk in x-dir, need to jump
-          // over to the next layer in z-direction
+          // increment: in regular case, just go to the next point in
+          // x-direction. If we are at the end of one chunk in x-dir, need to
+          // jump over to the next layer in z-direction
           switch (direction)
             {
             case 0:
@@ -5013,11 +4598,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
     {
       for (int i1=0; i1<n_blocks1; ++i1)
         {
-          // For the specialized loop used for the
-          // gradient computation in here, the 1D shape
-          // values read (sorted lexicographically, rows
-          // run over 1D dofs, columns over quadrature
-          // points):
+          // For the specialized loop used for the gradient computation in
+          // here, the 1D shape values read (sorted lexicographically, rows
+          // run over 1D dofs, columns over quadrature points):
           // Q2 --> [-2.549 -1  0.549 ]
           //        [ 3.098  0 -3.098 ]
           //        [-0.549  1  2.549 ]
@@ -5031,9 +4614,8 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
           //        [ 2.306 -2.066  2.667  2.76 -10.111 ]
           //        [-0.413  0.353 -0.333 -0.353  0.413 ]
           //
-          // In these matrices, we want to use avoid
-          // computations involving zeros and ones and
-          // in addition use the symmetry in entries to
+          // In these matrices, we want to use avoid computations involving
+          // zeros and ones and in addition use the symmetry in entries to
           // reduce the number of read operations.
           for (int col=0; col<n_cols; ++col)
             {
@@ -5124,11 +4706,10 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                 out[stride*n_cols] += res0;
             }
 
-          // increment: in regular case, just go to the
-          // next point in x-direction. for y-part in 3D
-          // and if we are at the end of one chunk in
-          // x-dir, need to jump over to the next layer
-          // in z-direction
+          // increment: in regular case, just go to the next point in
+          // x-direction. for y-part in 3D and if we are at the end of one
+          // chunk in x-dir, need to jump over to the next layer in
+          // z-direction
           switch (direction)
             {
             case 0:
@@ -5155,11 +4736,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 
 
 
-// Laplacian operator application. Very
-// similar to value application because the
-// same symmetry relations hold. However, it
-// is not possible to omit some values that
-// are zero for the values
+// Laplacian operator application. Very similar to value application because
+// the same symmetry relations hold. However, it is not possible to omit some
+// values that are zero for the values
 template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
          typename Number>
 template <int direction, bool dof_to_quad, bool add>
@@ -5285,10 +4864,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                 out[stride*n_cols] += res0;
             }
 
-          // increment: in regular case, just go to the
-          // next point in x-direction. If we are at the
-          // end of one chunk in x-dir, need to jump
-          // over to the next layer in z-direction
+          // increment: in regular case, just go to the next point in
+          // x-direction. If we are at the end of one chunk in x-dir, need to
+          // jump over to the next layer in z-direction
           switch (direction)
             {
             case 0:
@@ -5382,9 +4960,8 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
       this->values_quad_initialized = true;
 #endif
     }
-  // separate implementation here compared to
-  // the general case because the values are an
-  // identity operation
+  // separate implementation here compared to the general case because the
+  // values are an identity operation
   if (evaluate_grad == true)
     {
       for (unsigned int comp=0; comp<n_components; comp++)
@@ -5561,16 +5138,12 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
   const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
   const int stride    = Utilities::fixed_int_power<nn,direction>::value;
 
-  // This loop specializes the application of
-  // the tensor product loop for Gauss-Lobatto
-  // elements which are symmetric about 0.5 just
-  // as the general class of elements treated by
-  // FEEvaluation, have diagonal shape matrices
-  // for the values and have the following
-  // gradient matrices (notice the zeros on the
-  // diagonal in the interior points, which is
-  // due to the construction of Legendre
-  // polynomials):
+  // This loop specializes the application of the tensor product loop for
+  // Gauss-Lobatto elements which are symmetric about 0.5 just as the general
+  // class of elements treated by FEEvaluation, have diagonal shape matrices
+  // for the values and have the following gradient matrices (notice the zeros
+  // on the diagonal in the interior points, which is due to the construction
+  // of Legendre polynomials):
   // Q2 --> [-3 -1  1 ]
   //        [ 4  0 -4 ]
   //        [-1  1  3 ]
@@ -5641,8 +5214,7 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
                           val1 = this->data.shape_gradients[(nn-col-1)*mm+ind];
                         }
 
-                      // at inner points, the gradient is zero for
-                      // ind==col
+                      // at inner points, the gradient is zero for ind==col
                       in0 = in[stride*ind];
                       in1 = in[stride*(mm-1-ind)];
                       if (ind == col)
@@ -5712,11 +5284,10 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
                 out[stride*n_cols] += res0;
             }
 
-          // increment: in regular case, just go to the
-          // next point in x-direction. for y-part in 3D
-          // and if we are at the end of one chunk in
-          // x-dir, need to jump over to the next layer
-          // in z-direction
+          // increment: in regular case, just go to the next point in
+          // x-direction. for y-part in 3D and if we are at the end of one
+          // chunk in x-dir, need to jump over to the next layer in
+          // z-direction
           switch (direction)
             {
             case 0:

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.