]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a few out-of-bound array accesses. Refill some comments.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 28 Dec 2012 10:41:54 +0000 (10:41 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 28 Dec 2012 10:41:54 +0000 (10:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@27865 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/matrix_free/dof_info.h
deal.II/include/deal.II/matrix_free/dof_info.templates.h
deal.II/include/deal.II/matrix_free/mapping_info.h
deal.II/include/deal.II/matrix_free/mapping_info.templates.h
deal.II/include/deal.II/matrix_free/matrix_free.h
deal.II/include/deal.II/matrix_free/matrix_free.templates.h

index 36b008669165c0638026da08b99c005bcc7d55d6..90c7fccdabb9a907fc3f8f56580f7891f3d4344f 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2011 by the deal.II authors
+//    Copyright (C) 2011-2012 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -39,13 +39,13 @@ namespace internal
      * ConstraintMatrix.
      *
      * This class only stores index relations. The weights for hanging node
-     * constraints are stored in a different field. This is because a different
-     * field allows for the same compressed weight data on different DoFHandlers
-     * for vector-valued problems. There, the indices might be constrained
-     * differently on different components (e.g. Dirichlet conditions only on
-     * selected components), whereas the weights from hanging nodes are the same
-     * and need to be stored only once. The combination will be handled in the
-     * MatrixFree class.
+     * constraints are stored in a different field. This is because a
+     * different field allows for the same compressed weight data on different
+     * DoFHandlers for vector-valued problems. There, the indices might be
+     * constrained differently on different components (e.g. Dirichlet
+     * conditions only on selected components), whereas the weights from
+     * hanging nodes are the same and need to be stored only once. The
+     * combination will be handled in the MatrixFree class.
      *
      * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
      */
@@ -67,93 +67,78 @@ namespace internal
       void clear ();
 
       /**
-       * Returns a pointer to the first index in the
-       * DoF row @p row.
+       * Returns a pointer to the first index in the DoF row @p row.
        */
       const unsigned int *begin_indices (const unsigned int row) const;
 
       /**
-       * Returns a pointer to the one past the last
-       * DoF index in the row @p row.
+       * Returns a pointer to the one past the last DoF index in the row @p
+       * row.
        */
       const unsigned int *end_indices (const unsigned int row) const;
 
       /**
-       * Returns the number of entries in the
-       * indices field for the given row.
+       * Returns the number of entries in the indices field for the given row.
        */
       unsigned int row_length_indices (const unsigned int row) const;
 
       /**
-       * Returns a pointer to the first constraint
-       * indicator in the row @p row.
+       * Returns a pointer to the first constraint indicator in the row @p
+       * row.
        */
       const std::pair<unsigned short,unsigned short> *
       begin_indicators (const unsigned int row) const;
 
       /**
-       * Returns a pointer to the one past the last
-       * constraint indicator in the row @p row.
+       * Returns a pointer to the one past the last constraint indicator in
+       * the row @p row.
        */
       const std::pair<unsigned short,unsigned short> *
       end_indicators (const unsigned int row) const;
 
       /**
-       * Returns the number of entries in the
-       * constraint indicator field for the given
-       * row.
+       * Returns the number of entries in the constraint indicator field for
+       * the given row.
        */
       unsigned int row_length_indicators (const unsigned int row) const;
 
       /**
-       * Returns a pointer to the first index in the
-       * DoF row @p row for plain indices (i.e., the
-       * entries where constraints are not
-       * embedded).
+       * Returns a pointer to the first index in the DoF row @p row for plain
+       * indices (i.e., the entries where constraints are not embedded).
        */
       const unsigned int *begin_indices_plain (const unsigned int row) const;
 
       /**
-       * Returns a pointer to the one past the last
-       * DoF index in the row @p row (i.e., the
-       * entries where constraints are not
-       * embedded).
+       * Returns a pointer to the one past the last DoF index in the row @p
+       * row (i.e., the entries where constraints are not embedded).
        */
       const unsigned int *end_indices_plain (const unsigned int row) const;
 
       /**
-       * Returns the FE index for a given finite
-       * element degree. If not in hp mode, this
-       * function always returns index 0. If an
-       * index is not found in hp mode, it returns
-       * max_fe_degree, i.e., one index past the
-       * last valid one.
+       * Returns the FE index for a given finite element degree. If not in hp
+       * mode, this function always returns index 0. If an index is not found
+       * in hp mode, it returns max_fe_degree, i.e., one index past the last
+       * valid one.
        */
       unsigned int fe_index_from_degree (const unsigned int fe_degree) const;
 
 
       /**
-       * Returns the FE index for a given finite
-       * element degree. If not in hp mode or if the
-       * index is not found, this function always
-       * returns index 0. Hence, this function does
-       * not check whether the given degree is
+       * Returns the FE index for a given finite element degree. If not in hp
+       * mode or if the index is not found, this function always returns index
+       * 0. Hence, this function does not check whether the given degree is
        * actually present.
        */
       unsigned int
       fe_index_from_dofs_per_cell (const unsigned int dofs_per_cell) const;
 
       /**
-       * This internal method takes the local
-       * indices on a cell and fills them into this
-       * class. It resolves the constraints and
-       * distributes the results. Ghost indices,
-       * i.e., indices that are located on another
-       * processor, get a temporary number by this
-       * function, and will later be assigned the
-       * correct index after all the ghost indices
-       * have been collected by the call to @p
-       * assign_ghosts.
+       * This internal method takes the local indices on a cell and fills them
+       * into this class. It resolves the constraints and distributes the
+       * results. Ghost indices, i.e., indices that are located on another
+       * processor, get a temporary number by this function, and will later be
+       * assigned the correct index after all the ghost indices have been
+       * collected by the call to @p assign_ghosts.
        */
       void read_dof_indices (const std::vector<unsigned int> &local_indices,
                              const std::vector<unsigned int> &lexicographic_inv,
@@ -163,60 +148,48 @@ namespace internal
                              bool                            &cell_at_boundary);
 
       /**
-       * This method assigns the correct indices to
-       * ghost indices from the temporary numbering
-       * employed by the @p read_dof_indices
-       * function. The numbers are localized with
-       * respect to the MPI process, and ghosts
-       * start at the end of the locally owned
-       * range. This way, we get direct access to
-       * all vector entries.
+       * This method assigns the correct indices to ghost indices from the
+       * temporary numbering employed by the @p read_dof_indices function. The
+       * numbers are localized with respect to the MPI process, and ghosts
+       * start at the end of the locally owned range. This way, we get direct
+       * access to all vector entries.
        */
       void assign_ghosts(const std::vector<unsigned int> &boundary_cells);
 
       /**
-       * Reorganizes cells for serial
-       * (non-thread-parallelized) such that
-       * boundary cells are places in the
-       * middle. This way, computations and
-       * communication can be overlapped. Should
-       * only be called by one DoFInfo object when
-       * used on a system of several DoFHandlers.
+       * Reorganizes cells for serial (non-thread-parallelized) such that
+       * boundary cells are places in the middle. This way, computations and
+       * communication can be overlapped. Should only be called by one DoFInfo
+       * object when used on a system of several DoFHandlers.
        */
       void compute_renumber_serial (const std::vector<unsigned int> &boundary_cells,
                                     const SizeInfo                  &size_info,
                                     std::vector<unsigned int>       &renumbering);
 
       /**
-       * Reorganizes cells in the hp case without
-       * parallelism such that all cells with the
-       * same FE index are placed
-       * consecutively. Should only be called by one
-       * DoFInfo object when used on a system of
-       * several DoFHandlers.
+       * Reorganizes cells in the hp case without parallelism such that all
+       * cells with the same FE index are placed consecutively. Should only be
+       * called by one DoFInfo object when used on a system of several
+       * DoFHandlers.
        */
       void compute_renumber_hp_serial (SizeInfo                  &size_info,
                                        std::vector<unsigned int> &renumbering,
                                        std::vector<unsigned int> &irregular_cells);
 
       /**
-       * Computes the initial renumbering of cells
-       * such that all cells with ghosts are put
-       * first. This is the first step before
-       * building the thread graph and used to
-       * overlap computations and communication.
+       * Computes the initial renumbering of cells such that all cells with
+       * ghosts are put first. This is the first step before building the
+       * thread graph and used to overlap computations and communication.
        */
       void compute_renumber_parallel (const std::vector<unsigned int> &boundary_cells,
                                       SizeInfo                        &size_info,
                                       std::vector<unsigned int>       &renumbering);
 
       /**
-       * This method reorders the way cells are gone
-       * through based on a given renumbering of the
-       * cells. It also takes @p
-       * vectorization_length cells together and
-       * interprets them as one cell only, as is
-       * needed for vectorization.
+       * This method reorders the way cells are gone through based on a given
+       * renumbering of the cells. It also takes @p vectorization_length cells
+       * together and interprets them as one cell only, as is needed for
+       * vectorization.
        */
       void reorder_cells (const SizeInfo                  &size_info,
                           const std::vector<unsigned int> &renumbering,
@@ -225,37 +198,27 @@ namespace internal
                           const unsigned int               vectorization_length);
 
       /**
-       * This helper function determines a block
-       * size if the user decided not to force a
-       * block size through
-       * MatrixFree::AdditionalData. This is
-       * computed based on the number of hardware
-       * threads on the system
-       * (multithread_info.n_default_threads) and
-       * the number of macro cells that we should
-       * work on.
+       * This helper function determines a block size if the user decided not
+       * to force a block size through MatrixFree::AdditionalData. This is
+       * computed based on the number of hardware threads on the system
+       * (multithread_info.n_default_threads) and the number of macro cells
+       * that we should work on.
        */
       void guess_block_size (const SizeInfo &size_info,
                              TaskInfo       &task_info);
 
       /**
-       * This method goes through all cells that
-       * have been filled into @p dof_indices and
-       * finds out which cells can be worked on
-       * independently and which ones are
-       * neighboring and need to be done at
-       * different times when used in parallel.
+       * This method goes through all cells that have been filled into @p
+       * dof_indices and finds out which cells can be worked on independently
+       * and which ones are neighboring and need to be done at different times
+       * when used in parallel.
        *
-       * The strategy is based on a two-level
-       * approach. The outer level is subdivided
-       * into partitions similar to the type of
-       * neighbors in Cuthill-McKee, and the inner
-       * level is subdivided via colors (for chunks
-       * within the same color, can work
-       * independently). One task is represented by
-       * a chunk of cells. The cell chunks are
-       * formed before subdivision into partitions
-       * and colors.
+       * The strategy is based on a two-level approach. The outer level is
+       * subdivided into partitions similar to the type of neighbors in
+       * Cuthill-McKee, and the inner level is subdivided via colors (for
+       * chunks within the same color, can work independently). One task is
+       * represented by a chunk of cells. The cell chunks are formed before
+       * subdivision into partitions and colors.
        */
       void
       make_thread_graph_partition_color (SizeInfo                  &size_info,
@@ -265,24 +228,18 @@ namespace internal
                                          const bool                 hp_bool);
 
       /**
-       * This function goes through all cells that
-       * have been filled into @p dof_indices and
-       * finds out which cells can be worked on
-       * independently and which ones are
-       * neighboring and need to be done at
-       * different times when used in parallel.
+       * This function goes through all cells that have been filled into @p
+       * dof_indices and finds out which cells can be worked on independently
+       * and which ones are neighboring and need to be done at different times
+       * when used in parallel.
        *
-       * The strategy is based on a two-level
-       * approach. The outer level is subdivided
-       * into partitions similar to the type of
-       * neighbors in Cuthill-McKee, and the inner
-       * level is again subdivided into
-       * Cuthill-McKee-like partitions (partitions
-       * whose level differs by more than 2 can be
-       * worked on independently). One task is
-       * represented by a chunk of cells. The cell
-       * chunks are formed after subdivision into
-       * the two levels of partitions.
+       * The strategy is based on a two-level approach. The outer level is
+       * subdivided into partitions similar to the type of neighbors in
+       * Cuthill-McKee, and the inner level is again subdivided into
+       * Cuthill-McKee-like partitions (partitions whose level differs by more
+       * than 2 can be worked on independently). One task is represented by a
+       * chunk of cells. The cell chunks are formed after subdivision into the
+       * two levels of partitions.
        */
       void
       make_thread_graph_partition_partition (SizeInfo                  &size_info,
@@ -292,12 +249,10 @@ namespace internal
                                              const bool                 hp_bool);
 
       /**
-       * This function computes the connectivity of
-       * the currently stored indices and fills the
-       * structure into a sparsity pattern. The
-       * parameter block_size can be used to specify
-       * whether several cells should be treated as
-       * one.
+       * This function computes the connectivity of the currently stored
+       * indices and fills the structure into a sparsity pattern. The
+       * parameter block_size can be used to specify whether several cells
+       * should be treated as one.
        */
       void
       make_connectivity_graph (const SizeInfo                  &size_info,
@@ -308,29 +263,25 @@ namespace internal
                                CompressedSimpleSparsityPattern &connectivity) const;
 
       /**
-       * Renumbers the degrees of freedom to give
-       * good access for this class.
+       * Renumbers the degrees of freedom to give good access for this class.
        */
       void renumber_dofs (std::vector<unsigned int> &renumbering);
 
       /**
-       * Returns the memory consumption in bytes of
-       * this class.
+       * Returns the memory consumption in bytes of this class.
        */
       std::size_t memory_consumption() const;
 
       /**
-       * Prints a detailed summary of memory
-       * consumption in the different structures of
-       * this class to the given output stream.
+       * Prints a detailed summary of memory consumption in the different
+       * structures of this class to the given output stream.
        */
       template <typename STREAM>
       void print_memory_consumption(STREAM         &out,
                                     const SizeInfo &size_info) const;
 
       /**
-       * Prints a representation of the
-       * indices in the class to the given
+       * Prints a representation of the indices in the class to the given
        * output stream.
        */
       template <typename Number>
@@ -339,144 +290,113 @@ namespace internal
                   std::ostream                    &out) const;
 
       /**
-       * Stores the rowstart indices of the
-       * compressed row storage in the @p
-       * dof_indices and @p constraint_indicator
-       * fields. These two fields are always
-       * accessed together, so it is simpler to keep
-       * just one variable for them. This also
-       * obviates keeping two rowstart vectors in
-       * synch.
+       * Stores the rowstart indices of the compressed row storage in the @p
+       * dof_indices and @p constraint_indicator fields. These two fields are
+       * always accessed together, so it is simpler to keep just one variable
+       * for them. This also obviates keeping two rowstart vectors in synch.
        *
-       * In addition, the third field stores whether
-       * a particular cell has a certain structure
-       * in the indices, like indices for
-       * vector-valued problems or for cells where
-       * not all vector components are filled.
+       * In addition, the third field stores whether a particular cell has a
+       * certain structure in the indices, like indices for vector-valued
+       * problems or for cells where not all vector components are filled.
        */
       std::vector<std_cxx1x::tuple<unsigned int,
           unsigned int,
           unsigned int> > row_starts;
 
       /**
-       * Stores the indices of the degrees of
-       * freedom for each cell. This array also
-       * includes the indirect contributions from
-       * constraints, which are described by the @p
-       * constraint_indicator field. Because of
-       * variable lengths of rows, this would be a
-       * vector of a vector. However, we use one
-       * contiguous memory region and store the
+       * Stores the indices of the degrees of freedom for each cell. This
+       * array also includes the indirect contributions from constraints,
+       * which are described by the @p constraint_indicator field. Because of
+       * variable lengths of rows, this would be a vector of a
+       * vector. However, we use one contiguous memory region and store the
        * rowstart in the variable @p row_starts.
        */
       std::vector<unsigned int> dof_indices;
 
       /**
-       * This variable describes the position of
-       * constraints in terms of the local numbering
-       * of degrees of freedom on a cell. The first
-       * number stores the distance from one
-       * constrained degree of freedom to the
-       * next. This allows to identify the position
-       * of constrained DoFs as we loop through the
-       * local degrees of freedom of the cell when
-       * reading from or writing to a vector. The
-       * second number stores the index of the
-       * constraint weights, stored in the variable
-       * constraint_pool_data.
+       * This variable describes the position of constraints in terms of the
+       * local numbering of degrees of freedom on a cell. The first number
+       * stores the distance from one constrained degree of freedom to the
+       * next. This allows to identify the position of constrained DoFs as we
+       * loop through the local degrees of freedom of the cell when reading
+       * from or writing to a vector. The second number stores the index of
+       * the constraint weights, stored in the variable constraint_pool_data.
        */
       std::vector<std::pair<unsigned short,unsigned short> > constraint_indicator;
 
       /**
-       * This stores the parallel partitioning that
-       * can be used to set up vectors. The
-       * partitioner includes the description of the
-       * local range in the vector, and also
-       * includes how the ghosts look like. This
-       * enables initialization of vectors based on
-       * the DoFInfo field.
+       * This stores the parallel partitioning that can be used to set up
+       * vectors. The partitioner includes the description of the local range
+       * in the vector, and also includes how the ghosts look like. This
+       * enables initialization of vectors based on the DoFInfo field.
        */
       std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> vector_partitioner;
 
       /**
-       * This stores a (sorted) list of all locally
-       * owned degrees of freedom that are
-       * constrained.
+       * This stores a (sorted) list of all locally owned degrees of freedom
+       * that are constrained.
        */
       std::vector<unsigned int> constrained_dofs;
 
       /**
-       * Stores the rowstart indices of the
-       * compressed row storage in the @p
-       * dof_indices_plain
-       * fields.
+       * Stores the rowstart indices of the compressed row storage in the @p
+       * dof_indices_plain fields.
        */
       std::vector<unsigned int> row_starts_plain_indices;
 
       /**
-       * Stores the indices of the degrees of
-       * freedom for each cell. This array does not
-       * include the indirect contributions from
-       * constraints, which are included in @p
-       * dof_indices. Because of variable lengths of
-       * rows, this would be a vector of a
-       * vector. However, we use one contiguous
-       * memory region and store the rowstart in the
-       * variable @p row_starts_plain_indices.
+       * Stores the indices of the degrees of freedom for each cell. This
+       * array does not include the indirect contributions from constraints,
+       * which are included in @p dof_indices. Because of variable lengths of
+       * rows, this would be a vector of a vector. However, we use one
+       * contiguous memory region and store the rowstart in the variable @p
+       * row_starts_plain_indices.
        */
       std::vector<unsigned int> plain_dof_indices;
 
       /**
-       * Stores the number of components in the
-       * DoFHandler where the indices have been read
-       * from.
+       * Stores the number of components in the DoFHandler where the indices
+       * have been read from.
        */
       unsigned int n_components;
 
       /**
-       * Stores the number of degrees of freedom per
-       * cell.
+       * Stores the number of degrees of freedom per cell.
        */
       std::vector<unsigned int> dofs_per_cell;
 
       /**
-       * Stores the number of degrees of freedom per
-       * face.
+       * Stores the number of degrees of freedom per face.
        */
       std::vector<unsigned int> dofs_per_face;
 
       /**
-       * Informs on whether plain indices are
-       * cached.
+       * Informs on whether plain indices are cached.
        */
       bool store_plain_indices;
 
       /**
-       * Stores the index of the active
-       * finite element in the hp case.
+       * Stores the index of the active finite element in the hp case.
        */
       std::vector<unsigned int> cell_active_fe_index;
 
       /**
-       * Stores the maximum degree of
-       * different finite elements for the
-       * hp case.
+       * Stores the maximum degree of different finite elements for the hp
+       * case.
        */
       unsigned int max_fe_index;
 
       /**
-       * This variable stores the dofs per cell and
-       * the finite element degree associated for
-       * all fe indices in the underlying element
-       * for easier access to data in the hp case.
+       * This variable stores the dofs per cell and the finite element degree
+       * associated for all fe indices in the underlying element for easier
+       * access to data in the hp case.
        */
       std::vector<std::pair<unsigned int,unsigned int> > fe_index_conversion;
 
       /**
-       * Temporarily stores the numbers of ghosts
-       * during setup. Cleared when calling @p
-       * assign_ghosts. Then, all information is
-       * collected by the partitioner.
+       * Temporarily stores the numbers of ghosts during setup. Cleared when
+       * calling @p assign_ghosts. Then, all information is collected by the
+       * partitioner.
        */
       std::vector<unsigned int> ghost_dofs;
     };
@@ -493,7 +413,7 @@ namespace internal
       AssertIndexRange (row, row_starts.size()-1);
       const unsigned int index = std_cxx1x::get<0>(row_starts[row]);
       AssertIndexRange(index, dof_indices.size()+1);
-      return &dof_indices[index];
+      return dof_indices.empty() ? 0 : &dof_indices[0] + index;
     }
 
 
@@ -505,7 +425,7 @@ namespace internal
       AssertIndexRange (row, row_starts.size()-1);
       const unsigned int index = std_cxx1x::get<0>(row_starts[row+1]);
       AssertIndexRange(index, dof_indices.size()+1);
-      return &dof_indices[0] + index;
+      return dof_indices.empty() ? 0 : &dof_indices[0] + index;
     }
 
 
@@ -528,7 +448,7 @@ namespace internal
       AssertIndexRange (row, row_starts.size()-1);
       const unsigned int index = std_cxx1x::get<1>(row_starts[row]);
       AssertIndexRange (index, constraint_indicator.size()+1);
-      return &constraint_indicator[index];
+      return constraint_indicator.empty() ? 0 : &constraint_indicator[0] + index;
     }
 
 
@@ -540,7 +460,7 @@ namespace internal
       AssertIndexRange (row, row_starts.size()-1);
       const unsigned int index = std_cxx1x::get<1>(row_starts[row+1]);
       AssertIndexRange (index, constraint_indicator.size()+1);
-      return &constraint_indicator[0] + index;
+      return constraint_indicator.empty() ? 0 : &constraint_indicator[0] + index;
     }
 
 
@@ -560,8 +480,7 @@ namespace internal
     const unsigned int *
     DoFInfo::begin_indices_plain (const unsigned int row) const
     {
-      // if we have no constraints, should take the
-      // data from dof_indices
+      // if we have no constraints, should take the data from dof_indices
       if (row_length_indicators(row) == 0)
         {
           Assert (row_starts_plain_indices[row] == numbers::invalid_unsigned_int,
@@ -573,7 +492,7 @@ namespace internal
           AssertDimension (row_starts.size(), row_starts_plain_indices.size());
           const unsigned int index = row_starts_plain_indices[row];
           AssertIndexRange(index, plain_dof_indices.size()+1);
-          return &plain_dof_indices[index];
+          return plain_dof_indices.empty() ? 0 : &plain_dof_indices[0] + index;
         }
     }
 
index d38900c3bfc8bd85ece449df0a071005d2cdfa07..bf4003fedd86b7d77bd66d9609dab6ac8ac4663d 100644 (file)
@@ -35,10 +35,8 @@ namespace internal
     };
 
     /**
-     * A struct that takes entries describing
-     * a constraint and puts them into a
-     * sorted list where duplicates are
-     * filtered out
+     * A struct that takes entries describing a constraint and puts them into
+     * a sorted list where duplicates are filtered out
      */
     template <typename Number>
     struct ConstraintValues
@@ -46,18 +44,11 @@ namespace internal
       ConstraintValues();
 
       /**
-       * This function inserts some
-       * constrained entries to the
-       * collection of all
-       * values. It stores the
-       * (reordered) numbering of
-       * the dofs (according to the
-       * ordering that matches with
-       * the function) in
-       * new_indices, and returns
-       * the storage position the
-       * double array for access
-       * later on.
+       * This function inserts some constrained entries to the collection of
+       * all values. It stores the (reordered) numbering of the dofs
+       * (according to the ordering that matches with the function) in
+       * new_indices, and returns the storage position the double array for
+       * access later on.
        */
       unsigned short
       insert_entries (const std::vector<std::pair<unsigned int,double> > &entries);
@@ -89,27 +80,22 @@ namespace internal
                     ConstraintComparator());
           for (unsigned int j=0; j<constraint_entries.size(); j++)
             {
-              // copy the indices of the constraint entries after
-              // sorting.
+              // copy the indices of the constraint entries after sorting.
               constraint_indices[j] = constraint_entries[j].first;
 
-              // one_constraint takes the weights of the
-              // constraint
+              // one_constraint takes the weights of the constraint
               next_constraint.first[j] = constraint_entries[j].second;
             }
         }
       next_constraint.second = constraints.size();
 
-      // check whether or not constraint is already
-      // in pool. the initial implementation
-      // computed a hash value based on the truncated
-      // array (to given accuracy around 1e-13) in
-      // order to easily detect different arrays and
-      // then made a fine-grained check when the
-      // hash values were equal. this was quite
-      // lenghty and now we use a std::map with a
-      // user-defined comparator to compare floating
-      // point arrays to a tolerance 1e-13.
+      // check whether or not constraint is already in pool. the initial
+      // implementation computed a hash value based on the truncated array (to
+      // given accuracy around 1e-13) in order to easily detect different
+      // arrays and then made a fine-grained check when the hash values were
+      // equal. this was quite lenghty and now we use a std::map with a
+      // user-defined comparator to compare floating point arrays to a
+      // tolerance 1e-13.
       std::pair<typename std::map<std::vector<double>, unsigned int,
           FPArrayComparator<double> >::iterator,
           bool> it = constraints.insert(next_constraint);
@@ -119,10 +105,8 @@ namespace internal
       else
         insert_position = next_constraint.second;
 
-      // we want to store the result as a short
-      // variable, so we have to make sure that the
-      // result does not exceed the limits when
-      // casting.
+      // we want to store the result as a short variable, so we have to make
+      // sure that the result does not exceed the limits when casting.
       Assert(insert_position < (1<<(8*sizeof(unsigned short))),
              ExcInternalError());
       return static_cast<unsigned short>(insert_position);
@@ -208,19 +192,17 @@ namespace internal
           // dof is constrained
           if (entries_ptr != 0)
             {
-              // in case we want to access plain indices, we
-              // need to know about the location of
-              // constrained indices as well (all the other
-              // indices are collected by the cases below)
+              // in case we want to access plain indices, we need to know
+              // about the location of constrained indices as well (all the
+              // other indices are collected by the cases below)
               if (current_dof < first_owned || current_dof >= last_owned)
                 {
                   ghost_dofs.push_back (current_dof);
                   cell_at_boundary = true;
                 }
 
-              // check whether this dof is identity
-              // constrained to another dof. then we can
-              // simply insert that dof and there is no need
+              // check whether this dof is identity constrained to another
+              // dof. then we can simply insert that dof and there is no need
               // to actually resolve the constraint entries
               const std::vector<std::pair<unsigned int,double> >
               &entries = *entries_ptr;
@@ -239,10 +221,9 @@ namespace internal
               // reset constraint iterator for next round
               constraint_iterator.first = 0;
 
-              // add the local_to_global indices computed in
-              // the insert_entries function. transform the
-              // index to local index space or mark it as
-              // ghost if necessary
+              // add the local_to_global indices computed in the
+              // insert_entries function. transform the index to local index
+              // space or mark it as ghost if necessary
               if (n_entries > 0)
                 {
                   const std::vector<unsigned int> &constraint_indices =
@@ -256,15 +237,15 @@ namespace internal
                           dof_indices.push_back
                           (n_owned + ghost_dofs.size());
 
-                          // collect ghosts so that we can later
-                          // construct an IndexSet for them. also store
-                          // whether the current cell is on the boundary
+                          // collect ghosts so that we can later construct an
+                          // IndexSet for them. also store whether the current
+                          // cell is on the boundary
                           ghost_dofs.push_back(constraint_indices[j]);
                           cell_at_boundary = true;
                         }
                       else
-                        // not ghost, so transform to the local index
-                        // space directly
+                        // not ghost, so transform to the local index space
+                        // directly
                         dof_indices.push_back
                         (constraint_indices[j] - first_owned);
                     }
@@ -273,11 +254,9 @@ namespace internal
           else
             {
 no_constraint:
-              // Not constrained, we simply have to add the
-              // local index to the indices_local_to_global
-              // list and increment constraint
-              // iterator. transform to local index
-              // space/mark as ghost
+              // Not constrained, we simply have to add the local index to the
+              // indices_local_to_global list and increment constraint
+              // iterator. transform to local index space/mark as ghost
               if (n_mpi_procs > 1 &&
                   (current_dof < first_owned ||
                    current_dof >= last_owned))
@@ -291,8 +270,8 @@ no_constraint:
 
               dof_indices.push_back (current_dof);
 
-              // make sure constraint_iterator.first is
-              // always within the bounds of unsigned short
+              // make sure constraint_iterator.first is always within the
+              // bounds of unsigned short
               Assert (constraint_iterator.first <
                       (1<<(8*sizeof(unsigned short)))-1,
                       ExcInternalError());
@@ -303,9 +282,8 @@ no_constraint:
                                   unsigned int>
                                   (dof_indices.size(), constraint_indicator.size(), 0);
 
-      // now to the plain indices: in case we have
-      // constraints on this cell, store the indices
-      // without the constraints resolve once again
+      // now to the plain indices: in case we have constraints on this cell,
+      // store the indices without the constraints resolve once again
       if (store_plain_indices == true)
         {
           if (cell_number == 0)
@@ -358,11 +336,9 @@ no_constraint:
       IndexSet ghost_indices (vector_partitioner->size());
       if (n_ghosts > 0)
         {
-          // since we need to go back to the
-          // local_to_global indices and replace the
-          // temporary numbering of ghosts by the real
-          // number in the index set, we need to store
-          // these values
+          // since we need to go back to the local_to_global indices and
+          // replace the temporary numbering of ghosts by the real number in
+          // the index set, we need to store these values
           std::vector<std::pair<unsigned int,unsigned int> > ghost_origin(n_ghosts);
           for (unsigned int i=0; i<n_ghosts; ++i)
             {
@@ -389,9 +365,8 @@ no_constraint:
                                    ghost_origin.back().first+1);
           ghost_indices.compress();
 
-          // make sure that we got the correct local
-          // numbering of the ghost dofs. the ghost
-          // index set should store the same number
+          // make sure that we got the correct local numbering of the ghost
+          // dofs. the ghost index set should store the same number
           {
             AssertDimension (n_unique_ghosts, ghost_indices.n_elements());
             for (unsigned int i=0; i<n_ghosts; ++i)
@@ -400,12 +375,10 @@ no_constraint:
                       ExcInternalError());
           }
 
-          // apply correct numbering for ghost indices:
-          // We previously just enumerated them
-          // according to their appearance in the
-          // local_to_global structure. Above, we
-          // derived a relation between this enumeration
-          // and the actual number
+          // apply correct numbering for ghost indices: We previously just
+          // enumerated them according to their appearance in the
+          // local_to_global structure. Above, we derived a relation between
+          // this enumeration and the actual number
           const unsigned int n_boundary_cells = boundary_cells.size();
           for (unsigned int i=0; i<n_boundary_cells; ++i)
             {
@@ -444,10 +417,9 @@ no_constraint:
       std::vector<unsigned int> new_ghosts;
       ghost_dofs.swap(new_ghosts);
 
-      // set the ghost indices now. need to cast
-      // away constness here, but that is uncritical
-      // since we reset the Partitioner in the same
-      // initialize call as this call here.
+      // set the ghost indices now. need to cast away constness here, but that
+      // is uncritical since we reset the Partitioner in the same initialize
+      // call as this call here.
       Utilities::MPI::Partitioner *vec_part =
         const_cast<Utilities::MPI::Partitioner *>(vector_partitioner.get());
       vec_part->set_ghost_indices (ghost_indices);
@@ -617,8 +589,7 @@ no_constraint:
               for (unsigned int j=0; j<n_comp; ++j)
                 fe_indices[j]=cell_active_fe_index[renumbering[position_cell+j]];
 
-              // by construction, all cells should
-              // have the same fe index.
+              // by construction, all cells should have the same fe index.
               for (unsigned int j=1; j<n_comp; ++j)
                 Assert (fe_indices[j] == fe_indices[0], ExcInternalError());
 
@@ -645,13 +616,11 @@ no_constraint:
           new_plain_indices.reserve (plain_dof_indices.size());
         }
 
-      // copy the indices and the constraint
-      // indicators to the new data field: Store the
-      // indices in a way so that adjacent data
-      // fields in local vectors are adjacent, i.e.,
-      // first dof index 0 for all vectors, then dof
-      // index 1 for all vectors, and so on. This
-      // involves some extra resorting.
+      // copy the indices and the constraint indicators to the new data field:
+      // Store the indices in a way so that adjacent data fields in local
+      // vectors are adjacent, i.e., first dof index 0 for all vectors, then
+      // dof index 1 for all vectors, and so on. This involves some extra
+      // resorting.
       std::vector<const unsigned int *> glob_indices (vectorization_length);
       std::vector<const unsigned int *> plain_glob_indices (vectorization_length);
       std::vector<const std::pair<unsigned short,unsigned short>*>
@@ -707,14 +676,12 @@ no_constraint:
                     continue;
                   }
 
-                // otherwise, check if we are a constrained
-                // dof. The dof is not constrained if we are
-                // at the end of the row for the constraints
-                // (indi[j] == n_indi[j]) or if the local
-                // index[j] is smaller than the next position
-                // for a constraint. Then, just copy
-                // it. otherwise, copy all the entries that
-                // come with this dof
+                // otherwise, check if we are a constrained dof. The dof is
+                // not constrained if we are at the end of the row for the
+                // constraints (indi[j] == n_indi[j]) or if the local index[j]
+                // is smaller than the next position for a constraint. Then,
+                // just copy it. otherwise, copy all the entries that come
+                // with this dof
                 if (constr_ind[j] == constr_end[j] ||
                     index[j] < constr_ind[j]->first)
                   {
@@ -763,20 +730,17 @@ no_constraint:
       new_rowstart_plain.swap (row_starts_plain_indices);
 
 #ifdef DEBUG
-      // sanity check 1: all indices should be
-      // smaller than the number of dofs locally
-      // owned plus the number of ghosts
+      // sanity check 1: all indices should be smaller than the number of dofs
+      // locally owned plus the number of ghosts
       const unsigned int index_range = (vector_partitioner->local_range().second-
                                         vector_partitioner->local_range().first)
                                        + vector_partitioner->ghost_indices().n_elements();
       for (unsigned int i=0; i<dof_indices.size(); ++i)
         AssertIndexRange (dof_indices[i], index_range);
 
-      // sanity check 2: for the constraint
-      // indicators, the first index should be
-      // smaller than the number of indices in the
-      // row, and the second index should be smaller
-      // than the number of constraints in the
+      // sanity check 2: for the constraint indicators, the first index should
+      // be smaller than the number of indices in the row, and the second
+      // index should be smaller than the number of constraints in the
       // constraint pool.
       for (unsigned int row=0; row<size_info.n_macro_cells; ++row)
         {
@@ -791,9 +755,8 @@ no_constraint:
             }
         }
 
-      // sanity check 3: all non-boundary cells
-      // should have indices that only refer to the
-      // locally owned range
+      // sanity check 3: all non-boundary cells should have indices that only
+      // refer to the locally owned range
       const unsigned int local_size = (vector_partitioner->local_range().second-
                                        vector_partitioner->local_range().first);
       for (unsigned int row=0; row<size_info.boundary_cells_start; ++row)
@@ -819,19 +782,16 @@ no_constraint:
     void DoFInfo::guess_block_size (const SizeInfo &size_info,
                                     TaskInfo       &task_info)
     {
-      // user did not say a positive number, so we
-      // have to guess
+      // user did not say a positive number, so we have to guess
       if (task_info.block_size == 0)
         {
-          // we would like to have enough work to do, so
-          // as first guess, try to get 50 times as many
-          // chunks as we have threads on the
-          // system.
+          // we would like to have enough work to do, so as first guess, try
+          // to get 50 times as many chunks as we have threads on the system.
           task_info.block_size =
             size_info.n_macro_cells / (multithread_info.n_default_threads * 50);
 
-          // if there are too few degrees of freedom per
-          // cell, need to increase the block size
+          // if there are too few degrees of freedom per cell, need to
+          // increase the block size
           const unsigned int minimum_parallel_grain_size = 500;
           if (dofs_per_cell[0] * task_info.block_size <
               minimum_parallel_grain_size)
@@ -859,11 +819,9 @@ no_constraint:
 
       guess_block_size (size_info, task_info);
 
-      // set up partitions. if we just use coloring
-      // without partitions, do nothing here, assume
-      // all cells to belong to the zero partition
-      // (that we otherwise use for MPI boundary
-      // cells)
+      // set up partitions. if we just use coloring without partitions, do
+      // nothing here, assume all cells to belong to the zero partition (that
+      // we otherwise use for MPI boundary cells)
       unsigned int partition = 0, start_up = 0, counter = 0;
       unsigned int start_nonboundary = numbers::invalid_unsigned_int;
       bool work = true;
@@ -941,51 +899,45 @@ no_constraint:
             task_info.block_size_last = task_info.block_size;
         }
 
-      // assume that all FEs have the same
-      // connectivity graph, so take the zeroth FE
+      // assume that all FEs have the same connectivity graph, so take the
+      // zeroth FE
       task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)/
                            task_info.block_size;
       task_info.block_size_last = size_info.n_macro_cells-
                                   (task_info.block_size*(task_info.n_blocks-1));
 
-      // create the connectivity graph with internal
-      // blocking
+      // create the connectivity graph with internal blocking
       CompressedSimpleSparsityPattern connectivity;
       make_connectivity_graph (size_info, task_info, renumbering,irregular_cells,
                                true, connectivity);
 
       // Create cell-block  partitioning.
 
-      // For each block of cells, this variable
-      // saves to which partitions the block
-      // belongs. Initialize all to n_macro_cells to
-      // mark them as not yet assigned a partition.
+      // For each block of cells, this variable saves to which partitions the
+      // block belongs. Initialize all to n_macro_cells to mark them as not
+      // yet assigned a partition.
       std::vector<unsigned int> cell_partition(task_info.n_blocks,
                                                size_info.n_macro_cells);
       std::vector<unsigned int> neighbor_list;
       std::vector<unsigned int> neighbor_neighbor_list;
 
-      // In element j of this variable, one puts the
-      // old number of the block that should be the
-      // jth block in the new numeration.
+      // In element j of this variable, one puts the old number of the block
+      // that should be the jth block in the new numeration.
       std::vector<unsigned int> partition_list      (task_info.n_blocks,0);
       std::vector<unsigned int> partition_color_list(task_info.n_blocks,0);
 
-      // This vector points to the start of each
-      // partition.
+      // This vector points to the start of each partition.
       std::vector<unsigned int> partition_blocks (2,0);
       std::vector<unsigned int> cell_color(task_info.n_blocks,
                                            size_info.n_macro_cells);
       std::vector<bool> color_finder;
 
-      // this performs a classical breath-first
-      // search in the connectivity graph of the
-      // cell chunks
+      // this performs a classical breath-first search in the connectivity
+      // graph of the cell chunks
       while (work)
         {
-          // put all cells up to begin_inner_cells into
-          // first partition. if the numbers do not add
-          // up exactly, assign an additional block
+          // put all cells up to begin_inner_cells into first partition. if
+          // the numbers do not add up exactly, assign an additional block
           if (start_nonboundary>0 && start_up == start_nonboundary)
             {
               unsigned int n_blocks = ((start_nonboundary+task_info.block_size-1)
@@ -1000,8 +952,8 @@ no_constraint:
             }
           else
             {
-              // To start up, set the start_up cell to
-              // partition and list all its neighbors.
+              // To start up, set the start_up cell to partition and list all
+              // its neighbors.
               AssertIndexRange(start_up, cell_partition.size());
               cell_partition[start_up] = partition;
               neighbor_list.push_back(start_up);
@@ -1035,9 +987,8 @@ no_constraint:
               neighbor_neighbor_list.resize(0);
             }
 
-          // One has to check if the graph is not
-          // connected so we have to find another
-          // partition.
+          // One has to check if the graph is not connected so we have to find
+          // another partition.
           work = false;
           for (unsigned int j=start_up; j<task_info.n_blocks; ++j)
             if (cell_partition[j] == size_info.n_macro_cells)
@@ -1063,9 +1014,8 @@ no_constraint:
               unsigned int cell = partition_list[k];
               unsigned int no_neighbors = connectivity.row_length(cell);
 
-              // In the worst case, each neighbor has a
-              // different color. So we find at least one
-              // available color between 0 and no_neighbors.
+              // In the worst case, each neighbor has a different color. So we
+              // find at least one available color between 0 and no_neighbors.
               color_finder.resize(no_neighbors+1);
               for (unsigned int j=0; j<=no_neighbors; ++j)
                 color_finder[j]=true;
@@ -1074,24 +1024,22 @@ no_constraint:
               end      = connectivity.row_end(cell);
               for (; neighbor!=end ; ++neighbor)
                 {
-                  // Mark the color that a neighbor within the
-                  // partition has as taken
+                  // Mark the color that a neighbor within the partition has
+                  // as taken
                   if (cell_partition[*neighbor] == part &&
                       cell_color[*neighbor] <= no_neighbors)
                     color_finder[cell_color[*neighbor]] = false;
                 }
-              // Choose the smallest color that is not taken
-              // for the block
+              // Choose the smallest color that is not taken for the block
               cell_color[cell]=0;
               while (color_finder[cell_color[cell]] == false)
                 cell_color[cell]++;
               if (cell_color[cell] > max_color)
                 max_color = cell_color[cell];
             }
-          // Reorder within partition: First, all blocks
-          // that belong the 0 and then so on until
-          // those with color max (Note that the smaller
-          // the number the larger the partition)
+          // Reorder within partition: First, all blocks that belong the 0 and
+          // then so on until those with color max (Note that the smaller the
+          // number the larger the partition)
           for (unsigned int color=0; color<=max_color; color++)
             {
               task_info.partition_color_blocks_data.push_back(color_counter);
@@ -1113,8 +1061,7 @@ no_constraint:
 
       partition_list = renumbering;
 
-      // in debug mode, check that the partition
-      // color list is one-to-one
+      // in debug mode, check that the partition color list is one-to-one
 #ifdef DEBUG
       {
         std::vector<unsigned int> sorted_pc_list (partition_color_list);
@@ -1124,8 +1071,8 @@ no_constraint:
       }
 #endif
 
-      // set the start list for each block and
-      // compute the renumbering of cells
+      // set the start list for each block and compute the renumbering of
+      // cells
       std::vector<unsigned int> block_start(size_info.n_macro_cells+1);
       std::vector<unsigned int> irregular(size_info.n_macro_cells);
 
@@ -1202,8 +1149,8 @@ no_constraint:
 
       guess_block_size (size_info, task_info);
 
-      // assume that all FEs have the same
-      // connectivity graph, so take the zeroth FE
+      // assume that all FEs have the same connectivity graph, so take the
+      // zeroth FE
       task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)/
                            task_info.block_size;
       task_info.block_size_last = size_info.n_macro_cells-
@@ -1211,31 +1158,27 @@ no_constraint:
       task_info.position_short_block = task_info.n_blocks-1;
       unsigned int cluster_size = task_info.block_size*vectorization_length;
 
-      // create the connectivity graph without
-      // internal blocking
+      // create the connectivity graph without internal blocking
       CompressedSimpleSparsityPattern connectivity;
       make_connectivity_graph (size_info, task_info, renumbering,irregular_cells,
                                false, connectivity);
 
       // Create cell-block  partitioning.
 
-      // For each block of cells, this variable
-      // saves to which partitions the block
-      // belongs. Initialize all to n_macro_cells to
-      // mark them as not yet assigned a partition.
+      // For each block of cells, this variable saves to which partitions the
+      // block belongs. Initialize all to n_macro_cells to mark them as not
+      // yet assigned a partition.
       std::vector<unsigned int> cell_partition (size_info.n_active_cells,
                                                 size_info.n_active_cells);
       std::vector<unsigned int> neighbor_list;
       std::vector<unsigned int> neighbor_neighbor_list;
 
-      // In element j of this variable, one puts the
-      // old number of the block that should be the
-      // jth block in the new numeration.
+      // In element j of this variable, one puts the old number of the block
+      // that should be the jth block in the new numeration.
       std::vector<unsigned int> partition_list(size_info.n_active_cells,0);
       std::vector<unsigned int> partition_partition_list(size_info.n_active_cells,0);
 
-      // This vector points to the start of each
-      // partition.
+      // This vector points to the start of each partition.
       std::vector<unsigned int> partition_size(2,0);
 
       unsigned int partition = 0,start_up=0,counter=0;
@@ -1245,15 +1188,12 @@ no_constraint:
       bool work = true;
       unsigned int remainder = cluster_size;
 
-      // this performs a classical breath-first
-      // search in the connectivity graph of the
-      // cells under the restriction that the size
-      // of the partitions should be a multiple of
-      // the given block size
+      // this performs a classical breath-first search in the connectivity
+      // graph of the cells under the restriction that the size of the
+      // partitions should be a multiple of the given block size
       while (work)
         {
-          // put the cells with neighbors on remote MPI
-          // processes up front
+          // put the cells with neighbors on remote MPI processes up front
           if (start_nonboundary>0)
             {
               for (unsigned int cell=0; cell<start_nonboundary; ++cell)
@@ -1268,14 +1208,13 @@ no_constraint:
               if (remainder == cluster_size)
                 remainder = 0;
 
-              // adjust end of boundary cells to the
-              // remainder
+              // adjust end of boundary cells to the remainder
               size_info.boundary_cells_end += (remainder+vectorization_length-1)/vectorization_length;
             }
           else
             {
-              // To start up, set the start_up cell to
-              // partition and list all its neighbors.
+              // To start up, set the start_up cell to partition and list all
+              // its neighbors.
               cell_partition[start_up] = partition;
               neighbor_list.push_back(start_up);
               partition_list[counter++] = start_up;
@@ -1387,9 +1326,8 @@ no_constraint:
               neighbor_neighbor_list.resize(0);
             }
 not_connect:
-          // One has to check if the graph is not
-          // connected so we have to find another
-          // partition.
+          // One has to check if the graph is not connected so we have to find
+          // another partition.
           work = false;
           for (unsigned int j=start_up; j<size_info.n_active_cells; ++j)
             if (cell_partition[j] == size_info.n_active_cells)
@@ -1412,10 +1350,9 @@ not_connect:
       {
         // Create partitioning within partitions.
 
-        // For each block of cells, this variable
-        // saves to which partitions the block
-        // belongs. Initialize all to n_macro_cells to
-        // mark them as not yet assigned a partition.
+        // For each block of cells, this variable saves to which partitions
+        // the block belongs. Initialize all to n_macro_cells to mark them as
+        // not yet assigned a partition.
         std::vector<unsigned int> cell_partition_l2(size_info.n_active_cells,
                                                     size_info.n_active_cells);
         task_info.partition_color_blocks_row_index.resize(partition+1,0);
@@ -1445,8 +1382,8 @@ not_connect:
                           start_up = j;
                           work = true;
                           partition_counter = 1;
-                          // To start up, set the start_up cell to
-                          // partition and list all its neighbors.
+                          // To start up, set the start_up cell to partition
+                          // and list all its neighbors.
                           AssertIndexRange (start_up, partition_size[part+1]);
                           cell_partition_l2[partition_list[start_up]] =
                             partition_l2;
@@ -1489,8 +1426,8 @@ not_connect:
                     int index_before = neighbor_neighbor_list.size(),
                         index = index_before;
                     {
-                      // put the cells into separate lists for each
-                      // FE index within one partition-partition
+                      // put the cells into separate lists for each FE index
+                      // within one partition-partition
                       missing_macros = 0;
                       std::vector<unsigned int> remaining_per_macro_cell
                       (max_fe_index);
@@ -1509,8 +1446,7 @@ not_connect:
                                                     [cell]]].
                               push_back(partition_partition_list[cell]);
                             }
-                          // check how many more cells are needed in the
-                          // lists
+                          // check how many more cells are needed in the lists
                           for (unsigned int j=0; j<max_fe_index; j++)
                             {
                               remaining_per_macro_cell[j] =
@@ -1536,8 +1472,7 @@ not_connect:
                       missing_macros = task_info.block_size -
                                        (missing_macros%task_info.block_size);
 
-                      // now we realized that there are some cells
-                      // missing.
+                      // now we realized that there are some cells missing.
                       while (missing_macros>0 || filled == false)
                         {
                           if (index==0)
@@ -1557,9 +1492,9 @@ not_connect:
                           unsigned int additional = neighbor_neighbor_list
                                                     [index];
 
-                          // go through the neighbors of the last cell
-                          // in the current partition and check if we
-                          // find some to fill up with.
+                          // go through the neighbors of the last cell in the
+                          // current partition and check if we find some to
+                          // fill up with.
                           CompressedSimpleSparsityPattern::row_iterator
                           neighbor =
                             connectivity.row_begin(additional),
@@ -1576,9 +1511,9 @@ not_connect:
                                                  cell_active_fe_index[*neighbor];
 
                                   // Only add this cell if we need more macro
-                                  // cells in the current block or if there is a
-                                  // macro cell with the FE index that is not
-                                  // yet fully populated
+                                  // cells in the current block or if there is
+                                  // a macro cell with the FE index that is
+                                  // not yet fully populated
                                   if (missing_macros > 0 ||
                                       remaining_per_macro_cell[this_index] > 0)
                                     {
@@ -1616,10 +1551,9 @@ not_connect:
                         }
                       if (hp_bool == true)
                         {
-                          // set the renumbering according to their
-                          // active FE index within one
-                          // partition-partition which was implicitly
-                          // assumed above
+                          // set the renumbering according to their active FE
+                          // index within one partition-partition which was
+                          // implicitly assumed above
                           cell = counter - partition_counter;
                           for (unsigned int j=0; j<max_fe_index; j++)
                             {
@@ -1698,14 +1632,9 @@ not_connect:
     }
 
 
-    // internal sparsity pattern structure that is
-    // specialized to the situation in the setup
-    // of the cells the dofs belong to: allows
-    // only insertion of entries at the end
     namespace internal
     {
-      // rudimentary version of a vector that keeps
-      // entries always ordered
+      // rudimentary version of a vector that keeps entries always ordered
       class ordered_vector : public std::vector<unsigned int>
       {
       public:
@@ -1721,10 +1650,8 @@ not_connect:
         }
 
 
-        // insert a given entry. dat is a pointer
-        // within this vector (make sure that it
-        // really stays there, otherwise there will be
-        // an error)
+        // insert a given entry. dat is a pointer within this vector (the user
+        // needs to make sure that it really stays there)
         void insert (const unsigned int entry,
                      std::vector<unsigned int>::iterator &dat)
         {
@@ -1737,12 +1664,15 @@ not_connect:
           if (dat == end())
             {
               push_back(entry);
+              dat = end();
             }
           else if (*dat > entry)
             {
-              this->std::vector<unsigned int>::insert (dat, entry);
+              dat = this->std::vector<unsigned int>::insert (dat, entry);
+              ++dat;
             }
-          ++dat;
+          else
+            ++dat;
         }
       };
     }
@@ -1926,8 +1856,7 @@ not_connect:
 
     void DoFInfo::renumber_dofs (std::vector<unsigned int> &renumbering)
     {
-      // first renumber all locally owned degrees of
-      // freedom
+      // first renumber all locally owned degrees of freedom
       AssertDimension (vector_partitioner->local_size(),
                        vector_partitioner->size());
       const unsigned int local_size = vector_partitioner->local_size();
@@ -1956,9 +1885,8 @@ not_connect:
       for (unsigned int i=0; i<constrained_dofs.size(); ++i)
         new_constrained_dofs[i] = renumbering[constrained_dofs[i]];
 
-      // the new constrained DoFs should be sorted
-      // already as they are not contained in
-      // dof_indices and then get contiguous numbers
+      // the new constrained DoFs should be sorted already as they are not
+      // contained in dof_indices and then get contiguous numbers
 #ifdef DEBUG
       for (unsigned int i=1; i<new_constrained_dofs.size(); ++i)
         Assert (new_constrained_dofs[i] > new_constrained_dofs[i-1], ExcInternalError());
index 6a46b630fffd3a3930c7873472397fa2c7c21b0c..481fe26814a1c480f59f83757a7634e5febc4909 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2011 by the deal.II authors
+//    Copyright (C) 2011-2012 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -42,17 +42,14 @@ namespace internal
     struct MappingInfo
     {
       /**
-       * Determines how many bits of an unsigned int
-       * are used to distinguish the cell types
-       * (Cartesian, with constant Jacobian, or
-       * general)
+       * Determines how many bits of an unsigned int are used to distinguish
+       * the cell types (Cartesian, with constant Jacobian, or general)
        */
       static const std::size_t  n_cell_type_bits = 2;
 
       /**
-       * Determines how many types of different
-       * cells can be detected at most. Corresponds
-       * to the number of bits we reserved for it.
+       * Determines how many types of different cells can be detected at
+       * most. Corresponds to the number of bits we reserved for it.
        */
       static const unsigned int n_cell_types = 1U<<n_cell_type_bits;
 
@@ -62,16 +59,12 @@ namespace internal
       MappingInfo();
 
       /**
-       * Computes the information in the given
-       * cells. The cells are specified by the level
-       * and the index within the level (as given by
-       * CellIterator::level() and
-       * CellIterator::index(), in order to allow
-       * for different kinds of iterators,
-       * e.g. standard DoFHandler, multigrid, etc.)
-       * on a fixed Triangulation. In addition, a
-       * mapping and several quadrature formulas are
-       * given.
+       * Computes the information in the given cells. The cells are specified
+       * by the level and the index within the level (as given by
+       * CellIterator::level() and CellIterator::index(), in order to allow
+       * for different kinds of iterators, e.g. standard DoFHandler,
+       * multigrid, etc.)  on a fixed Triangulation. In addition, a mapping
+       * and several quadrature formulas are given.
        */
       void initialize (const dealii::Triangulation<dim>                &tria,
                        const std::vector<std::pair<unsigned int,unsigned int> > &cells,
@@ -81,24 +74,20 @@ namespace internal
                        const UpdateFlags                        update_flags);
 
       /**
-       * Helper function to determine which update
-       * flags must be set in the internal functions
-       * to initialize all data as requested by the
-       * user.
+       * Helper function to determine which update flags must be set in the
+       * internal functions to initialize all data as requested by the user.
        */
       UpdateFlags
       compute_update_flags (const UpdateFlags                        update_flags,
                             const std::vector<dealii::hp::QCollection<1> > &quad) const;
 
       /**
-       * Returns the type of a given cell as
-       * detected during initialization.
+       * Returns the type of a given cell as detected during initialization.
        */
       CellType get_cell_type (const unsigned int cell_chunk_no) const;
 
       /**
-       * Returns the type of a given cell as
-       * detected during initialization.
+       * Returns the type of a given cell as detected during initialization.
        */
       unsigned int get_cell_data_index (const unsigned int cell_chunk_no) const;
 
@@ -108,215 +97,171 @@ namespace internal
       void clear ();
 
       /**
-       * Returns the memory consumption of this
-       * class in bytes.
+       * Returns the memory consumption of this class in bytes.
        */
       std::size_t memory_consumption() const;
 
       /**
-       * Prints a detailed summary of memory
-       * consumption in the different structures of
-       * this class to the given output stream.
+       * Prints a detailed summary of memory consumption in the different
+       * structures of this class to the given output stream.
        */
       template <typename STREAM>
       void print_memory_consumption(STREAM         &out,
                                     const SizeInfo &size_info) const;
 
       /**
-       * Stores whether a cell is Cartesian, has
-       * constant transform data (Jacobians) or is
-       * general. cell_type % 4 gives this
-       * information (0: Cartesian, 1: constant
-       * Jacobian throughout cell, 2: general cell),
-       * and cell_type / 4 gives the index in the
-       * data field of where to find the information
-       * in the fields Jacobian and JxW values
-       * (except for quadrature points, for which
-       * the index runs as usual).
+       * Stores whether a cell is Cartesian, has constant transform data
+       * (Jacobians) or is general. cell_type % 4 gives this information (0:
+       * Cartesian, 1: constant Jacobian throughout cell, 2: general cell),
+       * and cell_type / 4 gives the index in the data field of where to find
+       * the information in the fields Jacobian and JxW values (except for
+       * quadrature points, for which the index runs as usual).
        */
       std::vector<unsigned int> cell_type;
 
       /**
-       * The first field stores the inverse Jacobian
-       * for Cartesian cells: There, it is a
-       * diagonal rank-2 tensor, so we actually just
-       * store a rank-1 tensor. It is the same on
-       * all cells, therefore we only store it once
-       * per cell, and use similarities from one
-       * cell to another, too (on structured meshes,
-       * there are usually many cells with the same
+       * The first field stores the inverse Jacobian for Cartesian cells:
+       * There, it is a diagonal rank-2 tensor, so we actually just store a
+       * rank-1 tensor. It is the same on all cells, therefore we only store
+       * it once per cell, and use similarities from one cell to another, too
+       * (on structured meshes, there are usually many cells with the same
        * Jacobian).
        *
-       * The second field stores the Jacobian
-       * determinant for Cartesian cells (without
-       * the quadrature weight, which depends on the
-       * quadrature point, whereas the determinant
-       * is the same on each quadrature point).
+       * The second field stores the Jacobian determinant for Cartesian cells
+       * (without the quadrature weight, which depends on the quadrature
+       * point, whereas the determinant is the same on each quadrature point).
        */
       AlignedVector<std::pair<Tensor<1,dim,VectorizedArray<Number> >,
                     VectorizedArray<Number> > > cartesian_data;
 
       /**
-       * The first field stores the Jacobian for
-       * non-Cartesian cells where all the Jacobians
-       * on the cell are the same (i.e., constant,
-       * which comes from a linear transformation
-       * from unit to real cell). Also use
-       * similarities from one cell to another (on
-       * structured meshes, there are usually many
-       * cells with the same Jacobian).
+       * The first field stores the Jacobian for non-Cartesian cells where all
+       * the Jacobians on the cell are the same (i.e., constant, which comes
+       * from a linear transformation from unit to real cell). Also use
+       * similarities from one cell to another (on structured meshes, there
+       * are usually many cells with the same Jacobian).
        *
-       * The second field stores the Jacobian
-       * determinant for non-Cartesian cells with
-       * constant Jacobian throughout the cell
-       * (without the quadrature weight, which
-       * depends on the quadrature point, whereas
-       * the determinant is the same on each
-       * quadrature point).
+       * The second field stores the Jacobian determinant for non-Cartesian
+       * cells with constant Jacobian throughout the cell (without the
+       * quadrature weight, which depends on the quadrature point, whereas the
+       * determinant is the same on each quadrature point).
        */
       AlignedVector<std::pair<Tensor<2,dim,VectorizedArray<Number> >,
                     VectorizedArray<Number> > > affine_data;
 
       /**
-       * Definition of a structure that stores data
-       * that depends on the quadrature formula (if
-       * we have more than one quadrature formula on
-       * a given problem, these fields will be
-       * different)
+       * Definition of a structure that stores data that depends on the
+       * quadrature formula (if we have more than one quadrature formula on a
+       * given problem, these fields will be different)
        */
       struct MappingInfoDependent
       {
         /**
-         * This field stores the row starts for the
-         * inverse Jacobian transformations,
-         * quadrature weights and second derivatives.
+         * This field stores the row starts for the inverse Jacobian
+         * transformations, quadrature weights and second derivatives.
          */
         std::vector<unsigned int> rowstart_jacobians;
 
         /**
-         * This field stores the inverse Jacobian
-         * transformation from unit to real cell,
-         * which is needed for most gradient
-         * transformations (corresponds to
-         * FEValues::inverse_jacobian) for general
-         * cells.
+         * This field stores the inverse Jacobian transformation from unit to
+         * real cell, which is needed for most gradient transformations
+         * (corresponds to FEValues::inverse_jacobian) for general cells.
          */
         AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians;
 
         /**
-         * This field stores the Jacobian
-         * determinant times the quadrature weights
-         * (JxW in deal.II speak) for general cells.
+         * This field stores the Jacobian determinant times the quadrature
+         * weights (JxW in deal.II speak) for general cells.
          */
         AlignedVector<VectorizedArray<Number> > JxW_values;
 
         /**
-         * Stores the diagonal part of the gradient of
-         * the inverse Jacobian transformation. The
-         * first index runs over the derivatives
-         * $\partial^2/\partial x_i^2$, the second
-         * over the space coordinate. Needed for
-         * computing the Laplacian of FE functions on
-         * the real cell. Uses a separate storage from
-         * the off-diagonal part $\partial^2/\partial
-         * x_i \partial x_j, i\neq j$ because that is
+         * Stores the diagonal part of the gradient of the inverse Jacobian
+         * transformation. The first index runs over the derivatives
+         * $\partial^2/\partial x_i^2$, the second over the space
+         * coordinate. Needed for computing the Laplacian of FE functions on
+         * the real cell. Uses a separate storage from the off-diagonal part
+         * $\partial^2/\partial x_i \partial x_j, i\neq j$ because that is
          * only needed for computing a full Hessian.
          */
         AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians_grad_diag;
 
         /**
-         * Stores the off-diagonal part of the
-         * gradient of the inverse Jacobian
-         * transformation. Because of symmetry, only
-         * the upper diagonal part is needed. The
-         * first index runs through the derivatives
-         * row-wise, i.e., $\partial^2/\partial x_1
-         * \partial x_2$ first, then
-         * $\partial^2/\partial x_1 \partial x_3$, and
-         * so on. The second index is the spatial
-         * coordinate. Not filled currently.
+         * Stores the off-diagonal part of the gradient of the inverse
+         * Jacobian transformation. Because of symmetry, only the upper
+         * diagonal part is needed. The first index runs through the
+         * derivatives row-wise, i.e., $\partial^2/\partial x_1 \partial x_2$
+         * first, then $\partial^2/\partial x_1 \partial x_3$, and so on. The
+         * second index is the spatial coordinate. Not filled currently.
          */
         AlignedVector<Tensor<1,(dim>1?dim*(dim-1)/2:1),
                       Tensor<1,dim,VectorizedArray<Number> > > > jacobians_grad_upper;
 
         /**
-         * Stores the row start for quadrature points
-         * in real coordinates for both types of
-         * cells. Note that Cartesian cells will have
-         * shorter fields (length is @p n_q_points_1d)
-         * than non-Cartesian cells (length is @p
-         * n_q_points).
+         * Stores the row start for quadrature points in real coordinates for
+         * both types of cells. Note that Cartesian cells will have shorter
+         * fields (length is @p n_q_points_1d) than non-Cartesian cells
+         * (length is @p n_q_points).
          */
         std::vector<unsigned int> rowstart_q_points;
 
         /**
-         * Stores the quadrature points in real
-         * coordinates for Cartesian cells (does not
-         * need to store the full data on all points)
+         * Stores the quadrature points in real coordinates for Cartesian
+         * cells (does not need to store the full data on all points)
          */
         AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
 
         /**
-         * The dim-dimensional quadrature formula
-         * underlying the problem (constructed from a
-         * 1D tensor product quadrature formula).
+         * The dim-dimensional quadrature formula underlying the problem
+         * (constructed from a 1D tensor product quadrature formula).
          */
         dealii::hp::QCollection<dim>    quadrature;
 
         /**
-         * The (dim-1)-dimensional quadrature formula
-         * corresponding to face evaluation
-         * (constructed from a 1D tensor product
-         * quadrature formula).
+         * The (dim-1)-dimensional quadrature formula corresponding to face
+         * evaluation (constructed from a 1D tensor product quadrature
+         * formula).
          */
         dealii::hp::QCollection<dim-1>  face_quadrature;
 
         /**
-         * The number of quadrature points for the
-         * current quadrature formula.
+         * The number of quadrature points for the current quadrature formula.
          */
         std::vector<unsigned int> n_q_points;
 
         /**
-         * The number of quadrature points for the
-         * current quadrature formula when applied to
-         * a face. Only set if the quadrature formula
-         * is derived from a tensor product, since it
-         * is not defined from the full quadrature
-         * formula otherwise.
+         * The number of quadrature points for the current quadrature formula
+         * when applied to a face. Only set if the quadrature formula is
+         * derived from a tensor product, since it is not defined from the
+         * full quadrature formula otherwise.
          */
         std::vector<unsigned int> n_q_points_face;
 
         /**
-         * The quadrature weights (vectorized data
-         * format) on the unit cell.
+         * The quadrature weights (vectorized data format) on the unit cell.
          */
         std::vector<AlignedVector<VectorizedArray<Number> > > quadrature_weights;
 
         /**
-         * This variable stores the number of
-         * quadrature points for all quadrature
-         * indices in the underlying element for
-         * easier access to data in the hp case.
+         * This variable stores the number of quadrature points for all
+         * quadrature indices in the underlying element for easier access to
+         * data in the hp case.
          */
         std::vector<unsigned int> quad_index_conversion;
 
         /**
-         * Returns the quadrature index for a given
-         * number of quadrature points. If not in hp
-         * mode or if the index is not found, this
-         * function always returns index 0. Hence,
-         * this function does not check whether the
-         * given degree is actually present.
+         * Returns the quadrature index for a given number of quadrature
+         * points. If not in hp mode or if the index is not found, this
+         * function always returns index 0. Hence, this function does not
+         * check whether the given degree is actually present.
          */
         unsigned int
         quad_index_from_n_q_points (const unsigned int n_q_points) const;
 
 
         /**
-         * Prints a detailed summary of memory
-         * consumption in the different structures of
-         * this class to the given output stream.
+         * Prints a detailed summary of memory consumption in the different
+         * structures of this class to the given output stream.
          */
         template <typename STREAM>
         void print_memory_consumption(STREAM         &out,
@@ -329,33 +274,28 @@ namespace internal
       };
 
       /**
-       * Contains all the stuff that depends on the
-       * quadrature formula
+       * Contains all the stuff that depends on the quadrature formula
 
        */
       std::vector<MappingInfoDependent> mapping_data_gen;
 
       /**
-       * Stores whether JxW values have been
-       * initialized
+       * Stores whether JxW values have been initialized
        */
       bool JxW_values_initialized;
 
       /**
-       * Stores whether we computed second
-       * derivatives.
+       * Stores whether we computed second derivatives.
        */
       bool second_derivatives_initialized;
 
       /**
-       * Stores whether we computed quadrature
-       * points.
+       * Stores whether we computed quadrature points.
        */
       bool quadrature_points_initialized;
 
       /**
-       * Internal temporary data used for the
-       * initialization.
+       * Internal temporary data used for the initialization.
        */
       struct CellData
       {
@@ -370,8 +310,7 @@ namespace internal
       };
 
       /**
-       * Helper function called internally during
-       * the initialize function.
+       * Helper function called internally during the initialize function.
        */
       void evaluate_on_cell (const dealii::Triangulation<dim> &tria,
                              const std::pair<unsigned int,unsigned int> *cells,
index 493e39b124a12f669c34356225729e3bc31fb4d3..647ed1476c1d20e0ee034ef560a52e7078f30e06 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2011 by the deal.II authors
+//    Copyright (C) 2011-2012 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -60,14 +60,11 @@ namespace internal
     compute_update_flags (const UpdateFlags update_flags,
                           const std::vector<dealii::hp::QCollection<1> > &quad) const
     {
-      // this class is build around the evaluation
-      // this class is build around the evaluation
-      // of inverse gradients, so compute them in
-      // any case
+      // this class is build around the evaluation this class is build around
+      // the evaluation of inverse gradients, so compute them in any case
       UpdateFlags new_flags = update_inverse_jacobians;
 
-      // if the user requested gradients, need
-      // inverse Jacobians
+      // if the user requested gradients, need inverse Jacobians
       if (update_flags & update_gradients || update_flags & update_inverse_jacobians)
         new_flags |= update_inverse_jacobians;
 
@@ -75,24 +72,19 @@ namespace internal
       if (update_flags & update_JxW_values)
         new_flags |= update_JxW_values;
 
-      // for Hessian information, need inverse
-      // Jacobians and the derivative of Jacobians
-      // (these two together will give use the
-      // gradients of the inverse Jacobians, which
-      // is what we need)
+      // for Hessian information, need inverse Jacobians and the derivative of
+      // Jacobians (these two together will give use the gradients of the
+      // inverse Jacobians, which is what we need)
       if (update_flags & update_hessians || update_flags & update_jacobian_grads)
         new_flags |= update_jacobian_grads;
 
       if (update_flags & update_quadrature_points)
         new_flags |= update_quadrature_points;
 
-      // there is one more thing: if we have a
-      // quadrature formula with only one quadrature
-      // point on the first component, but more
-      // points on later components, we need to have
-      // Jacobian gradients anyway in order to
-      // determine whether the Jacobian is constant
-      // throughout a cell
+      // there is one more thing: if we have a quadrature formula with only
+      // one quadrature point on the first component, but more points on later
+      // components, we need to have Jacobian gradients anyway in order to
+      // determine whether the Jacobian is constant throughout a cell
       bool formula_with_one_point = false;
       for (unsigned int i=0; i<quad[0].size(); ++i)
         if (quad[0][i].size() == 1)
@@ -147,11 +139,9 @@ end_set:
       mapping_data_gen.resize (n_quads);
       cell_type.resize (n_macro_cells);
 
-      // dummy FE that is used to set up an FEValues
-      // object. Do not need the actual finite
-      // element because we will only evaluate
-      // quantities for the mapping that are
-      // independent of the FE
+      // dummy FE that is used to set up an FEValues object. Do not need the
+      // actual finite element because we will only evaluate quantities for
+      // the mapping that are independent of the FE
       FE_Nothing<dim> dummy_fe;
       UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
 
@@ -162,25 +152,19 @@ end_set:
       if (update_flags & update_quadrature_points)
         quadrature_points_initialized = true;
 
-      // when we make comparisons about the size of
-      // Jacobians we need to know the approximate
-      // size of typical entries in Jacobians. We
-      // need to fix the Jacobian size once and for
-      // all. We choose the diameter of the first
-      // cell (on level zero, which is the best
-      // accuracy we can hope for, since diameters
-      // on finer levels are computed by differences
-      // of nearby cells). If the mesh extends over
-      // a certain domain, the precision of double
-      // values is essentially limited by this
-      // precision.
+      // when we make comparisons about the size of Jacobians we need to know
+      // the approximate size of typical entries in Jacobians. We need to fix
+      // the Jacobian size once and for all. We choose the diameter of the
+      // first cell (on level zero, which is the best accuracy we can hope
+      // for, since diameters on finer levels are computed by differences of
+      // nearby cells). If the mesh extends over a certain domain, the
+      // precision of double values is essentially limited by this precision.
       const double jacobian_size = internal::get_jacobian_size(tria);
 
-      // objects that hold the data for up to
-      // vectorization_length cells while we fill them up. Only
-      // after all vectorization_length cells have been
-      // processed, we can insert the data into the
-      // data structures of this class
+      // objects that hold the data for up to vectorization_length cells while
+      // we fill them up. Only after all vectorization_length cells have been
+      // processed, we can insert the data into the data structures of this
+      // class
       CellData data (jacobian_size);
 
       for (unsigned int my_q=0; my_q<n_quads; ++my_q)
@@ -218,11 +202,10 @@ end_set:
               if (n_hp_quads > 1)
                 current_data.quad_index_conversion[q] = n_q_points;
 
-              // To walk on the diagonal for lexicographic
-              // ordering, we have to jump one index ahead
-              // in each direction. For direction 0, this is
-              // just the next point, for direction 1, it
-              // means adding n_q_points_1d, and so on.
+              // To walk on the diagonal for lexicographic ordering, we have
+              // to jump one index ahead in each direction. For direction 0,
+              // this is just the next point, for direction 1, it means adding
+              // n_q_points_1d, and so on.
               step_size_cartesian[q] = 0;
               unsigned int factor = 1;
               for (unsigned int d=0; d<dim; ++d)
@@ -232,32 +215,26 @@ end_set:
                 }
             }
 
-          // if there are no cells, there is nothing to
-          // do
+          // if there are no cells, there is nothing to do
           if (cells.size() == 0)
             continue;
 
           Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
           Tensor<1,dim,VectorizedArray<Number> > tmp;
 
-          // encodes the cell types of the current
-          // cell. Since several cells must be
-          // considered together, this variable holds
-          // the individual info of the last chunk of
-          // cells
+          // encodes the cell types of the current cell. Since several cells
+          // must be considered together, this variable holds the individual
+          // info of the last chunk of cells
           CellType cell_t [vectorization_length],
                    cell_t_prev [vectorization_length];
           for (unsigned int j=0; j<vectorization_length; ++j)
             cell_t_prev[j] = undefined;
 
-          // fe_values object that is used to compute
-          // the mapping data. for the hp case there
-          // might be more than one finite
-          // element. since we manually select the
-          // active FE index and not via a
-          // hp::DoFHandler<dim>::active_cell_iterator,
-          // we need to manually select the correct
-          // finite element, so just hold a vector of
+          // fe_values object that is used to compute the mapping data. for
+          // the hp case there might be more than one finite element. since we
+          // manually select the active FE index and not via a
+          // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
+          // select the correct finite element, so just hold a vector of
           // FEValues
           std::vector<std_cxx1x::shared_ptr<FEValues<dim> > >
           fe_values (current_data.quadrature.size());
@@ -266,9 +243,8 @@ end_set:
             (update_flags & update_jacobian_grads ? update_jacobian_grads : update_default) |
             (update_flags & update_quadrature_points ? update_quadrature_points : update_default);
 
-          // resize the fields that have fixed size or
-          // for which we know something from an earlier
-          // loop
+          // resize the fields that have fixed size or for which we know
+          // something from an earlier loop
           current_data.rowstart_q_points.resize (n_macro_cells+1);
           if (my_q > 0)
             {
@@ -301,11 +277,10 @@ end_set:
           // loop over all cells
           for (unsigned int cell=0; cell<n_macro_cells; ++cell)
             {
-              // GENERAL OUTLINE: First generate the data in
-              // format "number" for vectorization_length cells, and
-              // then find the most general type of cell for
-              // appropriate vectorized formats. then fill
-              // this data in
+              // GENERAL OUTLINE: First generate the data in format "number"
+              // for vectorization_length cells, and then find the most
+              // general type of cell for appropriate vectorized formats. then
+              // fill this data in
               const unsigned int fe_index = active_fe_index.size() > 0 ?
                                             active_fe_index[cell] : 0;
               const unsigned int n_q_points = current_data.n_q_points[fe_index];
@@ -317,9 +292,8 @@ end_set:
               FEValues<dim> &fe_val = *fe_values[fe_index];
               data.resize (n_q_points);
 
-              // if the fe index has changed from the
-              // previous cell, set the old cell type to
-              // invalid (otherwise, we might detect
+              // if the fe index has changed from the previous cell, set the
+              // old cell type to invalid (otherwise, we might detect
               // similarity due to some cells further ahead)
               if (cell > 0 && active_fe_index.size() > 0 &&
                   active_fe_index[cell] != active_fe_index[cell-1])
@@ -327,20 +301,17 @@ end_set:
               evaluate_on_cell (tria, &cells[cell*vectorization_length],
                                 cell, my_q, cell_t_prev, cell_t, fe_val, data);
 
-              // now reorder the data into vectorized
-              // types. if we are here for the first time,
-              // we need to find out whether the Jacobian
-              // allows for some simplification (Cartesian,
-              // affine) taking vectorization_length cell together and
-              // we have to insert that data into the
-              // respective fields. Also, we have to
-              // compress different cell indicators into one
-              // structure.
+              // now reorder the data into vectorized types. if we are here
+              // for the first time, we need to find out whether the Jacobian
+              // allows for some simplification (Cartesian, affine) taking
+              // vectorization_length cell together and we have to insert that
+              // data into the respective fields. Also, we have to compress
+              // different cell indicators into one structure.
 
               if (my_q == 0)
                 {
-                  // find the most general cell type (most
-                  // general type is 2 (general cell))
+                  // find the most general cell type (most general type is 2
+                  // (general cell))
                   CellType most_general_type = cartesian;
                   for (unsigned int j=0; j<vectorization_length; ++j)
                     if (cell_t[j] > most_general_type)
@@ -348,18 +319,15 @@ end_set:
                   AssertIndexRange (most_general_type, 3);
                   unsigned int insert_position = numbers::invalid_unsigned_int;
 
-                  // Cartesian cell with diagonal Jacobian: only
-                  // insert the diagonal of the inverse and the
-                  // Jacobian determinant. We do this by using
-                  // an std::map that collects pointers to all
-                  // Cartesian Jacobians. We need a pointer in
-                  // the std::map because it cannot store data
-                  // based on VectorizedArray (alignment
-                  // issue). We circumvent the problem by
-                  // temporarily filling the next value into the
-                  // cartesian_data field and, in case we did an
-                  // insertion, the data is already in the
-                  // correct place.
+                  // Cartesian cell with diagonal Jacobian: only insert the
+                  // diagonal of the inverse and the Jacobian determinant. We
+                  // do this by using an std::map that collects pointers to
+                  // all Cartesian Jacobians. We need a pointer in the
+                  // std::map because it cannot store data based on
+                  // VectorizedArray (alignment issue). We circumvent the
+                  // problem by temporarily filling the next value into the
+                  // cartesian_data field and, in case we did an insertion,
+                  // the data is already in the correct place.
                   if (most_general_type == cartesian)
                     {
                       std::pair<Tensor<1,dim,VEC_ARRAY>,unsigned int> new_entry;
@@ -377,8 +345,8 @@ end_set:
                         insert_position = new_entry.second;
                     }
 
-                  // Constant Jacobian case. same strategy as
-                  // before, but with other data fields
+                  // Constant Jacobian case. same strategy as before, but with
+                  // other data fields
                   else if (most_general_type == affine)
                     {
                       std::pair<Tensor<2,dim,VEC_ARRAY>,unsigned int> new_entry;
@@ -397,13 +365,12 @@ end_set:
                         insert_position = new_entry.second;
                     }
 
-                  // general cell case: first resize the data
-                  // field to fit the new data. if we are here
-                  // the first time, assume that there are many
-                  // general cells to come, so reserve some
-                  // memory in order to not have too many
-                  // reallocations and memcpy's. The scheme used
-                  // here involves at most one reallocation.
+                  // general cell case: first resize the data field to fit the
+                  // new data. if we are here the first time, assume that
+                  // there are many general cells to come, so reserve some
+                  // memory in order to not have too many reallocations and
+                  // memcpy's. The scheme used here involves at most one
+                  // reallocation.
                   else
                     {
                       Assert (most_general_type == general, ExcInternalError());
@@ -429,10 +396,9 @@ end_set:
 
                 } // end if (my_q == 0)
 
-              // general cell case: now go through all
-              // quadrature points and collect the
-              // data. done for all different quadrature
-              // formulas, so do it outside the above loop.
+              // general cell case: now go through all quadrature points and
+              // collect the data. done for all different quadrature formulas,
+              // so do it outside the above loop.
               if (get_cell_type(cell) == general)
                 {
                   const unsigned int previous_size =
@@ -470,33 +436,30 @@ end_set:
                       current_data.jacobians.push_back (transpose(invert(jac)));
                       const Tensor<2,dim,VectorizedArray<Number> > &inv_jac = current_data.jacobians.back();
 
-                      // TODO: deal.II does not use abs on
-                      // determinants. Is there an assumption
-                      // somewhere that negative determinants are
-                      // not allowed?
+                      // TODO: deal.II does not use abs on determinants. Is
+                      // there an assumption somewhere that negative
+                      // determinants are not allowed?
                       if (update_flags & update_JxW_values)
                         current_data.JxW_values.push_back
                         (std::abs(det) * current_data.quadrature_weights[fe_index][q]);
 
                       if (update_flags & update_jacobian_grads)
                         {
-                          // for second derivatives on the real cell,
-                          // need also the gradient of the inverse
-                          // Jacobian J. This involves some calculus and
-                          // is done vectorized. This is very cheap
-                          // compared to what fe_values does (in early
-                          // 2011). If L is the gradient of the jacobian
-                          // on the unit cell, the gradient of the
-                          // inverse is given by (multidimensional
-                          // calculus) - J * (J * L) * J (the third J is
-                          // because we need to transform the gradient L
-                          // from the unit to the real cell, and then
-                          // apply the inverse Jacobian). Compare this
-                          // with 1D with j(x) = 1/k(phi(x)), where j =
-                          // phi' is the inverse of the jacobian and k
-                          // is the derivative of the jacobian on the
-                          // unit cell. Then j' = phi' k'/k^2 = j k'
-                          // j^2.
+                          // for second derivatives on the real cell, need
+                          // also the gradient of the inverse Jacobian J. This
+                          // involves some calculus and is done
+                          // vectorized. This is very cheap compared to what
+                          // fe_values does (in early 2011). If L is the
+                          // gradient of the jacobian on the unit cell, the
+                          // gradient of the inverse is given by
+                          // (multidimensional calculus) - J * (J * L) * J
+                          // (the third J is because we need to transform the
+                          // gradient L from the unit to the real cell, and
+                          // then apply the inverse Jacobian). Compare this
+                          // with 1D with j(x) = 1/k(phi(x)), where j = phi'
+                          // is the inverse of the jacobian and k is the
+                          // derivative of the jacobian on the unit cell. Then
+                          // j' = phi' k'/k^2 = j k' j^2.
 
                           // compute: jac_grad = J*grad_unit(J^-1)
                           for (unsigned int d=0; d<dim; ++d)
@@ -553,18 +516,15 @@ end_set:
 
               if (update_flags & update_quadrature_points)
                 {
-                  // eventually we turn to the quadrature points
-                  // that we can compress in case we have
-                  // Cartesian cells. we also need to reorder
-                  // them into arrays of vectorized data types.
-                  // first go through the cells and find out how
-                  // much memory we need to allocate for the
-                  // quadrature points. We store n_q_points for
-                  // all cells but Cartesian cells. For
-                  // Cartesian cells, only need to store the
-                  // values on a diagonal through the cell
-                  // (n_q_points_1d). This will give (with some
-                  // little indexing) the location of all
+                  // eventually we turn to the quadrature points that we can
+                  // compress in case we have Cartesian cells. we also need to
+                  // reorder them into arrays of vectorized data types.  first
+                  // go through the cells and find out how much memory we need
+                  // to allocate for the quadrature points. We store
+                  // n_q_points for all cells but Cartesian cells. For
+                  // Cartesian cells, only need to store the values on a
+                  // diagonal through the cell (n_q_points_1d). This will give
+                  // (with some little indexing) the location of all
                   // quadrature points.
                   const unsigned int old_size =
                     current_data.quadrature_points.size();
@@ -596,10 +556,9 @@ end_set:
           current_data.rowstart_q_points[n_macro_cells] =
             current_data.quadrature_points.size();
 
-          // finally, need to invert and
-          // transpose the Jacobians in the
-          // cartesian_data and affine_data
-          // fields and compute the JxW value.
+          // finally, fill the accumulated data for Cartesian and affine cells
+          //  into cartesian_data and affine_data, invert and transpose the
+          //  Jacobians, and compute the JxW value.
           if (my_q == 0)
             {
               cartesian_data.resize(cartesians.size());
@@ -656,8 +615,7 @@ end_set:
       const unsigned int n_q_points = fe_val.n_quadrature_points;
       const UpdateFlags update_flags = fe_val.get_update_flags();
 
-      // this should be the same value as
-      // used in HashValue::scaling (but we
+      // this should be the same value as used in HashValue::scaling (but we
       // not have that field here)
       const double zero_tolerance_double = data.jac_size *
                                            std::numeric_limits<double>::epsilon() * 1024.;
@@ -668,9 +626,8 @@ end_set:
           fe_val.reinit(cell_it);
           cell_t[j] = undefined;
 
-          // extract quadrature points and store them
-          // temporarily. if we have Cartesian cells, we
-          // can compress the indices
+          // extract quadrature points and store them temporarily. if we have
+          // Cartesian cells, we can compress the indices
           if (update_flags & update_quadrature_points)
             for (unsigned int q=0; q<n_q_points; ++q)
               {
@@ -679,17 +636,15 @@ end_set:
                   data.quadrature_points[q][d][j] = point[d];
               }
 
-          // if this is not the first quadrature formula
-          // and we already have determined that this
-          // cell is either Cartesian or with constant
+          // if this is not the first quadrature formula and we already have
+          // determined that this cell is either Cartesian or with constant
           // Jacobian, we have nothing more to do.
           if (my_q > 0 && (get_cell_type(cell) == cartesian
                            || get_cell_type(cell) == affine) )
             continue;
 
-          // first round: if the transformation is
-          // detected to be the same as on the old cell,
-          // we only need to copy over the data.
+          // first round: if the transformation is detected to be the same as
+          // on the old cell, we only need to copy over the data.
           if (fe_val.get_cell_similarity() == CellSimilarity::translation
               &&
               my_q == 0)
@@ -707,9 +662,8 @@ end_set:
 
           if (my_q == 0)
             {
-              // check whether the Jacobian is constant on
-              // this cell the first time we come around
-              // here
+              // check whether the Jacobian is constant on this cell the first
+              // time we come around here
               if (cell_t[j] == undefined)
                 {
                   bool jacobian_constant = true;
@@ -725,8 +679,8 @@ end_set:
                         break;
                     }
 
-                  // check whether the Jacobian is diagonal to
-                  // machine accuracy
+                  // check whether the Jacobian is diagonal to machine
+                  // accuracy
                   bool cell_cartesian = jacobian_constant;
                   for (unsigned int d=0; d<dim; ++d)
                     for (unsigned int e=0; e<dim; ++e)
@@ -779,8 +733,7 @@ end_set:
               // Cartesian cell
               if (cell_t[j] == cartesian)
                 {
-                  // set Jacobian into diagonal and clear
-                  // off-diagonal part
+                  // set Jacobian into diagonal and clear off-diagonal part
                   for (unsigned int d=0; d<dim; ++d)
                     {
                       data.const_jac[d][d][j] = jac_0[d][d];
@@ -808,30 +761,26 @@ end_set:
 
           // general cell case
 
-          // go through all quadrature points and fill
-          // in the data into the temporary data
-          // structures with slots for the vectorized
-          // data types
+          // go through all quadrature points and fill in the data into the
+          // temporary data structures with slots for the vectorized data
+          // types
           for (unsigned int q=0; q<n_q_points; ++q)
             {
 
-              // compress out very small numbers which are
-              // only noise. Then it is cleaner to use zero
-              // straight away (though it does not save any
-              // memory)
+              // compress out very small numbers which are only noise. Then it
+              // is cleaner to use zero straight away (though it does not save
+              // any memory)
               const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
               for (unsigned int d=0; d<dim; ++d)
                 for (unsigned int e=0; e<dim; ++e)
                   data.general_jac[q][d][e][j] =
                     std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
 
-              // need to do some calculus based on the
-              // gradient of the Jacobian, in order to find
-              // the gradient of the inverse Jacobian which
-              // is needed in user code. however, we would
-              // like to perform that on vectorized data
-              // types instead of doubles or floats. to this
-              // end, copy the gradients first
+              // need to do some calculus based on the gradient of the
+              // Jacobian, in order to find the gradient of the inverse
+              // Jacobian which is needed in user code. however, we would like
+              // to perform that on vectorized data types instead of doubles
+              // or floats. to this end, copy the gradients first
               if (update_flags & update_jacobian_grads)
                 {
                   const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
@@ -841,7 +790,8 @@ end_set:
                         data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
                 }
             }
-        } // end loop over all entries in vectorization (vectorization_length cells)
+        } // end loop over all entries in vectorization (vectorization_length
+          // cells)
 
       // set information for next cell
       for (unsigned int j=0; j<vectorization_length; ++j)
@@ -912,10 +862,8 @@ end_set:
     (STREAM         &out,
      const SizeInfo &size_info) const
     {
-      // print_memory_statistics involves
-      // global communication, so we can
-      // disable the check here only if no
-      // processor has any such data
+      // print_memory_statistics involves global communication, so we can
+      // disable the check here only if no processor has any such data
 #if DEAL_II_COMPILER_SUPPORTS_MPI
       unsigned int general_size_glob = 0, general_size_loc = jacobians.size();
       MPI_Allreduce (&general_size_loc, &general_size_glob, 1, MPI_UNSIGNED,
index 8351f2e8dc5d6e8b9c7c88db7b09d8139e6d2fa6..97b71f042584fa97845953baa0001e454db0fa33 100644 (file)
@@ -51,10 +51,10 @@ DEAL_II_NAMESPACE_OPEN
 
 /**
  * This class collects all the data that is stored for the matrix free
- * implementation. The storage scheme is tailored towards several
- * loops performed with the same data, i.e., typically doing many
- * matrix-vector products or residual computations on the same
- * mesh. The class is used in step-37 and step-48.
+ * implementation. The storage scheme is tailored towards several loops
+ * performed with the same data, i.e., typically doing many matrix-vector
+ * products or residual computations on the same mesh. The class is used in
+ * step-37 and step-48.
  *
  * This class does not implement any operations involving finite element basis
  * functions, i.e., regarding the operation performed on the cells. For these
@@ -103,44 +103,32 @@ class MatrixFree
 public:
 
   /**
-   * Collects the options for initialization of
-   * the MatrixFree class. The first
-   * parameter specifies the MPI communicator to
-   * be used, the second the parallelization
-   * options in shared memory (task-based
-   * parallelism, where one can choose between
-   * no parallelism and three schemes that avoid
-   * that cells with access to the same vector
-   * entries are accessed simultaneously), the
-   * third with the block size for task parallel
-   * scheduling, the fourth the update flags
-   * that should be stored by this class.
+   * Collects the options for initialization of the MatrixFree class. The
+   * first parameter specifies the MPI communicator to be used, the second the
+   * parallelization options in shared memory (task-based parallelism, where
+   * one can choose between no parallelism and three schemes that avoid that
+   * cells with access to the same vector entries are accessed
+   * simultaneously), the third with the block size for task parallel
+   * scheduling, the fourth the update flags that should be stored by this
+   * class.
    *
-   * The fifth parameter specifies the level in
-   * the triangulation from which the indices
-   * are to be used. If the level is set to
-   * numbers::invalid_unsigned_int, the active
-   * cells are traversed, and otherwise the
-   * cells in the given level. This option has
-   * no effect in case a DoFHandler or
-   * hp::DoFHandler is given.
+   * The fifth parameter specifies the level in the triangulation from which
+   * the indices are to be used. If the level is set to
+   * numbers::invalid_unsigned_int, the active cells are traversed, and
+   * otherwise the cells in the given level. This option has no effect in case
+   * a DoFHandler or hp::DoFHandler is given.
    *
-   * The parameter @p initialize_plain_indices
-   * indicates whether the DoFInfo class should
-   * also allow for access to vectors without
-   * resolving constraints.
+   * The parameter @p initialize_plain_indices indicates whether the DoFInfo
+   * class should also allow for access to vectors without resolving
+   * constraints.
    *
-   * The last two parameters allow the user to
-   * disable some of the initialization
-   * processes. For example, if only the
-   * scheduling that avoids touching the same
-   * vector/matrix indices simultaneously is to
-   * be found, the mapping needs not be
-   * initialized. Likewise, if the mapping has
-   * changed from one iteration to the next but
-   * the topology has not (like when using a
-   * deforming mesh with MappingQEulerian), it
-   * suffices to initialize the mapping only.
+   * The last two parameters allow the user to disable some of the
+   * initialization processes. For example, if only the scheduling that avoids
+   * touching the same vector/matrix indices simultaneously is to be found,
+   * the mapping needs not be initialized. Likewise, if the mapping has
+   * changed from one iteration to the next but the topology has not (like
+   * when using a deforming mesh with MappingQEulerian), it suffices to
+   * initialize the mapping only.
    */
   struct AdditionalData
   {
@@ -172,142 +160,100 @@ public:
     {};
 
     /**
-     * Sets the MPI communicator that the parallel
-     * layout of the operator should be based
-     * upon. Defaults to MPI_COMM_SELF, but should
-     * be set to a communicator similar to the one
-     * used for a distributed triangulation in
-     * order to inform this class over all cells
-     * that are present.
+     * Sets the MPI communicator that the parallel layout of the operator
+     * should be based upon. Defaults to MPI_COMM_SELF, but should be set to a
+     * communicator similar to the one used for a distributed triangulation in
+     * order to inform this class over all cells that are present.
      */
     MPI_Comm            mpi_communicator;
 
     /**
-     * Sets the scheme for task parallelism. There
-     * are four options available. If set to @p
-     * none, the operator application is done in
-     * serial without shared memory
-     * parallelism. If this class is used together
-     * with MPI and MPI is also used for
-     * parallelism within the nodes, this flag
-     * should be set to @p none. The default value
-     * is @p partition_partition, i.e. we actually
-     * use multithreading with the first option
-     * described below.
+     * Sets the scheme for task parallelism. There are four options
+     * available. If set to @p none, the operator application is done in
+     * serial without shared memory parallelism. If this class is used
+     * together with MPI and MPI is also used for parallelism within the
+     * nodes, this flag should be set to @p none. The default value is @p
+     * partition_partition, i.e. we actually use multithreading with the first
+     * option described below.
      *
-     * The first option @p partition_partition is
-     * to partition the cells on two levels in
-     * onion-skin-like partitions and forming
-     * chunks of tasks_block_size after the
-     * partitioning. The partitioning finds sets
-     * of independent cells that enable working in
-     * parallel without accessing the same vector
-     * entries at the same time.
+     * The first option @p partition_partition is to partition the cells on
+     * two levels in onion-skin-like partitions and forming chunks of
+     * tasks_block_size after the partitioning. The partitioning finds sets of
+     * independent cells that enable working in parallel without accessing the
+     * same vector entries at the same time.
      *
-     * The second option @p partition_color is to
-     * use a partition on the global level and
-     * color cells within the partitions (where
-     * all chunks within a color are
-     * independent). Here, the subdivision into
-     * chunks of cells is done before the
-     * partitioning, which might give worse
-     * partitions but better cache performance if
-     * degrees of freedom are not renumbered.
+     * The second option @p partition_color is to use a partition on the
+     * global level and color cells within the partitions (where all chunks
+     * within a color are independent). Here, the subdivision into chunks of
+     * cells is done before the partitioning, which might give worse
+     * partitions but better cache performance if degrees of freedom are not
+     * renumbered.
      *
-     * The third option @p color is to use a
-     * traditional algorithm of coloring on the
-     * global level. This scheme is a special case
-     * of the second option where only one
-     * partition is present. Note that for
-     * problems with hanging nodes, there are
-     * quite many colors (50 or more in 3D), which
-     * might degrade parallel performance (bad
-     * cache behavior, many synchronization
-     * points).
+     * The third option @p color is to use a traditional algorithm of coloring
+     * on the global level. This scheme is a special case of the second option
+     * where only one partition is present. Note that for problems with
+     * hanging nodes, there are quite many colors (50 or more in 3D), which
+     * might degrade parallel performance (bad cache behavior, many
+     * synchronization points).
      */
     TasksParallelScheme tasks_parallel_scheme;
 
     /**
-     * Sets the number of so-called macro cells
-     * that should form one partition. If zero
-     * size is given, the class tries to find a
-     * good size for the blocks based on the
-     * number of threads stored in
-     * multithread_info.n_default_threads and the
-     * number of cells present. Otherwise, the
-     * given number is used. If the given number
-     * is larger than one third of the number of
-     * total cells, this means no
-     * parallelism. Note that in the case
-     * vectorization is used, a macro cell
+     * Sets the number of so-called macro cells that should form one
+     * partition. If zero size is given, the class tries to find a good size
+     * for the blocks based on the number of threads stored in
+     * multithread_info.n_default_threads and the number of cells
+     * present. Otherwise, the given number is used. If the given number is
+     * larger than one third of the number of total cells, this means no
+     * parallelism. Note that in the case vectorization is used, a macro cell
      * consists of more than one physical cell.
      */
     unsigned int        tasks_block_size;
 
     /**
-     * This flag is used to determine which
-     * quantities should be cached. This class can
-     * cache data needed for gradient computations
-     * (inverse Jacobians), Jacobian determinants
-     * (JxW), quadrature points as well as data
-     * for Hessians (derivative of Jacobians). By
-     * default, only data for gradients and
-     * Jacobian determinants times quadrature
-     * weights, JxW, are cached. If quadrature
-     * points or second derivatives are needed,
-     * they must be specified by this field (even
-     * though second derivatives might still be
-     * evaluated on Cartesian cells without this
-     * option set here, since there the Jacobian
-     * describes the mapping completely).
+     * This flag is used to determine which quantities should be cached. This
+     * class can cache data needed for gradient computations (inverse
+     * Jacobians), Jacobian determinants (JxW), quadrature points as well as
+     * data for Hessians (derivative of Jacobians). By default, only data for
+     * gradients and Jacobian determinants times quadrature weights, JxW, are
+     * cached. If quadrature points or second derivatives are needed, they
+     * must be specified by this field (even though second derivatives might
+     * still be evaluated on Cartesian cells without this option set here,
+     * since there the Jacobian describes the mapping completely).
      */
     UpdateFlags         mapping_update_flags;
 
     /**
-     * If working on a MGDoFHandler, this option
-     * can be used to define whether we work on a
-     * certain level of the mesh, and not the
-     * active cells. If set to
-     * invalid_unsigned_int (which is the default
-     * value), the active cells are gone through,
-     * otherwise the level given by this
-     * parameter.
+     * If working on a MGDoFHandler, this option can be used to define whether
+     * we work on a certain level of the mesh, and not the active cells. If
+     * set to invalid_unsigned_int (which is the default value), the active
+     * cells are gone through, otherwise the level given by this parameter.
      */
     unsigned int        level_mg_handler;
 
     /**
-     * Controls whether to allow reading from
-     * vectors without resolving constraints,
-     * i.e., just read the local values of the
-     * vector. By default, this option is
-     * disabled, so if you want to use
-     * FEEvaluationBase::read_dof_values_plain,
-     * this flag needs to be set.
+     * Controls whether to allow reading from vectors without resolving
+     * constraints, i.e., just read the local values of the vector. By
+     * default, this option is disabled, so if you want to use
+     * FEEvaluationBase::read_dof_values_plain, this flag needs to be set.
      */
     bool                store_plain_indices;
 
     /**
-     * Option to control whether the indices
-     * stored in the DoFHandler should be read and
-     * the pattern for task parallelism should be
-     * set up in the initialize method of
-     * MatrixFree. Defaults to true. Can be
-     * disabled in case the mapping should be
-     * recomputed (e.g. when using a deforming
-     * mesh described through MappingEulerian) but
-     * the topology of cells has remained the
-     * same.
+     * Option to control whether the indices stored in the DoFHandler should
+     * be read and the pattern for task parallelism should be set up in the
+     * initialize method of MatrixFree. Defaults to true. Can be disabled in
+     * case the mapping should be recomputed (e.g. when using a deforming mesh
+     * described through MappingEulerian) but the topology of cells has
+     * remained the same.
      */
     bool                initialize_indices;
 
     /**
-     * Option to control whether the mapping
-     * information should be computed in the
-     * initialize method of
-     * MatrixFree. Defaults to true. Can be
-     * disabled when only some indices should be
-     * set up (e.g. when only a set of independent
-     * cells should be computed).
+     * Option to control whether the mapping information should be computed in
+     * the initialize method of MatrixFree. Defaults to true. Can be disabled
+     * when only some indices should be set up (e.g. when only a set of
+     * independent cells should be computed).
      */
     bool                initialize_mapping;
   };
@@ -317,8 +263,7 @@ public:
    */
   //@{
   /**
-   * Default empty constructor. Does
-   * nothing.
+   * Default empty constructor. Does nothing.
    */
   MatrixFree ();
 
@@ -328,34 +273,20 @@ public:
   ~MatrixFree();
 
   /**
-  * Extracts the information needed to
-  * perform loops over cells. The
-  * DoFHandler and ConstraintMatrix
-  * describe the layout of degrees of
-  * freedom, the DoFHandler and the
-  * mapping describe the transformations
-  * from unit to real cell, and the finite
-  * element underlying the DoFHandler
-  * together with the quadrature formula
-  * describe the local operations. Note
-  * that the finite element underlying the
-  * DoFHandler must either be scalar or
-  * contain several copies of the same
-  * element. Mixing several different
-  * elements into one FESystem is not
-  * allowed. In that case, use the
-  * initialization function with several
-  * DoFHandler arguments.
+  * Extracts the information needed to perform loops over cells. The
+  * DoFHandler and ConstraintMatrix describe the layout of degrees of freedom,
+  * the DoFHandler and the mapping describe the transformations from unit to
+  * real cell, and the finite element underlying the DoFHandler together with
+  * the quadrature formula describe the local operations. Note that the finite
+  * element underlying the DoFHandler must either be scalar or contain several
+  * copies of the same element. Mixing several different elements into one
+  * FESystem is not allowed. In that case, use the initialization function
+  * with several DoFHandler arguments.
   *
-  * The @p IndexSet @p locally_owned_dofs
-  * is used to specify the parallel
-  * partitioning with MPI. Usually, this
-  * needs not be specified, and the other
-  * initialization function without and @p
-  * IndexSet description can be used,
-  * which gets the partitioning
-  * information builtin into the
-  * DoFHandler.
+  * The @p IndexSet @p locally_owned_dofs is used to specify the parallel
+  * partitioning with MPI. Usually, this needs not be specified, and the other
+  * initialization function without and @p IndexSet description can be used,
+  * which gets the partitioning information builtin into the DoFHandler.
   */
   template <typename DH, typename Quadrature>
   void reinit (const Mapping<dim>     &mapping,
@@ -366,10 +297,8 @@ public:
                const AdditionalData    additional_data = AdditionalData());
 
   /**
-  * Initializes the data structures. Same
-  * as above, but with index set stored in
-  * the DoFHandler for describing the
-  * locally owned degrees of freedom.
+  * Initializes the data structures. Same as above, but with index set stored
+  * in the DoFHandler for describing the locally owned degrees of freedom.
   */
   template <typename DH, typename Quadrature>
   void reinit (const Mapping<dim>     &mapping,
@@ -379,8 +308,7 @@ public:
                const AdditionalData    additional_data = AdditionalData());
 
   /**
-  * Initializes the data structures. Same
-  * as above, but with mapping @p
+  * Initializes the data structures. Same as above, but with mapping @p
   * MappingQ1.
   */
   template <typename DH, typename Quadrature>
@@ -390,51 +318,30 @@ public:
                const AdditionalData    additional_data = AdditionalData());
 
   /**
-  * Extracts the information needed to
-  * perform loops over cells. The
-  * DoFHandler and ConstraintMatrix
-  * describe the layout of degrees of
-  * freedom, the DoFHandler and the
-  * mapping describe the transformations
-  * from unit to real cell, and the finite
-  * element underlying the DoFHandler
-  * together with the quadrature formula
-  * describe the local operations. As
-  * opposed to the scalar case treated
-  * with the other initialization
-  * functions, this function allows for
-  * problems with two or more different
-  * finite elements. The DoFHandlers to
-  * each element must be passed as
-  * pointers to the initialization
-  * function. Note that the finite element
-  * underlying an DoFHandler must either
-  * be scalar or contain several copies of
-  * the same element. Mixing several
-  * different elements into one @p
-  * FE_System is not allowed.
+  * Extracts the information needed to perform loops over cells. The
+  * DoFHandler and ConstraintMatrix describe the layout of degrees of freedom,
+  * the DoFHandler and the mapping describe the transformations from unit to
+  * real cell, and the finite element underlying the DoFHandler together with
+  * the quadrature formula describe the local operations. As opposed to the
+  * scalar case treated with the other initialization functions, this function
+  * allows for problems with two or more different finite elements. The
+  * DoFHandlers to each element must be passed as pointers to the
+  * initialization function. Note that the finite element underlying an
+  * DoFHandler must either be scalar or contain several copies of the same
+  * element. Mixing several different elements into one @p FE_System is not
+  * allowed.
   *
-  * This function also allows for using
-  * several quadrature formulas, e.g. when
-  * the description contains independent
-  * integrations of elements of different
-  * degrees. However, the number of
-  * different quadrature formulas can be
-  * sets independently from the number of
-  * DoFHandlers, when several elements are
-  * always integrated with the same
-  * quadrature formula.
+  * This function also allows for using several quadrature formulas, e.g. when
+  * the description contains independent integrations of elements of different
+  * degrees. However, the number of different quadrature formulas can be sets
+  * independently from the number of DoFHandlers, when several elements are
+  * always integrated with the same quadrature formula.
   *
-  * The @p IndexSet @p locally_owned_dofs
-  * is used to specify the parallel
-  * partitioning with MPI. Usually, this
-  * needs not be specified, and the other
-  * initialization function without and @p
-  * IndexSet description can be used,
-  * which gets the partitioning
-  * information from the DoFHandler. This
-  * is the most general initialization
-  * function.
+  * The @p IndexSet @p locally_owned_dofs is used to specify the parallel
+  * partitioning with MPI. Usually, this needs not be specified, and the other
+  * initialization function without and @p IndexSet description can be used,
+  * which gets the partitioning information from the DoFHandler. This is the
+  * most general initialization function.
   */
   template <typename DH, typename Quadrature>
   void reinit (const Mapping<dim>                         &mapping,
@@ -445,10 +352,8 @@ public:
                const AdditionalData                        additional_data = AdditionalData());
 
   /**
-  * Initializes the data structures. Same
-  * as before, but now the index set
-  * description of the locally owned range
-  * of degrees of freedom is taken from
+  * Initializes the data structures. Same as before, but now the index set
+  * description of the locally owned range of degrees of freedom is taken from
   * the DoFHandler.
   */
   template <typename DH, typename Quadrature>
@@ -459,8 +364,7 @@ public:
                const AdditionalData                        additional_data = AdditionalData());
 
   /**
-  * Initializes the data structures. Same
-  * as above, but with mapping @p
+  * Initializes the data structures. Same as above, but with mapping @p
   * MappingQ1.
   */
   template <typename DH, typename Quadrature>
@@ -470,16 +374,11 @@ public:
                const AdditionalData                        additional_data = AdditionalData());
 
   /**
-  * Initializes the data structures. Same
-  * as before, but now the index set
-  * description of the locally owned range
-  * of degrees of freedom is taken from
-  * the DoFHandler. Moreover, only a
-  * single quadrature formula is used, as
-  * might be necessary when several
-  * components in a vector-valued problem
-  * are integrated together based on the
-  * same quadrature formula.
+  * Initializes the data structures. Same as before, but now the index set
+  * description of the locally owned range of degrees of freedom is taken from
+  * the DoFHandler. Moreover, only a single quadrature formula is used, as
+  * might be necessary when several components in a vector-valued problem are
+  * integrated together based on the same quadrature formula.
   */
   template <typename DH, typename Quadrature>
   void reinit (const Mapping<dim>                         &mapping,
@@ -489,8 +388,7 @@ public:
                const AdditionalData                        additional_data = AdditionalData());
 
   /**
-  * Initializes the data structures. Same
-  * as above, but with mapping @p
+  * Initializes the data structures. Same as above, but with mapping @p
   * MappingQ1.
   */
   template <typename DH, typename Quadrature>
@@ -500,18 +398,15 @@ public:
                const AdditionalData                        additional_data = AdditionalData());
 
   /**
-   * Copy function. Creates a deep copy of all
-   * data structures. It is usually enough to
-   * keep the data for different operations
-   * once, so this function should not be needed
-   * very often.
+   * Copy function. Creates a deep copy of all data structures. It is usually
+   * enough to keep the data for different operations once, so this function
+   * should not be needed very often.
    */
   void copy_from (const MatrixFree<dim,Number> &matrix_free_base);
 
   /**
-   * Clears all data fields and brings the class
-   * into a condition similar to after having
-   * called the default constructor.
+   * Clears all data fields and brings the class into a condition similar to
+   * after having called the default constructor.
    */
   void clear();
 
@@ -522,29 +417,17 @@ public:
    */
   //@{
   /**
-   * This method runs the loop over all
-   * cells (in parallel) and performs
-   * the MPI data exchange on the source
-   * vector and destination vector. The
-   * first argument indicates a function
-   * object that has the following
-   * signature: <code>cell_operation
-   * (const MatrixFree<dim,Number>
-   * &, OutVector &, InVector &,
-   * std::pair<unsigned int,unsigned
-   * int> &)</code>, where the first
-   * argument passes the data of the
-   * calling class and the last argument
-   * defines the range of cells which
-   * should be worked on (typically more
-   * than one cell should be worked on
-   * in order to reduce overheads).  One
-   * can pass a pointer to an object in
-   * this place if it has an
-   * <code>operator()</code> with the
-   * correct set of arguments since such
-   * a pointer can be converted to the
-   * function object.
+   * This method runs the loop over all cells (in parallel) and performs the
+   * MPI data exchange on the source vector and destination vector. The first
+   * argument indicates a function object that has the following signature:
+   * <code>cell_operation (const MatrixFree<dim,Number> &, OutVector &,
+   * InVector &, std::pair<unsigned int,unsigned int> &)</code>, where the
+   * first argument passes the data of the calling class and the last argument
+   * defines the range of cells which should be worked on (typically more than
+   * one cell should be worked on in order to reduce overheads).  One can pass
+   * a pointer to an object in this place if it has an <code>operator()</code>
+   * with the correct set of arguments since such a pointer can be converted
+   * to the function object.
    */
   template <typename OutVector, typename InVector>
   void cell_loop (const std_cxx1x::function<void (const MatrixFree<dim,Number> &,
@@ -556,18 +439,13 @@ public:
                   const InVector &src) const;
 
   /**
-   * This is the second variant to run the loop
-   * over all cells, now providing a function
-   * pointer to a member function of class @p
-   * CLASS with the signature
-   * <code>cell_operation (const
-   * MatrixFree<dim,Number> &, OutVector &,
-   * InVector &, std::pair<unsigned int,unsigned
-   * int>&)const</code>. This method obviates the need
-   * to call std_cxx1x::bind to bind the class
-   * into the given function in case the local
-   * function needs to access data in the class
-   * (i.e., it is a non-static member function).
+   * This is the second variant to run the loop over all cells, now providing
+   * a function pointer to a member function of class @p CLASS with the
+   * signature <code>cell_operation (const MatrixFree<dim,Number> &, OutVector
+   * &, InVector &, std::pair<unsigned int,unsigned int>&)const</code>. This
+   * method obviates the need to call std_cxx1x::bind to bind the class into
+   * the given function in case the local function needs to access data in the
+   * class (i.e., it is a non-static member function).
    */
   template <typename CLASS, typename OutVector, typename InVector>
   void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
@@ -580,8 +458,7 @@ public:
                   const InVector &src) const;
 
   /**
-   * Same as above, but for class member
-   * functions which are non-const.
+   * Same as above, but for class member functions which are non-const.
    */
   template <typename CLASS, typename OutVector, typename InVector>
   void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
@@ -594,14 +471,11 @@ public:
                   const InVector &src) const;
 
   /**
-   * In the hp adaptive case, a subrange of
-   * cells as computed during the cell loop
-   * might contain elements of different
-   * degrees. Use this function to compute what
-   * the subrange for an individual finite
-   * element degree is. The finite element
-   * degree is associated to the vector
-   * component given in the function call.
+   * In the hp adaptive case, a subrange of cells as computed during the cell
+   * loop might contain elements of different degrees. Use this function to
+   * compute what the subrange for an individual finite element degree is. The
+   * finite element degree is associated to the vector component given in the
+   * function call.
    */
   std::pair<unsigned int,unsigned int>
   create_cell_subrange_hp (const std::pair<unsigned int,unsigned int> &range,
@@ -609,13 +483,10 @@ public:
                            const unsigned int vector_component = 0) const;
 
   /**
-   * In the hp adaptive case, a subrange of
-   * cells as computed during the cell loop
-   * might contain elements of different
-   * degrees. Use this function to compute what
-   * the subrange for a given index the hp
-   * finite element, as opposed to the finite
-   * element degree in the other function.
+   * In the hp adaptive case, a subrange of cells as computed during the cell
+   * loop might contain elements of different degrees. Use this function to
+   * compute what the subrange for a given index the hp finite element, as
+   * opposed to the finite element degree in the other function.
    */
   std::pair<unsigned int,unsigned int>
   create_cell_subrange_hp_by_index (const std::pair<unsigned int,unsigned int> &range,
@@ -629,15 +500,11 @@ public:
    */
   //@{
   /**
-   * Initialize function for a general
-   * vector. The length of the vector is equal
-   * to the total number of degrees in the
-   * DoFHandler. If the vector is of class
-   * parallel::distributed::Vector@<Number@>, the ghost
-   * entries are set accordingly. For
-   * vector-valued problems with several
-   * DoFHandlers underlying this class, the
-   * parameter @p vector_component defines which
+   * Initialize function for a general vector. The length of the vector is
+   * equal to the total number of degrees in the DoFHandler. If the vector is
+   * of class parallel::distributed::Vector@<Number@>, the ghost entries are
+   * set accordingly. For vector-valued problems with several DoFHandlers
+   * underlying this class, the parameter @p vector_component defines which
    * component is to be used.
    */
   template <typename VectorType>
@@ -645,15 +512,11 @@ public:
                              const unsigned int vector_component=0) const;
 
   /**
-   * Initialize function for a distributed
-   * vector. The length of the vector is equal
-   * to the total number of degrees in the
-   * DoFHandler. If the vector is of class
-   * parallel::distributed::Vector@<Number@>, the ghost
-   * entries are set accordingly. For
-   * vector-valued problems with several
-   * DoFHandlers underlying this class, the
-   * parameter @p vector_component defines which
+   * Initialize function for a distributed vector. The length of the vector is
+   * equal to the total number of degrees in the DoFHandler. If the vector is
+   * of class parallel::distributed::Vector@<Number@>, the ghost entries are
+   * set accordingly. For vector-valued problems with several DoFHandlers
+   * underlying this class, the parameter @p vector_component defines which
    * component is to be used.
    */
   template <typename Number2>
@@ -661,55 +524,42 @@ public:
                              const unsigned int vector_component=0) const;
 
   /**
-   * Returns the partitioner that represents the
-   * locally owned data and the ghost indices
-   * where access is needed to for the cell
-   * loop. The partitioner is constructed from
-   * the locally owned dofs and ghost dofs given
-   * by the respective fields. If you want to
-   * have specific information about these
-   * objects, you can query them with the
-   * respective access functions. If you just
-   * want to initialize a (parallel) vector, you
-   * should usually prefer this data structure
-   * as the data exchange information can be
-   * reused from one vector to another.
+   * Returns the partitioner that represents the locally owned data and the
+   * ghost indices where access is needed to for the cell loop. The
+   * partitioner is constructed from the locally owned dofs and ghost dofs
+   * given by the respective fields. If you want to have specific information
+   * about these objects, you can query them with the respective access
+   * functions. If you just want to initialize a (parallel) vector, you should
+   * usually prefer this data structure as the data exchange information can
+   * be reused from one vector to another.
    */
   const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &
   get_vector_partitioner (const unsigned int vector_component=0) const;
 
   /**
-   * Returns the set of cells that are
-   * oned by the processor.
+   * Returns the set of cells that are oned by the processor.
    */
   const IndexSet &
   get_locally_owned_set (const unsigned int fe_component = 0) const;
 
   /**
-   * Returns the set of ghost cells
-   * needed but not owned by the
-   * processor.
+   * Returns the set of ghost cells needed but not owned by the processor.
    */
   const IndexSet &
   get_ghost_set (const unsigned int fe_component = 0) const;
 
   /**
-   * Returns a list of all degrees of freedom
-   * that are constrained. The list is returned
-   * in local index space for the locally owned
-   * range of the vector, not in global
-   * numbers. In addition, it only returns the
-   * indices for degrees of freedom that are
-   * owned locally, not for ghosts.
+   * Returns a list of all degrees of freedom that are constrained. The list
+   * is returned in local index space for the locally owned range of the
+   * vector, not in global numbers. In addition, it only returns the indices
+   * for degrees of freedom that are owned locally, not for ghosts.
    */
   const std::vector<unsigned int> &
   get_constrained_dofs (const unsigned int fe_component = 0) const;
 
   /**
-   * Calls renumber_dofs function in dof
-   * info which renumbers the the
-   * degrees of freedom according to the
-   * ordering for parallelization.
+   * Calls renumber_dofs function in dof info which renumbers the the degrees
+   * of freedom according to the ordering for parallelization.
    */
   void renumber_dofs (std::vector<unsigned int> &renumbering,
                       const unsigned int vector_component = 0);
@@ -721,85 +571,61 @@ public:
    */
   //@{
   /**
-   * Returns the number of different DoFHandlers
-   * specified at initialization.
+   * Returns the number of different DoFHandlers specified at initialization.
    */
   unsigned int n_components () const;
 
   /**
-   * Returns the number of cells this structure
-   * is based on. If you are using a usual
-   * DoFHandler, it corresponds to the number of
-   * (locally owned) active cells. Note that
-   * most data structures in this class do not
-   * directly act on this number but rather on
-   * n_macro_cells() which gives the number of
-   * cells as seen when lumping several cells
-   * together with vectorization.
+   * Returns the number of cells this structure is based on. If you are using
+   * a usual DoFHandler, it corresponds to the number of (locally owned)
+   * active cells. Note that most data structures in this class do not
+   * directly act on this number but rather on n_macro_cells() which gives the
+   * number of cells as seen when lumping several cells together with
+   * vectorization.
    */
   unsigned int n_physical_cells () const;
 
   /**
-   * Returns the number of macro cells
-   * that this structure works on, i.e.,
-   * the number of cell chunks that are
-   * worked on after the application of
-   * vectorization which in general
-   * works on several cells at once. The
-   * cell range in @p cell_loop runs
-   * from zero to n_macro_cells()
-   * (exclusive), so this is the
-   * appropriate size if you want to
-   * store arrays of data for all cells
-   * to be worked on. This number is
-   * approximately
-   * n_physical_cells()/VectorizedArray::n_array_elements
-   * (depending on how many cell chunks
-   * that do not get filled up
-   * completely).
+   * Returns the number of macro cells that this structure works on, i.e., the
+   * number of cell chunks that are worked on after the application of
+   * vectorization which in general works on several cells at once. The cell
+   * range in @p cell_loop runs from zero to n_macro_cells() (exclusive), so
+   * this is the appropriate size if you want to store arrays of data for all
+   * cells to be worked on. This number is approximately
+   * n_physical_cells()/VectorizedArray::n_array_elements (depending on how
+   * many cell chunks that do not get filled up completely).
    */
   unsigned int n_macro_cells () const;
 
   /**
-   * In case this structure was built based on a
-   * DoFHandler, this returns the
-   * DoFHandler. Note that this function returns
-   * an exception in case the structure was
-   * based on MGDoFHandler and a level has been
-   * specified in InitializationOption.
+   * In case this structure was built based on a DoFHandler, this returns the
+   * DoFHandler. Note that this function returns an exception in case the
+   * structure was based on MGDoFHandler and a level has been specified in
+   * InitializationOption.
    */
   const DoFHandler<dim> &
   get_dof_handler (const unsigned int fe_component = 0) const;
 
   /**
-   * In case this structure was built based on a
-   * DoFHandler, this returns the
-   * DoFHandler. Note that this function returns
-   * an exception in case the structure was
-   * based on MGDoFHandler and a level has been
-   * specified in InitializationOption.
+   * In case this structure was built based on a DoFHandler, this returns the
+   * DoFHandler. Note that this function returns an exception in case the
+   * structure was based on MGDoFHandler and a level has been specified in
+   * InitializationOption.
    */
   const MGDoFHandler<dim> &
   get_mg_dof_handler (const unsigned int fe_component = 0) const;
 
   /**
-   * This returns the cell iterator in deal.II
-   * speak to a given cell in the renumbering of
-   * this structure. This function returns an
-   * exception in case the structure was
-   * constructed based on an MGDoFHandler with
-   * level specified, as these cells are in
-   * general not active.
+   * This returns the cell iterator in deal.II speak to a given cell in the
+   * renumbering of this structure. This function returns an exception in case
+   * the structure was constructed based on an MGDoFHandler with level
+   * specified, as these cells are in general not active.
    *
-   * Note that the cell iterators in deal.II go
-   * through cells differently to what the cell
-   * loop of this class does. This is because
-   * several cells are worked on together
-   * (vectorization), and since cells with
-   * neighbors on different MPI processors need
-   * to be accessed at a certain time when
-   * accessing remote data and overlapping
-   * communication with computation.
+   * Note that the cell iterators in deal.II go through cells differently to
+   * what the cell loop of this class does. This is because several cells are
+   * worked on together (vectorization), and since cells with neighbors on
+   * different MPI processors need to be accessed at a certain time when
+   * accessing remote data and overlapping communication with computation.
    */
   typename DoFHandler<dim>::active_cell_iterator
   get_cell_iterator (const unsigned int macro_cell_number,
@@ -807,22 +633,15 @@ public:
                      const unsigned int fe_component = 0) const;
 
   /**
-   * This returns the cell iterator in deal.II
-   * speak to a given cell in the renumbering of
-   * this structure. This function returns an
-   * exception in case the structure was
-   * constructed based on a DoFHandler and not
-   * MGDoFHandler.
+   * This returns the cell iterator in deal.II speak to a given cell in the
+   * renumbering of this structure. This function returns an exception in case
+   * the structure was constructed based on a DoFHandler and not MGDoFHandler.
    *
-   * Note that the cell iterators in deal.II go
-   * through cells differently to what the cell
-   * loop of this class does. This is because
-   * several cells are worked on together
-   * (vectorization), and since cells with
-   * neighbors on different MPI processors need
-   * to be accessed at a certain time when
-   * accessing remote data and overlapping
-   * communication with computation.
+   * Note that the cell iterators in deal.II go through cells differently to
+   * what the cell loop of this class does. This is because several cells are
+   * worked on together (vectorization), and since cells with neighbors on
+   * different MPI processors need to be accessed at a certain time when
+   * accessing remote data and overlapping communication with computation.
    */
   typename MGDoFHandler<dim>::cell_iterator
   get_mg_cell_iterator (const unsigned int macro_cell_number,
@@ -830,21 +649,15 @@ public:
                         const unsigned int fe_component = 0) const;
 
   /**
-   * This returns the cell iterator in deal.II
-   * speak to a given cell in the renumbering of
-   * this structure. This function returns an
-   * exception in case the structure was
-   * not constructed based on an hp::DoFHandler.
+   * This returns the cell iterator in deal.II speak to a given cell in the
+   * renumbering of this structure. This function returns an exception in case
+   * the structure was not constructed based on an hp::DoFHandler.
    *
-   * Note that the cell iterators in deal.II go
-   * through cells differently to what the cell
-   * loop of this class does. This is because
-   * several cells are worked on together
-   * (vectorization), and since cells with
-   * neighbors on different MPI processors need
-   * to be accessed at a certain time when
-   * accessing remote data and overlapping
-   * communication with computation.
+   * Note that the cell iterators in deal.II go through cells differently to
+   * what the cell loop of this class does. This is because several cells are
+   * worked on together (vectorization), and since cells with neighbors on
+   * different MPI processors need to be accessed at a certain time when
+   * accessing remote data and overlapping communication with computation.
    */
   typename hp::DoFHandler<dim>::active_cell_iterator
   get_hp_cell_iterator (const unsigned int macro_cell_number,
@@ -852,75 +665,54 @@ public:
                         const unsigned int fe_component = 0) const;
 
   /**
-   * Since this class uses vectorized
-   * data types with usually more than
-   * one value in the data field, a
-   * situation might occur when some
-   * components of the vector type do
-   * not correspond to an actual cell in
-   * the mesh. When using only this
-   * class, one usually does not need to
-   * bother about that fact since the
-   * values are padded with
-   * zeros. However, when this class is
-   * mixed with deal.II access to cells,
-   * care needs to be taken. This
-   * function returns @p true if not all
-   * @p vectorization_length cells for the given @p
-   * macro_cell are real cells. To find
-   * out how many cells are actually
-   * used, use the function @p
-   * n_components_filled.
+   * Since this class uses vectorized data types with usually more than one
+   * value in the data field, a situation might occur when some components of
+   * the vector type do not correspond to an actual cell in the mesh. When
+   * using only this class, one usually does not need to bother about that
+   * fact since the values are padded with zeros. However, when this class is
+   * mixed with deal.II access to cells, care needs to be taken. This function
+   * returns @p true if not all @p vectorization_length cells for the given @p
+   * macro_cell are real cells. To find out how many cells are actually used,
+   * use the function @p n_components_filled.
    */
   bool
   at_irregular_cell (const unsigned int macro_cell_number) const;
 
   /**
-   * Use this function to find out how
-   * many cells over the length of
-   * vectorization data types correspond
-   * to real cells in the mesh. For most
-   * given @p macro_cells, this is just
-   * @p vectorization_length many, but there might
-   * be one or a few meshes (where the
-   * numbers do not add up) where there
-   * are less such components filled,
-   * indicated by the function @p
+   * Use this function to find out how many cells over the length of
+   * vectorization data types correspond to real cells in the mesh. For most
+   * given @p macro_cells, this is just @p vectorization_length many, but
+   * there might be one or a few meshes (where the numbers do not add up)
+   * where there are less such components filled, indicated by the function @p
    * at_irregular_cell.
    */
   unsigned int
   n_components_filled (const unsigned int macro_cell_number) const;
 
   /**
-   * Returns the number of degrees of
-   * freedom per cell for a given hp
-   * index.
+   * Returns the number of degrees of freedom per cell for a given hp index.
    */
   unsigned int
   get_dofs_per_cell (const unsigned int fe_component = 0,
                      const unsigned int hp_active_fe_index = 0) const;
 
   /**
-   * Returns the number of quadrature
-   * points per cell for a given hp
-   * index.
+   * Returns the number of quadrature points per cell for a given hp index.
    */
   unsigned int
   get_n_q_points (const unsigned int quad_index = 0,
                   const unsigned int hp_active_fe_index = 0) const;
 
   /**
-   * Returns the number of degrees of
-   * freedom on each face of the cell
-   * for given hp index.
+   * Returns the number of degrees of freedom on each face of the cell for
+   * given hp index.
    */
   unsigned int
   get_dofs_per_face (const unsigned int fe_component = 0,
                      const unsigned int hp_active_fe_index = 0) const;
 
   /**
-   * Returns the number of quadrature
-   * points on each face of the cell for
+   * Returns the number of quadrature points on each face of the cell for
    * given hp index.
    */
   unsigned int
@@ -928,54 +720,47 @@ public:
                        const unsigned int hp_active_fe_index = 0) const;
 
   /**
-   * Returns the quadrature rule for
-   * given hp index.
+   * Returns the quadrature rule for given hp index.
    */
   const Quadrature<dim> &
   get_quadrature (const unsigned int quad_index = 0,
                   const unsigned int hp_active_fe_index = 0) const;
 
   /**
-   * Returns the quadrature rule for
-   * given hp index.
+   * Returns the quadrature rule for given hp index.
    */
   const Quadrature<dim-1> &
   get_face_quadrature (const unsigned int quad_index = 0,
                        const unsigned int hp_active_fe_index = 0) const;
 
   /**
-   * Queries whether or not the
-   * indexation has been set.
+   * Queries whether or not the indexation has been set.
    */
   bool indices_initialized () const;
 
   /**
-   * Queries whether or not the
-   * geometry-related information for
-   * the cells has been set.
+   * Queries whether or not the geometry-related information for the cells has
+   * been set.
    */
 
   bool mapping_initialized () const;
 
   /**
-   * Returns an approximation of the memory
-   * consumption of this class in bytes.
+   * Returns an approximation of the memory consumption of this class in
+   * bytes.
    */
   std::size_t memory_consumption() const;
 
   /**
-   * Prints a detailed summary of memory
-   * consumption in the different structures of
-   * this class to the given output stream.
+   * Prints a detailed summary of memory consumption in the different
+   * structures of this class to the given output stream.
    */
   template <typename STREAM>
   void print_memory_consumption(STREAM &out) const;
 
   /**
-   * Prints a summary of this class to the given
-   * output stream. It is focused on the
-   * indices, and does not print all the data
-   * stored.
+   * Prints a summary of this class to the given output stream. It is focused
+   * on the indices, and does not print all the data stored.
    */
   void print (std::ostream &out) const;
 
@@ -998,48 +783,39 @@ public:
   get_size_info () const;
 
   /*
-   * Returns geometry-dependent
-   * information on the cells.
+   * Returns geometry-dependent information on the cells.
    */
   const internal::MatrixFreeFunctions::MappingInfo<dim,Number> &
   get_mapping_info () const;
 
   /**
-   * Returns information on indexation
-   * degrees of freedom.
+   * Returns information on indexation degrees of freedom.
    */
   const internal::MatrixFreeFunctions::DoFInfo &
   get_dof_info (const unsigned int fe_component = 0) const;
 
   /**
-   * Returns the number of weights in the
-   * constraint pool.
+   * Returns the number of weights in the constraint pool.
    */
   unsigned int n_constraint_pool_entries() const;
 
   /**
-   * Returns a pointer to the first
-   * number in the constraint pool data
-   * with index @p pool_index (to
-   * be used together with @p
-   * constraint_pool_end()).
+   * Returns a pointer to the first number in the constraint pool data with
+   * index @p pool_index (to be used together with @p constraint_pool_end()).
    */
   const Number *
   constraint_pool_begin (const unsigned int pool_index) const;
 
   /**
-   * Returns a pointer to one past the
-   * last number in the constraint pool
-   * data with index @p pool_index (to
-   * be used together with @p
+   * Returns a pointer to one past the last number in the constraint pool data
+   * with index @p pool_index (to be used together with @p
    * constraint_pool_begin()).
    */
   const Number *
   constraint_pool_end (const unsigned int pool_index) const;
 
   /**
-   * Returns the unit cell information
-   * for given hp index.
+   * Returns the unit cell information for given hp index.
    */
   const internal::MatrixFreeFunctions::ShapeInfo<Number> &
   get_shape_info (const unsigned int fe_component = 0,
@@ -1052,8 +828,7 @@ public:
 private:
 
   /**
-   * This is the actual reinit function
-   * that sets up the indices for the
+   * This is the actual reinit function that sets up the indices for the
    * DoFHandler and MGDoFHandler case.
    */
   template <typename DoFHandler>
@@ -1065,8 +840,7 @@ private:
                         const AdditionalData               additional_data);
 
   /**
-   * Same as before but for hp::DoFHandler
-   * instead of generic DoFHandler type.
+   * Same as before but for hp::DoFHandler instead of generic DoFHandler type.
    */
   void internal_reinit (const Mapping<dim>               &mapping,
                         const std::vector<const hp::DoFHandler<dim>*> &dof_handler,
@@ -1076,45 +850,37 @@ private:
                         const AdditionalData              additional_data);
 
   /**
-   * Initializes the fields in DoFInfo
-   * together with the constraint pool
-   * that holds all different weights in
-   * the constraints (not part of
-   * DoFInfo because several DoFInfo
-   * classes can have the same weights
-   * which consequently only need to be
-   * stored once).
+   * Initializes the fields in DoFInfo together with the constraint pool that
+   * holds all different weights in the constraints (not part of DoFInfo
+   * because several DoFInfo classes can have the same weights which
+   * consequently only need to be stored once).
    */
   void
   initialize_indices (const std::vector<const ConstraintMatrix *> &constraint,
                       const std::vector<IndexSet> &locally_owned_set);
 
   /**
-   * Initializes the DoFHandlers based on a
-   * DoFHandler<dim> argument.
+   * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
    */
   void initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handlers,
                                 const unsigned int                         level);
 
   /**
-   * Initializes the DoFHandlers based on a
-   * DoFHandler<dim> argument.
+   * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
    */
   void initialize_dof_handlers (const std::vector<const MGDoFHandler<dim>*> &dof_handlers,
                                 const unsigned int                           level);
 
   /**
-   * Initializes the DoFHandlers based on a
-   * DoFHandler<dim> argument.
+   * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
    */
   void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
                                 const unsigned int                             level);
 
   /**
-   * This struct defines which DoFHandler has
-   * actually been given at construction, in
-   * order to define the correct behavior when
-   * querying the underlying DoFHandler.
+   * This struct defines which DoFHandler has actually been given at
+   * construction, in order to define the correct behavior when querying the
+   * underlying DoFHandler.
    */
   struct DoFHandlers
   {
@@ -1128,83 +894,67 @@ private:
   };
 
   /**
-   * Pointers to the DoFHandlers underlying the
-   * current problem.
+   * Pointers to the DoFHandlers underlying the current problem.
    */
   DoFHandlers dof_handlers;
 
   /**
-   * Contains the information about degrees of
-   * freedom on the individual cells and
-   * constraints.
+   * Contains the information about degrees of freedom on the individual cells
+   * and constraints.
    */
   std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
 
   /**
-   * Contains the weights for constraints stored
-   * in DoFInfo. Filled into a separate field
-   * since several vector components might share
-   * similar weights, which reduces memory
-   * consumption. Moreover, it obviates template
-   * arguments on DoFInfo and keeps it a plain
-   * field of indices only.
+   * Contains the weights for constraints stored in DoFInfo. Filled into a
+   * separate field since several vector components might share similar
+   * weights, which reduces memory consumption. Moreover, it obviates template
+   * arguments on DoFInfo and keeps it a plain field of indices only.
    */
   std::vector<Number> constraint_pool_data;
 
   /**
-   * Contains an indicator to the start
-   * of the ith index in the constraint
+   * Contains an indicator to the start of the ith index in the constraint
    * pool data.
    */
   std::vector<unsigned int> constraint_pool_row_index;
 
   /**
-   * Holds information on transformation of
-   * cells from reference cell to real cell that
-   * is needed for evaluating integrals.
+   * Holds information on transformation of cells from reference cell to real
+   * cell that is needed for evaluating integrals.
    */
   internal::MatrixFreeFunctions::MappingInfo<dim,Number> mapping_info;
 
   /**
-   * Contains shape value information on the
-   * unit cell.
+   * Contains shape value information on the unit cell.
    */
   Table<4,internal::MatrixFreeFunctions::ShapeInfo<Number> > shape_info;
 
   /**
-   * Describes how the cells are gone
-   * through. With the cell level (first index
-   * in this field) and the index within the
-   * level, one can reconstruct a deal.II
-   * cell iterator and use all the traditional
-   * things deal.II offers to do with cell
-   * iterators.
+   * Describes how the cells are gone through. With the cell level (first
+   * index in this field) and the index within the level, one can reconstruct
+   * a deal.II cell iterator and use all the traditional things deal.II offers
+   * to do with cell iterators.
    */
   std::vector<std::pair<unsigned int,unsigned int> > cell_level_index;
 
   /**
-   * Stores how many cells we have, how many
-   * cells that we see after applying
-   * vectorization (i.e., the number of macro
-   * cells), and MPI-related stuff.
+   * Stores how many cells we have, how many cells that we see after applying
+   * vectorization (i.e., the number of macro cells), and MPI-related stuff.
    */
   internal::MatrixFreeFunctions::SizeInfo size_info;
 
   /**
-   * Information regarding the shared memory
-   * parallelization.
+   * Information regarding the shared memory parallelization.
    */
   internal::MatrixFreeFunctions::TaskInfo task_info;
 
   /**
-   * Stores whether indices have been
-   * initialized.
+   * Stores whether indices have been initialized.
    */
   bool                               indices_are_initialized;
 
   /**
-   * Stores whether indices have been
-   * initialized.
+   * Stores whether indices have been initialized.
    */
   bool                               mapping_is_initialized;
 };
@@ -1352,7 +1102,8 @@ const Number *
 MatrixFree<dim,Number>::constraint_pool_begin (const unsigned int row) const
 {
   AssertIndexRange (row, constraint_pool_row_index.size()-1);
-  return &constraint_pool_data[0] + constraint_pool_row_index[row];
+  return constraint_pool_data.empty() ? 0 :
+    &constraint_pool_data[0] + constraint_pool_row_index[row];
 }
 
 
@@ -1363,7 +1114,8 @@ const Number *
 MatrixFree<dim,Number>::constraint_pool_end (const unsigned int row) const
 {
   AssertIndexRange (row, constraint_pool_row_index.size()-1);
-  return &constraint_pool_data[0] + constraint_pool_row_index[row+1];
+  return constraint_pool_data.empty() ? 0 :
+    &constraint_pool_data[0] + constraint_pool_row_index[row+1];
 }
 
 
@@ -1411,8 +1163,8 @@ MatrixFree<dim,Number>::create_cell_subrange_hp_by_index
     return range;
   else
     {
-      // the range over which we are searching must be ordered,
-      // otherwise we got a range that spans over too many cells
+      // the range over which we are searching must be ordered, otherwise we
+      // got a range that spans over too many cells
 #ifdef DEBUG
       for (unsigned int i=range.first+1; i<range.second; ++i)
         Assert (fe_indices[i] >= fe_indices[i-1],
@@ -1422,13 +1174,13 @@ MatrixFree<dim,Number>::create_cell_subrange_hp_by_index
 #endif
       std::pair<unsigned int,unsigned int> return_range;
       return_range.first =
-        std::lower_bound(&fe_indices[range.first], &fe_indices[range.second],
-                         fe_index)
+        std::lower_bound(&fe_indices[0] + range.first,
+                         &fe_indices[0] + range.second, fe_index)
         -&fe_indices[0] ;
       return_range.second =
-        std::lower_bound(&fe_indices[return_range.first],
-                         &fe_indices[range.second],
-                         fe_index+1)-&fe_indices[0];
+        std::lower_bound(&fe_indices[0] + return_range.first,
+                         &fe_indices[0] + range.second,
+                         fe_index + 1)-&fe_indices[0];
       Assert(return_range.first >= range.first &&
              return_range.second <= range.second, ExcInternalError());
       return return_range;
@@ -1473,9 +1225,8 @@ MatrixFree<dim,Number>::get_dof_handler (const unsigned int dof_index) const
   else
     {
       Assert (false, ExcNotImplemented());
-      // put pseudo return argument to avoid
-      // compiler error, but trigger a segfault in
-      // case this is only run in optimized mode
+      // put pseudo return argument to avoid compiler error, but trigger a
+      // segfault in case this is only run in optimized mode
       return *dof_handlers.dof_handler[numbers::invalid_unsigned_int];
     }
 }
@@ -1944,8 +1695,7 @@ reinit(const Mapping<dim>                         &mapping,
        const std::vector<Quad>                    &quad,
        const MatrixFree<dim,Number>::AdditionalData additional_data)
 {
-  // find out whether we use a hp Quadrature or
-  // a standard quadrature
+  // find out whether we use a hp Quadrature or a standard quadrature
   std::vector<hp::QCollection<1> > quad_hp;
   for (unsigned int q=0; q<quad.size(); ++q)
     quad_hp.push_back (hp::QCollection<1>(quad[q]));
@@ -1957,13 +1707,10 @@ reinit(const Mapping<dim>                         &mapping,
 
 // ------------------------------ implementation of cell_loop ---------------
 
-// internal helper functions that define how
-// to call MPI data exchange functions: for
-// generic vectors, do nothing at all. For
-// distributed vectors, can call
-// update_ghost_values_start function and so on. If we
-// have collections of vectors, just do the
-// individual functions of the components
+// internal helper functions that define how to call MPI data exchange
+// functions: for generic vectors, do nothing at all. For distributed vectors,
+// can call update_ghost_values_start function and so on. If we have
+// collections of vectors, just do the individual functions of the components
 namespace internal
 {
   template<typename VectorStruct>
@@ -2119,8 +1866,7 @@ namespace internal
 
 #if DEAL_II_USE_MT==1
 
-  // This defines the TBB data structures that
-  // are needed to schedule the
+  // This defines the TBB data structures that are needed to schedule the
   // partition-partition variant
 
   namespace partition
@@ -2390,14 +2136,12 @@ MatrixFree<dim, Number>::cell_loop
 {
 #if DEAL_II_USE_MT==1
 
-  // Use multithreading if so requested and if
-  // there is enough work to do in parallel (the
-  // code might hang if there are less than two
-  // chunks!)
+  // Use multithreading if so requested and if there is enough work to do in
+  // parallel (the code might hang if there are less than two chunks!)
   if (task_info.use_multithreading == true && task_info.n_blocks > 3)
     {
-      // to simplify the function calls, bind away
-      // all arguments except the cell range
+      // to simplify the function calls, bind away all arguments except the
+      // cell range
       typedef
       std_cxx1x::function<void (const std::pair<unsigned int,unsigned int> &range)>
       Worker;
@@ -2489,9 +2233,8 @@ MatrixFree<dim, Number>::cell_loop
           unsigned int evens = task_info.evens;
           unsigned int odds  = task_info.odds;
 
-          // check whether there is only one
-          // partition. if not, build up the tree of
-          // partitions
+          // check whether there is only one partition. if not, build up the
+          // tree of partitions
           if (odds > 0)
             {
               tbb::empty_task *root = new( tbb::task::allocate_root() ) tbb::empty_task;
@@ -2597,9 +2340,8 @@ MatrixFree<dim, Number>::cell_loop
               root->wait_for_all();
               root->destroy(*root);
             }
-          // case when we only have one partition: this
-          // is the usual coloring scheme, and we just
-          // schedule a parallel for loop for each color
+          // case when we only have one partition: this is the usual coloring
+          // scheme, and we just schedule a parallel for loop for each color
           else
             {
               Assert(evens==1,ExcInternalError());
@@ -2629,20 +2371,18 @@ MatrixFree<dim, Number>::cell_loop
 
       internal::update_ghost_values_start (src);
 
-      // First operate on cells where no ghost data is needed (inner
-      // cells)
+      // First operate on cells where no ghost data is needed (inner cells)
       {
         cell_range.first = 0;
         cell_range.second = size_info.boundary_cells_start;
         cell_operation (*this, dst, src, cell_range);
       }
 
-      // before starting operations on cells that contain ghost nodes
-      // (outer cells), wait for the MPI commands to finish
+      // before starting operations on cells that contain ghost nodes (outer
+      // cells), wait for the MPI commands to finish
       internal::update_ghost_values_finish(src);
 
-      // For the outer cells, do the same procedure as for inner
-      // cells.
+      // For the outer cells, do the same procedure as for inner cells.
       if (size_info.boundary_cells_end > size_info.boundary_cells_start)
         {
           cell_range.first = size_info.boundary_cells_start;
@@ -2652,8 +2392,7 @@ MatrixFree<dim, Number>::cell_loop
 
       internal::compress_start(dst);
 
-      // Finally operate on cells where no ghost data is needed (inner
-      // cells)
+      // Finally operate on cells where no ghost data is needed (inner cells)
       if (size_info.n_macro_cells > size_info.boundary_cells_end)
         {
           cell_range.first = size_info.boundary_cells_end;
@@ -2681,8 +2420,7 @@ MatrixFree<dim,Number>::cell_loop
  OutVector      &dst,
  const InVector &src) const
 {
-  // here, use std_cxx1x::bind to hand a
-  // function handler with the appropriate
+  // here, use std_cxx1x::bind to hand a function handler with the appropriate
   // argument to the other loop function
   std_cxx1x::function<void (const MatrixFree<dim,Number> &,
                             OutVector &,
@@ -2714,8 +2452,7 @@ MatrixFree<dim,Number>::cell_loop
  OutVector      &dst,
  const InVector &src) const
 {
-  // here, use std_cxx1x::bind to hand a
-  // function handler with the appropriate
+  // here, use std_cxx1x::bind to hand a function handler with the appropriate
   // argument to the other loop function
   std_cxx1x::function<void (const MatrixFree<dim,Number> &,
                             OutVector &,
index b7b2b4c0d4a661191833f3792a0395d721c3851a..7315162db3e1d560740bb3d240ba7c494f375b87 100644 (file)
@@ -136,9 +136,8 @@ internal_reinit(const Mapping<dim>                         &mapping,
       for (unsigned int no=0; no<dof_handler.size(); ++no)
         dof_info[no].store_plain_indices = additional_data.store_plain_indices;
 
-      // initialize the basic multithreading
-      // information that needs to be passed to the
-      // DoFInfo structure
+      // initialize the basic multithreading information that needs to be
+      // passed to the DoFInfo structure
 #if DEAL_II_USE_MT == 1
       if (additional_data.tasks_parallel_scheme != AdditionalData::none)
         {
@@ -155,18 +154,15 @@ internal_reinit(const Mapping<dim>                         &mapping,
 #endif
         task_info.use_multithreading = false;
 
-      // set dof_indices together with
-      // constraint_indicator and
-      // constraint_pool_data. It also reorders the way
-      // cells are gone through (to separate cells
-      // with overlap to other processors from
-      // others without).
+      // set dof_indices together with constraint_indicator and
+      // constraint_pool_data. It also reorders the way cells are gone through
+      // (to separate cells with overlap to other processors from others
+      // without).
       initialize_indices (constraint, locally_owned_set);
     }
 
-  // Reads out the FE information and stores the
-  // shape function values, gradients and
-  // Hessians for quadrature points.
+  // Reads out the FE information and stores the shape function values,
+  // gradients and Hessians for quadrature points.
   const unsigned int n_fe   = dof_handler.size();
   const unsigned int n_quad = quad.size();
   shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1));
@@ -180,16 +176,12 @@ internal_reinit(const Mapping<dim>                         &mapping,
         }
     }
 
-  // Evaluates transformations from unit to real
-  // cell, Jacobian determinants, quadrature
-  // points in real space, based on the ordering
-  // of the cells determined in @p
-  // extract_local_to_global_indices. The
-  // algorithm assumes that the active FE index
-  // for the transformations is given the active
-  // FE index in the zeroth DoFHandler. TODO:
-  // how do things look like in the more general
-  // case?
+  // Evaluates transformations from unit to real cell, Jacobian determinants,
+  // quadrature points in real space, based on the ordering of the cells
+  // determined in @p extract_local_to_global_indices. The algorithm assumes
+  // that the active FE index for the transformations is given the active FE
+  // index in the zeroth DoFHandler. TODO: how do things look like in the more
+  // general case?
   if (additional_data.initialize_mapping == true)
     {
       mapping_info.initialize (dof_handler[0]->get_tria(), cell_level_index,
@@ -239,9 +231,8 @@ internal_reinit(const Mapping<dim>                            &mapping,
       for (unsigned int no=0; no<dof_handler.size(); ++no)
         dof_info[no].store_plain_indices = additional_data.store_plain_indices;
 
-      // initialize the basic multithreading
-      // information that needs to be passed to the
-      // DoFInfo structure
+      // initialize the basic multithreading information that needs to be
+      // passed to the DoFInfo structure
 #if DEAL_II_USE_MT == 1
       if (additional_data.tasks_parallel_scheme != AdditionalData::none)
         {
@@ -258,18 +249,15 @@ internal_reinit(const Mapping<dim>                            &mapping,
 #endif
         task_info.use_multithreading = false;
 
-      // set dof_indices together with
-      // constraint_indicator and
-      // constraint_pool_data. It also reorders the way
-      // cells are gone through (to separate cells
-      // with overlap to other processors from
-      // others without).
+      // set dof_indices together with constraint_indicator and
+      // constraint_pool_data. It also reorders the way cells are gone through
+      // (to separate cells with overlap to other processors from others
+      // without).
       initialize_indices (constraint, locally_owned_set);
     }
 
-  // Reads out the FE information and stores the
-  // shape function values, gradients and
-  // Hessians for quadrature points.
+  // Reads out the FE information and stores the shape function values,
+  // gradients and Hessians for quadrature points.
   const unsigned int n_components = dof_handler.size();
   const unsigned int n_quad       = quad.size();
   unsigned int n_fe_in_collection = 0;
@@ -292,11 +280,9 @@ internal_reinit(const Mapping<dim>                            &mapping,
                                                  fe.base_element(0));
       }
 
-  // Evaluates transformations from unit to real
-  // cell, Jacobian determinants, quadrature
-  // points in real space, based on the ordering
-  // of the cells determined in @p
-  // extract_local_to_global_indices.
+  // Evaluates transformations from unit to real cell, Jacobian determinants,
+  // quadrature points in real space, based on the ordering of the cells
+  // determined in @p extract_local_to_global_indices.
   if (additional_data.initialize_mapping == true)
     {
       mapping_info.initialize (dof_handler[0]->get_tria(), cell_level_index,
@@ -347,12 +333,9 @@ initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handler,
 
   dof_info.resize (dof_handlers.n_dof_handlers);
 
-  // go through cells on zeroth level and then
-  // successively step down into children. This
-  // gives a z-ordering of the cells, which is
-  // beneficial when setting up neighboring
-  // relations between cells for thread
-  // parallelization
+  // go through cells on zeroth level and then successively step down into
+  // children. This gives a z-ordering of the cells, which is beneficial when
+  // setting up neighboring relations between cells for thread parallelization
   const unsigned int n_mpi_procs = size_info.n_procs;
   const unsigned int my_pid = size_info.my_pid;
 
@@ -383,18 +366,14 @@ initialize_dof_handlers (const std::vector<const MGDoFHandler<dim>*> &dof_handle
 
   dof_info.resize (dof_handlers.n_dof_handlers);
 
-  // go through cells on zeroth level and then
-  // successively step down into children. This
-  // gives a z-ordering of the cells, which is
-  // beneficial when setting up neighboring
-  // relations between cells for thread
-  // parallelization
+  // go through cells on zeroth level and then successively step down into
+  // children. This gives a z-ordering of the cells, which is beneficial when
+  // setting up neighboring relations between cells for thread parallelization
   const unsigned int n_mpi_procs = size_info.n_procs;
   const unsigned int my_pid = size_info.my_pid;
 
-  // if we have no level given, use the same as
-  // for the standard DoFHandler, otherwise we
-  // must loop through the respective level
+  // if we have no level given, use the same as for the standard DoFHandler,
+  // otherwise we must loop through the respective level
   const Triangulation<dim> &tria = dof_handlers.mg_dof_handler[0]->get_tria();
 
   if (level == numbers::invalid_unsigned_int)
@@ -433,18 +412,14 @@ initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_hand
 
   dof_info.resize (dof_handlers.n_dof_handlers);
 
-  // go through cells on zeroth level and then
-  // successively step down into children. This
-  // gives a z-ordering of the cells, which is
-  // beneficial when setting up neighboring
-  // relations between cells for thread
-  // parallelization
+  // go through cells on zeroth level and then successively step down into
+  // children. This gives a z-ordering of the cells, which is beneficial when
+  // setting up neighboring relations between cells for thread parallelization
   const unsigned int n_mpi_procs = size_info.n_procs;
   const unsigned int my_pid = size_info.my_pid;
 
-  // if we have no level given, use the same as
-  // for the standard DoFHandler, otherwise we
-  // must loop through the respective level
+  // if we have no level given, use the same as for the standard DoFHandler,
+  // otherwise we must loop through the respective level
   const Triangulation<dim> &tria = dof_handler[0]->get_tria();
 
   if (n_mpi_procs == 1)
@@ -521,29 +496,24 @@ void MatrixFree<dim,Number>::initialize_indices
                   ExcMessage ("MatrixFree only works for DoFHandler with one base element"));
           const unsigned int n_fe_components = fe.element_multiplicity (0);
 
-          // cache number of finite elements and
-          // dofs_per_cell
+          // cache number of finite elements and dofs_per_cell
           dof_info[no].dofs_per_cell.push_back (fe.dofs_per_cell);
           dof_info[no].dofs_per_face.push_back (fe.dofs_per_face);
           dof_info[no].n_components  = n_fe_components;
 
 
-          // get permutation that gives lexicographic
-          // renumbering of the cell dofs
-          // renumber (this is necessary for FE_Q, for
-          // example, since there the vertex DoFs come
-          // first, which is incompatible with the
-          // lexicographic ordering necessary to apply
-          // tensor products efficiently)
+          // get permutation that gives lexicographic renumbering of the cell
+          // dofs renumber (this is necessary for FE_Q, for example, since
+          // there the vertex DoFs come first, which is incompatible with the
+          // lexicographic ordering necessary to apply tensor products
+          // efficiently)
           const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
             dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>
             (&fe.base_element(0));
 
-          // This class currently only works for
-          // elements derived from
-          // FE_Poly<TensorProductPolynomials<dim>,dim,dim>.
-          // For any other element, the dynamic cast
-          // above will fail and give fe_poly == 0.
+          // This class currently only works for elements derived from
+          // FE_Poly<TensorProductPolynomials<dim>,dim,dim>.  For any other
+          // element, the dynamic cast above will fail and give fe_poly == 0.
           Assert (fe_poly != 0, ExcNotImplemented());
           if (n_fe_components == 1)
             {
@@ -586,8 +556,7 @@ void MatrixFree<dim,Number>::initialize_indices
       dof_info[no].dof_indices.reserve
       ((n_active_cells*dof_info[no].dofs_per_cell[0]*3)/2);
 
-      // cache the constrained indices for use in
-      // matrix-vector products
+      // cache the constrained indices for use in matrix-vector products
       {
         const unsigned int
         start_index = dof_info[no].vector_partitioner->local_range().first,
@@ -601,18 +570,16 @@ void MatrixFree<dim,Number>::initialize_indices
         ghost_dofs[no].reserve (locally_owned_set[no].n_elements()/10+1);
     }
 
-  // extract all the global indices associated
-  // with the computation, and form the ghost
-  // indices
+  // extract all the global indices associated with the computation, and form
+  // the ghost indices
   std::vector<unsigned int> boundary_cells;
   for (unsigned int counter = 0 ; counter < n_active_cells ; ++counter)
     {
       bool cell_at_boundary = false;
       for (unsigned int no=0; no<n_fe; ++no)
         {
-          // OK, read indices from standard DoFHandler
-          // or active indices in MGDoFHandler. That is
-          // the usual stuff
+          // OK, read indices from standard DoFHandler or active indices in
+          // MGDoFHandler. That is the usual stuff
           if (dof_handlers.active_dof_handler == DoFHandlers::usual ||
               (dof_handlers.active_dof_handler == DoFHandlers::multigrid &&
                dof_handlers.level == numbers::invalid_unsigned_int))
@@ -633,8 +600,7 @@ void MatrixFree<dim,Number>::initialize_indices
                                              constraint_values,
                                              cell_at_boundary);
             }
-          // ok, now we are requested to use a level in
-          // a MGDoFHandler
+          // ok, now we are requested to use a level in a MGDoFHandler
           else if (dof_handlers.active_dof_handler == DoFHandlers::multigrid &&
                    dof_handlers.level != numbers::invalid_unsigned_int)
             {
@@ -680,9 +646,8 @@ void MatrixFree<dim,Number>::initialize_indices
             }
         }
 
-      // if we found dofs on some FE component that
-      // belong to other processors, the cell is
-      // added to the boundary cells.
+      // if we found dofs on some FE component that belong to other
+      // processors, the cell is added to the boundary cells.
       if (cell_at_boundary == true)
         boundary_cells.push_back(counter);
     }
@@ -696,11 +661,9 @@ void MatrixFree<dim,Number>::initialize_indices
   for (unsigned int no=0; no<n_fe; ++no)
     dof_info[no].assign_ghosts (boundary_cells);
 
-  // reorganize the indices in order to overlap
-  // communication in MPI with computations:
-  // Place all cells with ghost indices into one
-  // chunk. Also reorder cells so that we can
-  // parallelize by threads
+  // reorganize the indices in order to overlap communication in MPI with
+  // computations: Place all cells with ghost indices into one chunk. Also
+  // reorder cells so that we can parallelize by threads
   std::vector<unsigned int> renumbering;
   if (task_info.use_multithreading == true)
     {
@@ -717,9 +680,8 @@ void MatrixFree<dim,Number>::initialize_indices
     }
   else
     {
-      // In case, we have an hp-dofhandler, we have
-      // to reorder the cell according to the
-      // polynomial degree on the cell.
+      // In case, we have an hp-dofhandler, we have to reorder the cell
+      // according to the polynomial degree on the cell.
       dof_info[0].compute_renumber_serial (boundary_cells, size_info,
                                            renumbering);
       if (dof_handlers.active_dof_handler == DoFHandlers::hp)
@@ -727,11 +689,9 @@ void MatrixFree<dim,Number>::initialize_indices
                                                 irregular_cells);
     }
 
-  // Finally perform the renumbering. We also
-  // want to group several cells together to one
-  // "macro-cell" for vectorization (where the
-  // arithmetic operations will then be done
-  // simultaneously).
+  // Finally perform the renumbering. We also want to group several cells
+  // together to one "macro-cell" for vectorization (where the arithmetic
+  // operations will then be done simultaneously).
 #ifdef DEBUG
   {
     std::vector<unsigned int> sorted_renumbering (renumbering);
@@ -754,13 +714,10 @@ void MatrixFree<dim,Number>::initialize_indices
           cell_level_index.push_back
           (cell_level_index_old[renumbering[position_cell+j]]);
 
-        // generate a cell and level index
-        // also when we have not filled up
-        // vectorization_length cells. This is needed for
-        // MappingInfo when the transformation
-        // data is initialized. We just set
-        // the value to the last valid cell in
-        // that case.
+        // generate a cell and level index also when we have not filled up
+        // vectorization_length cells. This is needed for MappingInfo when the
+        // transformation data is initialized. We just set the value to the
+        // last valid cell in that case.
         for (unsigned int j=n_comp; j<vectorization_length; ++j)
           cell_level_index.push_back
           (cell_level_index_old[renumbering[position_cell+n_comp-1]]);
@@ -770,8 +727,7 @@ void MatrixFree<dim,Number>::initialize_indices
     AssertDimension (cell_level_index.size(),size_info.n_macro_cells*vectorization_length);
   }
 
-  // set constraint pool from the std::map and
-  // reorder the indices
+  // set constraint pool from the std::map and reorder the indices
   typename std::map<std::vector<double>, unsigned int,
            internal::MatrixFreeFunctions::FPArrayComparator<double> >::iterator
            it = constraint_values.constraints.begin(),
@@ -995,13 +951,12 @@ namespace internal
       unsigned int n_max_boundary_cells = boundary_cells.size();
       unsigned int n_boundary_cells = n_max_boundary_cells;
 
-      // try to make the number of boundary cells
-      // divisible by the number of vectors in
-      // vectorization
+      // try to make the number of boundary cells divisible by the number of
+      // vectors in vectorization
+
       /*
-                                  // try to balance the number of cells before
-                                  // and after the boundary part on each
-                                  // processor. probably not worth it!
+      // try to balance the number of cells before and after the boundary part
+      // on each processor. probably not worth it!
       #if DEAL_II_COMPILER_SUPPORTS_MPI
       MPI_Allreduce (&n_boundary_cells, &n_max_boundary_cells, 1, MPI_UNSIGNED,
                      MPI_MAX, size_info.communicator);
@@ -1014,9 +969,8 @@ namespace internal
         (vectorization_length - n_boundary_cells%vectorization_length)%vectorization_length;
       if (fillup_needed > 0 && n_boundary_cells < n_active_cells)
         {
-          // fill additional cells into the list of
-          // boundary cells to get a balanced number. Go
-          // through the indices successively until we
+          // fill additional cells into the list of boundary cells to get a
+          // balanced number. Go through the indices successively until we
           // found enough indices
           std::vector<unsigned int> new_boundary_cells;
           new_boundary_cells.reserve (n_max_boundary_cells);
@@ -1026,8 +980,8 @@ namespace internal
             {
               if (next_free_slot < boundary_cells[bound_index])
                 {
-                  // check if there are enough cells to fill
-                  // with in the current slot
+                  // check if there are enough cells to fill with in the
+                  // current slot
                   if (next_free_slot + fillup_needed <= boundary_cells[bound_index])
                     {
                       for (unsigned int j=boundary_cells[bound_index]-fillup_needed;
@@ -1035,8 +989,8 @@ namespace internal
                         new_boundary_cells.push_back(j);
                       fillup_needed = 0;
                     }
-                  // ok, not enough indices, so just take them
-                  // all up to the next boundary cell
+                  // ok, not enough indices, so just take them all up to the
+                  // next boundary cell
                   else
                     {
                       for (unsigned int j=next_free_slot;
@@ -1062,10 +1016,8 @@ namespace internal
       std::sort (boundary_cells.begin(), boundary_cells.end());
       n_boundary_cells = boundary_cells.size();
 
-      // check that number of boundary cells
-      // is divisible by
-      // vectorization_length or that it
-      // contains all cells
+      // check that number of boundary cells is divisible by
+      // vectorization_length or that it contains all cells
       Assert (n_boundary_cells % vectorization_length == 0 ||
               n_boundary_cells == n_active_cells, ExcInternalError());
       n_macro_cells = (n_active_cells+vectorization_length-1)/vectorization_length;
@@ -1088,6 +1040,8 @@ namespace internal
 
 
 
+    /* ------------------------------------------------------------------ */
+
     template <typename Number>
     FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
       :

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.