Initialize restriction and prolongation matrix in FESystem and FE_Q only upon request...

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

diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h
index 97ba56192fc1049ed9e0f3cdb6b142908482461c..fbf458b9fe2a2bb2800e012544c6e5cf126f7fd7 100644 (file)
--- a/deal.II/include/deal.II/fe/fe.h
+++ b/deal.II/include/deal.II/fe/fe.h
@@ -343,14 +343,11 @@ class FiniteElement : public Subscriptor,
 {
 public:
   /**
-   * Base class for internal data.
-   * Adds data for second derivatives to
+   * Base class for internal data.  Adds data for second derivatives to
    * Mapping::InternalDataBase()
    *
-   * For information about the
-   * general purpose of this class,
-   * see the documentation of the
-   * base class.
+   * For information about the general purpose of this class, see the
+   * documentation of the base class.
    *
    * @author Guido Kanschat, 2001
    */
@@ -358,35 +355,24 @@ public:
   {
   public:
     /**
-     * Destructor. Needed to
-     * avoid memory leaks with
-     * difference quotients.
+     * Destructor. Needed to avoid memory leaks with difference quotients.
      */
     virtual ~InternalDataBase ();
 
     /**
-     * Initialize some pointers
-     * used in the computation of
-     * second derivatives by
-     * finite differencing of
-     * gradients.
+     * Initialize some pointers used in the computation of second derivatives
+     * by finite differencing of gradients.
      */
     void initialize_2nd (const FiniteElement<dim,spacedim> *element,
                          const Mapping<dim,spacedim>       &mapping,
                          const Quadrature<dim>    &quadrature);
 
     /**
-     * Storage for FEValues
-     * objects needed to
-     * approximate second
-     * derivatives.
+     * Storage for FEValues objects needed to approximate second derivatives.
      *
-     * The ordering is <i>p+hx</i>,
-     * <i>p+hy</i>, <i>p+hz</i>,
-     * <i>p-hx</i>, <i>p-hy</i>,
-     * <i>p-hz</i>, where unused
-     * entries in lower dimensions
-     * are missing.
+     * The ordering is <i>p+hx</i>, <i>p+hy</i>, <i>p+hz</i>, <i>p-hx</i>,
+     * <i>p-hy</i>, <i>p-hz</i>, where unused entries in lower dimensions are
+     * missing.
      */
     std::vector<FEValues<dim,spacedim>*> differences;
   };
@@ -400,74 +386,43 @@ public:
                  const std::vector<ComponentMask> &nonzero_components);
 
   /**
-   * Virtual destructor. Makes sure
-   * that pointers to this class
-   * are deleted properly.
+   * Virtual destructor. Makes sure that pointers to this class are deleted
+   * properly.
    */
   virtual ~FiniteElement ();
 
   /**
-   * Return a string that uniquely
-   * identifies a finite
-   * element. The general
-   * convention is that this is the
-   * class name, followed by the
-   * dimension in angle
-   * brackets, and the polynomial
-   * degree and whatever else is
-   * necessary in parentheses. For
-   * example, <tt>FE_Q<2>(3)</tt> is the
-   * value returned for a cubic
-   * element in 2d.
+   * Return a string that uniquely identifies a finite element. The general
+   * convention is that this is the class name, followed by the dimension in
+   * angle brackets, and the polynomial degree and whatever else is necessary
+   * in parentheses. For example, <tt>FE_Q<2>(3)</tt> is the value returned
+   * for a cubic element in 2d.
    *
-   * Systems of elements have their
-   * own naming convention, see the
-   * FESystem class.
+   * Systems of elements have their own naming convention, see the FESystem
+   * class.
    */
   virtual std::string get_name () const = 0;
 
   /**
-   * This operator returns a
-   * reference to the present
-   * object if the argument given
-   * equals to zero. While this
-   * does not seem particularly
-   * useful, it is helpful in
-   * writing code that works with
-   * both ::DoFHandler and the hp
-   * version hp::DoFHandler, since
-   * one can then write code like
-   * this:
+   * This operator returns a reference to the present object if the argument
+   * given equals to zero. While this does not seem particularly useful, it is
+   * helpful in writing code that works with both ::DoFHandler and the hp
+   * version hp::DoFHandler, since one can then write code like this:
    * @code
    *   dofs_per_cell
    *     = dof_handler->get_fe()[cell->active_fe_index()].dofs_per_cell;
    * @endcode
    *
-   * This code doesn't work in both
-   * situations without the present
-   * operator because
-   * DoFHandler::get_fe() returns a
-   * finite element, whereas
-   * hp::DoFHandler::get_fe()
-   * returns a collection of finite
-   * elements that doesn't offer a
-   * <code>dofs_per_cell</code>
-   * member variable: one first has
-   * to select which finite element
-   * to work on, which is done
-   * using the
-   * operator[]. Fortunately,
-   * <code>cell-@>active_fe_index()</code>
-   * also works for non-hp classes
-   * and simply returns zero in
-   * that case. The present
-   * operator[] accepts this zero
-   * argument, by returning the
-   * finite element with index zero
-   * within its collection (that,
-   * of course, consists only of
-   * the present finite element
-   * anyway).
+   * This code doesn't work in both situations without the present operator
+   * because DoFHandler::get_fe() returns a finite element, whereas
+   * hp::DoFHandler::get_fe() returns a collection of finite elements that
+   * doesn't offer a <code>dofs_per_cell</code> member variable: one first has
+   * to select which finite element to work on, which is done using the
+   * operator[]. Fortunately, <code>cell-@>active_fe_index()</code> also works
+   * for non-hp classes and simply returns zero in that case. The present
+   * operator[] accepts this zero argument, by returning the finite element
+   * with index zero within its collection (that, of course, consists only of
+   * the present finite element anyway).
    */
   const FiniteElement<dim,spacedim> &operator[] (const unsigned int fe_index) const;
 
@@ -477,139 +432,75 @@ public:
    */
 
   /**
-   * Return the value of the
-   * @p ith shape function at the
-   * point @p p. @p p is a point
-   * on the reference element. If
-   * the finite element is
-   * vector-valued, then return the
-   * value of the only non-zero
-   * component of the vector value
-   * of this shape function. If the
-   * shape function has more than
-   * one non-zero component (which
-   * we refer to with the term
-   * non-primitive), then derived
-   * classes implementing this
-   * function should throw an
-   * exception of type
-   * ExcShapeFunctionNotPrimitive. In
-   * that case, use the
-   * shape_value_component()
-   * function.
+   * Return the value of the @p ith shape function at the point @p p. @p p is
+   * a point on the reference element. If the finite element is vector-valued,
+   * then return the value of the only non-zero component of the vector value
+   * of this shape function. If the shape function has more than one non-zero
+   * component (which we refer to with the term non-primitive), then derived
+   * classes implementing this function should throw an exception of type
+   * ExcShapeFunctionNotPrimitive. In that case, use the
+   * shape_value_component() function.
    *
-   * An
-   * ExcUnitShapeValuesDoNotExist
-   * is thrown if the shape values
-   * of the FiniteElement under
-   * consideration depends on the
-   * shape of the cell in real
-   * space.
+   * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+   * FiniteElement under consideration depends on the shape of the cell in
+   * real space.
    */
   virtual double shape_value (const unsigned int  i,
                               const Point<dim>   &p) const;
 
   /**
-   * Just like for shape_value(),
-   * but this function will be
-   * called when the shape function
-   * has more than one non-zero
-   * vector component. In that
-   * case, this function should
-   * return the value of the
-   * @p component-th vector
-   * component of the @p ith shape
-   * function at point @p p.
+   * Just like for shape_value(), but this function will be called when the
+   * shape function has more than one non-zero vector component. In that case,
+   * this function should return the value of the @p component-th vector
+   * component of the @p ith shape function at point @p p.
    */
   virtual double shape_value_component (const unsigned int i,
                                         const Point<dim>   &p,
                                         const unsigned int component) const;
 
   /**
-   * Return the gradient of the
-   * @p ith shape function at the
-   * point @p p. @p p is a point
-   * on the reference element, and
-   * likewise the gradient is the
-   * gradient on the unit cell with
-   * respect to unit cell
-   * coordinates. If
-   * the finite element is
-   * vector-valued, then return the
-   * value of the only non-zero
-   * component of the vector value
-   * of this shape function. If the
-   * shape function has more than
-   * one non-zero component (which
-   * we refer to with the term
-   * non-primitive), then derived
-   * classes implementing this
-   * function should throw an
-   * exception of type
-   * ExcShapeFunctionNotPrimitive. In
-   * that case, use the
-   * shape_grad_component()
-   * function.
+   * Return the gradient of the @p ith shape function at the point @p p. @p p
+   * is a point on the reference element, and likewise the gradient is the
+   * gradient on the unit cell with respect to unit cell coordinates. If the
+   * finite element is vector-valued, then return the value of the only
+   * non-zero component of the vector value of this shape function. If the
+   * shape function has more than one non-zero component (which we refer to
+   * with the term non-primitive), then derived classes implementing this
+   * function should throw an exception of type
+   * ExcShapeFunctionNotPrimitive. In that case, use the
+   * shape_grad_component() function.
    *
-   * An
-   * ExcUnitShapeValuesDoNotExist
-   * is thrown if the shape values
-   * of the FiniteElement under
-   * consideration depends on the
-   * shape of the cell in real
-   * space.
+   * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+   * FiniteElement under consideration depends on the shape of the cell in
+   * real space.
    */
   virtual Tensor<1,dim> shape_grad (const unsigned int  i,
                                     const Point<dim>   &p) const;
 
   /**
-   * Just like for shape_grad(),
-   * but this function will be
-   * called when the shape function
-   * has more than one non-zero
-   * vector component. In that
-   * case, this function should
-   * return the gradient of the
-   * @p component-th vector
-   * component of the @p ith shape
-   * function at point @p p.
+   * Just like for shape_grad(), but this function will be called when the
+   * shape function has more than one non-zero vector component. In that case,
+   * this function should return the gradient of the @p component-th vector
+   * component of the @p ith shape function at point @p p.
    */
   virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
                                               const Point<dim>   &p,
                                               const unsigned int component) const;
 
   /**
-   * Return the tensor of second
-   * derivatives of the @p ith
-   * shape function at point @p p
-   * on the unit cell. The
-   * derivatives are derivatives on
-   * the unit cell with respect to
-   * unit cell coordinates. If
-   * the finite element is
-   * vector-valued, then return the
-   * value of the only non-zero
-   * component of the vector value
-   * of this shape function. If the
-   * shape function has more than
-   * one non-zero component (which
-   * we refer to with the term
-   * non-primitive), then derived
-   * classes implementing this
-   * function should throw an
-   * exception of type
-   * ExcShapeFunctionNotPrimitive. In
-   * that case, use the
-   * shape_grad_grad_component()
-   * function.
+   * Return the tensor of second derivatives of the @p ith shape function at
+   * point @p p on the unit cell. The derivatives are derivatives on the unit
+   * cell with respect to unit cell coordinates. If the finite element is
+   * vector-valued, then return the value of the only non-zero component of
+   * the vector value of this shape function. If the shape function has more
+   * than one non-zero component (which we refer to with the term
+   * non-primitive), then derived classes implementing this function should
+   * throw an exception of type ExcShapeFunctionNotPrimitive. In that case,
+   * use the shape_grad_grad_component() function.
    *
-   * An
-   * ExcUnitShapeValuesDoNotExist
-   * is thrown if the shape values
-   * of the FiniteElement under
-   * consideration depends on the
-   * shape of the cell in real
-   * space.
+   * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+   * FiniteElement under consideration depends on the shape of the cell in
+   * real space.
    */
   virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
                                          const Point<dim>   &p) const;
@@ -630,21 +521,14 @@ public:
                                                    const Point<dim>   &p,
                                                    const unsigned int component) const;
   /**
-   * Check for non-zero values on a
-   * face in order to optimize out
-   * matrix elements.
+   * Check for non-zero values on a face in order to optimize out matrix
+   * elements.
    *
-   * This function returns
-   * @p true, if the shape
-   * function @p shape_index has
-   * non-zero values on the face
-   * @p face_index.
+   * This function returns @p true, if the shape function @p shape_index has
+   * non-zero values on the face @p face_index.
    *
-   * A default implementation is
-   * provided in this basis class
-   * which always returns @p
-   * true. This is the safe way to
-   * go.
+   * A default implementation is provided in this basis class which always
+   * returns @p true. This is the safe way to go.
    */
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
@@ -656,374 +540,214 @@ public:
    */
 
   /**
-   * Projection from a fine grid
-   * space onto a coarse grid
-   * space. If this projection
-   * operator is associated with a
-   * matrix @p P, then the
-   * restriction of this matrix
-   * @p P_i to a single child cell
-   * is returned here.
-   *
-   * The matrix @p P is the
-   * concatenation or the sum of
-   * the cell matrices @p P_i,
-   * depending on the
-   * #restriction_is_additive_flags. This
-   * distinguishes interpolation
-   * (concatenation) and projection
-   * with respect to scalar
+   * Projection from a fine grid space onto a coarse grid space. If this
+   * projection operator is associated with a matrix @p P, then the
+   * restriction of this matrix @p P_i to a single child cell is returned
+   * here.
+   *
+   * The matrix @p P is the concatenation or the sum of the cell matrices @p
+   * P_i, depending on the #restriction_is_additive_flags. This distinguishes
+   * interpolation (concatenation) and projection with respect to scalar
    * products (summation).
    *
-   * Row and column indices are
-   * related to coarse grid and
-   * fine grid spaces,
-   * respectively, consistent with
-   * the definition of the
-   * associated operator.
-   *
-   * If projection matrices are not
-   * implemented in the derived
-   * finite element class, this
-   * function aborts with
-   * ExcProjectionVoid. You can
-   * check whether this is the case
-   * by calling the
-   * restriction_is_implemented()
-   * or the
-   * isotropic_restriction_is_implemented()
-   * function.
+   * Row and column indices are related to coarse grid and fine grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * If projection matrices are not implemented in the derived finite element
+   * class, this function aborts with ExcProjectionVoid. You can check whether
+   * this is the case by calling the restriction_is_implemented() or the
+   * isotropic_restriction_is_implemented() function.
    */
-  const FullMatrix<double> &
+  virtual const FullMatrix<double> &
   get_restriction_matrix (const unsigned int child,
                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
 
   /**
    * Embedding matrix between grids.
    *
-   * The identity operator from a
-   * coarse grid space into a fine
-   * grid space is associated with
-   * a matrix @p P. The
-   * restriction of this matrix @p P_i to
-   * a single child cell is
-   * returned here.
-   *
-   * The matrix @p P is the
-   * concatenation, not the sum of
-   * the cell matrices
-   * @p P_i. That is, if the same
-   * non-zero entry <tt>j,k</tt> exists
-   * in in two different child
-   * matrices @p P_i, the value
-   * should be the same in both
-   * matrices and it is copied into
-   * the matrix @p P only once.
-   *
-   * Row and column indices are
-   * related to fine grid and
-   * coarse grid spaces,
-   * respectively, consistent with
-   * the definition of the
-   * associated operator.
-   *
-   * These matrices are used by
-   * routines assembling the
-   * prolongation matrix for
-   * multi-level methods.  Upon
-   * assembling the transfer matrix
-   * between cells using this
-   * matrix array, zero elements in
-   * the prolongation matrix are
-   * discarded and will not fill up
-   * the transfer matrix.
-   *
-   * If projection matrices are not
-   * implemented in the derived
-   * finite element class, this
-   * function aborts with
-   * ExcEmbeddingVoid. You can
-   * check whether this is the case
-   * by calling the
-   * prolongation_is_implemented()
-   * or the
-   * isotropic_prolongation_is_implemented()
-   * function.
+   * The identity operator from a coarse grid space into a fine grid space is
+   * associated with a matrix @p P. The restriction of this matrix @p P_i to a
+   * single child cell is returned here.
+   *
+   * The matrix @p P is the concatenation, not the sum of the cell matrices @p
+   * P_i. That is, if the same non-zero entry <tt>j,k</tt> exists in in two
+   * different child matrices @p P_i, the value should be the same in both
+   * matrices and it is copied into the matrix @p P only once.
+   *
+   * Row and column indices are related to fine grid and coarse grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * These matrices are used by routines assembling the prolongation matrix
+   * for multi-level methods.  Upon assembling the transfer matrix between
+   * cells using this matrix array, zero elements in the prolongation matrix
+   * are discarded and will not fill up the transfer matrix.
+   *
+   * If projection matrices are not implemented in the derived finite element
+   * class, this function aborts with ExcEmbeddingVoid. You can check whether
+   * this is the case by calling the prolongation_is_implemented() or the
+   * isotropic_prolongation_is_implemented() function.
    */
-  const FullMatrix<double> &
+  virtual const FullMatrix<double> &
   get_prolongation_matrix (const unsigned int child,
                            const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
 
   /**
-   * Return whether this element implements
-   * its prolongation matrices. The return
-   * value also indicates whether a call to
-   * the get_prolongation_matrix()
-   * function will generate an error or
-   * not.
-   *
-   * Note, that this function
-   * returns <code>true</code> only
-   * if the prolongation matrices of
-   * the isotropic and all
-   * anisotropic refinement cases
-   * are implemented. If you are
-   * interested in the prolongation
-   * matrices for isotropic
-   * refinement only, use the
-   * isotropic_prolongation_is_implemented
-   * function instead.
-   *
-   * This function is mostly here in order
-   * to allow us to write more efficient
-   * test programs which we run on all
-   * kinds of weird elements, and for which
-   * we simply need to exclude certain
-   * tests in case something is not
-   * implemented. It will in general
-   * probably not be a great help in
-   * applications, since there is not much
-   * one can do if one needs these features
-   * and they are not implemented. This
-   * function could be used to check
-   * whether a call to
-   * <tt>get_prolongation_matrix()</tt> will
-   * succeed; however, one then still needs
-   * to cope with the lack of information
-   * this just expresses.
+   * Return whether this element implements its prolongation matrices. The
+   * return value also indicates whether a call to the
+   * get_prolongation_matrix() function will generate an error or not.
+   *
+   * Note, that this function returns <code>true</code> only if the
+   * prolongation matrices of the isotropic and all anisotropic refinement
+   * cases are implemented. If you are interested in the prolongation matrices
+   * for isotropic refinement only, use the
+   * isotropic_prolongation_is_implemented function instead.
+   *
+   * This function is mostly here in order to allow us to write more efficient
+   * test programs which we run on all kinds of weird elements, and for which
+   * we simply need to exclude certain tests in case something is not
+   * implemented. It will in general probably not be a great help in
+   * applications, since there is not much one can do if one needs these
+   * features and they are not implemented. This function could be used to
+   * check whether a call to <tt>get_prolongation_matrix()</tt> will succeed;
+   * however, one then still needs to cope with the lack of information this
+   * just expresses.
    */
   bool prolongation_is_implemented () const;
 
   /**
-   * Return whether this element implements
-   * its prolongation matrices for isotropic
-   * children. The return value also
-   * indicates whether a call to the @p
-   * get_prolongation_matrix function will
-   * generate an error or not.
-   *
-   * This function is mostly here in order
-   * to allow us to write more efficient
-   * test programs which we run on all
-   * kinds of weird elements, and for which
-   * we simply need to exclude certain
-   * tests in case something is not
-   * implemented. It will in general
-   * probably not be a great help in
-   * applications, since there is not much
-   * one can do if one needs these features
-   * and they are not implemented. This
-   * function could be used to check
-   * whether a call to
-   * <tt>get_prolongation_matrix()</tt> will
-   * succeed; however, one then still needs
-   * to cope with the lack of information
-   * this just expresses.
+   * Return whether this element implements its prolongation matrices for
+   * isotropic children. The return value also indicates whether a call to the
+   * @p get_prolongation_matrix function will generate an error or not.
+   *
+   * This function is mostly here in order to allow us to write more efficient
+   * test programs which we run on all kinds of weird elements, and for which
+   * we simply need to exclude certain tests in case something is not
+   * implemented. It will in general probably not be a great help in
+   * applications, since there is not much one can do if one needs these
+   * features and they are not implemented. This function could be used to
+   * check whether a call to <tt>get_prolongation_matrix()</tt> will succeed;
+   * however, one then still needs to cope with the lack of information this
+   * just expresses.
    */
   bool isotropic_prolongation_is_implemented () const;
 
   /**
-   * Return whether this element implements
-   * its restriction matrices. The return
-   * value also indicates whether a call to
-   * the get_restriction_matrix()
-   * function will generate an error or
-   * not.
-   *
-   * Note, that this function
-   * returns <code>true</code> only
-   * if the restriction matrices of
-   * the isotropic and all
-   * anisotropic refinement cases
-   * are implemented. If you are
-   * interested in the restriction
-   * matrices for isotropic
-   * refinement only, use the
-   * isotropic_restriction_is_implemented
-   * function instead.
-   *
-   * This function is mostly here in order
-   * to allow us to write more efficient
-   * test programs which we run on all
-   * kinds of weird elements, and for which
-   * we simply need to exclude certain
-   * tests in case something is not
-   * implemented. It will in general
-   * probably not be a great help in
-   * applications, since there is not much
-   * one can do if one needs these features
-   * and they are not implemented. This
-   * function could be used to check
-   * whether a call to
-   * <tt>get_restriction_matrix()</tt> will
-   * succeed; however, one then still needs
-   * to cope with the lack of information
-   * this just expresses.
+   * Return whether this element implements its restriction matrices. The
+   * return value also indicates whether a call to the
+   * get_restriction_matrix() function will generate an error or not.
+   *
+   * Note, that this function returns <code>true</code> only if the
+   * restriction matrices of the isotropic and all anisotropic refinement
+   * cases are implemented. If you are interested in the restriction matrices
+   * for isotropic refinement only, use the
+   * isotropic_restriction_is_implemented function instead.
+   *
+   * This function is mostly here in order to allow us to write more efficient
+   * test programs which we run on all kinds of weird elements, and for which
+   * we simply need to exclude certain tests in case something is not
+   * implemented. It will in general probably not be a great help in
+   * applications, since there is not much one can do if one needs these
+   * features and they are not implemented. This function could be used to
+   * check whether a call to <tt>get_restriction_matrix()</tt> will succeed;
+   * however, one then still needs to cope with the lack of information this
+   * just expresses.
    */
   bool restriction_is_implemented () const;
 
   /**
-   * Return whether this element implements
-   * its restriction matrices for isotropic
-   * children. The return value also
-   * indicates whether a call to the @p
-   * get_restriction_matrix function will
-   * generate an error or not.
-   *
-   * This function is mostly here in order
-   * to allow us to write more efficient
-   * test programs which we run on all
-   * kinds of weird elements, and for which
-   * we simply need to exclude certain
-   * tests in case something is not
-   * implemented. It will in general
-   * probably not be a great help in
-   * applications, since there is not much
-   * one can do if one needs these features
-   * and they are not implemented. This
-   * function could be used to check
-   * whether a call to
-   * <tt>get_restriction_matrix()</tt> will
-   * succeed; however, one then still needs
-   * to cope with the lack of information
-   * this just expresses.
+   * Return whether this element implements its restriction matrices for
+   * isotropic children. The return value also indicates whether a call to the
+   * @p get_restriction_matrix function will generate an error or not.
+   *
+   * This function is mostly here in order to allow us to write more efficient
+   * test programs which we run on all kinds of weird elements, and for which
+   * we simply need to exclude certain tests in case something is not
+   * implemented. It will in general probably not be a great help in
+   * applications, since there is not much one can do if one needs these
+   * features and they are not implemented. This function could be used to
+   * check whether a call to <tt>get_restriction_matrix()</tt> will succeed;
+   * however, one then still needs to cope with the lack of information this
+   * just expresses.
    */
   bool isotropic_restriction_is_implemented () const;
 
 
   /**
-   * Access the
-   * #restriction_is_additive_flags
-   * field. See there for more
+   * Access the #restriction_is_additive_flags field. See there for more
    * information on its contents.
    *
-   * The index must be between zero
-   * and the number of shape
-   * functions of this element.
+   * The index must be between zero and the number of shape functions of this
+   * element.
    */
   bool restriction_is_additive (const unsigned int index) const;
 
   /**
-   * Return a read only reference to
-   * the matrix which describes the
-   * constraints at the interface
-   * between a refined and an
-   * unrefined cell.
-   *
-   * The matrix is obviously empty
-   * in only one dimension,
-   * since there are no constraints
-   * then.
-   *
-   * Note that some finite elements
-   * do not (yet) implement hanging
-   * node constraints. If this is
-   * the case, then this function
-   * will generate an exception,
-   * since no useful return value
-   * can be generated. If you
-   * should have a way to live with
-   * this, then you might want to
-   * use the
-   * constraints_are_implemented()
-   * function to check up front
-   * whether this function will
-   * succeed or generate the
-   * exception.
+   * Return a read only reference to the matrix which describes the
+   * constraints at the interface between a refined and an unrefined cell.
+   *
+   * The matrix is obviously empty in only one dimension, since there are no
+   * constraints then.
+   *
+   * Note that some finite elements do not (yet) implement hanging node
+   * constraints. If this is the case, then this function will generate an
+   * exception, since no useful return value can be generated. If you should
+   * have a way to live with this, then you might want to use the
+   * constraints_are_implemented() function to check up front whether this
+   * function will succeed or generate the exception.
    */
   const FullMatrix<double> &constraints (const dealii::internal::SubfaceCase<dim> &subface_case=dealii::internal::SubfaceCase<dim>::case_isotropic) const;
 
   /**
-   * Return whether this element
-   * implements its hanging node
-   * constraints. The return value
-   * also indicates whether a call
-   * to the constraints() function
+   * Return whether this element implements its hanging node constraints. The
+   * return value also indicates whether a call to the constraints() function
    * will generate an error or not.
    *
-   * This function is mostly here
-   * in order to allow us to write
-   * more efficient test programs
-   * which we run on all kinds of
-   * weird elements, and for which
-   * we simply need to exclude
-   * certain tests in case hanging
-   * node constraints are not
-   * implemented. It will in
-   * general probably not be a
-   * great help in applications,
-   * since there is not much one
-   * can do if one needs hanging
-   * node constraints and they are
-   * not implemented. This function
-   * could be used to check whether
-   * a call to <tt>constraints()</tt>
-   * will succeed; however, one
-   * then still needs to cope with
-   * the lack of information this
+   * This function is mostly here in order to allow us to write more efficient
+   * test programs which we run on all kinds of weird elements, and for which
+   * we simply need to exclude certain tests in case hanging node constraints
+   * are not implemented. It will in general probably not be a great help in
+   * applications, since there is not much one can do if one needs hanging
+   * node constraints and they are not implemented. This function could be
+   * used to check whether a call to <tt>constraints()</tt> will succeed;
+   * however, one then still needs to cope with the lack of information this
    * just expresses.
    */
   bool constraints_are_implemented (const dealii::internal::SubfaceCase<dim> &subface_case=dealii::internal::SubfaceCase<dim>::case_isotropic) const;
 
 
   /**
-   * Return whether this element
-   * implements its hanging node
-   * constraints in the new way,
-   * which has to be used to make
-   * elements "hp compatible".
-   * That means, the element properly
-   * implements the
-   * get_face_interpolation_matrix
-   * and get_subface_interpolation_matrix
-   * methods. Therefore the return
-   * value also indicates whether a call
-   * to the get_face_interpolation_matrix()
-   * method and the get_subface_interpolation_matrix()
-   * method will generate an error or not.
-   *
-   * Currently the main purpose of this
-   * function is to allow the
-   * make_hanging_node_constraints method
-   * to decide whether the new procedures,
-   * which are supposed to work in the hp
-   * framework can be used, or if the old
-   * well verified but not hp capable
-   * functions should be used.  Once the
-   * transition to the new scheme for
-   * computing the interface constraints is
-   * complete, this function will be
-   * superfluous and will probably go away.
-   *
-   * Derived classes should implement this
-   * function accordingly. The default
-   * assumption is that a finite element
-   * does not provide hp capable face
-   * interpolation, and the default
-   * implementation therefore returns @p
-   * false.
+   * Return whether this element implements its hanging node constraints in
+   * the new way, which has to be used to make elements "hp compatible".  That
+   * means, the element properly implements the get_face_interpolation_matrix
+   * and get_subface_interpolation_matrix methods. Therefore the return value
+   * also indicates whether a call to the get_face_interpolation_matrix()
+   * method and the get_subface_interpolation_matrix() method will generate an
+   * error or not.
+   *
+   * Currently the main purpose of this function is to allow the
+   * make_hanging_node_constraints method to decide whether the new
+   * procedures, which are supposed to work in the hp framework can be used,
+   * or if the old well verified but not hp capable functions should be used.
+   * Once the transition to the new scheme for computing the interface
+   * constraints is complete, this function will be superfluous and will
+   * probably go away.
+   *
+   * Derived classes should implement this function accordingly. The default
+   * assumption is that a finite element does not provide hp capable face
+   * interpolation, and the default implementation therefore returns @p false.
    */
   virtual bool hp_constraints_are_implemented () const;
 
 
   /**
-   * Return the matrix
-   * interpolating from the given
-   * finite element to the present
-   * one. The size of the matrix is
-   * then #dofs_per_cell times
+   * Return the matrix interpolating from the given finite element to the
+   * present one. The size of the matrix is then #dofs_per_cell times
    * <tt>source.#dofs_per_cell</tt>.
    *
-   * Derived elements will have to
-   * implement this function. They
-   * may only provide interpolation
-   * matrices for certain source
-   * finite elements, for example
-   * those from the same family. If
-   * they don't implement
-   * interpolation from a given
-   * element, then they must throw
-   * an exception of type
+   * Derived elements will have to implement this function. They may only
+   * provide interpolation matrices for certain source finite elements, for
+   * example those from the same family. If they don't implement interpolation
+   * from a given element, then they must throw an exception of type
    * ExcInterpolationNotImplemented.
    */
   virtual void
@@ -1038,24 +762,14 @@ public:
 
 
   /**
-   * Return the matrix
-   * interpolating from a face of
-   * of one element to the face of
-   * the neighboring element.
-   * The size of the matrix is
-   * then <tt>source.#dofs_per_face</tt> times
-   * <tt>this->#dofs_per_face</tt>.
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.#dofs_per_face</tt> times <tt>this->#dofs_per_face</tt>.
    *
-   * Derived elements will have to
-   * implement this function. They
-   * may only provide interpolation
-   * matrices for certain source
-   * finite elements, for example
-   * those from the same family. If
-   * they don't implement
-   * interpolation from a given
-   * element, then they must throw
-   * an exception of type
+   * Derived elements will have to implement this function. They may only
+   * provide interpolation matrices for certain source finite elements, for
+   * example those from the same family. If they don't implement interpolation
+   * from a given element, then they must throw an exception of type
    * ExcInterpolationNotImplemented.
    */
   virtual void
@@ -1064,24 +778,14 @@ public:
 
 
   /**
-   * Return the matrix
-   * interpolating from a face of
-   * of one element to the subface of
-   * the neighboring element.
-   * The size of the matrix is
-   * then <tt>source.#dofs_per_face</tt> times
-   * <tt>this->#dofs_per_face</tt>.
-   *
-   * Derived elements will have to
-   * implement this function. They
-   * may only provide interpolation
-   * matrices for certain source
-   * finite elements, for example
-   * those from the same family. If
-   * they don't implement
-   * interpolation from a given
-   * element, then they must throw
-   * an exception of type
+   * Return the matrix interpolating from a face of of one element to the
+   * subface of the neighboring element.  The size of the matrix is then
+   * <tt>source.#dofs_per_face</tt> times <tt>this->#dofs_per_face</tt>.
+   *
+   * Derived elements will have to implement this function. They may only
+   * provide interpolation matrices for certain source finite elements, for
+   * example those from the same family. If they don't implement interpolation
+   * from a given element, then they must throw an exception of type
    * ExcInterpolationNotImplemented.
    */
   virtual void
@@ -1134,15 +838,11 @@ public:
   hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
 
   /**
-   * Return whether this element dominates
-   * the one given as argument when they
-   * meet at a common face,
-   * whether it is the other way around,
-   * whether neither dominates, or if
-   * either could dominate.
+   * Return whether this element dominates the one given as argument when they
+   * meet at a common face, whether it is the other way around, whether
+   * neither dominates, or if either could dominate.
    *
-   * For a definition of domination, see
-   * FiniteElementBase::Domination and in
+   * For a definition of domination, see FiniteElementBase::Domination and in
    * particular the @ref hp_paper "hp paper".
    */
   virtual
@@ -1152,21 +852,13 @@ public:
   //@}
 
   /**
-   * Comparison operator. We also
-   * check for equality of the
-   * constraint matrix, which is
-   * quite an expensive operation.
-   * Do therefore use this function
-   * with care, if possible only
-   * for debugging purposes.
+   * Comparison operator. We also check for equality of the constraint matrix,
+   * which is quite an expensive operation.  Do therefore use this function
+   * with care, if possible only for debugging purposes.
    *
-   * Since this function is not
-   * that important, we avoid an
-   * implementational question
-   * about comparing arrays and do
-   * not compare the matrix arrays
-   * #restriction and
-   * #prolongation.
+   * Since this function is not that important, we avoid an implementational
+   * question about comparing arrays and do not compare the matrix arrays
+   * #restriction and #prolongation.
    */
   bool operator == (const FiniteElement<dim,spacedim> &) const;
 
@@ -1175,116 +867,65 @@ public:
    * @{
    */
   /**
-   * Compute vector component and
-   * index of this shape function
-   * within the shape functions
-   * corresponding to this
-   * component from the index of a
-   * shape function within this
-   * finite element.
+   * Compute vector component and index of this shape function within the
+   * shape functions corresponding to this component from the index of a shape
+   * function within this finite element.
    *
-   * If the element is scalar, then
-   * the component is always zero,
-   * and the index within this
-   * component is equal to the
-   * overall index.
-   *
-   * If the shape function
-   * referenced has more than one
-   * non-zero component, then it
-   * cannot be associated with one
-   * vector component, and an
-   * exception of type
-   * ExcShapeFunctionNotPrimitive
-   * will be raised.
-   *
-   * Note that if the element is
-   * composed of other (base)
-   * elements, and a base element
-   * has more than one component
-   * but all its shape functions
-   * are primitive (i.e. are
-   * non-zero in only one
-   * component), then this mapping
-   * contains valid
-   * information. However, the
-   * index of a shape function of
-   * this element within one
-   * component (i.e. the second
-   * number of the respective entry
-   * of this array) does not
-   * indicate the index of the
-   * respective shape function
-   * within the base element (since
-   * that has more than one
-   * vector-component). For this
-   * information, refer to the
-   * #system_to_base_table field
-   * and the
-   * system_to_base_index()
-   * function.
+   * If the element is scalar, then the component is always zero, and the
+   * index within this component is equal to the overall index.
+   *
+   * If the shape function referenced has more than one non-zero component,
+   * then it cannot be associated with one vector component, and an exception
+   * of type ExcShapeFunctionNotPrimitive will be raised.
+   *
+   * Note that if the element is composed of other (base) elements, and a base
+   * element has more than one component but all its shape functions are
+   * primitive (i.e. are non-zero in only one component), then this mapping
+   * contains valid information. However, the index of a shape function of
+   * this element within one component (i.e. the second number of the
+   * respective entry of this array) does not indicate the index of the
+   * respective shape function within the base element (since that has more
+   * than one vector-component). For this information, refer to the
+   * #system_to_base_table field and the system_to_base_index() function.
    *
-   * The use of this function is
-   * explained extensively in the
-   * step-8 and @ref
-   * step_20 "step-20" tutorial
-   * programs as well as in the
-   * @ref vector_valued module.
+   * The use of this function is explained extensively in the step-8 and @ref
+   * step_20 "step-20" tutorial programs as well as in the @ref vector_valued
+   * module.
    */
   std::pair<unsigned int, unsigned int>
   system_to_component_index (const unsigned int index) const;
 
   /**
-   * Compute the shape function for
-   * the given vector component and
-   * index.
+   * Compute the shape function for the given vector component and index.
    *
-   * If the element is scalar, then
-   * the component must be zero,
-   * and the index within this
-   * component is equal to the
-   * overall index.
+   * If the element is scalar, then the component must be zero, and the index
+   * within this component is equal to the overall index.
    *
-   * This is the opposite operation
-   * from the system_to_component_index()
+   * This is the opposite operation from the system_to_component_index()
    * function.
    */
   unsigned int component_to_system_index(const unsigned int component,
                                          const unsigned int index) const;
 
   /**
-   * Same as
-   * system_to_component_index(),
-   * but do it for shape functions
-   * and their indices on a
-   * face. The range of allowed
-   * indices is therefore
+   * Same as system_to_component_index(), but do it for shape functions and
+   * their indices on a face. The range of allowed indices is therefore
    * 0..#dofs_per_face.
    *
-   * You will rarely need this
-   * function in application
-   * programs, since almost all
-   * application codes only need to
-   * deal with cell indices, not
-   * face indices. The function is
-   * mainly there for use inside
-   * the library.
+   * You will rarely need this function in application programs, since almost
+   * all application codes only need to deal with cell indices, not face
+   * indices. The function is mainly there for use inside the library.
    */
   std::pair<unsigned int, unsigned int>
   face_system_to_component_index (const unsigned int index) const;
 
   /**
-   * For faces with non-standard
-   * face_orientation in 3D, the dofs on
-   * faces (quads) have to be permuted in
-   * order to be combined with the correct
-   * shape functions. Given a local dof @p
-   * index on a quad, return the local index,
-   * if the face has non-standard
-   * face_orientation, face_flip or
-   * face_rotation. In 2D and 1D there is no
-   * need for permutation and consequently
-   * an exception is thrown.
+   * For faces with non-standard face_orientation in 3D, the dofs on faces
+   * (quads) have to be permuted in order to be combined with the correct
+   * shape functions. Given a local dof @p index on a quad, return the local
+   * index, if the face has non-standard face_orientation, face_flip or
+   * face_rotation. In 2D and 1D there is no need for permutation and
+   * consequently an exception is thrown.
    */
   unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index,
                                                            const bool face_orientation,
@@ -1292,268 +933,166 @@ public:
                                                            const bool face_rotation) const;
 
   /**
-   * For lines with non-standard
-   * line_orientation in 3D, the dofs on
-   * lines have to be permuted in order to be
-   * combined with the correct shape
-   * functions. Given a local dof @p index on
-   * a line, return the local index, if the
-   * line has non-standard
-   * line_orientation. In 2D and 1D there is
-   * no need for permutation, so the given
-   * index is simply returned.
+   * For lines with non-standard line_orientation in 3D, the dofs on lines
+   * have to be permuted in order to be combined with the correct shape
+   * functions. Given a local dof @p index on a line, return the local index,
+   * if the line has non-standard line_orientation. In 2D and 1D there is no
+   * need for permutation, so the given index is simply returned.
    */
   unsigned int adjust_line_dof_index_for_line_orientation (const unsigned int index,
                                                            const bool line_orientation) const;
 
   /**
-   * Return in which of the vector
-   * components of this finite
-   * element the @p ith shape
-   * function is non-zero. The
-   * length of the returned array
-   * is equal to the number of
-   * vector components of this
-   * element.
+   * Return in which of the vector components of this finite element the @p
+   * ith shape function is non-zero. The length of the returned array is equal
+   * to the number of vector components of this element.
+   *
+   * For most finite element spaces, the result of this function will be a
+   * vector with exactly one element being @p true, since for most spaces the
+   * individual vector components are independent. In that case, the component
+   * with the single zero is also the first element of what
+   * system_to_component_index() returns.
    *
-   * For most finite element
-   * spaces, the result of this
-   * function will be a vector with
-   * exactly one element being
-   * @p true, since for most
-   * spaces the individual vector
-   * components are independent. In
-   * that case, the component with
-   * the single zero is also the
-   * first element of what
-   * system_to_component_index()
-   * returns.
-   *
-   * Only for those spaces that couple the
-   * components, for example to make a
-   * shape function divergence free, will
-   * there be more than one @p true entry.
-   * Elements for which this is true are
-   * called non-primitive (see
-  * @ref GlossPrimitive).
+   * Only for those spaces that couple the components, for example to make a
+   * shape function divergence free, will there be more than one @p true
+   * entry.  Elements for which this is true are called non-primitive (see
+   * @ref GlossPrimitive).
    */
   const ComponentMask &
   get_nonzero_components (const unsigned int i) const;
 
   /**
-   * Return in how many vector
-   * components the @p ith shape
-   * function is non-zero. This
-   * value equals the number of
-   * entries equal to @p true in
-   * the result of the
-   * get_nonzero_components()
-   * function.
+   * Return in how many vector components the @p ith shape function is
+   * non-zero. This value equals the number of entries equal to @p true in the
+   * result of the get_nonzero_components() function.
    *
-   * For most finite element
-   * spaces, the result will be
-   * equal to one. It is not equal
-   * to one only for those ansatz
-   * spaces for which vector-valued
-   * shape functions couple the
-   * individual components, for
-   * example in order to make them
-   * divergence-free.
+   * For most finite element spaces, the result will be equal to one. It is
+   * not equal to one only for those ansatz spaces for which vector-valued
+   * shape functions couple the individual components, for example in order to
+   * make them divergence-free.
    */
   unsigned int
   n_nonzero_components (const unsigned int i) const;
 
   /**
-   * Return whether the @p ith
-   * shape function is primitive in
-   * the sense that the shape
-   * function is non-zero in only
-   * one vector
-   * component. Non-primitive shape
-   * functions would then, for
-   * example, be those of
-   * divergence free ansatz spaces,
-   * in which the individual vector
+   * Return whether the @p ith shape function is primitive in the sense that
+   * the shape function is non-zero in only one vector
+   * component. Non-primitive shape functions would then, for example, be
+   * those of divergence free ansatz spaces, in which the individual vector
    * components are coupled.
    *
-   * The result of the function is
-   * @p true if and only if the
-   * result of
-   * <tt>n_nonzero_components(i)</tt> is
-   * equal to one.
+   * The result of the function is @p true if and only if the result of
+   * <tt>n_nonzero_components(i)</tt> is equal to one.
    */
   bool
   is_primitive (const unsigned int i) const;
 
   /**
-   * Import function that is overloaded
-   * by the one above and would otherwise
+   * Import function that is overloaded by the one above and would otherwise
    * be hidden.
    */
   using FiniteElementData<dim>::is_primitive;
 
   /**
-   * Number of base elements in a
-   * mixed discretization.
-   *
-   * Note that even for vector
-   * valued finite elements, the
-   * number of components needs not
-   * coincide with the number of
-   * base elements, since they may
-   * be reused. For example, if you
-   * create a FESystem with
-   * three identical finite element
-   * classes by using the
-   * constructor that takes one
-   * finite element and a
-   * multiplicity, then the number
-   * of base elements is still one,
-   * although the number of
-   * components of the finite
-   * element is equal to the
-   * multiplicity.
+   * Number of base elements in a mixed discretization.
+   *
+   * Note that even for vector valued finite elements, the number of
+   * components needs not coincide with the number of base elements, since
+   * they may be reused. For example, if you create a FESystem with three
+   * identical finite element classes by using the constructor that takes one
+   * finite element and a multiplicity, then the number of base elements is
+   * still one, although the number of components of the finite element is
+   * equal to the multiplicity.
    */
   unsigned int n_base_elements () const;
 
   /**
-   * Access to base element
-   * objects. If the element is
-   * atomic, then
-   * <code>base_element(0)</code> is
-   * @p this.
+   * Access to base element objects. If the element is atomic, then
+   * <code>base_element(0)</code> is @p this.
    */
   virtual
   const FiniteElement<dim,spacedim> &
   base_element (const unsigned int index) const;
 
   /**
-   * This index denotes how often
-   * the base element @p index is
-   * used in a composed element. If
-   * the element is atomic, then
-   * the result is always equal to
-   * one. See the documentation for
-   * the n_base_elements()
-   * function for more details.
+   * This index denotes how often the base element @p index is used in a
+   * composed element. If the element is atomic, then the result is always
+   * equal to one. See the documentation for the n_base_elements() function
+   * for more details.
    */
   unsigned int
   element_multiplicity (const unsigned int index) const;
 
   /**
-   * Return for shape function
-   * @p index the base element it
-   * belongs to, the number of the
-   * copy of this base element
-   * (which is between zero and the
-   * multiplicity of this element),
-   * and the index of this shape
-   * function within this base
-   * element.
+   * Return for shape function @p index the base element it belongs to, the
+   * number of the copy of this base element (which is between zero and the
+   * multiplicity of this element), and the index of this shape function
+   * within this base element.
    *
-   * If the element is not composed of
-   * others, then base and instance
-   * are always zero, and the index
-   * is equal to the number of the
-   * shape function. If the element
-   * is composed of single
-   * instances of other elements
-   * (i.e. all with multiplicity
-   * one) all of which are scalar,
-   * then base values and dof
-   * indices within this element
-   * are equal to the
-   * #system_to_component_table. It
-   * differs only in case the
-   * element is composed of other
-   * elements and at least one of
-   * them is vector-valued itself.
-   *
-   * This function returns valid
-   * values also in the case of
-   * vector-valued
-   * (i.e. non-primitive) shape
-   * functions, in contrast to the
-   * system_to_component_index()
-   * function.
+   * If the element is not composed of others, then base and instance are
+   * always zero, and the index is equal to the number of the shape
+   * function. If the element is composed of single instances of other
+   * elements (i.e. all with multiplicity one) all of which are scalar, then
+   * base values and dof indices within this element are equal to the
+   * #system_to_component_table. It differs only in case the element is
+   * composed of other elements and at least one of them is vector-valued
+   * itself.
+   *
+   * This function returns valid values also in the case of vector-valued
+   * (i.e. non-primitive) shape functions, in contrast to the
+   * system_to_component_index() function.
    */
   std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
   system_to_base_index (const unsigned int index) const;
 
   /**
-   * Same as
-   * system_to_base_index(), but
-   * for degrees of freedom located
-   * on a face. The range of allowed
-   * indices is therefore
-   * 0..#dofs_per_face.
+   * Same as system_to_base_index(), but for degrees of freedom located on a
+   * face. The range of allowed indices is therefore 0..#dofs_per_face.
    *
-   * You will rarely need this
-   * function in application
-   * programs, since almost all
-   * application codes only need to
-   * deal with cell indices, not
-   * face indices. The function is
-   * mainly there for use inside
-   * the library.
+   * You will rarely need this function in application programs, since almost
+   * all application codes only need to deal with cell indices, not face
+   * indices. The function is mainly there for use inside the library.
    */
   std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
   face_system_to_base_index (const unsigned int index) const;
 
   /**
-   * Given a base element number,
-   * return the first block of a
-   * BlockVector it would generate.
+   * Given a base element number, return the first block of a BlockVector it
+   * would generate.
    */
   unsigned int first_block_of_base (const unsigned int b) const;
 
   /**
-   * For each vector component,
-   * return which base
-   * element implements this
-   * component and which vector
-   * component in this base element
-   * this is. This information is
-   * only of interest for
-   * vector-valued finite elements
-   * which are composed of several
-   * sub-elements. In that case,
-   * one may want to obtain
-   * information about the element
-   * implementing a certain vector
-   * component, which can be done
-   * using this function and the
-   * FESystem::base_element()
-   * function.
+   * For each vector component, return which base element implements this
+   * component and which vector component in this base element this is. This
+   * information is only of interest for vector-valued finite elements which
+   * are composed of several sub-elements. In that case, one may want to
+   * obtain information about the element implementing a certain vector
+   * component, which can be done using this function and the
+   * FESystem::base_element() function.
    *
-   * If this is a scalar finite
-   * element, then the return value
-   * is always equal to a pair of
-   * zeros.
+   * If this is a scalar finite element, then the return value is always equal
+   * to a pair of zeros.
    */
   std::pair<unsigned int, unsigned int>
   component_to_base_index (const unsigned int component) const;
 
 
   /**
-   * Return the base element for
-   * this block and the number of
-   * the copy of the base element.
+   * Return the base element for this block and the number of the copy of the
+   * base element.
    */
   std::pair<unsigned int,unsigned int>
   block_to_base_index (const unsigned int block) const;
 
   /**
-   * The vector block and the index
-   * inside the block for this
-   * shape function.
+   * The vector block and the index inside the block for this shape function.
    */
   std::pair<unsigned int,unsigned int>
   system_to_block_index (const unsigned int component) const;
 
   /**
-   * The vector block for this
-   * component.
+   * The vector block for this component.
    */
   unsigned int
   component_to_block_index (const unsigned int component) const;
@@ -1566,57 +1105,53 @@ public:
    */
 
   /**
-   * Return a component mask with as many elements as this
-   * object has vector components and of which exactly the
-   * one component is true that corresponds to the given
-   * argument. See @ref GlossComponentMask "the glossary"
-   *  for more information.
+   * Return a component mask with as many elements as this object has vector
+   * components and of which exactly the one component is true that
+   * corresponds to the given argument. See @ref GlossComponentMask "the
+   * glossary" for more information.
    *
-   * @param scalar An object that represents a single scalar
-   * vector component of this finite element.
-   * @return A component mask that is false in all components
-   * except for the one that corresponds to the argument.
+   * @param scalar An object that represents a single scalar vector component
+   * of this finite element.
+   * @return A component mask that is false in all components except for the
+   * one that corresponds to the argument.
    */
   ComponentMask
   component_mask (const FEValuesExtractors::Scalar &scalar) const;
 
   /**
-   * Return a component mask with as many elements as this
-   * object has vector components and of which exactly the
-   * <code>dim</code> components are true that correspond to the given
-   * argument. See @ref GlossComponentMask "the glossary"
-   * for more information.
+   * Return a component mask with as many elements as this object has vector
+   * components and of which exactly the <code>dim</code> components are true
+   * that correspond to the given argument. See @ref GlossComponentMask "the
+   * glossary" for more information.
    *
-   * @param vector An object that represents dim
-   * vector components of this finite element.
-   * @return A component mask that is false in all components
-   * except for the ones that corresponds to the argument.
+   * @param vector An object that represents dim vector components of this
+   * finite element.
+   * @return A component mask that is false in all components except for the
+   * ones that corresponds to the argument.
    */
   ComponentMask
   component_mask (const FEValuesExtractors::Vector &vector) const;
 
   /**
-   * Return a component mask with as many elements as this
-   * object has vector components and of which exactly the
-   * <code>dim*(dim+1)/2</code> components are true that
-   * correspond to the given argument. See @ref GlossComponentMask "the glossary"
-   * for more information.
+   * Return a component mask with as many elements as this object has vector
+   * components and of which exactly the <code>dim*(dim+1)/2</code> components
+   * are true that correspond to the given argument. See @ref
+   * GlossComponentMask "the glossary" for more information.
    *
-   * @param sym_tensor An object that represents dim*(dim+1)/2
-   * components of this finite element that are jointly to be
-   * interpreted as forming a symmetric tensor.
-   * @return A component mask that is false in all components
-   * except for the ones that corresponds to the argument.
+   * @param sym_tensor An object that represents dim*(dim+1)/2 components of
+   * this finite element that are jointly to be interpreted as forming a
+   * symmetric tensor.  @return A component mask that is false in all
+   * components except for the ones that corresponds to the argument.
    */
   ComponentMask
   component_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
 
   /**
    * Given a block mask (see @ref GlossBlockMask "this glossary entry"),
-   * produce a component mask (see @ref GlossComponentMask "this glossary entry")
-   * that represents the components that correspond to the blocks selected in
-   * the input argument. This is essentially a conversion operator from
-   * BlockMask to ComponentMask.
+   * produce a component mask (see @ref GlossComponentMask "this glossary
+   * entry") that represents the components that correspond to the blocks
+   * selected in the input argument. This is essentially a conversion operator
+   * from BlockMask to ComponentMask.
    *
    * @param block_mask The mask that selects individual blocks of the finite
    * element
@@ -1627,84 +1162,80 @@ public:
   component_mask (const BlockMask &block_mask) const;
 
   /**
-   * Return a block mask with as many elements as this
-   * object has blocks and of which exactly the
-   * one component is true that corresponds to the given
-   * argument. See @ref GlossBlockMask "the glossary"
-   * for more information.
+   * Return a block mask with as many elements as this object has blocks and
+   * of which exactly the one component is true that corresponds to the given
+   * argument. See @ref GlossBlockMask "the glossary" for more information.
    *
-   * @note This function will only succeed if the scalar referenced
-   * by the argument encompasses a complete block. In other words,
-   * if, for example, you pass an extractor for the single
-   * $x$ velocity and this object represents an FE_RaviartThomas
-   * object, then the single scalar object you selected is part
-   * of a larger block and consequently there is no block mask that
-   * would represent it. The function will then produce an exception.
+   * @note This function will only succeed if the scalar referenced by the
+   * argument encompasses a complete block. In other words, if, for example,
+   * you pass an extractor for the single $x$ velocity and this object
+   * represents an FE_RaviartThomas object, then the single scalar object you
+   * selected is part of a larger block and consequently there is no block
+   * mask that would represent it. The function will then produce an
+   * exception.
    *
-   * @param scalar An object that represents a single scalar
-   * vector component of this finite element.
-   * @return A component mask that is false in all components
-   * except for the one that corresponds to the argument.
+   * @param scalar An object that represents a single scalar vector component
+   * of this finite element.
+   * @return A component mask that is false in all components except for the
+   * one that corresponds to the argument.
    */
   BlockMask
   block_mask (const FEValuesExtractors::Scalar &scalar) const;
 
   /**
-   * Return a component mask with as many elements as this
-   * object has vector components and of which exactly the
-   * <code>dim</code> components are true that correspond to the given
-   * argument. See @ref GlossBlockMask "the glossary"
-   * for more information.
+   * Return a component mask with as many elements as this object has vector
+   * components and of which exactly the <code>dim</code> components are true
+   * that correspond to the given argument. See @ref GlossBlockMask "the
+   * glossary" for more information.
    *
    * @note The same caveat applies as to the version of the function above:
-   * The extractor object passed as argument must be so that it corresponds
-   * to full blocks and does not split blocks of this element.
+   * The extractor object passed as argument must be so that it corresponds to
+   * full blocks and does not split blocks of this element.
    *
-   * @param vector An object that represents dim
-   * vector components of this finite element.
-   * @return A component mask that is false in all components
-   * except for the ones that corresponds to the argument.
+   * @param vector An object that represents dim vector components of this
+   * finite element.
+   * @return A component mask that is false in all components except for the
+   * ones that corresponds to the argument.
    */
   BlockMask
   block_mask (const FEValuesExtractors::Vector &vector) const;
 
   /**
-   * Return a component mask with as many elements as this
-   * object has vector components and of which exactly the
-   * <code>dim*(dim+1)/2</code> components are true that
-   * correspond to the given argument. See @ref GlossBlockMask "the glossary"
-   * for more information.
+   * Return a component mask with as many elements as this object has vector
+   * components and of which exactly the <code>dim*(dim+1)/2</code> components
+   * are true that correspond to the given argument. See @ref GlossBlockMask
+   * "the glossary" for more information.
    *
    * @note The same caveat applies as to the version of the function above:
-   * The extractor object passed as argument must be so that it corresponds
-   * to full blocks and does not split blocks of this element.
+   * The extractor object passed as argument must be so that it corresponds to
+   * full blocks and does not split blocks of this element.
    *
-   * @param sym_tensor An object that represents dim*(dim+1)/2
-   * components of this finite element that are jointly to be
-   * interpreted as forming a symmetric tensor.
-   * @return A component mask that is false in all components
-   * except for the ones that corresponds to the argument.
+   * @param sym_tensor An object that represents dim*(dim+1)/2 components of
+   * this finite element that are jointly to be interpreted as forming a
+   * symmetric tensor.
+   * @return A component mask that is false in all components except for the
+   * ones that corresponds to the argument.
    */
   BlockMask
   block_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
 
   /**
-   * Given a component mask (see @ref GlossComponentMask "this glossary entry"),
-   * produce a block mask (see @ref GlossBlockMask "this glossary entry")
-   * that represents the blocks that correspond to the components selected in
-   * the input argument. This is essentially a conversion operator from
-   * ComponentMask to BlockMask.
+   * Given a component mask (see @ref GlossComponentMask "this glossary
+   * entry"), produce a block mask (see @ref GlossBlockMask "this glossary
+   * entry") that represents the blocks that correspond to the components
+   * selected in the input argument. This is essentially a conversion operator
+   * from ComponentMask to BlockMask.
    *
-   * @note This function will only succeed if the components referenced
-   * by the argument encompasses complete blocks. In other words,
-   * if, for example, you pass an component mask for the single
-   * $x$ velocity and this object represents an FE_RaviartThomas
-   * object, then the single component you selected is part
-   * of a larger block and consequently there is no block mask that
-   * would represent it. The function will then produce an exception.
+   * @note This function will only succeed if the components referenced by the
+   * argument encompasses complete blocks. In other words, if, for example,
+   * you pass an component mask for the single $x$ velocity and this object
+   * represents an FE_RaviartThomas object, then the single component you
+   * selected is part of a larger block and consequently there is no block
+   * mask that would represent it. The function will then produce an
+   * exception.
    *
-   * @param component_mask The mask that selects individual components of the finite
-   * element
+   * @param component_mask The mask that selects individual components of the
+   * finite element
    * @return A mask that selects those blocks corresponding to the selected
    * blocks of the input argument.
    */
@@ -1719,199 +1250,118 @@ public:
    */
 
   /**
-   * Return the support points of
-   * the trial functions on the
-   * unit cell, if the derived
-   * finite element defines some.
-   * Finite elements that allow
-   * some kind of interpolation
-   * operation usually have support
-   * points. On the other hand,
-   * elements that define their
-   * degrees of freedom by, for
-   * example, moments on faces, or
-   * as derivatives, don't have
-   * support points. In that case,
-   * the returned field is empty.
-   *
-   * If the finite element defines
-   * support points, then their
-   * number equals the number of
-   * degrees of freedom of the
-   * element.  The order of points
-   * in the array matches that
-   * returned by the
-   * <tt>cell->get_dof_indices</tt>
+   * Return the support points of the trial functions on the unit cell, if the
+   * derived finite element defines some.  Finite elements that allow some
+   * kind of interpolation operation usually have support points. On the other
+   * hand, elements that define their degrees of freedom by, for example,
+   * moments on faces, or as derivatives, don't have support points. In that
+   * case, the returned field is empty.
+   *
+   * If the finite element defines support points, then their number equals
+   * the number of degrees of freedom of the element.  The order of points in
+   * the array matches that returned by the <tt>cell->get_dof_indices</tt>
    * function.
    *
-   * See the class documentation
-   * for details on support points.
+   * See the class documentation for details on support points.
    */
   const std::vector<Point<dim> > &
   get_unit_support_points () const;
 
   /**
-   * Return whether a finite
-   * element has defined support
-   * points. If the result is true,
-   * then a call to the
-   * get_unit_support_points()
-   * yields a non-empty array.
+   * Return whether a finite element has defined support points. If the result
+   * is true, then a call to the get_unit_support_points() yields a non-empty
+   * array.
    *
-   * The result may be false if an
-   * element is not defined by
-   * interpolating shape functions,
-   * for example by P-elements on
-   * quadrilaterals. It will
-   * usually only be true if the
-   * element constructs its shape
-   * functions by the requirement
-   * that they be one at a certain
-   * point and zero at all the
-   * points associated with the
-   * other shape functions.
-   *
-   * In composed elements (i.e. for
-   * the FESystem class, the
-   * result will be true if all all
-   * the base elements have defined
-   * support points.
+   * The result may be false if an element is not defined by interpolating
+   * shape functions, for example by P-elements on quadrilaterals. It will
+   * usually only be true if the element constructs its shape functions by the
+   * requirement that they be one at a certain point and zero at all the
+   * points associated with the other shape functions.
+   *
+   * In composed elements (i.e. for the FESystem class, the result will be
+   * true if all all the base elements have defined support points.
    */
   bool has_support_points () const;
 
   /**
-   * Return the position of the
-   * support point of the
-   * @p indexth shape function. If
-   * it does not exist, raise an
-   * exception.
+   * Return the position of the support point of the @p indexth shape
+   * function. If it does not exist, raise an exception.
    *
-   * The default implementation
-   * simply returns the respective
-   * element from the array you get
-   * from
-   * get_unit_support_points(),
-   * but derived elements may
-   * overload this function. In
-   * particular, note that the
-   * FESystem class overloads
-   * it so that it can return the
-   * support points of individual
-   * base elements, if not all the
-   * base elements define support
-   * points. In this way, you can
-   * still ask for certain support
-   * points, even if
-   * get_unit_support_points()
-   * only returns an empty array.
+   * The default implementation simply returns the respective element from the
+   * array you get from get_unit_support_points(), but derived elements may
+   * overload this function. In particular, note that the FESystem class
+   * overloads it so that it can return the support points of individual base
+   * elements, if not all the base elements define support points. In this
+   * way, you can still ask for certain support points, even if
+   * get_unit_support_points() only returns an empty array.
    */
   virtual
   Point<dim>
   unit_support_point (const unsigned int index) const;
 
   /**
-   * Return the support points of
-   * the trial functions on the
-   * unit face, if the derived
-   * finite element defines some.
-   * Finite elements that allow
-   * some kind of interpolation
-   * operation usually have support
-   * points. On the other hand,
-   * elements that define their
-   * degrees of freedom by, for
-   * example, moments on faces, or
-   * as derivatives, don't have
-   * support points. In that case,
-   * the returned field is empty
-   *
-   * Note that elements that have
-   * support points need not
-   * necessarily have some on the
-   * faces, even if the
-   * interpolation points are
-   * located physically on a
-   * face. For example, the
-   * discontinuous elements have
-   * interpolation points on the
-   * vertices, and for higher
-   * degree elements also on the
-   * faces, but they are not
-   * defined to be on faces since
-   * in that case degrees of
-   * freedom from both sides of a
-   * face (or from all adjacent
-   * elements to a vertex) would be
-   * identified with each other,
-   * which is not what we would
-   * like to have). Logically,
-   * these degrees of freedom are
-   * therefore defined to belong to
-   * the cell, rather than the face
-   * or vertex. In that case, the
-   * returned element would
-   * therefore have length zero.
-   *
-   * If the finite element defines
-   * support points, then their
-   * number equals the number of
-   * degrees of freedom on the face
-   * (#dofs_per_face). The order
-   * of points in the array matches
-   * that returned by the
-   * <tt>cell->get_dof_indices</tt>
-   * function.
-   *
-   * See the class documentation
-   * for details on support points.
+   * Return the support points of the trial functions on the unit face, if the
+   * derived finite element defines some.  Finite elements that allow some
+   * kind of interpolation operation usually have support points. On the other
+   * hand, elements that define their degrees of freedom by, for example,
+   * moments on faces, or as derivatives, don't have support points. In that
+   * case, the returned field is empty
+   *
+   * Note that elements that have support points need not necessarily have
+   * some on the faces, even if the interpolation points are located
+   * physically on a face. For example, the discontinuous elements have
+   * interpolation points on the vertices, and for higher degree elements also
+   * on the faces, but they are not defined to be on faces since in that case
+   * degrees of freedom from both sides of a face (or from all adjacent
+   * elements to a vertex) would be identified with each other, which is not
+   * what we would like to have). Logically, these degrees of freedom are
+   * therefore defined to belong to the cell, rather than the face or
+   * vertex. In that case, the returned element would therefore have length
+   * zero.
+   *
+   * If the finite element defines support points, then their number equals
+   * the number of degrees of freedom on the face (#dofs_per_face). The order
+   * of points in the array matches that returned by the
+   * <tt>cell->get_dof_indices</tt> function.
+   *
+   * See the class documentation for details on support points.
    */
   const std::vector<Point<dim-1> > &
   get_unit_face_support_points () const;
 
   /**
-   * Return whether a finite
-   * element has defined support
-   * points on faces. If the result
-   * is true, then a call to the
-   * get_unit_face_support_points()
+   * Return whether a finite element has defined support points on faces. If
+   * the result is true, then a call to the get_unit_face_support_points()
    * yields a non-empty array.
    *
-   * For more information, see the
-   * documentation for the
-   * has_support_points()
+   * For more information, see the documentation for the has_support_points()
    * function.
    */
   bool has_face_support_points () const;
 
   /**
-   * The function corresponding to
-   * the unit_support_point()
-   * function, but for faces. See
-   * there for more information.
+   * The function corresponding to the unit_support_point() function, but for
+   * faces. See there for more information.
    */
   virtual
   Point<dim-1>
   unit_face_support_point (const unsigned int index) const;
 
   /**
-   * Return a support point vector
-   * for generalized interpolation.
-  *
-  * See the @ref GlossGeneralizedSupport "glossary entry on generalized support points"
-  * for more information.
+   * Return a support point vector for generalized interpolation.
+   *
+   * See the @ref GlossGeneralizedSupport "glossary entry on generalized
+   * support points" for more information.
    */
   const std::vector<Point<dim> > &
   get_generalized_support_points () const;
 
   /**
-   * Returns <tt>true</tt> if the
-   * class provides nonempty
-   * vectors either from
-   * get_unit_support_points() or
-   * get_generalized_support_points().
-  *
-  * See the @ref GlossGeneralizedSupport "glossary entry on generalized support points"
-  * for more information.
+   * Returns <tt>true</tt> if the class provides nonempty vectors either from
+   * get_unit_support_points() or get_generalized_support_points().
+   *
+   * See the @ref GlossGeneralizedSupport "glossary entry on generalized
+   * support points" for more information.
    */
   bool has_generalized_support_points () const;
 
@@ -1922,61 +1372,43 @@ public:
   get_generalized_face_support_points () const;
 
   /**
-   * Return whether a finite
-   * element has defined
-   * generalized support
-   * points on faces. If the result
-   * is true, then a call to the
-   * get_generalized_face_support_points
-   * yields a non-empty array.
+   * Return whether a finite element has defined generalized support points on
+   * faces. If the result is true, then a call to the
+   * get_generalized_face_support_points yields a non-empty array.
    *
-   * For more information, see the
-   * documentation for the
-   * has_support_points()
+   * For more information, see the documentation for the has_support_points()
    * function.
    */
   bool has_generalized_face_support_points () const;
 
   /**
-   * Interpolate a set of scalar
-   * values, computed in the
-   * generalized support points.
+   * Interpolate a set of scalar values, computed in the generalized support
+   * points.
    *
-   * @note This function is
-   * implemented in
-   * FiniteElement for the case
-   * that the element has support
-   * points. In this case, the
-   * resulting coefficients are
-   * just the values in the suport
-   * points. All other elements
-   * must reimplement it.
+   * @note This function is implemented in FiniteElement for the case that the
+   * element has support points. In this case, the resulting coefficients are
+   * just the values in the suport points. All other elements must reimplement
+   * it.
    */
   virtual void interpolate(std::vector<double>       &local_dofs,
                            const std::vector<double> &values) const;
 
   /**
-   * Interpolate a set of vector
-   * values, computed in the
-   * generalized support points.
+   * Interpolate a set of vector values, computed in the generalized support
+   * points.
    *
-   * Since a finite element often
-   * only interpolates part of a
-   * vector, <tt>offset</tt> is
-   * used to determine the first
-   * component of the vector to be
-   * interpolated. Maybe consider
-   * changing your data structures
-   * to use the next function.
+   * Since a finite element often only interpolates part of a vector,
+   * <tt>offset</tt> is used to determine the first component of the vector to
+   * be interpolated. Maybe consider changing your data structures to use the
+   * next function.
    */
   virtual void interpolate(std::vector<double>                &local_dofs,
                            const std::vector<Vector<double> > &values,
                            unsigned int offset = 0) const;
 
   /**
-   * Interpolate a set of vector
-   * values, computed in the
-   * generalized support points.
+   * Interpolate a set of vector values, computed in the generalized support
+   * points.
    */
   virtual void interpolate(
     std::vector<double> &local_dofs,
@@ -1985,15 +1417,12 @@ public:
   //@}
 
   /**
-   * Determine an estimate for the
-   * memory consumption (in bytes)
-   * of this object.
+   * Determine an estimate for the memory consumption (in bytes) of this
+   * object.
    *
-   * This function is made virtual,
-   * since finite element objects
-   * are usually accessed through
-   * pointers to their base class,
-   * rather than the class itself.
+   * This function is made virtual, since finite element objects are usually
+   * accessed through pointers to their base class, rather than the class
+   * itself.
    */
   virtual std::size_t memory_consumption () const;
   /**
@@ -2023,29 +1452,24 @@ public:
   DeclException0 (ExcUnitShapeValuesDoNotExist);
 
   /**
-   * Attempt to access support
-   * points of a finite element
-   * which is not Lagrangian.
+   * Attempt to access support points of a finite element which is not
+   * Lagrangian.
    *
    * @ingroup Exceptions
    */
   DeclException0 (ExcFEHasNoSupportPoints);
 
   /**
-   * Attempt to access embedding
-   * matrices of a finite element
-   * which did not implement these
-   * matrices.
+   * Attempt to access embedding matrices of a finite element which did not
+   * implement these matrices.
    *
    * @ingroup Exceptions
    */
   DeclException0 (ExcEmbeddingVoid);
 
   /**
-   * Attempt to access restriction
-   * matrices of a finite element
-   * which did not implement these
-   * matrices.
+   * Attempt to access restriction matrices of a finite element which did not
+   * implement these matrices.
    *
    * Exception
    * @ingroup Exceptions
@@ -2053,10 +1477,8 @@ public:
   DeclException0 (ExcProjectionVoid);
 
   /**
-   * Attempt to access constraint
-   * matrices of a finite element
-   * which did not implement these
-   * matrices.
+   * Attempt to access constraint matrices of a finite element which did not
+   * implement these matrices.
    *
    * Exception
    * @ingroup Exceptions
@@ -2102,235 +1524,140 @@ public:
 protected:
 
   /**
-   * Reinit the vectors of
-   * restriction and prolongation
-   * matrices to the right sizes:
-   * For every refinement case,
-   * except for
-   * RefinementCase::no_refinement,
-   * and for every child of that
-   * refinement case the space of
-   * one restriction and
-   * prolongation matrix is
-   * allocated, see the
-   * documentation of the
-   * restriction and prolongation
-   * vectors for more detail on the
-   * actual vector sizes.
-   *
-   * @param
-   * isotropic_restriction_only
-   * only the restriction matrices
-   * required for isotropic
-   * refinement are reinited to the
-   * right size.
-   * @param
-   * isotropic_prolongation_only
-   * only the prolongation matrices
-   * required for isotropic
-   * refinement are reinited to the
-   * right size.
+   * Reinit the vectors of restriction and prolongation matrices to the right
+   * sizes: For every refinement case, except for
+   * RefinementCase::no_refinement, and for every child of that refinement
+   * case the space of one restriction and prolongation matrix is allocated,
+   * see the documentation of the restriction and prolongation vectors for
+   * more detail on the actual vector sizes.
+   *
+   * @param isotropic_restriction_only only the restriction matrices required
+   * for isotropic refinement are reinited to the right size.
+   * @param isotropic_prolongation_only only the prolongation matrices
+   * required for isotropic refinement are reinited to the right size.
    */
   void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false,
                                                     const bool isotropic_prolongation_only=false);
 
   /**
-   * Vector of projection
-   * matrices. See
-   * get_restriction_matrix()
-   * above. The constructor
-   * initializes these matrices to
-   * zero dimensions, which can be
-   * changed by derived classes
-   * implementing them.
-   *
-   * Note, that
-   * <code>restriction[refinement_case-1][child]</code>
-   * includes the restriction
-   * matrix of child
-   * <code>child</code> for the
-   * RefinementCase
-   * <code>refinement_case</code>. Here,
-   * we use
-   * <code>refinement_case-1</code>
-   * instead of
-   * <code>refinement_case</code>
-   * as for
-   * RefinementCase::no_refinement(=0)
-   * there are no restriction
-   * matrices available.
+   * Vector of projection matrices. See get_restriction_matrix() above. The
+   * constructor initializes these matrices to zero dimensions, which can be
+   * changed by derived classes implementing them.
+   *
+   * Note, that <code>restriction[refinement_case-1][child]</code> includes
+   * the restriction matrix of child <code>child</code> for the RefinementCase
+   * <code>refinement_case</code>. Here, we use <code>refinement_case-1</code>
+   * instead of <code>refinement_case</code> as for
+   * RefinementCase::no_refinement(=0) there are no restriction matrices
+   * available.
    */
   std::vector<std::vector<FullMatrix<double> > > restriction;
 
   /**
-   * Vector of embedding
-   * matrices. See
-   * <tt>get_prolongation_matrix()</tt>
-   * above. The constructor
-   * initializes these matrices to
-   * zero dimensions, which can be
-   * changed by derived classes
-   * implementing them.
-   *
-   * Note, that
-   * <code>prolongation[refinement_case-1][child]</code>
-   * includes the prolongation
-   * matrix of child
-   * <code>child</code> for the
-   * RefinementCase
-   * <code>refinement_case</code>. Here,
-   * we use
-   * <code>refinement_case-1</code>
-   * instead of
-   * <code>refinement_case</code>
-   * as for
-   * RefinementCase::no_refinement(=0)
-   * there are no prolongation
-   * matrices available.
+   * Vector of embedding matrices. See <tt>get_prolongation_matrix()</tt>
+   * above. The constructor initializes these matrices to zero dimensions,
+   * which can be changed by derived classes implementing them.
+   *
+   * Note, that <code>prolongation[refinement_case-1][child]</code> includes
+   * the prolongation matrix of child <code>child</code> for the
+   * RefinementCase <code>refinement_case</code>. Here, we use
+   * <code>refinement_case-1</code> instead of <code>refinement_case</code> as
+   * for RefinementCase::no_refinement(=0) there are no prolongation matrices
+   * available.
    */
   std::vector<std::vector<FullMatrix<double> > > prolongation;
 
   /**
-   * Specify the constraints which
-   * the dofs on the two sides of a
-   * cell interface underly if the
-   * line connects two cells of
-   * which one is refined once.
+   * Specify the constraints which the dofs on the two sides of a cell
+   * interface underly if the line connects two cells of which one is refined
+   * once.
    *
-   * For further details see the
-   * general description of the
-   * derived class.
+   * For further details see the general description of the derived class.
    *
-   * This field is obviously
-   * useless in one dimension
-   * and has there a zero size.
+   * This field is obviously useless in one dimension and has there a zero
+   * size.
    */
   FullMatrix<double> interface_constraints;
 
   /**
-   * List of support points on the
-   * unit cell, in case the finite
-   * element has any. The
-   * constructor leaves this field
-   * empty, derived classes may
-   * write in some contents.
-   *
-   * Finite elements that allow
-   * some kind of interpolation
-   * operation usually have support
-   * points. On the other hand,
-   * elements that define their
-   * degrees of freedom by, for
-   * example, moments on faces, or
-   * as derivatives, don't have
-   * support points. In that case,
-   * this field remains empty.
+   * List of support points on the unit cell, in case the finite element has
+   * any. The constructor leaves this field empty, derived classes may write
+   * in some contents.
+   *
+   * Finite elements that allow some kind of interpolation operation usually
+   * have support points. On the other hand, elements that define their
+   * degrees of freedom by, for example, moments on faces, or as derivatives,
+   * don't have support points. In that case, this field remains empty.
    */
   std::vector<Point<dim> > unit_support_points;
 
   /**
-   * Same for the faces. See the
-   * description of the
-   * get_unit_face_support_points()
-   * function for a discussion of
-   * what contributes a face
-   * support point.
+   * Same for the faces. See the description of the
+   * get_unit_face_support_points() function for a discussion of what
+   * contributes a face support point.
    */
   std::vector<Point<dim-1> > unit_face_support_points;
 
   /**
-   * Support points used for
-   * interpolation functions of
-   * non-Lagrangian elements.
+   * Support points used for interpolation functions of non-Lagrangian
+   * elements.
    */
   std::vector<Point<dim> > generalized_support_points;
 
   /**
-   * Face support points used for
-   * interpolation functions of
-   * non-Lagrangian elements.
+   * Face support points used for interpolation functions of non-Lagrangian
+   * elements.
    */
   std::vector<Point<dim-1> > generalized_face_support_points;
 
   /**
-   * For faces with non-standard
-   * face_orientation in 3D, the dofs on
-   * faces (quads) have to be permuted in
-   * order to be combined with the correct
-   * shape functions. Given a local dof @p
-   * index on a quad, return the shift in the
-   * local index, if the face has
-   * non-standard face_orientation,
-   * i.e. <code>old_index + shift =
-   * new_index</code>. In 2D and 1D there is
-   * no need for permutation so the vector is
-   * empty. In 3D it has the size of <code>
-   * #dofs_per_quad * 8 </code>, where 8 is
-   * the number of orientations, a face can
-   * be in (all combinations of the three
-   * bool flags face_orientation, face_flip
-   * and face_rotation).
-   *
-   * The standard implementation fills this
-   * with zeros, i.e. no permuatation at
-   * all. Derived finite element classes have
-   * to fill this Table with the correct
-   * values.
+   * For faces with non-standard face_orientation in 3D, the dofs on faces
+   * (quads) have to be permuted in order to be combined with the correct
+   * shape functions. Given a local dof @p index on a quad, return the shift
+   * in the local index, if the face has non-standard face_orientation,
+   * i.e. <code>old_index + shift = new_index</code>. In 2D and 1D there is no
+   * need for permutation so the vector is empty. In 3D it has the size of
+   * <code> #dofs_per_quad * 8 </code>, where 8 is the number of orientations,
+   * a face can be in (all combinations of the three bool flags
+   * face_orientation, face_flip and face_rotation).
+   *
+   * The standard implementation fills this with zeros, i.e. no permuatation
+   * at all. Derived finite element classes have to fill this Table with the
+   * correct values.
    */
   Table<2,int> adjust_quad_dof_index_for_face_orientation_table;
 
   /**
-   * For lines with non-standard
-   * line_orientation in 3D, the dofs on
-   * lines have to be permuted in
-   * order to be combined with the correct
-   * shape functions. Given a local dof @p
-   * index on a line, return the shift in the
-   * local index, if the line has
-   * non-standard line_orientation,
-   * i.e. <code>old_index + shift =
-   * new_index</code>. In 2D and 1D there is
-   * no need for permutation so the vector is
-   * empty. In 3D it has the size of
+   * For lines with non-standard line_orientation in 3D, the dofs on lines
+   * have to be permuted in order to be combined with the correct shape
+   * functions. Given a local dof @p index on a line, return the shift in the
+   * local index, if the line has non-standard line_orientation,
+   * i.e. <code>old_index + shift = new_index</code>. In 2D and 1D there is no
+   * need for permutation so the vector is empty. In 3D it has the size of
    * #dofs_per_line.
    *
-   * The standard implementation fills this
-   * with zeros, i.e. no permutation at
-   * all. Derived finite element classes have
-   * to fill this vector with the correct
-   * values.
+   * The standard implementation fills this with zeros, i.e. no permutation at
+   * all. Derived finite element classes have to fill this vector with the
+   * correct values.
    */
   std::vector<int> adjust_line_dof_index_for_line_orientation_table;
 
   /**
-   * Return the size of interface
-   * constraint matrices. Since
-   * this is needed in every
-   * derived finite element class
-   * when initializing their size,
-   * it is placed into this
-   * function, to avoid having to
-   * recompute the
-   * dimension-dependent size of
-   * these matrices each time.
-   *
-   * Note that some elements do not
-   * implement the interface
-   * constraints for certain
-   * polynomial degrees. In this
-   * case, this function still
-   * returns the size these
-   * matrices should have when
-   * implemented, but the actual
-   * matrices are empty.
+   * Return the size of interface constraint matrices. Since this is needed in
+   * every derived finite element class when initializing their size, it is
+   * placed into this function, to avoid having to recompute the
+   * dimension-dependent size of these matrices each time.
+   *
+   * Note that some elements do not implement the interface constraints for
+   * certain polynomial degrees. In this case, this function still returns the
+   * size these matrices should have when implemented, but the actual matrices
+   * are empty.
    */
   TableIndices<2>
   interface_constraints_size () const;
 
   /**
-   * Compute second derivatives by
-   * finite differences of
-   * gradients.
+   * Compute second derivatives by finite differences of gradients.
    */
   void compute_2nd (const Mapping<dim,spacedim>                      &mapping,
                     const typename Triangulation<dim,spacedim>::cell_iterator    &cell,
@@ -2340,174 +1667,107 @@ protected:
                     FEValuesData<dim,spacedim>                       &data) const;
 
   /**
-   * Given the pattern of nonzero
-   * components for each shape
-   * function, compute for each
-   * entry how many components are
-   * non-zero for each shape
-   * function. This function is
-   * used in the constructor of
-   * this class.
+   * Given the pattern of nonzero components for each shape function, compute
+   * for each entry how many components are non-zero for each shape
+   * function. This function is used in the constructor of this class.
    */
   static
   std::vector<unsigned int>
   compute_n_nonzero_components (const std::vector<ComponentMask> &nonzero_components);
 
   /**
-   * Determine the values a finite
-   * element should compute on
-   * initialization of data for
-   * FEValues.
-   *
-   * Given a set of flags
-   * indicating what quantities are
-   * requested from a FEValues
-   * object, update_once() and
-   * update_each() compute which
-   * values must really be
-   * computed. Then, the
-   * <tt>fill_*_values</tt> functions
-   * are called with the result of
-   * these.
-   *
-   * Furthermore, values must be
-   * computed either on the unit
-   * cell or on the physical
-   * cell. For instance, the
-   * function values of FE_Q do
-   * only depend on the quadrature
-   * points on the unit
-   * cell. Therefore, this flags
-   * will be returned by
-   * update_once(). The gradients
-   * require computation of the
-   * covariant transformation
-   * matrix. Therefore,
-   * @p update_covariant_transformation
-   * and @p update_gradients will
-   * be returned by
-   * update_each().
-   *
-   * For an example see the same
-   * function in the derived class
-   * FE_Q.
+   * Determine the values a finite element should compute on initialization of
+   * data for FEValues.
+   *
+   * Given a set of flags indicating what quantities are requested from a
+   * FEValues object, update_once() and update_each() compute which values
+   * must really be computed. Then, the <tt>fill_*_values</tt> functions are
+   * called with the result of these.
+   *
+   * Furthermore, values must be computed either on the unit cell or on the
+   * physical cell. For instance, the function values of FE_Q do only depend
+   * on the quadrature points on the unit cell. Therefore, this flags will be
+   * returned by update_once(). The gradients require computation of the
+   * covariant transformation matrix. Therefore, @p
+   * update_covariant_transformation and @p update_gradients will be returned
+   * by update_each().
+   *
+   * For an example see the same function in the derived class FE_Q.
    */
   virtual UpdateFlags update_once (const UpdateFlags flags) const = 0;
 
   /**
-   * Complementary function for
-   * update_once().
+   * Complementary function for update_once().
    *
-   * While update_once() returns
-   * the values to be computed on
-   * the unit cell for yielding the
-   * required data, this function
-   * determines the values that
-   * must be recomputed on each
-   * cell.
+   * While update_once() returns the values to be computed on the unit cell
+   * for yielding the required data, this function determines the values that
+   * must be recomputed on each cell.
    *
-   * Refer to update_once() for
-   * more details.
+   * Refer to update_once() for more details.
    */
   virtual UpdateFlags update_each (const UpdateFlags flags) const = 0;
 
   /**
-   * A sort of virtual copy
-   * constructor. Some places in
-   * the library, for example the
-   * constructors of FESystem as
-   * well as the hp::FECollection
-   * class, need to make copies of
-   * finite elements without
-   * knowing their exact type. They
-   * do so through this function.
+   * A sort of virtual copy constructor. Some places in the library, for
+   * example the constructors of FESystem as well as the hp::FECollection
+   * class, need to make copies of finite elements without knowing their exact
+   * type. They do so through this function.
    */
   virtual FiniteElement<dim,spacedim> *clone() const = 0;
 
 private:
   /**
-   * Store what
-   * system_to_component_index()
-   * will return.
+   * Store what system_to_component_index() will return.
    */
   std::vector< std::pair<unsigned int, unsigned int> > system_to_component_table;
 
   /**
-   * Map between linear dofs and
-   * component dofs on face. This
-   * is filled with default values
-   * in the constructor, but
-   * derived classes will have to
-   * overwrite the information if
-   * necessary.
+   * Map between linear dofs and component dofs on face. This is filled with
+   * default values in the constructor, but derived classes will have to
+   * overwrite the information if necessary.
    *
-   * By component, we mean the
-   * vector component, not the base
-   * element. The information thus
-   * makes only sense if a shape
-   * function is non-zero in only
+   * By component, we mean the vector component, not the base element. The
+   * information thus makes only sense if a shape function is non-zero in only
    * one component.
    */
   std::vector< std::pair<unsigned int, unsigned int> > face_system_to_component_table;
 
   /**
-   * For each shape function, store
-   * to which base element and
-   * which instance of this base
-   * element (in case its
-   * multiplicity is greater than
-   * one) it belongs, and its index
-   * within this base element. If
-   * the element is not composed of
-   * others, then base and instance
-   * are always zero, and the index
-   * is equal to the number of the
-   * shape function. If the element
-   * is composed of single
-   * instances of other elements
-   * (i.e. all with multiplicity
-   * one) all of which are scalar,
-   * then base values and dof
-   * indices within this element
-   * are equal to the
-   * #system_to_component_table. It
-   * differs only in case the
-   * element is composed of other
-   * elements and at least one of
-   * them is vector-valued itself.
-   *
-   * This array has valid values
-   * also in the case of
-   * vector-valued
-   * (i.e. non-primitive) shape
-   * functions, in contrast to the
+   * For each shape function, store to which base element and which instance
+   * of this base element (in case its multiplicity is greater than one) it
+   * belongs, and its index within this base element. If the element is not
+   * composed of others, then base and instance are always zero, and the index
+   * is equal to the number of the shape function. If the element is composed
+   * of single instances of other elements (i.e. all with multiplicity one)
+   * all of which are scalar, then base values and dof indices within this
+   * element are equal to the #system_to_component_table. It differs only in
+   * case the element is composed of other elements and at least one of them
+   * is vector-valued itself.
+   *
+   * This array has valid values also in the case of vector-valued
+   * (i.e. non-primitive) shape functions, in contrast to the
    * #system_to_component_table.
    */
   std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
   system_to_base_table;
 
   /**
-   * Likewise for the indices on
-   * faces.
+   * Likewise for the indices on faces.
    */
   std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
   face_system_to_base_table;
 
   /**
-   * For each base element, store
-   * the number of blocks generated
-   * by the base and the first block in a block
-   * vector it will generate.
+   * For each base element, store the number of blocks generated by the base
+   * and the first block in a block vector it will generate.
    */
   BlockIndices base_to_block_indices;
 
   /**
-   * The base element establishing
-   * a component.
+   * The base element establishing a component.
    *
-   * For each component number
-   * <tt>c</tt>, the entries have
-   * the following meaning:
+   * For each component number <tt>c</tt>, the entries have the following
+   * meaning:
    * <dl>
    * <dt><tt>table[c].first.first</tt></dt>
    * <dd>Number of the base element for <tt>c</tt>.</dd>
@@ -2517,147 +1777,82 @@ private:
    * <dd>Multiple of the base element for <tt>c</tt>.</dd>
    * </dl>
    *
-   * This variable is set to the
-   * correct size by the
-   * constructor of this class, but
-   * needs to be initialized by
-   * derived classes, unless its
-   * size is one and the only entry
-   * is a zero, which is the case
-   * for scalar elements. In that
-   * case, the initialization by
-   * the base class is sufficient.
+   * This variable is set to the correct size by the constructor of this
+   * class, but needs to be initialized by derived classes, unless its size is
+   * one and the only entry is a zero, which is the case for scalar
+   * elements. In that case, the initialization by the base class is
+   * sufficient.
    */
   std::vector<std::pair<std::pair<unsigned int, unsigned int>, unsigned int> >
   component_to_base_table;
 
   /**
-   * Projection matrices are
-   * concatenated or summed up.
-   *
-   * This flags decides on how the
-   * projection matrices of the
-   * children of the same father
-   * are put together to one
-   * operator. The possible modes
-   * are concatenation and
-   * summation.
-   *
-   * If the projection is defined
-   * by an interpolation operator,
-   * the child matrices are
-   * concatenated, i.e. values
-   * belonging to the same node
-   * functional are identified and
-   * enter the interpolated value
-   * only once. In this case, the
-   * flag must be @p false.
-   *
-   * For projections with respect
-   * to scalar products, the child
-   * matrices must be summed up to
-   * build the complete matrix. The
-   * flag should be @p true.
-   *
-   * For examples of use of these
-   * flags, see the places in the
-   * library where it is queried.
-   *
-   * There is one flag per shape
-   * function, indicating whether
-   * it belongs to the class of
-   * shape functions that are
-   * additive in the restriction or
-   * not.
-   *
-   * Note that in previous versions
-   * of the library, there was one
-   * flag per vector component of
-   * the element. This is based on
-   * the fact that all the shape
-   * functions that belong to the
-   * same vector component must
-   * necessarily behave in the same
-   * way, to make things
-   * reasonable. However, the
-   * problem is that it is
-   * sometimes impossible to query
-   * this flag in the vector-valued
-   * case: this used to be done
-   * with the
-   * #system_to_component_index
-   * function that returns which
-   * vector component a shape
-   * function is associated
-   * with. The point is that since
-   * we now support shape functions
-   * that are associated with more
-   * than one vector component (for
-   * example the shape functions of
-   * Raviart-Thomas, or Nedelec
-   * elements), that function can
-   * no more be used, so it can be
-   * difficult to find out which
-   * for vector component we would
-   * like to query the
+   * Projection matrices are concatenated or summed up.
+   *
+   * This flags decides on how the projection matrices of the children of the
+   * same father are put together to one operator. The possible modes are
+   * concatenation and summation.
+   *
+   * If the projection is defined by an interpolation operator, the child
+   * matrices are concatenated, i.e. values belonging to the same node
+   * functional are identified and enter the interpolated value only once. In
+   * this case, the flag must be @p false.
+   *
+   * For projections with respect to scalar products, the child matrices must
+   * be summed up to build the complete matrix. The flag should be @p true.
+   *
+   * For examples of use of these flags, see the places in the library where
+   * it is queried.
+   *
+   * There is one flag per shape function, indicating whether it belongs to
+   * the class of shape functions that are additive in the restriction or not.
+   *
+   * Note that in previous versions of the library, there was one flag per
+   * vector component of the element. This is based on the fact that all the
+   * shape functions that belong to the same vector component must necessarily
+   * behave in the same way, to make things reasonable. However, the problem
+   * is that it is sometimes impossible to query this flag in the
+   * vector-valued case: this used to be done with the
+   * #system_to_component_index function that returns which vector component a
+   * shape function is associated with. The point is that since we now support
+   * shape functions that are associated with more than one vector component
+   * (for example the shape functions of Raviart-Thomas, or Nedelec elements),
+   * that function can no more be used, so it can be difficult to find out
+   * which for vector component we would like to query the
    * restriction-is-additive flags.
    */
   const std::vector<bool> restriction_is_additive_flags;
 
   /**
-   * For each shape function, give
-   * a vector of bools (with size
-   * equal to the number of vector
-   * components which this finite
-   * element has) indicating in
-   * which component each of these
-   * shape functions is non-zero.
+   * For each shape function, give a vector of bools (with size equal to the
+   * number of vector components which this finite element has) indicating in
+   * which component each of these shape functions is non-zero.
    *
-   * For primitive elements, there
-   * is only one non-zero
-   * component.
+   * For primitive elements, there is only one non-zero component.
    */
   const std::vector<ComponentMask> nonzero_components;
 
   /**
-   * This array holds how many
-   * values in the respective entry
-   * of the #nonzero_components
-   * element are non-zero. The
-   * array is thus a short-cut to
-   * allow faster access to this
-   * information than if we had to
-   * count the non-zero entries
-   * upon each request for this
-   * information. The field is
-   * initialized in the constructor
-   * of this class.
+   * This array holds how many values in the respective entry of the
+   * #nonzero_components element are non-zero. The array is thus a short-cut
+   * to allow faster access to this information than if we had to count the
+   * non-zero entries upon each request for this information. The field is
+   * initialized in the constructor of this class.
    */
   const std::vector<unsigned int> n_nonzero_components_table;
 
   /**
-   * Second derivatives of shapes
-   * functions are not computed
-   * analytically, but by finite
-   * differences of the
-   * gradients. This static
-   * variable denotes the step
-   * length to be used for
-   * that. It's value is set to
-   * 1e-6.
+   * Second derivatives of shapes functions are not computed analytically, but
+   * by finite differences of the gradients. This static variable denotes the
+   * step length to be used for that. It's value is set to 1e-6.
    */
   static const double fd_step_length;
 
   /**
-   * Prepare internal data
-   * structures and fill in values
-   * independent of the
-   * cell. Returns a pointer to an
-   * object of which the caller of
-   * this function then has to
-   * assume ownership (which
-   * includes destruction when it
-   * is no more needed).
+   * Prepare internal data structures and fill in values independent of the
+   * cell. Returns a pointer to an object of which the caller of this function
+   * then has to assume ownership (which includes destruction when it is no
+   * more needed).
    */
   virtual typename Mapping<dim,spacedim>::InternalDataBase *
   get_data (const UpdateFlags      flags,
@@ -2665,16 +1860,10 @@ private:
             const Quadrature<dim> &quadrature) const = 0;
 
   /**
-   * Prepare internal data
-   * structure for transformation
-   * of faces and fill in values
-   * independent of the
-   * cell. Returns a pointer to an
-   * object of which the caller of
-   * this function then has to
-   * assume ownership (which
-   * includes destruction when it
-   * is no more needed).
+   * Prepare internal data structure for transformation of faces and fill in
+   * values independent of the cell. Returns a pointer to an object of which
+   * the caller of this function then has to assume ownership (which includes
+   * destruction when it is no more needed).
    */
   virtual typename Mapping<dim,spacedim>::InternalDataBase *
   get_face_data (const UpdateFlags        flags,
@@ -2682,16 +1871,10 @@ private:
                  const Quadrature<dim-1> &quadrature) const;
 
   /**
-   * Prepare internal data
-   * structure for transformation
-   * of children of faces and fill
-   * in values independent of the
-   * cell. Returns a pointer to an
-   * object of which the caller of
-   * this function then has to
-   * assume ownership (which
-   * includes destruction when it
-   * is no more needed).
+   * Prepare internal data structure for transformation of children of faces
+   * and fill in values independent of the cell. Returns a pointer to an
+   * object of which the caller of this function then has to assume ownership
+   * (which includes destruction when it is no more needed).
    */
   virtual typename Mapping<dim,spacedim>::InternalDataBase *
   get_subface_data (const UpdateFlags        flags,
@@ -2699,15 +1882,11 @@ private:
                     const Quadrature<dim-1> &quadrature) const;
 
   /**
-   * Fill the fields of
-   * FEValues. This function
-   * performs all the operations
-   * needed to compute the data of an
-   * FEValues object.
+   * Fill the fields of FEValues. This function performs all the operations
+   * needed to compute the data of an FEValues object.
    *
-   * The same function in
-   * @p mapping must have been
-   * called for the same cell first!
+   * The same function in @p mapping must have been called for the same cell
+   * first!
    */
   virtual void
   fill_fe_values (const Mapping<dim,spacedim>                               &mapping,
@@ -2719,15 +1898,11 @@ private:
                   CellSimilarity::Similarity                           &cell_similarity) const = 0;
 
   /**
-   * Fill the fields of
-   * FEFaceValues. This function
-   * performs all the operations
-   * needed to compute the data of an
-   * FEFaceValues object.
+   * Fill the fields of FEFaceValues. This function performs all the
+   * operations needed to compute the data of an FEFaceValues object.
    *
-   * The same function in
-   * @p mapping must have been
-   * called for the same cell first!
+   * The same function in @p mapping must have been called for the same cell
+   * first!
    */
   virtual void
   fill_fe_face_values (const Mapping<dim,spacedim>                   &mapping,
@@ -2739,15 +1914,11 @@ private:
                        FEValuesData<dim,spacedim>                    &data) const = 0;
 
   /**
-   * Fill the fields of
-   * FESubfaceValues. This function
-   * performs all the operations
-   * needed to compute the data of an
-   * FESubfaceValues object.
+   * Fill the fields of FESubfaceValues. This function performs all the
+   * operations needed to compute the data of an FESubfaceValues object.
    *
-   * The same function in
-   * @p mapping must have been
-   * called for the same cell first!
+   * The same function in @p mapping must have been called for the same cell
+   * first!
    */
   virtual void
   fill_fe_subface_values (const Mapping<dim,spacedim>                   &mapping,

diff --git a/deal.II/include/deal.II/fe/fe_q_base.h b/deal.II/include/deal.II/fe/fe_q_base.h
index 9ecbfd0eb946dad10ecc138d647b56ef92b45109..e3fe259c17bf4b36950ab408dae64956bf4c555f 100644 (file)
--- a/deal.II/include/deal.II/fe/fe_q_base.h
+++ b/deal.II/include/deal.II/fe/fe_q_base.h
@@ -14,6 +14,7 @@
 
 #include <deal.II/base/config.h>
 #include <deal.II/fe/fe_poly.h>
+#include <deal.II/base/thread_management.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -95,6 +96,62 @@ public:
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
+  /**
+   * Projection from a fine grid space onto a coarse grid space. Overrides the
+   * respective method in FiniteElement, implementing lazy evaluation
+   * (initialize when requested).
+   *
+   * If this projection operator is associated with a matrix @p P, then the
+   * restriction of this matrix @p P_i to a single child cell is returned
+   * here.
+   *
+   * The matrix @p P is the concatenation or the sum of the cell matrices @p
+   * P_i, depending on the #restriction_is_additive_flags. This distinguishes
+   * interpolation (concatenation) and projection with respect to scalar
+   * products (summation).
+   *
+   * Row and column indices are related to coarse grid and fine grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * If projection matrices are not implemented in the derived finite element
+   * class, this function aborts with ExcProjectionVoid. You can check whether
+   * this is the case by calling the restriction_is_implemented() or the
+   * isotropic_restriction_is_implemented() function.
+   */
+  virtual const FullMatrix<double> &
+  get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
+  /**
+   * Embedding matrix between grids. Overrides the respective method in
+   * FiniteElement, implementing lazy evaluation (initialize when queried).
+   *
+   * The identity operator from a coarse grid space into a fine grid space is
+   * associated with a matrix @p P. The restriction of this matrix @p P_i to a
+   * single child cell is returned here.
+   *
+   * The matrix @p P is the concatenation, not the sum of the cell matrices @p
+   * P_i. That is, if the same non-zero entry <tt>j,k</tt> exists in in two
+   * different child matrices @p P_i, the value should be the same in both
+   * matrices and it is copied into the matrix @p P only once.
+   *
+   * Row and column indices are related to fine grid and coarse grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * These matrices are used by routines assembling the prolongation matrix
+   * for multi-level methods.  Upon assembling the transfer matrix between
+   * cells using this matrix array, zero elements in the prolongation matrix
+   * are discarded and will not fill up the transfer matrix.
+   *
+   * If projection matrices are not implemented in the derived finite element
+   * class, this function aborts with ExcEmbeddingVoid. You can check whether
+   * this is the case by calling the prolongation_is_implemented() or the
+   * isotropic_prolongation_is_implemented() function.
+   */
+  virtual const FullMatrix<double> &
+  get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
   /**
    * @name Functions to support hp
    * @{
@@ -182,16 +239,6 @@ protected:
    */
   void initialize_constraints (const std::vector<Point<1> > &points);
 
-  /**
-   * Initialize the embedding matrices. Called from initialize().
-   */
-  void initialize_embedding ();
-
-  /**
-   * Initialize the restriction matrices. Called from initialize().
-   */
-  void initialize_restriction ();
-
   /**
   * Initialize the @p unit_support_points field of the FiniteElement
   * class. Called from initialize().
@@ -222,6 +269,12 @@ protected:
    * Declare implementation friend.
    */
   friend struct FE_Q_Base<POLY,dim,spacedim>::Implementation;
+
+private:
+  /*
+   * Mutex for protecting initialization of restriction and embedding matrix.
+   */
+  mutable Threads::Mutex mutex;
 };
 
 

diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h
index e5f408b544d6ef8617b35509df17d9053170faf1..650dbc5a43599b726f8a45385115af4939052c98 100644 (file)
--- a/deal.II/include/deal.II/fe/fe_system.h
+++ b/deal.II/include/deal.II/fe/fe_system.h
@@ -17,6 +17,7 @@
 
 
 #include <deal.II/base/config.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/fe/fe.h>
 #include <vector>
 #include <utility>
@@ -152,32 +153,23 @@ class FESystem : public FiniteElement<dim,spacedim>
 public:
 
   /**
-   * Constructor. Take a finite element type
-   * and the number of elements you want to
-   * group together using this class.
+   * Constructor. Take a finite element type and the number of elements you
+   * want to group together using this class.
    *
-   * In fact, the object @p fe is not used,
-   * apart from getting the number of dofs
-   * per vertex, line, etc for that finite
-   * element class. The objects creates its
-   * own copy of the finite element object
-   * at construction time (but after
-   * the initialization of the base class
-   * @p FiniteElement, which is why we need
-   * a valid finite element object passed
-   * to the constructor).
+   * In fact, the object @p fe is not used, apart from getting the number of
+   * dofs per vertex, line, etc for that finite element class. The objects
+   * creates its own copy of the finite element object at construction time
+   * (but after the initialization of the base class @p FiniteElement, which
+   * is why we need a valid finite element object passed to the constructor).
    *
-   * Obviously, the template finite element
-   * class needs to be of the same dimension
-   * as is this object.
+   * Obviously, the template finite element class needs to be of the same
+   * dimension as is this object.
    */
   FESystem (const FiniteElement<dim,spacedim> &fe,
             const unsigned int n_elements);
 
   /**
-   * Constructor for mixed
-   * discretizations with two
-   * base elements.
+   * Constructor for mixed discretizations with two base elements.
    *
    * See the other constructor.
    */
@@ -185,9 +177,7 @@ public:
             const FiniteElement<dim,spacedim> &fe2, const unsigned int n2);
 
   /**
-   * Constructor for mixed
-   * discretizations with three
-   * base elements.
+   * Constructor for mixed discretizations with three base elements.
    *
    * See the other constructor.
    */
@@ -196,9 +186,7 @@ public:
             const FiniteElement<dim,spacedim> &fe3, const unsigned int n3);
 
   /**
-   * Constructor for mixed
-   * discretizations with four
-   * base elements.
+   * Constructor for mixed discretizations with four base elements.
    *
    * See the other constructor.
    */
@@ -208,9 +196,7 @@ public:
             const FiniteElement<dim,spacedim> &fe4, const unsigned int n4);
 
   /**
-   * Constructor for mixed
-   * discretizations with five
-   * base elements.
+   * Constructor for mixed discretizations with five base elements.
    *
    * See the other constructor.
    */
@@ -221,15 +207,9 @@ public:
             const FiniteElement<dim,spacedim> &fe5, const unsigned int n5);
 
   /**
-   * Same as above but for any
-   * number of base
-   * elements. Pointers to the base
-   * elements and their
-   * multiplicities are passed as
-   * vectors to this
-   * constructor. The length of
-   * these vectors is assumed to be
-   * equal.
+   * Same as above but for any number of base elements. Pointers to the base
+   * elements and their multiplicities are passed as vectors to this
+   * constructor. The length of these vectors is assumed to be equal.
    */
 
   FESystem (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
@@ -241,172 +221,95 @@ public:
   virtual ~FESystem ();
 
   /**
-   * Return a string that uniquely
-   * identifies a finite
-   * element. This element returns
-   * a string that is composed of
-   * the strings
-   * @p name1...@p nameN returned
-   * by the basis elements. From
-   * these, we create a sequence
-   * <tt>FESystem<dim>[name1^m1-name2^m2-...-nameN^mN]</tt>,
-   * where @p mi are the
-   * multiplicities of the basis
-   * elements. If a multiplicity is
-   * equal to one, then the
-   * superscript is omitted.
+   * Return a string that uniquely identifies a finite element. This element
+   * returns a string that is composed of the strings @p name1...@p nameN
+   * returned by the basis elements. From these, we create a sequence
+   * <tt>FESystem<dim>[name1^m1-name2^m2-...-nameN^mN]</tt>, where @p mi are
+   * the multiplicities of the basis elements. If a multiplicity is equal to
+   * one, then the superscript is omitted.
    */
   virtual std::string get_name () const;
 
   /**
-   * Return the value of the
-   * @p ith shape function at the
-   * point @p p.  @p p is a point
-   * on the reference element. Since
-   * this finite element is always
-   * vector-valued, we return the
-   * value of the only non-zero
-   * component of the vector value
-   * of this shape function. If the
-   * shape function has more than
-   * one non-zero component (which
-   * we refer to with the term
-   * non-primitive), then throw an
-   * exception of type
-   * @p ExcShapeFunctionNotPrimitive.
+   * Return the value of the @p ith shape function at the point @p p.  @p p is
+   * a point on the reference element. Since this finite element is always
+   * vector-valued, we return the value of the only non-zero component of the
+   * vector value of this shape function. If the shape function has more than
+   * one non-zero component (which we refer to with the term non-primitive),
+   * then throw an exception of type @p ExcShapeFunctionNotPrimitive.
    *
-   * An
-   * @p ExcUnitShapeValuesDoNotExist
-   * is thrown if the shape values
-   * of the @p FiniteElement
-   * (corresponding to the @p ith
-   * shape function) depend on the
-   * shape of the cell in real
-   * space.
+   * An @p ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+   * @p FiniteElement (corresponding to the @p ith shape function) depend on
+   * the shape of the cell in real space.
    */
   virtual double shape_value (const unsigned int i,
                               const Point<dim> &p) const;
 
   /**
-   * Return the value of the
-   * @p componentth vector
-   * component of the @p ith shape
-   * function at the point
-   * @p p. See the
-   * FiniteElement base
-   * class for more information
-   * about the semantics of this
-   * function.
+   * Return the value of the @p componentth vector component of the @p ith
+   * shape function at the point @p p. See the FiniteElement base class for
+   * more information about the semantics of this function.
    *
-   * Since this element is vector
-   * valued in general, it relays
-   * the computation of these
-   * values to the base elements.
+   * Since this element is vector valued in general, it relays the computation
+   * of these values to the base elements.
    */
   virtual double shape_value_component (const unsigned int i,
                                         const Point<dim> &p,
                                         const unsigned int component) const;
 
   /**
-   * Return the gradient of the
-   * @p ith shape function at the
-   * point @p p. @p p is a point
-   * on the reference element, and
-   * likewise the gradient is the
-   * gradient on the unit cell with
-   * respect to unit cell
-   * coordinates. Since
-   * this finite element is always
-   * vector-valued, we return the
-   * value of the only non-zero
-   * component of the vector value
-   * of this shape function. If the
-   * shape function has more than
-   * one non-zero component (which
-   * we refer to with the term
-   * non-primitive), then throw an
-   * exception of type
-   * @p ExcShapeFunctionNotPrimitive.
+   * Return the gradient of the @p ith shape function at the point @p p. @p p
+   * is a point on the reference element, and likewise the gradient is the
+   * gradient on the unit cell with respect to unit cell coordinates. Since
+   * this finite element is always vector-valued, we return the value of the
+   * only non-zero component of the vector value of this shape function. If
+   * the shape function has more than one non-zero component (which we refer
+   * to with the term non-primitive), then throw an exception of type @p
+   * ExcShapeFunctionNotPrimitive.
    *
-   * An
-   * @p ExcUnitShapeValuesDoNotExist
-   * is thrown if the shape values
-   * of the @p FiniteElement
-   * (corresponding to the @p ith
-   * shape function) depend on the
-   * shape of the cell in real
-   * space.
+   * An @p ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+   * @p FiniteElement (corresponding to the @p ith shape function) depend on
+   * the shape of the cell in real space.
    */
   virtual Tensor<1,dim> shape_grad (const unsigned int  i,
                                     const Point<dim>   &p) const;
 
   /**
-   * Return the gradient of the
-   * @p componentth vector
-   * component of the @p ith shape
-   * function at the point
-   * @p p. See the
-   * FiniteElement base
-   * class for more information
-   * about the semantics of this
-   * function.
+   * Return the gradient of the @p componentth vector component of the @p ith
+   * shape function at the point @p p. See the FiniteElement base class for
+   * more information about the semantics of this function.
    *
-   * Since this element is vector
-   * valued in general, it relays
-   * the computation of these
-   * values to the base elements.
+   * Since this element is vector valued in general, it relays the computation
+   * of these values to the base elements.
    */
   virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
                                               const Point<dim> &p,
                                               const unsigned int component) const;
 
   /**
-   * Return the tensor of second
-   * derivatives of the @p ith
-   * shape function at point @p p
-   * on the unit cell. The
-   * derivatives are derivatives on
-   * the unit cell with respect to
-   * unit cell coordinates. Since
-   * this finite element is always
-   * vector-valued, we return the
-   * value of the only non-zero
-   * component of the vector value
-   * of this shape function. If the
-   * shape function has more than
-   * one non-zero component (which
-   * we refer to with the term
-   * non-primitive), then throw an
-   * exception of type
-   * @p ExcShapeFunctionNotPrimitive.
+   * Return the tensor of second derivatives of the @p ith shape function at
+   * point @p p on the unit cell. The derivatives are derivatives on the unit
+   * cell with respect to unit cell coordinates. Since this finite element is
+   * always vector-valued, we return the value of the only non-zero component
+   * of the vector value of this shape function. If the shape function has
+   * more than one non-zero component (which we refer to with the term
+   * non-primitive), then throw an exception of type @p
+   * ExcShapeFunctionNotPrimitive.
    *
-   * An
-   * @p ExcUnitShapeValuesDoNotExist
-   * is thrown if the shape values
-   * of the @p FiniteElement
-   * (corresponding to the @p ith
-   * shape function) depend on the
-   * shape of the cell in real
-   * space.
+   * An @p ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+   * @p FiniteElement (corresponding to the @p ith shape function) depend on
+   * the shape of the cell in real space.
    */
   virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
                                          const Point<dim> &p) const;
 
   /**
-   * Return the second derivatives
-   * of the @p componentth vector
-   * component of the @p ith shape
-   * function at the point
-   * @p p. See the
-   * FiniteElement base
-   * class for more information
-   * about the semantics of this
-   * function.
+   * Return the second derivatives of the @p componentth vector component of
+   * the @p ith shape function at the point @p p. See the FiniteElement base
+   * class for more information about the semantics of this function.
    *
-   * Since this element is vector
-   * valued in general, it relays
-   * the computation of these
-   * values to the base elements.
+   * Since this element is vector valued in general, it relays the computation
+   * of these values to the base elements.
    */
   virtual
   Tensor<2,dim>
@@ -415,76 +318,105 @@ public:
                              const unsigned int component) const;
 
   /**
-   * Return the matrix
-   * interpolating from the given
-   * finite element to the present
-   * one. The size of the matrix is
-   * then @p dofs_per_cell times
+   * Return the matrix interpolating from the given finite element to the
+   * present one. The size of the matrix is then @p dofs_per_cell times
    * <tt>source.dofs_per_cell</tt>.
    *
-   * These matrices are available
-   * if source and destination
-   * element are both @p FESystem
-   * elements, have the same number
-   * of base elements with same
-   * element multiplicity, and if
-   * these base elements also
-   * implement their
-   * @p get_interpolation_matrix
-   * functions. Otherwise, an
-   * exception of type
-   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented
-   * is thrown.
+   * These matrices are available if source and destination element are both
+   * @p FESystem elements, have the same number of base elements with same
+   * element multiplicity, and if these base elements also implement their @p
+   * get_interpolation_matrix functions. Otherwise, an exception of type
+   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
    */
   virtual void
   get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                             FullMatrix<double>           &matrix) const;
 
   /**
-   * Access to a composing
-   * element. The index needs to be
-   * smaller than the number of
-   * base elements. Note that the
-   * number of base elements may in
-   * turn be smaller than the
-   * number of components of the
-   * system element, if the
-   * multiplicities are greater
-   * than one.
+   * Access to a composing element. The index needs to be smaller than the
+   * number of base elements. Note that the number of base elements may in
+   * turn be smaller than the number of components of the system element, if
+   * the multiplicities are greater than one.
    */
   virtual const FiniteElement<dim,spacedim> &
   base_element (const unsigned int index) const;
 
   /**
-   * Check for non-zero values on a
-   * face.
+   * Check for non-zero values on a face.
    *
-   * This function returns
-   * @p true, if the shape
-   * function @p shape_index has
-   * non-zero values on the face
-   * @p face_index.
+   * This function returns @p true, if the shape function @p shape_index has
+   * non-zero values on the face @p face_index.
    *
-   * Implementation of the
-   * interface in
-   * FiniteElement
+   * Implementation of the interface in FiniteElement
    */
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
   /**
-   * Implementation of the
-   * respective function in the
-   * base class.
+   * Projection from a fine grid space onto a coarse grid space. Overrides the
+   * respective method in FiniteElement, implementing lazy evaluation
+   * (initialize when requested).
+   *
+   * If this projection operator is associated with a matrix @p P, then the
+   * restriction of this matrix @p P_i to a single child cell is returned
+   * here.
+   *
+   * The matrix @p P is the concatenation or the sum of the cell matrices @p
+   * P_i, depending on the #restriction_is_additive_flags. This distinguishes
+   * interpolation (concatenation) and projection with respect to scalar
+   * products (summation).
+   *
+   * Row and column indices are related to coarse grid and fine grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * If projection matrices are not implemented in the derived finite element
+   * class, this function aborts with ExcProjectionVoid. You can check whether
+   * this is the case by calling the restriction_is_implemented() or the
+   * isotropic_restriction_is_implemented() function.
+   */
+  virtual const FullMatrix<double> &
+  get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
+  /**
+   * Embedding matrix between grids. Overrides the respective method in
+   * FiniteElement, implementing lazy evaluation (initialize when queried).
+   *
+   * The identity operator from a coarse grid space into a fine grid space is
+   * associated with a matrix @p P. The restriction of this matrix @p P_i to a
+   * single child cell is returned here.
+   *
+   * The matrix @p P is the concatenation, not the sum of the cell matrices @p
+   * P_i. That is, if the same non-zero entry <tt>j,k</tt> exists in in two
+   * different child matrices @p P_i, the value should be the same in both
+   * matrices and it is copied into the matrix @p P only once.
+   *
+   * Row and column indices are related to fine grid and coarse grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * These matrices are used by routines assembling the prolongation matrix
+   * for multi-level methods.  Upon assembling the transfer matrix between
+   * cells using this matrix array, zero elements in the prolongation matrix
+   * are discarded and will not fill up the transfer matrix.
+   *
+   * If projection matrices are not implemented in the derived finite element
+   * class, this function aborts with ExcEmbeddingVoid. You can check whether
+   * this is the case by calling the prolongation_is_implemented() or the
+   * isotropic_prolongation_is_implemented() function.
+   */
+  virtual const FullMatrix<double> &
+  get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
+  /**
+   * Implementation of the respective function in the base class.
    */
   virtual
   Point<dim>
   unit_support_point (const unsigned int index) const;
 
   /**
-   * Implementation of the
-   * respective function in the
-   * base class.
+   * Implementation of the respective function in the base class.
    */
   virtual
   Point<dim-1>
@@ -496,39 +428,25 @@ public:
    */
 
   /**
-   * Return whether this element
-   * implements its hanging node
-   * constraints in the new way,
-   * which has to be used to make
-   * elements "hp compatible".
+   * Return whether this element implements its hanging node constraints in
+   * the new way, which has to be used to make elements "hp compatible".
    *
-   * This function returns @p true iff all
-   * its base elements return @p true for
-   * this function.
+   * This function returns @p true iff all its base elements return @p true
+   * for this function.
    */
   virtual bool hp_constraints_are_implemented () const;
 
   /**
-   * Return the matrix
-   * interpolating from a face of
-   * of one element to the face of
-   * the neighboring element.
-   * The size of the matrix is
-   * then <tt>source.dofs_per_face</tt> times
-   * <tt>this->dofs_per_face</tt>.
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
    *
-   * Base elements of this element will
-   * have to implement this function. They
-   * may only provide interpolation
-   * matrices for certain source finite
-   * elements, for example those from the
-   * same family. If they don't implement
-   * interpolation from a given element,
-   * then they must throw an exception of
-   * type
-   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented,
-   * which will get propagated out from
-   * this element.
+   * Base elements of this element will have to implement this function. They
+   * may only provide interpolation matrices for certain source finite
+   * elements, for example those from the same family. If they don't implement
+   * interpolation from a given element, then they must throw an exception of
+   * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented, which
+   * will get propagated out from this element.
    */
   virtual void
   get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
@@ -536,26 +454,16 @@ public:
 
 
   /**
-   * Return the matrix
-   * interpolating from a face of
-   * of one element to the subface of
-   * the neighboring element.
-   * The size of the matrix is
-   * then <tt>source.dofs_per_face</tt> times
-   * <tt>this->dofs_per_face</tt>.
+   * Return the matrix interpolating from a face of of one element to the
+   * subface of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>.
    *
-   * Base elements of this element will
-   * have to implement this function. They
-   * may only provide interpolation
-   * matrices for certain source finite
-   * elements, for example those from the
-   * same family. If they don't implement
-   * interpolation from a given element,
-   * then they must throw an exception of
-   * type
-   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented,
-   * which will get propagated out from
-   * this element.
+   * Base elements of this element will have to implement this function. They
+   * may only provide interpolation matrices for certain source finite
+   * elements, for example those from the same family. If they don't implement
+   * interpolation from a given element, then they must throw an exception of
+   * type FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented, which
+   * will get propagated out from this element.
    */
   virtual void
   get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
@@ -600,15 +508,11 @@ public:
   hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
 
   /**
-   * Return whether this element dominates
-   * the one given as argument when they
-   * meet at a common face,
-   * whether it is the other way around,
-   * whether neither dominates, or if
-   * either could dominate.
+   * Return whether this element dominates the one given as argument when they
+   * meet at a common face, whether it is the other way around, whether
+   * neither dominates, or if either could dominate.
    *
-   * For a definition of domination, see
-   * FiniteElementBase::Domination and in
+   * For a definition of domination, see FiniteElementBase::Domination and in
    * particular the @ref hp_paper "hp paper".
    */
   virtual
@@ -617,37 +521,30 @@ public:
   //@}
 
   /**
-   * Determine an estimate for the
-   * memory consumption (in bytes)
-   * of this object.
+   * Determine an estimate for the memory consumption (in bytes) of this
+   * object.
    *
-   * This function is made virtual,
-   * since finite element objects
-   * are usually accessed through
-   * pointers to their base class,
-   * rather than the class itself.
+   * This function is made virtual, since finite element objects are usually
+   * accessed through pointers to their base class, rather than the class
+   * itself.
    */
   virtual std::size_t memory_consumption () const;
 
 protected:
   /**
-   * Compute flags for initial
-   * update only.
+   * Compute flags for initial update only.
    */
   virtual UpdateFlags update_once (const UpdateFlags flags) const;
 
   /**
-   * Compute flags for update on
-   * each cell.
+   * Compute flags for update on each cell.
    */
   virtual UpdateFlags update_each (const UpdateFlags flags) const;
 
   /**
-   * @p clone function instead of
-   * a copy constructor.
+   * @p clone function instead of a copy constructor.
    *
-   * This function is needed by the
-   * constructors of @p FESystem.
+   * This function is needed by the constructors of @p FESystem.
    */
   virtual FiniteElement<dim,spacedim> *clone() const;
 
@@ -667,15 +564,10 @@ protected:
                     const Quadrature<dim-1> &quadrature) const ;
 
   /**
-   * Implementation of the same
-   * function in
-   * FiniteElement.
+   * Implementation of the same function in FiniteElement.
    *
-   * Passes on control to
-   * @p compute_fill that does the
-   * work for all three
-   * <tt>fill_fe*_values</tt>
-   * functions.
+   * Passes on control to @p compute_fill that does the work for all three
+   * <tt>fill_fe*_values</tt> functions.
    */
   virtual void
   fill_fe_values (const Mapping<dim,spacedim>                      &mapping,
@@ -687,13 +579,9 @@ protected:
                   CellSimilarity::Similarity                  &cell_similarity) const;
 
   /**
-   * Implementation of the same
-   * function in
-   * FiniteElement.
+   * Implementation of the same function in FiniteElement.
    *
-   * Passes on control to
-   * @p compute_fill that does the
-   * work for all three
+   * Passes on control to @p compute_fill that does the work for all three
    * <tt>fill_fe*_values</tt> functions.
    */
   virtual void
@@ -706,13 +594,9 @@ protected:
                        FEValuesData<dim,spacedim>                    &data) const ;
 
   /**
-   * Implementation of the same
-   * function in
-   * FiniteElement.
+   * Implementation of the same function in FiniteElement.
    *
-   * Passes on control to
-   * @p compute_fill that does the
-   * work for all three
+   * Passes on control to @p compute_fill that does the work for all three
    * <tt>fill_fe*_values</tt> functions.
    */
   virtual void
@@ -727,26 +611,14 @@ protected:
 
 
   /**
-   * Do the work for the three
-   * <tt>fill_fe*_values</tt> functions.
+   * Do the work for the three <tt>fill_fe*_values</tt> functions.
    *
-   * Calls (among other things)
-   * <tt>fill_fe_([sub]face)_values</tt>
-   * of the base elements. Calls
-   * @p fill_fe_values if
-   * <tt>face_no==invalid_face_no</tt>
-   * and
-   * <tt>sub_no==invalid_face_no</tt>;
-   * calls @p fill_fe_face_values
-   * if
-   * <tt>face_no==invalid_face_no</tt>
-   * and
-   * <tt>sub_no!=invalid_face_no</tt>;
-   * and calls
-   * @p fill_fe_subface_values if
-   * <tt>face_no!=invalid_face_no</tt>
-   * and
-   * <tt>sub_no!=invalid_face_no</tt>.
+   * Calls (among other things) <tt>fill_fe_([sub]face)_values</tt> of the
+   * base elements. Calls @p fill_fe_values if
+   * <tt>face_no==invalid_face_no</tt> and <tt>sub_no==invalid_face_no</tt>;
+   * calls @p fill_fe_face_values if <tt>face_no==invalid_face_no</tt> and
+   * <tt>sub_no!=invalid_face_no</tt>; and calls @p fill_fe_subface_values if
+   * <tt>face_no!=invalid_face_no</tt> and <tt>sub_no!=invalid_face_no</tt>.
    */
   template <int dim_1>
   void compute_fill (const Mapping<dim,spacedim>                      &mapping,
@@ -762,9 +634,7 @@ protected:
 private:
 
   /**
-   * Value to indicate that a given
-   * face or subface number is
-   * invalid.
+   * Value to indicate that a given face or subface number is invalid.
    */
   static const unsigned int invalid_face_number = numbers::invalid_unsigned_int;
 
@@ -788,42 +658,28 @@ private:
 
 
   /**
-   * Initialize the
-   * @p unit_support_points field
-   * of the FiniteElement
-   * class. Called from the
-   * constructor.
+   * Initialize the @p unit_support_points field of the FiniteElement
+   * class. Called from the constructor.
    */
   void initialize_unit_support_points ();
 
   /**
-   * Initialize the
-   * @p unit_face_support_points field
-   * of the FiniteElement
-   * class. Called from the
-   * constructor.
+   * Initialize the @p unit_face_support_points field of the FiniteElement
+   * class. Called from the constructor.
    */
   void initialize_unit_face_support_points ();
 
   /**
-   * Initialize the
-   * @p adjust_quad_dof_index_for_face_orientation_table field
-   * of the FiniteElement
-   * class. Called from the
-   * constructor.
+   * Initialize the @p adjust_quad_dof_index_for_face_orientation_table field
+   * of the FiniteElement class. Called from the constructor.
    */
   void initialize_quad_dof_index_permutation ();
 
   /**
-   * Helper function used in the constructor:
-   * take a @p FiniteElementData object
-   * and return an object of the same type
-   * with the number of degrees of
-   * freedom per vertex, line, etc.
-   * multiplied by @p n. Don't touch the
-   * number of functions for the
-   * transformation from unit to real
-   * cell.
+   * Helper function used in the constructor: take a @p FiniteElementData
+   * object and return an object of the same type with the number of degrees
+   * of freedom per vertex, line, etc.  multiplied by @p n. Don't touch the
+   * number of functions for the transformation from unit to real cell.
    */
   static FiniteElementData<dim>
   multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
@@ -838,8 +694,7 @@ private:
                         const unsigned int            N5=0);
 
   /**
-   * Same as above but for
-   * any number of sub-elements.
+   * Same as above but for any number of sub-elements.
    */
   static FiniteElementData<dim>
   multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
@@ -848,15 +703,10 @@ private:
 
 
   /**
-   * Helper function used in the
-   * constructor: takes a
-   * @p FiniteElement object and
-   * returns an boolean vector
-   * including the
-   * @p restriction_is_additive_flags
-   * of the mixed element
-   * consisting of @p N elements
-   * of the sub-element @p fe.
+   * Helper function used in the constructor: takes a @p FiniteElement object
+   * and returns an boolean vector including the @p
+   * restriction_is_additive_flags of the mixed element consisting of @p N
+   * elements of the sub-element @p fe.
    */
   static std::vector<bool>
   compute_restriction_is_additive_flags (
@@ -872,11 +722,8 @@ private:
     const unsigned int        N5=0);
 
   /**
-   * Compute the named flags for a
-   * list of finite elements with
-   * multiplicities given in the
-   * second argument. This function
-   * is called from all the above
+   * Compute the named flags for a list of finite elements with multiplicities
+   * given in the second argument. This function is called from all the above
    * functions.
    */
   static std::vector<bool>
@@ -886,9 +733,7 @@ private:
 
 
   /**
-   * Compute the non-zero vector
-   * components of a composed
-   * finite element.
+   * Compute the non-zero vector components of a composed finite element.
    */
   static std::vector<ComponentMask>
   compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
@@ -903,26 +748,18 @@ private:
                               const unsigned int        N5=0);
 
   /**
-   * Compute the nonzero components
-   * of a list of finite elements
-   * with multiplicities given in
-   * the second argument. This
-   * function is called from all
-   * the above functions.
+   * Compute the nonzero components of a list of finite elements with
+   * multiplicities given in the second argument. This function is called from
+   * all the above functions.
   */
   static std::vector<ComponentMask>
   compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
                               const std::vector<unsigned int>              &multiplicities);
 
   /**
-   * This function is simply
-   * singled out of the
-   * constructors since there are
-   * several of them. It sets up
-   * the index table for the system
-   * as well as @p restriction and
-   * @p prolongation
-   * matrices.
+   * This function is simply singled out of the constructors since there are
+   * several of them. It sets up the index table for the system as well as @p
+   * restriction and @p prolongation matrices.
    */
   void initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
                    const std::vector<unsigned int> &multiplicities);
@@ -943,176 +780,116 @@ private:
   void build_interface_constraints ();
 
   /**
-   * A function that computes the
-   * hp_vertex_dof_identities(),
-   * hp_line_dof_identities(), or
-   * hp_quad_dof_identities(), depending on
-   * the value of the template parameter.
+   * A function that computes the hp_vertex_dof_identities(),
+   * hp_line_dof_identities(), or hp_quad_dof_identities(), depending on the
+   * value of the template parameter.
    */
   template <int structdim>
   std::vector<std::pair<unsigned int, unsigned int> >
   hp_object_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
 
   /**
-   * Usually: Fields of
-   * cell-independent data.
+   * Usually: Fields of cell-independent data.
    *
-   * However, here, this class does
-   * not itself store the data but
-   * only pointers to
-   * @p InternalData objects for
-   * each of the base elements.
+   * However, here, this class does not itself store the data but only
+   * pointers to @p InternalData objects for each of the base elements.
    */
   class InternalData : public FiniteElement<dim,spacedim>::InternalDataBase
   {
   public:
     /**
-     * Constructor. Is called by
-     * the @p get_data
-     * function. Sets the size of
-     * the @p base_fe_datas
-     * vector to
-     * @p n_base_elements.
+     * Constructor. Is called by the @p get_data function. Sets the size of
+     * the @p base_fe_datas vector to @p n_base_elements.
      */
     InternalData (const unsigned int n_base_elements);
 
     /**
-     * Destructor. Deletes all
-     * @p InternalDatas whose
-     * pointers are stored by the
-     * @p base_fe_datas
-     * vector.
+     * Destructor. Deletes all @p InternalDatas whose pointers are stored by
+     * the @p base_fe_datas vector.
      */
     ~InternalData();
 
     /**
-     * Flag indicating whether
-     * second derivatives shall
-     * be computed.
+     * Flag indicating whether second derivatives shall be computed.
      */
     bool compute_hessians;
 
     /**
-     * Gives write-access to the
-     * pointer to a
-     * @p InternalData of the
-     * @p base_noth base
-     * element.
+     * Gives write-access to the pointer to a @p InternalData of the @p
+     * base_noth base element.
      */
     void set_fe_data(const unsigned int                        base_no,
                      typename FiniteElement<dim,spacedim>::InternalDataBase *);
 
     /**
-     * Gives read-access to the
-     * pointer to a
-     * @p InternalData of the
-     * @p base_noth base element.
+     * Gives read-access to the pointer to a @p InternalData of the @p
+     * base_noth base element.
      */
     typename FiniteElement<dim,spacedim>::InternalDataBase &
     get_fe_data (const unsigned int base_no) const;
 
 
     /**
-     * Gives write-access to the
-     * pointer to a
-     * @p FEValuesData for the
-     * @p base_noth base
-     * element.
+     * Gives write-access to the pointer to a @p FEValuesData for the @p
+     * base_noth base element.
      */
     void set_fe_values_data (const unsigned int base_no,
                              FEValuesData<dim,spacedim> *);
 
     /**
-     * Gives read-access to the
-     * pointer to a
-     * @p FEValuesData for the
-     * @p base_noth base element.
+     * Gives read-access to the pointer to a @p FEValuesData for the @p
+     * base_noth base element.
      */
     FEValuesData<dim,spacedim> &get_fe_values_data (const unsigned int base_no) const;
 
     /**
-     * Deletes the
-     * @p FEValuesData the
-     * <tt>fe_datas[base_no]</tt>
-     * pointer is pointing
-     * to. Sets
-     * <tt>fe_datas[base_no]</tt> to
-     * zero.
+     * Deletes the @p FEValuesData the <tt>fe_datas[base_no]</tt> pointer is
+     * pointing to. Sets <tt>fe_datas[base_no]</tt> to zero.
      *
-     * This function is used to
-     * delete @p FEValuesData
-     * that are needed only on
-     * the first cell but not any
-     * more afterwards.  This is
-     * the case for
-     * e.g. Lagrangian elements
-     * (see e.g. @p FE_Q
-     * classes).
+     * This function is used to delete @p FEValuesData that are needed only on
+     * the first cell but not any more afterwards.  This is the case for
+     * e.g. Lagrangian elements (see e.g. @p FE_Q classes).
      */
     void delete_fe_values_data (const unsigned int base_no);
 
     /**
-     * Set the @p first_cell
-     * flag to @p false. Used by
-     * the @p FEValues class to
-     * indicate that we have
-     * already done the work on
-     * the first cell.
+     * Set the @p first_cell flag to @p false. Used by the @p FEValues class
+     * to indicate that we have already done the work on the first cell.
      *
-     * In addition to calling the
-     * respective function of the
-     * base class, this function
-     * also calls the functions
-     * of the sub-data objects.
+     * In addition to calling the respective function of the base class, this
+     * function also calls the functions of the sub-data objects.
      */
     virtual void clear_first_cell ();
 
   private:
 
     /**
-     * Pointers to
-     * @p InternalData objects
-     * for each of the base
-     * elements. They are
-     * accessed to by the
-     * @p set_ and
-     * @p get_fe_data functions.
+     * Pointers to @p InternalData objects for each of the base elements. They
+     * are accessed to by the @p set_ and @p get_fe_data functions.
      *
-     * The size of this vector is
-     * set to @p n_base_elements
-     * by the InternalData
-     * constructor.  It is filled
-     * by the @p get_data
-     * function. Note that since
-     * the data for each instance
-     * of a base class is
-     * necessarily the same, we
-     * only need as many of these
-     * objects as there are base
-     * elements, irrespective of
-     * their multiplicity.
+     * The size of this vector is set to @p n_base_elements by the
+     * InternalData constructor.  It is filled by the @p get_data
+     * function. Note that since the data for each instance of a base class is
+     * necessarily the same, we only need as many of these objects as there
+     * are base elements, irrespective of their multiplicity.
      */
     typename std::vector<typename FiniteElement<dim,spacedim>::InternalDataBase *> base_fe_datas;
 
     /**
-     * Pointers to the
-     * @p FEValuesData objects
-     * that are given to the
-     * @p fill_fe_values
-     * function of the base
-     * elements. They are
-     * accessed to by the
-     * @p set_ and
-     * @p get_fe_values_data
-     * functions.
+     * Pointers to the @p FEValuesData objects that are given to the @p
+     * fill_fe_values function of the base elements. They are accessed to by
+     * the @p set_ and @p get_fe_values_data functions.
      *
-     * The size of this vector is
-     * set to @p n_base_elements
-     * by the InternalData
-     * constructor.
+     * The size of this vector is set to @p n_base_elements by the
+     * InternalData constructor.
      */
     std::vector<FEValuesData<dim,spacedim> *> base_fe_values_datas;
   };
+
+  /*
+   * Mutex for protecting initialization of restriction and embedding matrix.
+   */
+  mutable Threads::Mutex mutex;
 };
 
 

diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc
index c5504e0cb8cb137f089ac140394971c956e2020d..1a9e9f9fa419a6c7d79e007a567664eb287d10fb 100644 (file)
--- a/deal.II/source/fe/fe.cc
+++ b/deal.II/source/fe/fe.cc
@@ -180,11 +180,9 @@ FiniteElement<dim,spacedim>::FiniteElement (
               ExcInternalError());
     };
 
-  // initialize some tables in the
-  // default way, i.e. if there is
-  // only one (vector-)component; if
-  // the element is not primitive,
-  // leave these tables empty.
+  // initialize some tables in the default way, i.e. if there is only one
+  // (vector-)component; if the element is not primitive, leave these tables
+  // empty.
   if (this->is_primitive())
     {
       system_to_component_table.resize(this->dofs_per_cell);
@@ -200,15 +198,11 @@ FiniteElement<dim,spacedim>::FiniteElement (
           face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
         }
     }
-  // Fill with default value; may be
-  // changed by constructor of
-  // derived class.
+  // Fill with default value; may be changed by constructor of derived class.
   base_to_block_indices.reinit(1,1);
 
-  // initialize the restriction and
-  // prolongation matrices. the default
-  // contructur of FullMatrix<dim> initializes
-  // them with size zero
+  // initialize the restriction and prolongation matrices. the default
+  // contructur of FullMatrix<dim> initializes them with size zero
   prolongation.resize(RefinementCase<dim>::isotropic_refinement);
   restriction.resize(RefinementCase<dim>::isotropic_refinement);
   for (unsigned int ref=RefinementCase<dim>::cut_x;
@@ -337,11 +331,9 @@ FiniteElement<dim,spacedim>::get_restriction_matrix (const unsigned int child,
           ExcMessage("Restriction matrices are only available for refined cells!"));
   Assert (child<GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
           ExcIndexRange(child,0,GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case))));
-  // we use refinement_case-1 here. the -1 takes care
-  // of the origin of the vector, as for
-  // RefinementCase<dim>::no_refinement (=0) there is no
-  // data available and so the vector indices
-  // are shifted
+  // we use refinement_case-1 here. the -1 takes care of the origin of the
+  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
+  // available and so the vector indices are shifted
   Assert (restriction[refinement_case-1][child].n() == this->dofs_per_cell, ExcProjectionVoid());
   return restriction[refinement_case-1][child];
 }
@@ -614,6 +606,8 @@ FiniteElement<dim,spacedim>::prolongation_is_implemented () const
     for (unsigned int c=0;
          c<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++c)
       {
+        // make sure also the lazily initialized matrices are created
+        get_prolongation_matrix(c, RefinementCase<dim>(ref_case));
         Assert ((prolongation[ref_case-1][c].m() == this->dofs_per_cell) ||
                 (prolongation[ref_case-1][c].m() == 0),
                 ExcInternalError());
@@ -638,6 +632,8 @@ FiniteElement<dim,spacedim>::restriction_is_implemented () const
     for (unsigned int c=0;
          c<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++c)
       {
+        // make sure also the lazily initialized matrices are created
+        get_restriction_matrix(c, RefinementCase<dim>(ref_case));
         Assert ((restriction[ref_case-1][c].m() == this->dofs_per_cell) ||
                 (restriction[ref_case-1][c].m() == 0),
                 ExcInternalError());
@@ -662,6 +658,8 @@ FiniteElement<dim,spacedim>::isotropic_prolongation_is_implemented () const
   for (unsigned int c=0;
        c<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++c)
     {
+      // make sure also the lazily initialized matrices are created
+      get_prolongation_matrix(c, RefinementCase<dim>(ref_case));
       Assert ((prolongation[ref_case-1][c].m() == this->dofs_per_cell) ||
               (prolongation[ref_case-1][c].m() == 0),
               ExcInternalError());
@@ -686,6 +684,8 @@ FiniteElement<dim,spacedim>::isotropic_restriction_is_implemented () const
   for (unsigned int c=0;
        c<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++c)
     {
+      // make sure also the lazily initialized matrices are created
+      get_restriction_matrix(c, RefinementCase<dim>(ref_case));
       Assert ((restriction[ref_case-1][c].m() == this->dofs_per_cell) ||
               (restriction[ref_case-1][c].m() == 0),
               ExcInternalError());

diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc
index b99c8b5d8a65fb1b681064b852d5a5cb605fb70f..9705535b4014e37bcedc71a75f038b100b7cc4d4 100644 (file)
--- a/deal.II/source/fe/fe_q_base.cc
+++ b/deal.II/source/fe/fe_q_base.cc
@@ -58,13 +58,13 @@ namespace FE_Q_Helper
 
 
 
-    // in initialize_embedding() and initialize_restriction(), want to undo
+    // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
     // tensorization on inner loops for performance reasons. this clears a
     // dim-array
     template <int dim>
     inline
     void
-    zero_indices (unsigned int indices[dim])
+    zero_indices (unsigned int (&indices)[dim])
     {
       for (unsigned int d=0; d<dim; ++d)
         indices[d] = 0;
@@ -72,14 +72,14 @@ namespace FE_Q_Helper
 
 
 
-    // in initialize_embedding() and initialize_restriction(), want to undo
+    // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
     // tensorization on inner loops for performance reasons. this increments
     // tensor product indices
     template <int dim>
     inline
     void
-    increment_indices (unsigned int       indices[dim],
-                       const unsigned int dofs1d)
+    increment_indices (unsigned int       (&indices)[dim],
+                       const unsigned int   dofs1d)
     {
       ++indices[0];
       for (int d=0; d<dim-1; ++d)
@@ -440,11 +440,9 @@ FE_Q_Base<POLY,dim,spacedim>::initialize (const std::vector<Point<1> > &points)
   // reinit constraints
   initialize_constraints (points);
 
-  // Reinit the vectors of restriction and prolongation matrices to the right
-  // sizes and compute the matrices
-  this->reinit_restriction_and_prolongation_matrices();
-  initialize_embedding();
-  initialize_restriction();
+  // do not initialize embedding and restriction here. these matrices are
+  // initialized on demand in get_restriction_matrix and
+  // get_prolongation_matrix
 
   this->initialize_quad_dof_index_permutation();
 }
@@ -1009,267 +1007,292 @@ FE_Q_Base<POLY,dim,spacedim>::initialize_constraints (const std::vector<Point<1>
 
 
 template <class POLY, int dim, int spacedim>
-void
-FE_Q_Base<POLY,dim,spacedim>::initialize_embedding ()
+const FullMatrix<double>&
+FE_Q_Base<POLY,dim,spacedim>
+::get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case) const
 {
-  // distinguish q/q_dg0 case: only treat Q dofs first
-  const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+          ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+          ExcMessage("Prolongation matrices are only available for refined cells!"));
+  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+  // needs initialization?
+  Threads::Mutex::ScopedLock(this->mutex);
+  if (this->prolongation[refinement_case-1][child].n() == 0)
+    {
+      // distinguish q/q_dg0 case: only treat Q dofs first
+      const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
 
-  // compute the interpolation matrices in much the same way as we do for the
-  // constraints. it's actually simpler here, since we don't have this weird
-  // renumbering stuff going on. The trick is again that we the interpolation
-  // matrix is formed by a permutation of the indices of the cell matrix. The
-  // value eps is used a threshold to decide when certain evaluations of the
-  // Lagrange polynomials are zero or one.
-  const double eps = 1e-15*this->degree*dim;
+      // compute the interpolation matrices in much the same way as we do for
+      // the constraints. it's actually simpler here, since we don't have this
+      // weird renumbering stuff going on. The trick is again that we the
+      // interpolation matrix is formed by a permutation of the indices of the
+      // cell matrix. The value eps is used a threshold to decide when certain
+      // evaluations of the Lagrange polynomials are zero or one.
+      const double eps = 1e-15*this->degree*dim;
 
 #ifdef DEBUG
-  // in DEBUG mode, check that the evaluation of support points in the current
-  // numbering gives the identity operation
-  for (unsigned int i=0; i<q_dofs_per_cell; ++i)
-    {
-      Assert (std::fabs (1.-this->poly_space.compute_value
-                         (i, this->unit_support_points[i])) < eps,
-              ExcInternalError());
-      for (unsigned int j=0; j<q_dofs_per_cell; ++j)
-        if (j!=i)
-          Assert (std::fabs (this->poly_space.compute_value
-                             (i, this->unit_support_points[j])) < eps,
+      // in DEBUG mode, check that the evaluation of support points in the
+      // current numbering gives the identity operation
+      for (unsigned int i=0; i<q_dofs_per_cell; ++i)
+        {
+          Assert (std::fabs (1.-this->poly_space.compute_value
+                             (i, this->unit_support_points[i])) < eps,
                   ExcInternalError());
-    }
+          for (unsigned int j=0; j<q_dofs_per_cell; ++j)
+            if (j!=i)
+              Assert (std::fabs (this->poly_space.compute_value
+                                 (i, this->unit_support_points[j])) < eps,
+                      ExcInternalError());
+        }
 #endif
 
-  // to efficiently evaluate the polynomial at the subcell, make use of the
-  // tensor product structure of this element and only evaluate 1D information
-  // from the polynomial. This makes the cost of this function almost
-  // negligible also for high order elements
-  const unsigned int dofs1d = this->degree+1;
-  std::vector<Table<2,double> >
-  subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d));
-  const std::vector<unsigned int> &index_map_inverse =
-    this->poly_space.get_numbering_inverse();
-
-  // helper value: step size how to walk through diagonal and how many points
-  // we have left apart from the first dimension
-  unsigned int step_size_diag = 0;
-  {
-    unsigned int factor = 1;
-    for (unsigned int d=0; d<dim; ++d)
-      {
-        step_size_diag += factor;
-        factor *= dofs1d;
-      }
-  }
+      // to efficiently evaluate the polynomial at the subcell, make use of
+      // the tensor product structure of this element and only evaluate 1D
+      // information from the polynomial. This makes the cost of this function
+      // almost negligible also for high order elements
+      const unsigned int dofs1d = this->degree+1;
+      std::vector<Table<2,double> >
+        subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d));
+      const std::vector<unsigned int> &index_map_inverse =
+        this->poly_space.get_numbering_inverse();
 
-  // next evaluate the functions for the different refinement cases.
-  for (unsigned int ref=0; ref<RefinementCase<dim>::isotropic_refinement; ++ref)
-    for (unsigned int child=0; child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref+1)); ++child)
+      // helper value: step size how to walk through diagonal and how many
+      // points we have left apart from the first dimension
+      unsigned int step_size_diag = 0;
       {
-        // go through the points in diagonal to capture variation in all
-        // directions simultaneously
-        for (unsigned int j=0; j<dofs1d; ++j)
+        unsigned int factor = 1;
+        for (unsigned int d=0; d<dim; ++d)
           {
-            const unsigned int diag_comp = index_map_inverse[j*step_size_diag];
-            const Point<dim> p_subcell = this->unit_support_points[diag_comp];
-            const Point<dim> p_cell =
-              GeometryInfo<dim>::child_to_cell_coordinates (p_subcell, child,
-                                                            RefinementCase<dim>(ref+1));
-            for (unsigned int i=0; i<dofs1d; ++i)
-              for (unsigned int d=0; d<dim; ++d)
-                {
-                  // evaluate along line where only x is different from zero
-                  Point<dim> point;
-                  point[0] = p_cell[d];
-                  const double cell_value =
-                    this->poly_space.compute_value(index_map_inverse[i], point);
-
-                  // cut off values that are too small. note that we have here
-                  // Lagrange interpolation functions, so they should be zero
-                  // at almost all points, and one at the others, at least on
-                  // the subcells. so set them to their exact values
-                  //
-                  // the actual cut-off value is somewhat fuzzy, but it works
-                  // for 2e-13*degree*dim (see above), which is kind of
-                  // reasonable given that we compute the values of the
-                  // polynomials via an degree-step recursion and then
-                  // multiply the 1d-values. this gives us a linear growth in
-                  // degree*dim, times a small constant.
-                  //
-                  // the embedding matrix is given by applying the inverse of
-                  // the subcell matrix on the cell_interpolation
-                  // matrix. since the subcell matrix is actually only a
-                  // permutation vector, all we need to do is to switch the
-                  // rows we write the data into. moreover, cut off very small
-                  // values here
-                  if (std::fabs(cell_value) < eps)
-                    subcell_evaluations[d](j,i) = 0;
-                  else
-                    subcell_evaluations[d](j,i) = cell_value;
-                }
+            step_size_diag += factor;
+            factor *= dofs1d;
           }
+      }
 
-        // now expand from 1D info. block innermost dimension (x_0) in order
-        // to avoid difficult checks at innermost loop
-        unsigned int j_indices[dim];
-        FE_Q_Helper::zero_indices<dim> (j_indices);
-        for (unsigned int j=0; j<q_dofs_per_cell; j+=dofs1d)
-          {
-            unsigned int i_indices[dim];
-            FE_Q_Helper::zero_indices<dim> (i_indices);
-            for (unsigned int i=0; i<q_dofs_per_cell; i+=dofs1d)
-              {
-                double val_extra_dim = 1.;
-                for (unsigned int d=1; d<dim; ++d)
-                  val_extra_dim *= subcell_evaluations[d](j_indices[d-1],
-                                                          i_indices[d-1]);
-
-                // innermost sum where we actually compute. the same as
-                // this->prolongation[ref][child](j,i) =
-                // this->poly_space.compute_value (i, p_cell);
-                for (unsigned int jj=0; jj<dofs1d; ++jj)
-                  {
-                    const unsigned int j_ind = index_map_inverse[j+jj];
-                    for (unsigned int ii=0; ii<dofs1d; ++ii)
-                      this->prolongation[ref][child](j_ind,index_map_inverse[i+ii])
-                        = val_extra_dim * subcell_evaluations[0](jj,ii);
-                  }
+      FullMatrix<double>& prolongate = const_cast<FullMatrix<double>&>
+        (this->prolongation[refinement_case-1][child]);
+      prolongate.reinit(this->dofs_per_cell, this->dofs_per_cell);
 
-                // update indices that denote the tensor product position. a
-                // bit fuzzy and therefore not done for innermost x_0
-                // direction
-                FE_Q_Helper::increment_indices<dim> (i_indices, dofs1d);
+      // go through the points in diagonal to capture variation in all
+      // directions simultaneously
+      for (unsigned int j=0; j<dofs1d; ++j)
+        {
+          const unsigned int diag_comp = index_map_inverse[j*step_size_diag];
+          const Point<dim> p_subcell = this->unit_support_points[diag_comp];
+          const Point<dim> p_cell =
+            GeometryInfo<dim>::child_to_cell_coordinates (p_subcell, child,
+                                                          refinement_case);
+          for (unsigned int i=0; i<dofs1d; ++i)
+            for (unsigned int d=0; d<dim; ++d)
+              {
+                // evaluate along line where only x is different from zero
+                Point<dim> point;
+                point[0] = p_cell[d];
+                const double cell_value =
+                  this->poly_space.compute_value(index_map_inverse[i], point);
+
+                // cut off values that are too small. note that we have here
+                // Lagrange interpolation functions, so they should be zero at
+                // almost all points, and one at the others, at least on the
+                // subcells. so set them to their exact values
+                //
+                // the actual cut-off value is somewhat fuzzy, but it works
+                // for 2e-13*degree*dim (see above), which is kind of
+                // reasonable given that we compute the values of the
+                // polynomials via an degree-step recursion and then multiply
+                // the 1d-values. this gives us a linear growth in degree*dim,
+                // times a small constant.
+                //
+                // the embedding matrix is given by applying the inverse of
+                // the subcell matrix on the cell_interpolation matrix. since
+                // the subcell matrix is actually only a permutation vector,
+                // all we need to do is to switch the rows we write the data
+                // into. moreover, cut off very small values here
+                if (std::fabs(cell_value) < eps)
+                  subcell_evaluations[d](j,i) = 0;
+                else
+                  subcell_evaluations[d](j,i) = cell_value;
               }
-            Assert (i_indices[dim-1] == 1, ExcInternalError());
-            FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
-          }
+        }
+
+      // now expand from 1D info. block innermost dimension (x_0) in order to
+      // avoid difficult checks at innermost loop
+      unsigned int j_indices[dim];
+      FE_Q_Helper::zero_indices<dim> (j_indices);
+      for (unsigned int j=0; j<q_dofs_per_cell; j+=dofs1d)
+        {
+          unsigned int i_indices[dim];
+          FE_Q_Helper::zero_indices<dim> (i_indices);
+          for (unsigned int i=0; i<q_dofs_per_cell; i+=dofs1d)
+            {
+              double val_extra_dim = 1.;
+              for (unsigned int d=1; d<dim; ++d)
+                val_extra_dim *= subcell_evaluations[d](j_indices[d-1],
+                                                        i_indices[d-1]);
+
+              // innermost sum where we actually compute. the same as
+              // prolongate(j,i) = this->poly_space.compute_value (i, p_cell)
+              for (unsigned int jj=0; jj<dofs1d; ++jj)
+                {
+                  const unsigned int j_ind = index_map_inverse[j+jj];
+                  for (unsigned int ii=0; ii<dofs1d; ++ii)
+                    prolongate(j_ind,index_map_inverse[i+ii])
+                      = val_extra_dim * subcell_evaluations[0](jj,ii);
+                }
+
+              // update indices that denote the tensor product position. a bit
+              // fuzzy and therefore not done for innermost x_0 direction
+              FE_Q_Helper::increment_indices<dim> (i_indices, dofs1d);
+            }
+          Assert (i_indices[dim-1] == 1, ExcInternalError());
+          FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+        }
 
-        // the discontinuous node is simply mapped on the discontinuous node
-        // on the child element
-        if (q_dofs_per_cell < this->dofs_per_cell)
-          this->prolongation[ref][child](q_dofs_per_cell,q_dofs_per_cell) = 1.;
+      // the discontinuous node is simply mapped on the discontinuous node on
+      // the child element
+      if (q_dofs_per_cell < this->dofs_per_cell)
+        prolongate(q_dofs_per_cell,q_dofs_per_cell) = 1.;
 
-        // and make sure that the row sum is 1. this must be so since for this
-        // element, the shape functions add up to one
+      // and make sure that the row sum is 1. this must be so since for this
+      // element, the shape functions add up to one
 #ifdef DEBUG
-        for (unsigned int row=0; row<this->dofs_per_cell; ++row)
-          {
-            double sum = 0;
-            for (unsigned int col=0; col<this->dofs_per_cell; ++col)
-              sum += this->prolongation[ref][child](row,col);
-            Assert (std::fabs(sum-1.) < eps, ExcInternalError());
-          }
+      for (unsigned int row=0; row<this->dofs_per_cell; ++row)
+        {
+          double sum = 0;
+          for (unsigned int col=0; col<this->dofs_per_cell; ++col)
+            sum += prolongate(row,col);
+          Assert (std::fabs(sum-1.) < eps, ExcInternalError());
+        }
 #endif
-      }
+    }
+
+  // finally return the matrix
+  return this->prolongation[refinement_case-1][child];
 }
 
 
 
 template <class POLY, int dim, int spacedim>
-void
-FE_Q_Base<POLY,dim,spacedim>::initialize_restriction ()
+const FullMatrix<double>&
+FE_Q_Base<POLY,dim,spacedim>
+::get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case) const
 {
-  // distinguish q/q_dg0 case
-  const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+          ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+          ExcMessage("Restriction matrices are only available for refined cells!"));
+  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+  // needs initialization?
+  Threads::Mutex::ScopedLock(this->mutex);
+  if (this->restriction[refinement_case-1][child].n() == 0)
+    {
+      FullMatrix<double> &restriction =
+        const_cast<FullMatrix<double> &>(this->restriction[refinement_case-1][child]);
+
+      // distinguish q/q_dg0 case
+      const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+
+      // for these Lagrange interpolation polynomials, construction of the
+      // restriction matrices is relatively simple. the reason is that the
+      // interpolation points on the mother cell are (except for the case with
+      // arbitrary nonequidistant nodes) always also interpolation points for
+      // some shape function on one or the other child, because we have chosen
+      // equidistant Lagrange interpolation points for the polynomials
+      //
+      // so the only thing we have to find out is: for each shape function on
+      // the mother cell, which is the child cell (possibly more than one) on
+      // which it is located, and which is the corresponding shape function
+      // there. rather than doing it for the shape functions on the mother
+      // cell, we take the interpolation points there
+      //
+      // note that the interpolation point of a shape function can be on the
+      // boundary between subcells. in that case, restriction from children to
+      // mother may produce two or more entries for a dof on the mother
+      // cell. however, this doesn't hurt: since the element is continuous,
+      // the contribution from each child should yield the same result, and
+      // since the element is non-additive we just overwrite one value
+      // (compute on one child) by the same value (compute on a later child),
+      // so we don't have to care about this
+
+      const double eps = 1e-15*this->degree*dim;
+      const std::vector<unsigned int> &index_map_inverse =
+        this->poly_space.get_numbering_inverse();
 
-  // for these Lagrange interpolation polynomials, construction of the
-  // restriction matrices is relatively simple. the reason is that the
-  // interpolation points on the mother cell are (except for the case with
-  // arbitrary nonequidistant nodes) always also interpolation points for some
-  // shape function on one or the other child, because we have chosen
-  // equidistant Lagrange interpolation points for the polynomials
-  //
-  // so the only thing we have to find out is: for each shape function on the
-  // mother cell, which is the child cell (possibly more than one) on which it
-  // is located, and which is the corresponding shape function there. rather
-  // than doing it for the shape functions on the mother cell, we take the
-  // interpolation points there
-  //
-  // note that the interpolation point of a shape function can be on the
-  // boundary between subcells. in that case, restriction from children to
-  // mother may produce two or more entries for a dof on the mother
-  // cell. however, this doesn't hurt: since the element is continuous, the
-  // contribution from each child should yield the same result, and since the
-  // element is non-additive we just overwrite one value (compute on one
-  // child) by the same value (compute on a later child), so we don't have to
-  // care about this
-
-  const double eps = 1e-15*this->degree*dim;
-  const std::vector<unsigned int> &index_map_inverse =
-    this->poly_space.get_numbering_inverse();
+      const unsigned int dofs1d = this->degree+1;
+      std::vector<Tensor<1,dim> > evaluations1d (dofs1d);
 
-  //  recreate 1D polynomials for faster evaluation of polynomial
-  const unsigned int dofs1d = this->degree+1;
-  std::vector<Tensor<1,dim> > evaluations1d (dofs1d);
+      restriction.reinit(this->dofs_per_cell, this->dofs_per_cell);
 
-  for (unsigned int i=0; i<q_dofs_per_cell; ++i)
-    {
-      unsigned int mother_dof = index_map_inverse[i];
-      const Point<dim> p_cell = this->unit_support_points[mother_dof];
-
-      // then find the children on which the interpolation point is located
-      for (unsigned int ref=RefinementCase<dim>::cut_x; ref<=RefinementCase<dim>::isotropic_refinement; ++ref)
-        for (unsigned int child=0; child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref)); ++child)
-          {
-            // check whether this interpolation point is inside this child
-            // cell
-            const Point<dim> p_subcell
-              = GeometryInfo<dim>::cell_to_child_coordinates (p_cell, child, RefinementCase<dim>(ref));
-            if (GeometryInfo<dim>::is_inside_unit_cell (p_subcell))
-              {
-                // same logic as in initialize_embedding to evaluate the
-                // polynomial faster than from the tensor product: since we
-                // evaluate all polynomials, it is much faster to just compute
-                // the 1D values for all polynomials before and then get the
-                // dim-data.
-                for (unsigned int j=0; j<dofs1d; ++j)
-                  for (unsigned int d=0; d<dim; ++d)
-                    {
-                      Point<dim> point;
-                      point[0] = p_subcell[d];
-                      evaluations1d[j][d] =
-                        this->poly_space.compute_value(index_map_inverse[j], point);
-                    }
-                unsigned int j_indices[dim];
-                FE_Q_Helper::zero_indices<dim> (j_indices);
-                double sum_check = 0;
-                for (unsigned int j = 0; j<q_dofs_per_cell; j += dofs1d)
+      for (unsigned int i=0; i<q_dofs_per_cell; ++i)
+        {
+          unsigned int mother_dof = index_map_inverse[i];
+          const Point<dim> p_cell = this->unit_support_points[mother_dof];
+
+          // check whether this interpolation point is inside this child cell
+          const Point<dim> p_subcell
+            = GeometryInfo<dim>::cell_to_child_coordinates (p_cell, child,
+                                                            refinement_case);
+          if (GeometryInfo<dim>::is_inside_unit_cell (p_subcell))
+            {
+              // same logic as in initialize_embedding to evaluate the
+              // polynomial faster than from the tensor product: since we
+              // evaluate all polynomials, it is much faster to just compute
+              // the 1D values for all polynomials before and then get the
+              // dim-data.
+              for (unsigned int j=0; j<dofs1d; ++j)
+                for (unsigned int d=0; d<dim; ++d)
                   {
-                    double val_extra_dim = 1.;
-                    for (unsigned int d=1; d<dim; ++d)
-                      val_extra_dim *= evaluations1d[j_indices[d-1]][d];
-                    for (unsigned int jj=0; jj<dofs1d; ++jj)
-                      {
-
-                        // find the child shape function(s) corresponding to
-                        // this point. Usually this is just one function;
-                        // however, when we use FE_Q on arbitrary nodes a
-                        // parent support point might not be a child support
-                        // point, and then we will get more than one nonzero
-                        // value per row. Still, the values should sum up to 1
-                        const double val
-                          = val_extra_dim * evaluations1d[jj][0];
-                        const unsigned int child_dof =
-                          index_map_inverse[j+jj];
-                        if (std::fabs (val-1.) < eps)
-                          this->restriction[ref-1][child](mother_dof,child_dof)=1.;
-                        else if (std::fabs(val) > eps)
-                          this->restriction[ref-1][child](mother_dof,child_dof)=val;
-                        sum_check += val;
-                      }
-                    FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+                    Point<dim> point;
+                    point[0] = p_subcell[d];
+                    evaluations1d[j][d] =
+                      this->poly_space.compute_value(index_map_inverse[j], point);
                   }
-                Assert (std::fabs(sum_check-1) < eps,
-                        ExcInternalError());
-              }
+              unsigned int j_indices[dim];
+              FE_Q_Helper::zero_indices<dim> (j_indices);
+              double sum_check = 0;
+              for (unsigned int j = 0; j<q_dofs_per_cell; j += dofs1d)
+                {
+                  double val_extra_dim = 1.;
+                  for (unsigned int d=1; d<dim; ++d)
+                    val_extra_dim *= evaluations1d[j_indices[d-1]][d];
+                  for (unsigned int jj=0; jj<dofs1d; ++jj)
+                    {
 
-            // part for FE_Q_DG0
-            if (q_dofs_per_cell < this->dofs_per_cell)
-              this->restriction[ref-1][child]
-                (this->dofs_per_cell-1,this->dofs_per_cell-1) =
-                1./(double)GeometryInfo<dim>::n_children
-                (RefinementCase<dim>(ref));
-          }
+                      // find the child shape function(s) corresponding to
+                      // this point. Usually this is just one function;
+                      // however, when we use FE_Q on arbitrary nodes a parent
+                      // support point might not be a child support point, and
+                      // then we will get more than one nonzero value per
+                      // row. Still, the values should sum up to 1
+                      const double val
+                        = val_extra_dim * evaluations1d[jj][0];
+                      const unsigned int child_dof =
+                        index_map_inverse[j+jj];
+                      if (std::fabs (val-1.) < eps)
+                        restriction(mother_dof,child_dof)=1.;
+                      else if (std::fabs(val) > eps)
+                        restriction(mother_dof,child_dof)=val;
+                      sum_check += val;
+                    }
+                  FE_Q_Helper::increment_indices<dim> (j_indices, dofs1d);
+                }
+              Assert (std::fabs(sum_check-1) < eps,
+                      ExcInternalError());
+            }
+
+          // part for FE_Q_DG0
+          if (q_dofs_per_cell < this->dofs_per_cell)
+            restriction(this->dofs_per_cell-1,this->dofs_per_cell-1) =
+              1./GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case));
+        }
     }
+  return this->restriction[refinement_case-1][child];
 }
 
 

diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc
index 6bb7de848964e7983b09014e97d34f68f45d661c..c3ecf565909a313eb87c5153f0760b1da5e3c161 100644 (file)
--- a/deal.II/source/fe/fe_system.cc
+++ b/deal.II/source/fe/fe_system.cc
@@ -54,8 +54,7 @@ FESystem<dim,spacedim>::InternalData::InternalData(const unsigned int n_base_ele
 template <int dim, int spacedim>
 FESystem<dim,spacedim>::InternalData::~InternalData()
 {
-  // delete pointers and set them to
-  // zero to avoid inadvertent use
+  // delete pointers and set them to zero to avoid inadvertent use
   for (unsigned int i=0; i<base_fe_datas.size(); ++i)
     if (base_fe_datas[i])
       {
@@ -453,24 +452,17 @@ FESystem<dim,spacedim>::shape_grad_component (const unsigned int i,
   Assert (component < this->n_components(),
           ExcIndexRange (component, 0, this->n_components()));
 
-  // if this value is supposed to be
-  // zero, then return right away...
+  // if this value is supposed to be zero, then return right away...
   if (this->nonzero_components[i][component] == false)
     return Tensor<1,dim>();
 
-  // ...otherwise: first find out to
-  // which of the base elements this
-  // desired component belongs, and
-  // which component within this base
-  // element it is
+  // ...otherwise: first find out to which of the base elements this desired
+  // component belongs, and which component within this base element it is
   const unsigned int base              = this->component_to_base_index(component).first;
   const unsigned int component_in_base = this->component_to_base_index(component).second;
 
-  // then get value from base
-  // element. note that that will
-  // throw an error should the
-  // respective shape function not be
-  // primitive; thus, there is no
+  // then get value from base element. note that that will throw an error
+  // should the respective shape function not be primitive; thus, there is no
   // need to check this here
   return (base_element(base).
           shape_grad_component(this->system_to_base_table[i].second,
@@ -505,24 +497,17 @@ FESystem<dim,spacedim>::shape_grad_grad_component (const unsigned int i,
   Assert (component < this->n_components(),
           ExcIndexRange (component, 0, this->n_components()));
 
-  // if this value is supposed to be
-  // zero, then return right away...
+  // if this value is supposed to be zero, then return right away...
   if (this->nonzero_components[i][component] == false)
     return Tensor<2,dim>();
 
-  // ...otherwise: first find out to
-  // which of the base elements this
-  // desired component belongs, and
-  // which component within this base
-  // element it is
+  // ...otherwise: first find out to which of the base elements this desired
+  // component belongs, and which component within this base element it is
   const unsigned int base              = this->component_to_base_index(component).first;
   const unsigned int component_in_base = this->component_to_base_index(component).second;
 
-  // then get value from base
-  // element. note that that will
-  // throw an error should the
-  // respective shape function not be
-  // primitive; thus, there is no
+  // then get value from base element. note that that will throw an error
+  // should the respective shape function not be primitive; thus, there is no
   // need to check this here
   return (base_element(base).
           shape_grad_grad_component(this->system_to_base_table[i].second,
@@ -545,12 +530,10 @@ FESystem<dim,spacedim>::get_interpolation_matrix (
           ExcDimensionMismatch (interpolation_matrix.m(),
                                 x_source_fe.dofs_per_cell));
 
-  // there are certain conditions
-  // that the two elements have to
-  // satisfy so that this can work.
+  // there are certain conditions that the two elements have to satisfy so
+  // that this can work.
   //
-  // condition 1: the other element
-  // must also be a system element
+  // condition 1: the other element must also be a system element
 
   typedef FiniteElement<dim,spacedim> FEL;
   AssertThrow ((x_source_fe.get_name().find ("FESystem<") == 0)
@@ -559,38 +542,29 @@ FESystem<dim,spacedim>::get_interpolation_matrix (
                typename FEL::
                ExcInterpolationNotImplemented());
 
-  // ok, source is a system element,
-  // so we may be able to do the work
+  // ok, source is a system element, so we may be able to do the work
   const FESystem<dim,spacedim> &source_fe
     = dynamic_cast<const FESystem<dim,spacedim>&>(x_source_fe);
 
-  // condition 2: same number of
-  // basis elements
+  // condition 2: same number of basis elements
   AssertThrow (this->n_base_elements() == source_fe.n_base_elements(),
                typename FEL::
                ExcInterpolationNotImplemented());
 
-  // condition 3: same number of
-  // basis elements
+  // condition 3: same number of basis elements
   for (unsigned int i=0; i< this->n_base_elements(); ++i)
     AssertThrow (this->element_multiplicity(i) ==
                  source_fe.element_multiplicity(i),
                  typename FEL::
                  ExcInterpolationNotImplemented());
 
-  // ok, so let's try whether it
-  // works:
-
-  // first let's see whether all the
-  // basis elements actually generate
-  // their interpolation matrices. if
-  // we get past the following loop,
-  // then apparently none of the
-  // called base elements threw an
-  // exception, so we're fine
-  // continuing and assembling the
-  // one big matrix from the small
-  // ones of the base elements
+  // ok, so let's try whether it works:
+
+  // first let's see whether all the basis elements actually generate their
+  // interpolation matrices. if we get past the following loop, then
+  // apparently none of the called base elements threw an exception, so we're
+  // fine continuing and assembling the one big matrix from the small ones of
+  // the base elements
   std::vector<FullMatrix<double> > base_matrices (this->n_base_elements());
   for (unsigned int i=0; i<this->n_base_elements(); ++i)
     {
@@ -600,11 +574,9 @@ FESystem<dim,spacedim>::get_interpolation_matrix (
                                                 base_matrices[i]);
     }
 
-  // first clear big matrix, to make
-  // sure that entries that would
-  // couple different bases (or
-  // multiplicity indices) are really
-  // zero. then assign entries
+  // first clear big matrix, to make sure that entries that would couple
+  // different bases (or multiplicity indices) are really zero. then assign
+  // entries
   interpolation_matrix = 0;
   for (unsigned int i=0; i<this->dofs_per_cell; ++i)
     for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
@@ -618,6 +590,169 @@ FESystem<dim,spacedim>::get_interpolation_matrix (
 
 
 
+template <int dim, int spacedim>
+const FullMatrix<double>&
+FESystem<dim,spacedim>
+::get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case) const
+{
+  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+          ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+          ExcMessage("Restriction matrices are only available for refined cells!"));
+  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+  // needs initialization?
+  Threads::Mutex::ScopedLock(this->mutex);
+  if (this->restriction[refinement_case-1][child].n() == 0)
+    {
+      // Check if some of the matrices of the base elements are void.
+      bool do_restriction = true;
+
+      // shortcut for accessing local restrictions further down
+      std::vector<const FullMatrix<double> *>
+        base_matrices(this->n_base_elements());
+
+      for (unsigned int i=0; i<this->n_base_elements(); ++i)
+        {
+          base_matrices[i] =
+            &base_element(i).get_restriction_matrix(child, refinement_case);
+          if (base_matrices[i]->n() != base_element(i).dofs_per_cell)
+            do_restriction = false;
+        }
+      Assert(do_restriction,
+             (typename FiniteElement<dim,spacedim>::ExcProjectionVoid()));
+
+      // if we did not encounter void matrices, initialize the matrix sizes
+      if (do_restriction)
+        {
+          FullMatrix<double> &restriction =
+            const_cast<FullMatrix<double>&>(this->restriction[refinement_case-1]
+                                            [child]);
+          restriction.reinit(this->dofs_per_cell, this->dofs_per_cell);
+
+          // distribute the matrices of the base finite elements to the
+          // matrices of this object. for this, loop over all degrees of
+          // freedom and take the respective entry of the underlying base
+          // element.
+          //
+          // note that we by definition of a base element, they are
+          // independent, i.e. do not couple. only DoFs that belong to the
+          // same instance of a base element may couple
+          for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+            for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+              {
+                // first find out to which base element indices i and j
+                // belong, and which instance thereof in case the base element
+                // has a multiplicity greater than one. if they should not
+                // happen to belong to the same instance of a base element,
+                // then they cannot couple, so go on with the next index
+                if (this->system_to_base_table[i].first !=
+                    this->system_to_base_table[j].first)
+                  continue;
+
+                // so get the common base element and the indices therein:
+                const unsigned int
+                  base = this->system_to_base_table[i].first.first;
+
+                const unsigned int
+                  base_index_i = this->system_to_base_table[i].second,
+                  base_index_j = this->system_to_base_table[j].second;
+
+                // if we are sure that DoFs i and j may couple, then copy
+                // entries of the matrices:
+                restriction(i,j) = (*base_matrices[base])(base_index_i,base_index_j);
+              }
+        }
+    }
+
+  return this->restriction[refinement_case-1][child];
+}
+
+
+
+template <int dim, int spacedim>
+const FullMatrix<double>&
+FESystem<dim,spacedim>
+::get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case) const
+{
+  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+          ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+          ExcMessage("Restriction matrices are only available for refined cells!"));
+  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+  // needs initialization?
+  Threads::Mutex::ScopedLock(this->mutex);
+  if (this->prolongation[refinement_case-1][child].n() == 0)
+    {
+      // Check if some of the matrices of the base elements are void.
+      bool do_prolongation = true;
+
+      // shortcut for accessing local prolongations further down
+      std::vector<const FullMatrix<double> *>
+        base_matrices(this->n_base_elements());
+
+      for (unsigned int i=0; i<this->n_base_elements(); ++i)
+        {
+          base_matrices[i] =
+            &base_element(i).get_prolongation_matrix(child, refinement_case);
+          if (base_matrices[i]->n() != base_element(i).dofs_per_cell)
+            do_prolongation = false;
+        }
+      Assert(do_prolongation,
+             (typename FiniteElement<dim,spacedim>::ExcEmbeddingVoid()));
+
+
+      if (do_prolongation)
+        {
+          FullMatrix<double> &prolongate =
+            const_cast<FullMatrix<double> &>(this->prolongation[refinement_case-1][child]);
+          prolongate.reinit(this->dofs_per_cell, this->dofs_per_cell);
+
+          // distribute the matrices of the base finite elements to the
+          // matrices of this object. for this, loop over all degrees of
+          // freedom and take the respective entry of the underlying base
+          // element.
+          //
+          // note that we by definition of a base element, they are
+          // independent, i.e. do not couple. only DoFs that belong to the
+          // same instance of a base element may couple
+          for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+            for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+              {
+                // first find out to which base element indices i and j
+                // belong, and which instance thereof in case the base element
+                // has a multiplicity greater than one. if they should not
+                // happen to belong to the same instance of a base element,
+                // then they cannot couple, so go on with the next index
+                if (this->system_to_base_table[i].first !=
+                    this->system_to_base_table[j].first)
+                  continue;
+
+                // so get the common base element and the indices therein:
+                const unsigned int
+                  base = this->system_to_base_table[i].first.first;
+
+                const unsigned int
+                  base_index_i = this->system_to_base_table[i].second,
+                  base_index_j = this->system_to_base_table[j].second;
+
+                // if we are sure that DoFs i and j may couple, then copy
+                // entries of the matrices:
+                prolongate(i,j) = (*base_matrices[base])(base_index_i,base_index_j);
+              }
+        }
+    }
+
+  return this->prolongation[refinement_case-1][child];
+}
+
+
+
 //---------------------------------------------------------------------------
 // Data field initialization
 //---------------------------------------------------------------------------
@@ -941,25 +1076,14 @@ FESystem<dim,spacedim>::compute_fill (
     {
       if (fe_data.is_first_cell())
         {
-          // Initialize the
-          // FEValuesDatas for the
-          // base elements.
-          // Originally this was the
-          // task of
-          // FEValues::FEValues() but
-          // the latter initializes
-          // the FEValuesDatas only
-          // of the FESystem, not of
-          // the FEValuesDatas needed
-          // by the base elements
-          // (and: how should it know
-          // of their existence,
-          // after all).
+          // Initialize the FEValuesDatas for the base elements.  Originally
+          // this was the task of FEValues::FEValues() but the latter
+          // initializes the FEValuesDatas only of the FESystem, not of the
+          // FEValuesDatas needed by the base elements (and: how should it
+          // know of their existence, after all).
           for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
             {
-              // Pointer needed to get
-              // the update flags of the
-              // base element
+              // Pointer needed to get the update flags of the base element
               typename Mapping<dim,spacedim>::InternalDataBase &
               base_fe_data = fe_data.get_fe_data(base_no);
 
@@ -967,28 +1091,20 @@ FESystem<dim,spacedim>::compute_fill (
               const UpdateFlags base_update_flags
                 = mapping_data.update_flags | base_fe_data.update_flags;
 
-              // Initialize the FEValuesDatas
-              // for the base elements.
+              // Initialize the FEValuesDatas for the base elements.
               FEValuesData<dim,spacedim> &base_data=fe_data.get_fe_values_data(base_no);
               const FiniteElement<dim,spacedim> &base_fe=base_element(base_no);
               base_data.initialize (n_q_points, base_fe, base_update_flags);
             }
         }
 
-      // fill_fe_face_values needs
-      // argument Quadrature<dim-1>
-      // for both cases
-      // dim_1==dim-1 and
-      // dim_1=dim. Hence the
-      // following workaround
+      // fill_fe_face_values needs argument Quadrature<dim-1> for both cases
+      // dim_1==dim-1 and dim_1=dim. Hence the following workaround
       const Quadrature<dim>   *cell_quadrature = 0;
       const Quadrature<dim-1> *face_quadrature = 0;
 
-      // static cast to the
-      // common base class of
-      // quadrature being either
-      // Quadrature<dim> or
-      // Quadrature<dim-1>:
+      // static cast to the common base class of quadrature being either
+      // Quadrature<dim> or Quadrature<dim-1>:
       const Subscriptor *quadrature_base_pointer = &quadrature;
 
       if (face_no==invalid_face_number)
@@ -1010,8 +1126,7 @@ FESystem<dim,spacedim>::compute_fill (
             = static_cast<const Quadrature<dim-1> *>(quadrature_base_pointer);
         }
 
-      // let base elements update the
-      // necessary data
+      // let base elements update the necessary data
       for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
         {
           const FiniteElement<dim,spacedim> &
@@ -1021,34 +1136,21 @@ FESystem<dim,spacedim>::compute_fill (
           FEValuesData<dim,spacedim> &
           base_data    = fe_data.get_fe_values_data(base_no);
 
-          //TODO: Think about a smarter alternative
-          // Copy quadrature points. These
-          // are required for computing the
-          // determinant in the FEPolyTensor
-          // class. The determinant is
-          // one ingredient of the Piola
-          // transformation, which is applied
-          // to correctly map the RT space from
-          // the reference element to the global
-          // coordinate system.
+          //TODO: Think about a smarter alternative Copy quadrature
+          // points. These are required for computing the determinant in the
+          // FEPolyTensor class. The determinant is one ingredient of the
+          // Piola transformation, which is applied to correctly map the RT
+          // space from the reference element to the global coordinate system.
           if (cell_similarity != CellSimilarity::translation)
             base_data.JxW_values = data.JxW_values;
 
 
-          // Make sure that in the
-          // case of fill_fe_values
-          // the data is only copied
-          // from base_data to data
-          // if base_data is
-          // changed. therefore use
-          // fe_fe_data.current_update_flags()
+          // Make sure that in the case of fill_fe_values the data is only
+          // copied from base_data to data if base_data is changed. therefore
+          // use fe_fe_data.current_update_flags()
           //
-          // for the case of
-          // fill_fe_(sub)face_values
-          // the data needs to be
-          // copied from base_data to
-          // data on each face,
-          // therefore use
+          // for the case of fill_fe_(sub)face_values the data needs to be
+          // copied from base_data to data on each face, therefore use
           // base_fe_data.update_flags.
           if (face_no==invalid_face_number)
             base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_data,
@@ -1060,45 +1162,25 @@ FESystem<dim,spacedim>::compute_fill (
             base_fe.fill_fe_subface_values(mapping, cell, face_no, sub_no,
                                            *face_quadrature, mapping_data, base_fe_data, base_data);
 
-          // now data has been
-          // generated, so copy
-          // it. we used to work by
-          // looping over all base
-          // elements (i.e. this
-          // outer loop), then over
-          // multiplicity, then over
-          // the shape functions from
-          // that base element, but
-          // that requires that we
-          // can infer the global
-          // number of a shape
-          // function from its number
-          // in the base element. for
-          // that we had the
-          // component_to_system_table.
+          // now data has been generated, so copy it. we used to work by
+          // looping over all base elements (i.e. this outer loop), then over
+          // multiplicity, then over the shape functions from that base
+          // element, but that requires that we can infer the global number of
+          // a shape function from its number in the base element. for that we
+          // had the component_to_system_table.
           //
-          // however, this does of
-          // course no longer work
-          // since we have
-          // non-primitive
-          // elements. so we go the
-          // other way round: loop
-          // over all shape functions
-          // of the composed element,
-          // and here only treat
-          // those shape functions
-          // that belong to a given
-          // base element
+          // however, this does of course no longer work since we have
+          // non-primitive elements. so we go the other way round: loop over
+          // all shape functions of the composed element, and here only treat
+          // those shape functions that belong to a given base element
 //TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style
           const UpdateFlags base_flags(dim_1==dim ?
                                        base_fe_data.current_update_flags() :
                                        base_fe_data.update_flags);
 
-          // if the current cell is just a
-          // translation of the previous
-          // one, the underlying data has
-          // not changed, and we don't even
-          // need to enter this section
+          // if the current cell is just a translation of the previous one,
+          // the underlying data has not changed, and we don't even need to
+          // enter this section
           if (cell_similarity != CellSimilarity::translation)
             for (unsigned int system_index=0; system_index<this->dofs_per_cell;
                  ++system_index)
@@ -1108,22 +1190,12 @@ FESystem<dim,spacedim>::compute_fill (
                   base_index = this->system_to_base_table[system_index].second;
                   Assert (base_index<base_fe.dofs_per_cell, ExcInternalError());
 
-                  // now copy. if the
-                  // shape function is
-                  // primitive, then
-                  // there is only one
-                  // value to be copied,
-                  // but for
-                  // non-primitive
-                  // elements, there
-                  // might be more values
-                  // to be copied
+                  // now copy. if the shape function is primitive, then there
+                  // is only one value to be copied, but for non-primitive
+                  // elements, there might be more values to be copied
                   //
-                  // so, find out from
-                  // which index to take
-                  // this one value, and
-                  // to which index to
-                  // put
+                  // so, find out from which index to take this one value, and
+                  // to which index to put
                   unsigned int out_index = 0;
                   for (unsigned int i=0; i<system_index; ++i)
                     out_index += this->n_nonzero_components(i);
@@ -1131,9 +1203,7 @@ FESystem<dim,spacedim>::compute_fill (
                   for (unsigned int i=0; i<base_index; ++i)
                     in_index += base_fe.n_nonzero_components(i);
 
-                  // then loop over the
-                  // number of components
-                  // to be copied
+                  // then loop over the number of components to be copied
                   Assert (this->n_nonzero_components(system_index) ==
                           base_fe.n_nonzero_components(base_index),
                           ExcInternalError());
@@ -1149,15 +1219,8 @@ FESystem<dim,spacedim>::compute_fill (
                           data.shape_gradients[out_index+s][q]=
                             base_data.shape_gradients[in_index+s][q];
 
-                      // _we_ handle
-                      // computation of
-                      // second
-                      // derivatives,
-                      // so the base
-                      // elements
-                      // should not
-                      // have computed
-                      // them!
+                      // _we_ handle computation of second derivatives, so the
+                      // base elements should not have computed them!
                       Assert (!(base_flags & update_hessians),
                               ExcInternalError());
                     };
@@ -1197,9 +1260,8 @@ template <int dim, int spacedim>
 void
 FESystem<dim,spacedim>::build_cell_tables()
 {
-  // If the system is not primitive,
-  // these have not been initialized
-  // by FiniteElement
+  // If the system is not primitive, these have not been initialized by
+  // FiniteElement
   this->system_to_component_table.resize(this->dofs_per_cell);
   this->face_system_to_component_table.resize(this->dofs_per_face);
 
@@ -1217,20 +1279,16 @@ FESystem<dim,spacedim>::build_cell_tables()
   Assert (total_index == this->component_to_base_table.size(),
           ExcInternalError());
 
-  // Initialize index tables.
-  // Multi-component base elements
-  // have to be thought of. For
-  // non-primitive shape functions,
-  // have a special invalid index.
+  // Initialize index tables.  Multi-component base elements have to be
+  // thought of. For non-primitive shape functions, have a special invalid
+  // index.
   const std::pair<unsigned int, unsigned int>
   non_primitive_index (numbers::invalid_unsigned_int,
                        numbers::invalid_unsigned_int);
 
-  // First enumerate vertex indices,
-  // where we first enumerate all
-  // indices on the first vertex in
-  // the order of the base elements,
-  // then of the second vertex, etc
+  // First enumerate vertex indices, where we first enumerate all indices on
+  // the first vertex in the order of the base elements, then of the second
+  // vertex, etc
   total_index = 0;
   for (unsigned int vertex_number=0;
        vertex_number<GeometryInfo<dim>::vertices_per_cell;
@@ -1388,14 +1446,10 @@ template <int dim, int spacedim>
 void
 FESystem<dim,spacedim>::build_face_tables()
 {
-  // Initialize index tables. do this
-  // in the same way as done for the
-  // cell tables, except that we now
-  // loop over the objects of faces
-
-  // For non-primitive shape
-  // functions, have a special
-  // invalid index
+  // Initialize index tables. do this in the same way as done for the cell
+  // tables, except that we now loop over the objects of faces
+
+  // For non-primitive shape functions, have a special invalid index
   const std::pair<unsigned int, unsigned int>
   non_primitive_index (numbers::invalid_unsigned_int,
                        numbers::invalid_unsigned_int);
@@ -1414,30 +1468,13 @@ FESystem<dim,spacedim>::build_face_tables()
                local_index < base_element(base).dofs_per_vertex;
                ++local_index, ++total_index)
             {
-              // get (cell) index of
-              // this shape function
-              // inside the base
-              // element to see
-              // whether the shape
-              // function is
-              // primitive (assume
-              // that all shape
-              // functions on
-              // vertices share the
-              // same primitivity
-              // property; assume
-              // likewise for all
-              // shape functions
-              // located on lines,
-              // quads, etc. this
-              // way, we can ask for
-              // primitivity of only
-              // _one_ shape
-              // function, which is
-              // taken as
-              // representative for
-              // all others located
-              // on the same type of
+              // get (cell) index of this shape function inside the base
+              // element to see whether the shape function is primitive
+              // (assume that all shape functions on vertices share the same
+              // primitivity property; assume likewise for all shape functions
+              // located on lines, quads, etc. this way, we can ask for
+              // primitivity of only _one_ shape function, which is taken as
+              // representative for all others located on the same type of
               // object):
               const unsigned int index_in_base
                 = (base_element(base).dofs_per_vertex*vertex_number +
@@ -1480,9 +1517,7 @@ FESystem<dim,spacedim>::build_face_tables()
                  local_index < base_element(base).dofs_per_line;
                  ++local_index, ++total_index)
               {
-                // do everything
-                // alike for this
-                // type of object
+                // do everything alike for this type of object
                 const unsigned int index_in_base
                   = (base_element(base).dofs_per_line*line_number +
                      local_index +
@@ -1526,9 +1561,7 @@ FESystem<dim,spacedim>::build_face_tables()
                  local_index < base_element(base).dofs_per_quad;
                  ++local_index, ++total_index)
               {
-                // do everything
-                // alike for this
-                // type of object
+                // do everything alike for this type of object
                 const unsigned int index_in_base
                   = (base_element(base).dofs_per_quad*quad_number +
                      local_index +
@@ -1569,14 +1602,10 @@ FESystem<dim,spacedim>::build_face_tables()
 template <int dim, int spacedim>
 void FESystem<dim,spacedim>::build_interface_constraints ()
 {
-  // check whether all base elements
-  // implement their interface
-  // constraint matrices. if this is
-  // not the case, then leave the
-  // interface costraints of this
-  // composed element empty as well;
-  // however, the rest of the element
-  // is usable
+  // check whether all base elements implement their interface constraint
+  // matrices. if this is not the case, then leave the interface costraints of
+  // this composed element empty as well; however, the rest of the element is
+  // usable
   for (unsigned int base=0; base<this->n_base_elements(); ++base)
     if (base_element(base).constraints_are_implemented() == false)
       return;
@@ -1584,58 +1613,35 @@ void FESystem<dim,spacedim>::build_interface_constraints ()
   this->interface_constraints.
   TableBase<2,double>::reinit (this->interface_constraints_size());
 
-  // the layout of the constraints
-  // matrix is described in the
-  // FiniteElement class. you may
-  // want to look there first before
-  // trying to understand the
-  // following, especially the
-  // mapping of the @p{m} index.
+  // the layout of the constraints matrix is described in the FiniteElement
+  // class. you may want to look there first before trying to understand the
+  // following, especially the mapping of the @p{m} index.
   //
-  // in order to map it to the
-  // fe-system class, we have to know
-  // which base element a degree of
-  // freedom within a vertex, line,
-  // etc belongs to. this can be
-  // accomplished by the
-  // system_to_component_index
-  // function in conjunction with the
-  // numbers
-  // first_{line,quad,...}_index
+  // in order to map it to the fe-system class, we have to know which base
+  // element a degree of freedom within a vertex, line, etc belongs to. this
+  // can be accomplished by the system_to_component_index function in
+  // conjunction with the numbers first_{line,quad,...}_index
   for (unsigned int n=0; n<this->interface_constraints.n(); ++n)
     for (unsigned int m=0; m<this->interface_constraints.m(); ++m)
       {
-        // for the pair (n,m) find
-        // out which base element
-        // they belong to and the
-        // number therein
+        // for the pair (n,m) find out which base element they belong to and
+        // the number therein
         //
-        // first for the n
-        // index. this is simple
-        // since the n indices are in
-        // the same order as they are
-        // usually on a face. note
-        // that for the data type,
-        // first value in pair is
-        // (base element,instance of
-        // base element), second is
-        // index within this instance
+        // first for the n index. this is simple since the n indices are in
+        // the same order as they are usually on a face. note that for the
+        // data type, first value in pair is (base element,instance of base
+        // element), second is index within this instance
         const std::pair<std::pair<unsigned int,unsigned int>, unsigned int> n_index
           = this->face_system_to_base_table[n];
 
-        // likewise for the m
-        // index. this is more
-        // complicated due to the
-        // strange ordering we have
-        // for the dofs on the
-        // refined faces.
+        // likewise for the m index. this is more complicated due to the
+        // strange ordering we have for the dofs on the refined faces.
         std::pair<std::pair<unsigned int,unsigned int>, unsigned int> m_index;
         switch (dim)
           {
           case 1:
           {
-            // we should never get here!
-            // (in 1d, the constraints matrix
+            // we should never get here!  (in 1d, the constraints matrix
             // should be of size zero)
             Assert (false, ExcInternalError());
             break;
@@ -1643,21 +1649,14 @@ void FESystem<dim,spacedim>::build_interface_constraints ()
 
           case 2:
           {
-            // the indices m=0..d_v-1 are
-            // from the center vertex.
-            // their order is the same
-            // as for the first vertex
-            // of the whole cell, so we
-            // can use the
-            // system_to_base_table
-            // variable (using the
-            // face_s_t_base_t function would
-            // yield the same)
+            // the indices m=0..d_v-1 are from the center vertex.  their order
+            // is the same as for the first vertex of the whole cell, so we
+            // can use the system_to_base_table variable (using the
+            // face_s_t_base_t function would yield the same)
             if (m < this->dofs_per_vertex)
               m_index = this->system_to_base_table[m];
             else
-              // then come the two sets of
-              // line indices
+              // then come the two sets of line indices
               {
                 const unsigned int index_in_line
                   = (m-this->dofs_per_vertex) % this->dofs_per_line;
@@ -1665,49 +1664,21 @@ void FESystem<dim,spacedim>::build_interface_constraints ()
                   = (m-this->dofs_per_vertex) / this->dofs_per_line;
                 Assert (sub_line < 2, ExcInternalError());
 
-                // from this
-                // information, try
-                // to get base
-                // element and
-                // instance of base
-                // element. we do
-                // so by
-                // constructing the
-                // corresponding
-                // face index of m
-                // in the present
-                // element, then
-                // use
+                // from this information, try to get base element and instance
+                // of base element. we do so by constructing the corresponding
+                // face index of m in the present element, then use
                 // face_system_to_base_table
                 const unsigned int tmp1 = 2*this->dofs_per_vertex+index_in_line;
                 m_index.first = this->face_system_to_base_table[tmp1].first;
 
-                // what we are
-                // still missing is
-                // the index of m
-                // within the base
-                // elements
-                // interface_constraints
-                // table
+                // what we are still missing is the index of m within the base
+                // elements interface_constraints table
                 //
-                // here, the second
-                // value of
-                // face_system_to_base_table
-                // can help: it
-                // denotes the face
-                // index of that
-                // shape function
-                // within the base
-                // element. since
-                // we know that it
-                // is a line dof,
-                // we can construct
-                // the rest: tmp2
-                // will denote the
-                // index of this
-                // shape function
-                // among the line
-                // shape functions:
+                // here, the second value of face_system_to_base_table can
+                // help: it denotes the face index of that shape function
+                // within the base element. since we know that it is a line
+                // dof, we can construct the rest: tmp2 will denote the index
+                // of this shape function among the line shape functions:
                 Assert (this->face_system_to_base_table[tmp1].second >=
                         2*base_element(m_index.first.first).dofs_per_vertex,
                         ExcInternalError());
@@ -1724,31 +1695,18 @@ void FESystem<dim,spacedim>::build_interface_constraints ()
 
           case 3:
           {
-            // same way as above,
-            // although a little
-            // more complicated...
-
-            // the indices
-            // m=0..5*d_v-1 are
-            // from the center and
-            // the four subline
-            // vertices.  their
-            // order is the same as
-            // for the first vertex
-            // of the whole cell,
-            // so we can use the
-            // simple arithmetic
+            // same way as above, although a little more complicated...
+
+            // the indices m=0..5*d_v-1 are from the center and the four
+            // subline vertices.  their order is the same as for the first
+            // vertex of the whole cell, so we can use the simple arithmetic
             if (m < 5*this->dofs_per_vertex)
               m_index = this->system_to_base_table[m];
             else
-              // then come the 12 sets of
-              // line indices
+              // then come the 12 sets of line indices
               if (m < 5*this->dofs_per_vertex + 12*this->dofs_per_line)
                 {
-                  // for the
-                  // meaning of all
-                  // this, see the
-                  // 2d part
+                  // for the meaning of all this, see the 2d part
                   const unsigned int index_in_line
                     = (m-5*this->dofs_per_vertex) % this->dofs_per_line;
                   const unsigned int sub_line
@@ -1772,10 +1730,7 @@ void FESystem<dim,spacedim>::build_interface_constraints ()
               else
                 // on one of the four sub-quads
                 {
-                  // for the
-                  // meaning of all
-                  // this, see the
-                  // 2d part
+                  // for the meaning of all this, see the 2d part
                   const unsigned int index_in_quad
                     = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) %
                       this->dofs_per_quad;
@@ -1815,14 +1770,9 @@ void FESystem<dim,spacedim>::build_interface_constraints ()
             Assert (false, ExcNotImplemented());
           };
 
-        // now that we gathered all
-        // information: use it to
-        // build the matrix. note
-        // that if n and m belong to
-        // different base elements or
-        // instances, then there
-        // definitely will be no
-        // coupling
+        // now that we gathered all information: use it to build the
+        // matrix. note that if n and m belong to different base elements or
+        // instances, then there definitely will be no coupling
         if (n_index.first == m_index.first)
           this->interface_constraints(m,n)
             = (base_element(n_index.first.first).constraints()(m_index.second,
@@ -1869,106 +1819,13 @@ void FESystem<dim,spacedim>::initialize (const std::vector<const FiniteElement<d
   build_cell_tables();
   build_face_tables();
 
-  for (unsigned int ref_case=RefinementCase<dim>::cut_x;
-       ref_case<RefinementCase<dim>::isotropic_refinement+1;
-       ++ref_case)
-    {
-      // Check if some of the matrices of
-      // the base elements are void.
-      // repeat this check for each RefineCase
-      bool do_restriction = true;
-      bool do_prolongation = true;
-
-      for (unsigned int i=0; i<this->n_base_elements(); ++i)
-        {
-          if (base_element(i).restriction[ref_case-1][0].n() != base_element(i).dofs_per_cell)
-            do_restriction = false;
-          if (base_element(i).prolongation[ref_case-1][0].n() != base_element(i).dofs_per_cell)
-            do_prolongation = false;
-        }
+  // restriction and prolongation matrices are build on demand
 
-      // if we did not encounter void
-      // matrices, initialize the
-      // respective matrix sizes
-      if (do_restriction)
-        for (unsigned int i=0; i<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++i)
-          this->restriction[ref_case-1][i].reinit(this->dofs_per_cell,
-                                                  this->dofs_per_cell);
-      if (do_prolongation)
-        for (unsigned int i=0; i<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++i)
-          this->prolongation[ref_case-1][i].reinit(this->dofs_per_cell,
-                                                   this->dofs_per_cell);
-
-      // distribute the matrices of the
-      // base finite elements to the
-      // matrices of this object. for
-      // this, loop over all degrees of
-      // freedom and take the respective
-      // entry of the underlying base
-      // element.
-      //
-      // note that we by definition of a
-      // base element, they are
-      // independent, i.e. do not
-      // couple. only DoFs that belong to
-      // the same instance of a base
-      // element may couple
-      for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-        for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-          {
-            // first find out to which
-            // base element indices i and
-            // j belong, and which
-            // instance thereof in case
-            // the base element has a
-            // multiplicity greater than
-            // one. if they should not
-            // happen to belong to the
-            // same instance of a base
-            // element, then they cannot
-            // couple, so go on with the
-            // next index
-            if (this->system_to_base_table[i].first !=
-                this->system_to_base_table[j].first)
-              continue;
-
-            // so get the common base
-            // element and the indices
-            // therein:
-            const unsigned int
-            base = this->system_to_base_table[i].first.first;
-
-            const unsigned int
-            base_index_i = this->system_to_base_table[i].second,
-            base_index_j = this->system_to_base_table[j].second;
-
-            // if we are sure that DoFs i
-            // and j may couple, then
-            // copy entries of the
-            // matrices:
-            for (unsigned int child=0; child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case)); ++child)
-              {
-                if (do_restriction)
-                  this->restriction[ref_case-1][child] (i,j)
-                    = (base_element(base)
-                       .get_restriction_matrix(child, RefinementCase<dim>(ref_case))(
-                         base_index_i,base_index_j));
-
-                if (do_prolongation)
-                  this->prolongation[ref_case-1][child] (i,j)
-                    = (base_element(base)
-                       .get_prolongation_matrix(child, RefinementCase<dim>(ref_case))(
-                         base_index_i,base_index_j));
-              };
-          };
-    }
-  // now set up the interface constraints.
-  // this is kind'o hairy, so don't try
+  // now set up the interface constraints.  this is kind'o hairy, so don't try
   // to do it dimension independent
   build_interface_constraints ();
 
-  // finally fill in support points
-  // on cell and face
+  // finally fill in support points on cell and face
   initialize_unit_support_points ();
   initialize_unit_face_support_points ();
 
@@ -1985,14 +1842,9 @@ void
 FESystem<dim,spacedim>::
 initialize_unit_support_points ()
 {
-  // if one of the base elements
-  // has no support points, then
-  // it makes no sense to define
-  // support points for the
-  // composed element, so return
-  // an empty array to
-  // demonstrate that
-  // fact
+  // if one of the base elements has no support points, then it makes no sense
+  // to define support points for the composed element, so return an empty
+  // array to demonstrate that fact
   for (unsigned int base_el=0; base_el<this->n_base_elements(); ++base_el)
     if (!base_element(base_el).has_support_points())
       {
@@ -2000,9 +1852,7 @@ initialize_unit_support_points ()
         return;
       };
 
-  // generate unit support points
-  // from unit support points of sub
-  // elements
+  // generate unit support points from unit support points of sub elements
   this->unit_support_points.resize(this->dofs_per_cell);
 
   for (unsigned int i=0; i<this->dofs_per_cell; ++i)
@@ -2028,22 +1878,15 @@ initialize_unit_face_support_points ()
   if (dim == 1)
     return;
 
-  // if one of the base elements has no
-  // support points, then it makes no sense
-  // to define support points for the
-  // composed element. In that case, return
-  // an empty array to demonstrate that fact
-  // (note that we ask whether the base
-  // element has no support points at all,
-  // not only none on the face!)
+  // if one of the base elements has no support points, then it makes no sense
+  // to define support points for the composed element. In that case, return
+  // an empty array to demonstrate that fact (note that we ask whether the
+  // base element has no support points at all, not only none on the face!)
   //
-  // on the other hand, if there is an
-  // element that simply has no degrees of
-  // freedom on the face at all, then we
-  // don't care whether it has support points
-  // or not. this is, for example, the case
-  // for the stable Stokes element Q(p)^dim
-  // \times DGP(p-1).
+  // on the other hand, if there is an element that simply has no degrees of
+  // freedom on the face at all, then we don't care whether it has support
+  // points or not. this is, for example, the case for the stable Stokes
+  // element Q(p)^dim \times DGP(p-1).
   for (unsigned int base_el=0; base_el<this->n_base_elements(); ++base_el)
     if (!base_element(base_el).has_support_points()
         &&
@@ -2054,8 +1897,7 @@ initialize_unit_face_support_points ()
       }
 
 
-  // generate unit face support points
-  // from unit support points of sub
+  // generate unit face support points from unit support points of sub
   // elements
   this->unit_face_support_points.resize(this->dofs_per_face);
 
@@ -2082,14 +1924,13 @@ FESystem<dim,spacedim>::initialize_quad_dof_index_permutation ()
   if (dim < 3)
     return;
 
-  // the array into which we want to write
-  // should have the correct size already.
+  // the array into which we want to write should have the correct size
+  // already.
   Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()==
           8*this->dofs_per_quad, ExcInternalError());
 
-  // to obtain the shifts for this composed
-  // element, copy the shift information of
-  // the base elements
+  // to obtain the shifts for this composed element, copy the shift
+  // information of the base elements
   unsigned int index = 0;
   for (unsigned int b=0; b<this->n_base_elements(); ++b)
     {
@@ -2107,8 +1948,7 @@ FESystem<dim,spacedim>::initialize_quad_dof_index_permutation ()
   Assert (index == this->dofs_per_quad,
           ExcInternalError());
 
-  // aditionally compose the permutation
-  // information for lines
+  // aditionally compose the permutation information for lines
   Assert (this->adjust_line_dof_index_for_line_orientation_table.size()==
           this->dofs_per_line, ExcInternalError());
   index = 0;
@@ -2164,31 +2004,23 @@ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
           ExcDimensionMismatch (interpolation_matrix.m(),
                                 x_source_fe.dofs_per_face));
 
-  // since dofs for each base are
-  // independent, we only have to stack
-  // things up from base element to base
-  // element
+  // since dofs for each base are independent, we only have to stack things up
+  // from base element to base element
   //
-  // the problem is that we have to work with
-  // two FEs (this and fe_other). only deal
-  // with the case that both are FESystems
-  // and that they both have the same number
-  // of bases (counting multiplicity) each of
-  // which match in their number of
-  // components. this covers
-  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
-  // FESystem(FE_Q(r),2,FE_Q(s),1), but not
+  // the problem is that we have to work with two FEs (this and
+  // fe_other). only deal with the case that both are FESystems and that they
+  // both have the same number of bases (counting multiplicity) each of which
+  // match in their number of components. this covers
+  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
   // FESystem(FE_Q(p),1,FE_Q(q),2) vs
   // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
   const FESystem<dim,spacedim> *fe_other_system
     = dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe);
 
-  // clear matrix, since we will not get to
-  // set all elements
+  // clear matrix, since we will not get to set all elements
   interpolation_matrix = 0;
 
-  // loop over all the base elements of
-  // this and the other element, counting
+  // loop over all the base elements of this and the other element, counting
   // their multiplicities
   unsigned int base_index       = 0,
                base_index_other = 0;
@@ -2212,11 +2044,9 @@ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
       base.get_face_interpolation_matrix (base_other,
                                           base_to_base_interpolation);
 
-      // now translate entries. we'd like to
-      // have something like
-      // face_base_to_system_index, but that
-      // doesn't exist. rather, all we have
-      // is the reverse. well, use that then
+      // now translate entries. we'd like to have something like
+      // face_base_to_system_index, but that doesn't exist. rather, all we
+      // have is the reverse. well, use that then
       for (unsigned int i=0; i<this->dofs_per_face; ++i)
         if (this->face_system_to_base_index(i).first
             ==
@@ -2229,12 +2059,9 @@ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
                 = base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second,
                                              this->face_system_to_base_index(i).second);
 
-      // advance to the next base element
-      // for this and the other
-      // fe_system; see if we can simply
-      // advance the multiplicity by one,
-      // or if have to move on to the
-      // next base element
+      // advance to the next base element for this and the other fe_system;
+      // see if we can simply advance the multiplicity by one, or if have to
+      // move on to the next base element
       ++multiplicity;
       if (multiplicity == this->element_multiplicity(base_index))
         {
@@ -2249,10 +2076,8 @@ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
           ++base_index_other;
         }
 
-      // see if we have reached the end
-      // of the present element. if so,
-      // we should have reached the end
-      // of the other one as well
+      // see if we have reached the end of the present element. if so, we
+      // should have reached the end of the other one as well
       if (base_index == this->n_base_elements())
         {
           Assert (base_index_other == fe_other_system->n_base_elements(),
@@ -2260,10 +2085,8 @@ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
           break;
         }
 
-      // if we haven't reached the end of
-      // this element, we shouldn't have
-      // reached the end of the other one
-      // either
+      // if we haven't reached the end of this element, we shouldn't have
+      // reached the end of the other one either
       Assert (base_index_other != fe_other_system->n_base_elements(),
               ExcInternalError());
     }
@@ -2292,31 +2115,23 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
           ExcDimensionMismatch (interpolation_matrix.m(),
                                 x_source_fe.dofs_per_face));
 
-  // since dofs for each base are
-  // independent, we only have to stack
-  // things up from base element to base
-  // element
+  // since dofs for each base are independent, we only have to stack things up
+  // from base element to base element
   //
-  // the problem is that we have to work with
-  // two FEs (this and fe_other). only deal
-  // with the case that both are FESystems
-  // and that they both have the same number
-  // of bases (counting multiplicity) each of
-  // which match in their number of
-  // components. this covers
-  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
-  // FESystem(FE_Q(r),2,FE_Q(s),1), but not
+  // the problem is that we have to work with two FEs (this and
+  // fe_other). only deal with the case that both are FESystems and that they
+  // both have the same number of bases (counting multiplicity) each of which
+  // match in their number of components. this covers
+  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
   // FESystem(FE_Q(p),1,FE_Q(q),2) vs
   // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
   const FESystem<dim,spacedim> *fe_other_system
     = dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe);
 
-  // clear matrix, since we will not get to
-  // set all elements
+  // clear matrix, since we will not get to set all elements
   interpolation_matrix = 0;
 
-  // loop over all the base elements of
-  // this and the other element, counting
+  // loop over all the base elements of this and the other element, counting
   // their multiplicities
   unsigned int base_index       = 0,
                base_index_other = 0;
@@ -2341,11 +2156,9 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
                                              subface,
                                              base_to_base_interpolation);
 
-      // now translate entries. we'd like to
-      // have something like
-      // face_base_to_system_index, but that
-      // doesn't exist. rather, all we have
-      // is the reverse. well, use that then
+      // now translate entries. we'd like to have something like
+      // face_base_to_system_index, but that doesn't exist. rather, all we
+      // have is the reverse. well, use that then
       for (unsigned int i=0; i<this->dofs_per_face; ++i)
         if (this->face_system_to_base_index(i).first
             ==
@@ -2358,12 +2171,9 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
                 = base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second,
                                              this->face_system_to_base_index(i).second);
 
-      // advance to the next base element
-      // for this and the other
-      // fe_system; see if we can simply
-      // advance the multiplicity by one,
-      // or if have to move on to the
-      // next base element
+      // advance to the next base element for this and the other fe_system;
+      // see if we can simply advance the multiplicity by one, or if have to
+      // move on to the next base element
       ++multiplicity;
       if (multiplicity == this->element_multiplicity(base_index))
         {
@@ -2378,10 +2188,8 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
           ++base_index_other;
         }
 
-      // see if we have reached the end
-      // of the present element. if so,
-      // we should have reached the end
-      // of the other one as well
+      // see if we have reached the end of the present element. if so, we
+      // should have reached the end of the other one as well
       if (base_index == this->n_base_elements())
         {
           Assert (base_index_other == fe_other_system->n_base_elements(),
@@ -2389,10 +2197,8 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
           break;
         }
 
-      // if we haven't reached the end of
-      // this element, we shouldn't have
-      // reached the end of the other one
-      // either
+      // if we haven't reached the end of this element, we shouldn't have
+      // reached the end of the other one either
       Assert (base_index_other != fe_other_system->n_base_elements(),
               ExcInternalError());
     }
@@ -2405,38 +2211,29 @@ template <int structdim>
 std::vector<std::pair<unsigned int, unsigned int> >
 FESystem<dim,spacedim>::hp_object_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
 {
-  // since dofs on each subobject (vertex,
-  // line, ...) are ordered such that first
-  // come all from the first base element all
-  // multiplicities, then second base element
-  // all multiplicities, etc., we simply have
-  // to stack all the identities after each
-  // other
+  // since dofs on each subobject (vertex, line, ...) are ordered such that
+  // first come all from the first base element all multiplicities, then
+  // second base element all multiplicities, etc., we simply have to stack all
+  // the identities after each other
   //
-  // the problem is that we have to work with
-  // two FEs (this and fe_other). only deal
-  // with the case that both are FESystems
-  // and that they both have the same number
-  // of bases (counting multiplicity) each of
-  // which match in their number of
-  // components. this covers
-  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
-  // FESystem(FE_Q(r),2,FE_Q(s),1), but not
+  // the problem is that we have to work with two FEs (this and
+  // fe_other). only deal with the case that both are FESystems and that they
+  // both have the same number of bases (counting multiplicity) each of which
+  // match in their number of components. this covers
+  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
   // FESystem(FE_Q(p),1,FE_Q(q),2) vs
   // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
   if (const FESystem<dim,spacedim> *fe_other_system
       = dynamic_cast<const FESystem<dim,spacedim>*>(&fe_other))
     {
-      // loop over all the base elements of
-      // this and the other element, counting
-      // their multiplicities
+      // loop over all the base elements of this and the other element,
+      // counting their multiplicities
       unsigned int base_index       = 0,
                    base_index_other = 0;
       unsigned int multiplicity       = 0,
                    multiplicity_other = 0;
 
-      // we also need to keep track of the
-      // number of dofs already treated for
+      // we also need to keep track of the number of dofs already treated for
       // each of the elements
       unsigned int dof_offset       = 0,
                    dof_offset_other = 0;
@@ -2452,10 +2249,8 @@ FESystem<dim,spacedim>::hp_object_dof_identities (const FiniteElement<dim,spaced
           Assert (base.n_components() == base_other.n_components(),
                   ExcNotImplemented());
 
-          // now translate the identities
-          // returned by the base elements to
-          // the indices of this system
-          // element
+          // now translate the identities returned by the base elements to the
+          // indices of this system element
           std::vector<std::pair<unsigned int, unsigned int> > base_identities;
           switch (structdim)
             {
@@ -2477,17 +2272,13 @@ FESystem<dim,spacedim>::hp_object_dof_identities (const FiniteElement<dim,spaced
             (std::make_pair (base_identities[i].first + dof_offset,
                              base_identities[i].second + dof_offset_other));
 
-          // record the dofs treated above as
-          // already taken care of
+          // record the dofs treated above as already taken care of
           dof_offset       += base.template n_dofs_per_object<structdim>();
           dof_offset_other += base_other.template n_dofs_per_object<structdim>();
 
-          // advance to the next base element
-          // for this and the other
-          // fe_system; see if we can simply
-          // advance the multiplicity by one,
-          // or if have to move on to the
-          // next base element
+          // advance to the next base element for this and the other
+          // fe_system; see if we can simply advance the multiplicity by one,
+          // or if have to move on to the next base element
           ++multiplicity;
           if (multiplicity == this->element_multiplicity(base_index))
             {
@@ -2502,10 +2293,8 @@ FESystem<dim,spacedim>::hp_object_dof_identities (const FiniteElement<dim,spaced
               ++base_index_other;
             }
 
-          // see if we have reached the end
-          // of the present element. if so,
-          // we should have reached the end
-          // of the other one as well
+          // see if we have reached the end of the present element. if so, we
+          // should have reached the end of the other one as well
           if (base_index == this->n_base_elements())
             {
               Assert (base_index_other == fe_other_system->n_base_elements(),
@@ -2513,10 +2302,8 @@ FESystem<dim,spacedim>::hp_object_dof_identities (const FiniteElement<dim,spaced
               break;
             }
 
-          // if we haven't reached the end of
-          // this element, we shouldn't have
-          // reached the end of the other one
-          // either
+          // if we haven't reached the end of this element, we shouldn't have
+          // reached the end of the other one either
           Assert (base_index_other != fe_other_system->n_base_elements(),
                   ExcInternalError());
         }
@@ -2562,9 +2349,8 @@ FiniteElementDomination::Domination
 FESystem<dim,spacedim>::
 compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
 {
-  // at present all we can do is to compare
-  // with other FESystems that have the same
-  // number of components and bases
+  // at present all we can do is to compare with other FESystems that have the
+  // same number of components and bases
   if (const FESystem<dim,spacedim> *fe_sys_other
       = dynamic_cast<const FESystem<dim,spacedim>*>(&fe_other))
     {
@@ -2576,8 +2362,7 @@ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
       FiniteElementDomination::Domination
       domination = FiniteElementDomination::no_requirements;
 
-      // loop over all base elements and do
-      // some sanity checks
+      // loop over all base elements and do some sanity checks
       for (unsigned int b=0; b<this->n_base_elements(); ++b)
         {
           Assert (this->base_element(b).n_components() ==
@@ -2587,8 +2372,7 @@ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
                   fe_sys_other->element_multiplicity(b),
                   ExcNotImplemented());
 
-          // for this pair of base elements,
-          // check who dominates and combine
+          // for this pair of base elements, check who dominates and combine
           // with previous result
           const FiniteElementDomination::Domination
           base_domination = (this->base_element(b)
@@ -2596,8 +2380,7 @@ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
           domination = domination & base_domination;
         }
 
-      // if we've gotten here, then we've
-      // either found a winner or either
+      // if we've gotten here, then we've either found a winner or either
       // element is fine being dominated
       Assert (domination !=
               FiniteElementDomination::neither_element_dominates,
@@ -2756,34 +2539,23 @@ compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim
 {
   Assert (fes.size() == multiplicities.size(), ExcInternalError());
 
-  // first count the number of dofs
-  // and components that will emerge
-  // from the given FEs
+  // first count the number of dofs and components that will emerge from the
+  // given FEs
   unsigned int n_shape_functions = 0;
   for (unsigned int i=0; i<fes.size(); ++i)
     if (multiplicities[i]>0) // check needed as fe might be NULL
       n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
 
-  // generate the array that will
-  // hold the output
+  // generate the array that will hold the output
   std::vector<bool> retval (n_shape_functions, false);
 
-  // finally go through all the shape
-  // functions of the base elements,
-  // and copy their flags. this
-  // somehow copies the code in
-  // build_cell_table, which is not
-  // nice as it uses too much
-  // implicit knowledge about the
-  // layout of the individual bases
-  // in the composed FE, but there
-  // seems no way around...
+  // finally go through all the shape functions of the base elements, and copy
+  // their flags. this somehow copies the code in build_cell_table, which is
+  // not nice as it uses too much implicit knowledge about the layout of the
+  // individual bases in the composed FE, but there seems no way around...
   //
-  // for each shape function, copy
-  // the flags from the base element
-  // to this one, taking into account
-  // multiplicities, and other
-  // complications
+  // for each shape function, copy the flags from the base element to this
+  // one, taking into account multiplicities, and other complications
   unsigned int total_index = 0;
   for (unsigned int vertex_number=0;
        vertex_number<GeometryInfo<dim>::vertices_per_cell;
@@ -2925,9 +2697,8 @@ compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*
 {
   Assert (fes.size() == multiplicities.size(), ExcInternalError());
 
-  // first count the number of dofs
-  // and components that will emerge
-  // from the given FEs
+  // first count the number of dofs and components that will emerge from the
+  // given FEs
   unsigned int n_shape_functions = 0;
   for (unsigned int i=0; i<fes.size(); ++i)
     if (multiplicities[i]>0) //needed because fe might be NULL
@@ -2938,28 +2709,18 @@ compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*
     if (multiplicities[i]>0) //needed because fe might be NULL
       n_components += fes[i]->n_components() * multiplicities[i];
 
-  // generate the array that will
-  // hold the output
+  // generate the array that will hold the output
   std::vector<std::vector<bool> >
   retval (n_shape_functions, std::vector<bool> (n_components, false));
 
-  // finally go through all the shape
-  // functions of the base elements,
-  // and copy their flags. this
-  // somehow copies the code in
-  // build_cell_table, which is not
-  // nice as it uses too much
-  // implicit knowledge about the
-  // layout of the individual bases
-  // in the composed FE, but there
-  // seems no way around...
+  // finally go through all the shape functions of the base elements, and copy
+  // their flags. this somehow copies the code in build_cell_table, which is
+  // not nice as it uses too much implicit knowledge about the layout of the
+  // individual bases in the composed FE, but there seems no way around...
   //
-  // for each shape function, copy
-  // the non-zero flags from the base
-  // element to this one, taking into
-  // account multiplicities, multiple
-  // components in base elements, and
-  // other complications
+  // for each shape function, copy the non-zero flags from the base element to
+  // this one, taking into account multiplicities, multiple components in base
+  // elements, and other complications
   unsigned int total_index = 0;
   for (unsigned int vertex_number=0;
        vertex_number<GeometryInfo<dim>::vertices_per_cell;
@@ -3143,14 +2904,12 @@ FESystem<dim,spacedim>::unit_support_point (const unsigned int index) const
           (this->unit_support_points.size() == 0),
           typename FEL::ExcFEHasNoSupportPoints ());
 
-  // let's see whether we have the
-  // information pre-computed
+  // let's see whether we have the information pre-computed
   if (this->unit_support_points.size() != 0)
     return this->unit_support_points[index];
   else
-    // no. ask the base element
-    // whether it would like to
-    // provide this information
+    // no. ask the base element whether it would like to provide this
+    // information
     return (base_element(this->system_to_base_index(index).first.first)
             .unit_support_point(this->system_to_base_index(index).second));
 }
@@ -3168,14 +2927,12 @@ FESystem<dim,spacedim>::unit_face_support_point (const unsigned int index) const
           (this->unit_face_support_points.size() == 0),
           typename FEL::ExcFEHasNoSupportPoints ());
 
-  // let's see whether we have the
-  // information pre-computed
+  // let's see whether we have the information pre-computed
   if (this->unit_face_support_points.size() != 0)
     return this->unit_face_support_points[index];
   else
-    // no. ask the base element
-    // whether it would like to
-    // provide this information
+    // no. ask the base element whether it would like to provide this
+    // information
     return (base_element(this->face_system_to_base_index(index).first.first)
             .unit_face_support_point(this->face_system_to_base_index(index).second));
 }
@@ -3186,11 +2943,8 @@ template <int dim, int spacedim>
 std::size_t
 FESystem<dim,spacedim>::memory_consumption () const
 {
-  // neglect size of data stored in
-  // @p{base_elements} due to some
-  // problems with the
-  // compiler. should be neglectable
-  // after all, considering the size
+  // neglect size of data stored in @p{base_elements} due to some problems
+  // with the compiler. should be neglectable after all, considering the size
   // of the data of the subelements
   std::size_t mem = (FiniteElement<dim,spacedim>::memory_consumption () +
                      sizeof (base_elements));