From: kronbichler Date: Fri, 28 Dec 2012 10:41:54 +0000 (+0000) Subject: Fix a few out-of-bound array accesses. Refill some comments. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a7a1861460cf965cb58105e8e9184528339fa938;p=dealii-svn.git Fix a few out-of-bound array accesses. Refill some comments. git-svn-id: https://svn.dealii.org/trunk@27865 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/matrix_free/dof_info.h b/deal.II/include/deal.II/matrix_free/dof_info.h index 36b0086691..90c7fccdab 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.h @@ -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 * 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 * 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 &local_indices, const std::vector &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 &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 &boundary_cells, const SizeInfo &size_info, std::vector &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 &renumbering, std::vector &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 &boundary_cells, SizeInfo &size_info, std::vector &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 &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 &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 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 @@ -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 > 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 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 > 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 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 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 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 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 dofs_per_cell; /** - * Stores the number of degrees of freedom per - * face. + * Stores the number of degrees of freedom per face. */ std::vector 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 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 > 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 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; } } diff --git a/deal.II/include/deal.II/matrix_free/dof_info.templates.h b/deal.II/include/deal.II/matrix_free/dof_info.templates.h index d38900c3bf..bf4003fedd 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.templates.h @@ -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 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 > &entries); @@ -89,27 +80,22 @@ namespace internal ConstraintComparator()); for (unsigned int j=0; j, unsigned int, FPArrayComparator >::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(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 > &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 &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 > ghost_origin(n_ghosts); for (unsigned int i=0; i 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(vector_partitioner.get()); vec_part->set_ghost_indices (ghost_indices); @@ -617,8 +589,7 @@ no_constraint: for (unsigned int j=0; j glob_indices (vectorization_length); std::vector plain_glob_indices (vectorization_length); std::vector*> @@ -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; ilocal_range().second- vector_partitioner->local_range().first); for (unsigned int row=0; row cell_partition(task_info.n_blocks, size_info.n_macro_cells); std::vector neighbor_list; std::vector 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 partition_list (task_info.n_blocks,0); std::vector 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 partition_blocks (2,0); std::vector cell_color(task_info.n_blocks, size_info.n_macro_cells); std::vector 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 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 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 block_start(size_info.n_macro_cells+1); std::vector 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 cell_partition (size_info.n_active_cells, size_info.n_active_cells); std::vector neighbor_list; std::vector 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 partition_list(size_info.n_active_cells,0); std::vector 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 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 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 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; j0 || 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 { 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::iterator &dat) { @@ -1737,12 +1664,15 @@ not_connect: if (dat == end()) { push_back(entry); + dat = end(); } else if (*dat > entry) { - this->std::vector::insert (dat, entry); + dat = this->std::vector::insert (dat, entry); + ++dat; } - ++dat; + else + ++dat; } }; } @@ -1926,8 +1856,7 @@ not_connect: void DoFInfo::renumber_dofs (std::vector &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 new_constrained_dofs[i-1], ExcInternalError()); diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.h b/deal.II/include/deal.II/matrix_free/mapping_info.h index 6a46b630ff..481fe26814 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.h @@ -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< &tria, const std::vector > &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 > &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 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 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 >, VectorizedArray > > 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 >, VectorizedArray > > 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 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 > > 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 > 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 > > 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. */ AlignedVector1?dim*(dim-1)/2:1), Tensor<1,dim,VectorizedArray > > > 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 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 > > 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 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 face_quadrature; /** - * The number of quadrature points for the - * current quadrature formula. + * The number of quadrature points for the current quadrature formula. */ std::vector 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 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 > > 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 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 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 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 &tria, const std::pair *cells, diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h index 493e39b124..647ed1476c 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h @@ -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 > &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 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 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 > jac_grad, grad_jac_inv; Tensor<1,dim,VectorizedArray > 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::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::active_cell_iterator, we need to manually + // select the correct finite element, so just hold a vector of // FEValues std::vector > > 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 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 &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 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,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,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 > &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::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 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 &jac = fe_val.jacobian(q); for (unsigned int d=0; d &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 void reinit (const Mapping &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 void reinit (const Mapping &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 @@ -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 void reinit (const Mapping &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 @@ -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 @@ -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 void reinit (const Mapping &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 @@ -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 &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: cell_operation - * (const MatrixFree - * &, OutVector &, InVector &, - * std::pair &), 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 - * operator() 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: + * cell_operation (const MatrixFree &, OutVector &, + * InVector &, std::pair &), 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 operator() + * with the correct set of arguments since such a pointer can be converted + * to the function object. */ template void cell_loop (const std_cxx1x::function &, @@ -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 - * cell_operation (const - * MatrixFree &, OutVector &, - * InVector &, std::pair&)const. 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 cell_operation (const MatrixFree &, OutVector + * &, InVector &, std::pair&)const. 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 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 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 create_cell_subrange_hp (const std::pair &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 create_cell_subrange_hp_by_index (const std::pair &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@, 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@, 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 @@ -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@, 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@, 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 @@ -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 & 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 & 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 &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 & 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 & 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::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::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::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 & 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 & 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 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 & 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 & 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 @@ -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 &mapping, const std::vector*> &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 &constraint, const std::vector &locally_owned_set); /** - * Initializes the DoFHandlers based on a - * DoFHandler argument. + * Initializes the DoFHandlers based on a DoFHandler argument. */ void initialize_dof_handlers (const std::vector*> &dof_handlers, const unsigned int level); /** - * Initializes the DoFHandlers based on a - * DoFHandler argument. + * Initializes the DoFHandlers based on a DoFHandler argument. */ void initialize_dof_handlers (const std::vector*> &dof_handlers, const unsigned int level); /** - * Initializes the DoFHandlers based on a - * DoFHandler argument. + * Initializes the DoFHandlers based on a DoFHandler argument. */ void initialize_dof_handlers (const std::vector*> &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 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 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 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 mapping_info; /** - * Contains shape value information on the - * unit cell. + * Contains shape value information on the unit cell. */ Table<4,internal::MatrixFreeFunctions::ShapeInfo > 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 > 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::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::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::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= fe_indices[i-1], @@ -1422,13 +1174,13 @@ MatrixFree::create_cell_subrange_hp_by_index #endif std::pair 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::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 &mapping, const std::vector &quad, const MatrixFree::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 > quad_hp; for (unsigned int q=0; q(quad[q])); @@ -1957,13 +1707,10 @@ reinit(const Mapping &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 @@ -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::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 &range)> Worker; @@ -2489,9 +2233,8 @@ MatrixFree::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::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::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::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::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 &, OutVector &, @@ -2714,8 +2452,7 @@ MatrixFree::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 &, OutVector &, diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h index b7b2b4c0d4..7315162db3 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h @@ -136,9 +136,8 @@ internal_reinit(const Mapping &mapping, for (unsigned int no=0; no &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 &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 &mapping, for (unsigned int no=0; no &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 &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*> &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*> &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 &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*> &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 &tria = dof_handler[0]->get_tria(); if (n_mpi_procs == 1) @@ -521,29 +496,24 @@ void MatrixFree::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,dim,dim> *fe_poly = dynamic_cast,dim,dim>*> (&fe.base_element(0)); - // This class currently only works for - // elements derived from - // FE_Poly,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,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::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::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 boundary_cells; for (unsigned int counter = 0 ; counter < n_active_cells ; ++counter) { bool cell_at_boundary = false; for (unsigned int no=0; no::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::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::initialize_indices for (unsigned int no=0; no renumbering; if (task_info.use_multithreading == true) { @@ -717,9 +680,8 @@ void MatrixFree::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::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 sorted_renumbering (renumbering); @@ -754,13 +714,10 @@ void MatrixFree::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::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, unsigned int, internal::MatrixFreeFunctions::FPArrayComparator >::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 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 FPArrayComparator::FPArrayComparator (const Number scaling) :