// $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
* 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
*/
void clear ();
/**
- * Returns a pointer to the first index in the
- * DoF row @p row.
+ * Returns a pointer to the first index in the DoF row @p row.
*/
const unsigned int *begin_indices (const unsigned int row) const;
/**
- * Returns a pointer to the one past the last
- * DoF index in the row @p row.
+ * Returns a pointer to the one past the last DoF index in the row @p
+ * row.
*/
const unsigned int *end_indices (const unsigned int row) const;
/**
- * Returns the number of entries in the
- * indices field for the given row.
+ * Returns the number of entries in the indices field for the given row.
*/
unsigned int row_length_indices (const unsigned int row) const;
/**
- * Returns a pointer to the first constraint
- * indicator in the row @p row.
+ * Returns a pointer to the first constraint indicator in the row @p
+ * row.
*/
const std::pair<unsigned short,unsigned short> *
begin_indicators (const unsigned int row) const;
/**
- * Returns a pointer to the one past the last
- * constraint indicator in the row @p row.
+ * Returns a pointer to the one past the last constraint indicator in
+ * the row @p row.
*/
const std::pair<unsigned short,unsigned short> *
end_indicators (const unsigned int row) const;
/**
- * Returns the number of entries in the
- * constraint indicator field for the given
- * row.
+ * Returns the number of entries in the constraint indicator field for
+ * the given row.
*/
unsigned int row_length_indicators (const unsigned int row) const;
/**
- * Returns a pointer to the first index in the
- * DoF row @p row for plain indices (i.e., the
- * entries where constraints are not
- * embedded).
+ * Returns a pointer to the first index in the DoF row @p row for plain
+ * indices (i.e., the entries where constraints are not embedded).
*/
const unsigned int *begin_indices_plain (const unsigned int row) const;
/**
- * Returns a pointer to the one past the last
- * DoF index in the row @p row (i.e., the
- * entries where constraints are not
- * embedded).
+ * Returns a pointer to the one past the last DoF index in the row @p
+ * row (i.e., the entries where constraints are not embedded).
*/
const unsigned int *end_indices_plain (const unsigned int row) const;
/**
- * Returns the FE index for a given finite
- * element degree. If not in hp mode, this
- * function always returns index 0. If an
- * index is not found in hp mode, it returns
- * max_fe_degree, i.e., one index past the
- * last valid one.
+ * Returns the FE index for a given finite element degree. If not in hp
+ * mode, this function always returns index 0. If an index is not found
+ * in hp mode, it returns max_fe_degree, i.e., one index past the last
+ * valid one.
*/
unsigned int fe_index_from_degree (const unsigned int fe_degree) const;
/**
- * Returns the FE index for a given finite
- * element degree. If not in hp mode or if the
- * index is not found, this function always
- * returns index 0. Hence, this function does
- * not check whether the given degree is
+ * Returns the FE index for a given finite element degree. If not in hp
+ * mode or if the index is not found, this function always returns index
+ * 0. Hence, this function does not check whether the given degree is
* actually present.
*/
unsigned int
fe_index_from_dofs_per_cell (const unsigned int dofs_per_cell) const;
/**
- * This internal method takes the local
- * indices on a cell and fills them into this
- * class. It resolves the constraints and
- * distributes the results. Ghost indices,
- * i.e., indices that are located on another
- * processor, get a temporary number by this
- * function, and will later be assigned the
- * correct index after all the ghost indices
- * have been collected by the call to @p
- * assign_ghosts.
+ * This internal method takes the local indices on a cell and fills them
+ * into this class. It resolves the constraints and distributes the
+ * results. Ghost indices, i.e., indices that are located on another
+ * processor, get a temporary number by this function, and will later be
+ * assigned the correct index after all the ghost indices have been
+ * collected by the call to @p assign_ghosts.
*/
void read_dof_indices (const std::vector<unsigned int> &local_indices,
const std::vector<unsigned int> &lexicographic_inv,
bool &cell_at_boundary);
/**
- * This method assigns the correct indices to
- * ghost indices from the temporary numbering
- * employed by the @p read_dof_indices
- * function. The numbers are localized with
- * respect to the MPI process, and ghosts
- * start at the end of the locally owned
- * range. This way, we get direct access to
- * all vector entries.
+ * This method assigns the correct indices to ghost indices from the
+ * temporary numbering employed by the @p read_dof_indices function. The
+ * numbers are localized with respect to the MPI process, and ghosts
+ * start at the end of the locally owned range. This way, we get direct
+ * access to all vector entries.
*/
void assign_ghosts(const std::vector<unsigned int> &boundary_cells);
/**
- * Reorganizes cells for serial
- * (non-thread-parallelized) such that
- * boundary cells are places in the
- * middle. This way, computations and
- * communication can be overlapped. Should
- * only be called by one DoFInfo object when
- * used on a system of several DoFHandlers.
+ * Reorganizes cells for serial (non-thread-parallelized) such that
+ * boundary cells are places in the middle. This way, computations and
+ * communication can be overlapped. Should only be called by one DoFInfo
+ * object when used on a system of several DoFHandlers.
*/
void compute_renumber_serial (const std::vector<unsigned int> &boundary_cells,
const SizeInfo &size_info,
std::vector<unsigned int> &renumbering);
/**
- * Reorganizes cells in the hp case without
- * parallelism such that all cells with the
- * same FE index are placed
- * consecutively. Should only be called by one
- * DoFInfo object when used on a system of
- * several DoFHandlers.
+ * Reorganizes cells in the hp case without parallelism such that all
+ * cells with the same FE index are placed consecutively. Should only be
+ * called by one DoFInfo object when used on a system of several
+ * DoFHandlers.
*/
void compute_renumber_hp_serial (SizeInfo &size_info,
std::vector<unsigned int> &renumbering,
std::vector<unsigned int> &irregular_cells);
/**
- * Computes the initial renumbering of cells
- * such that all cells with ghosts are put
- * first. This is the first step before
- * building the thread graph and used to
- * overlap computations and communication.
+ * Computes the initial renumbering of cells such that all cells with
+ * ghosts are put first. This is the first step before building the
+ * thread graph and used to overlap computations and communication.
*/
void compute_renumber_parallel (const std::vector<unsigned int> &boundary_cells,
SizeInfo &size_info,
std::vector<unsigned int> &renumbering);
/**
- * This method reorders the way cells are gone
- * through based on a given renumbering of the
- * cells. It also takes @p
- * vectorization_length cells together and
- * interprets them as one cell only, as is
- * needed for vectorization.
+ * This method reorders the way cells are gone through based on a given
+ * renumbering of the cells. It also takes @p vectorization_length cells
+ * together and interprets them as one cell only, as is needed for
+ * vectorization.
*/
void reorder_cells (const SizeInfo &size_info,
const std::vector<unsigned int> &renumbering,
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,
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,
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,
CompressedSimpleSparsityPattern &connectivity) const;
/**
- * Renumbers the degrees of freedom to give
- * good access for this class.
+ * Renumbers the degrees of freedom to give good access for this class.
*/
void renumber_dofs (std::vector<unsigned int> &renumbering);
/**
- * Returns the memory consumption in bytes of
- * this class.
+ * Returns the memory consumption in bytes of this class.
*/
std::size_t memory_consumption() const;
/**
- * Prints a detailed summary of memory
- * consumption in the different structures of
- * this class to the given output stream.
+ * Prints a detailed summary of memory consumption in the different
+ * structures of this class to the given output stream.
*/
template <typename STREAM>
void print_memory_consumption(STREAM &out,
const SizeInfo &size_info) const;
/**
- * Prints a representation of the
- * indices in the class to the given
+ * Prints a representation of the indices in the class to the given
* output stream.
*/
template <typename Number>
std::ostream &out) const;
/**
- * Stores the rowstart indices of the
- * compressed row storage in the @p
- * dof_indices and @p constraint_indicator
- * fields. These two fields are always
- * accessed together, so it is simpler to keep
- * just one variable for them. This also
- * obviates keeping two rowstart vectors in
- * synch.
+ * Stores the rowstart indices of the compressed row storage in the @p
+ * dof_indices and @p constraint_indicator fields. These two fields are
+ * always accessed together, so it is simpler to keep just one variable
+ * for them. This also obviates keeping two rowstart vectors in synch.
*
- * In addition, the third field stores whether
- * a particular cell has a certain structure
- * in the indices, like indices for
- * vector-valued problems or for cells where
- * not all vector components are filled.
+ * In addition, the third field stores whether a particular cell has a
+ * certain structure in the indices, like indices for vector-valued
+ * problems or for cells where not all vector components are filled.
*/
std::vector<std_cxx1x::tuple<unsigned int,
unsigned int,
unsigned int> > row_starts;
/**
- * Stores the indices of the degrees of
- * freedom for each cell. This array also
- * includes the indirect contributions from
- * constraints, which are described by the @p
- * constraint_indicator field. Because of
- * variable lengths of rows, this would be a
- * vector of a vector. However, we use one
- * contiguous memory region and store the
+ * Stores the indices of the degrees of freedom for each cell. This
+ * array also includes the indirect contributions from constraints,
+ * which are described by the @p constraint_indicator field. Because of
+ * variable lengths of rows, this would be a vector of a
+ * vector. However, we use one contiguous memory region and store the
* rowstart in the variable @p row_starts.
*/
std::vector<unsigned int> dof_indices;
/**
- * This variable describes the position of
- * constraints in terms of the local numbering
- * of degrees of freedom on a cell. The first
- * number stores the distance from one
- * constrained degree of freedom to the
- * next. This allows to identify the position
- * of constrained DoFs as we loop through the
- * local degrees of freedom of the cell when
- * reading from or writing to a vector. The
- * second number stores the index of the
- * constraint weights, stored in the variable
- * constraint_pool_data.
+ * This variable describes the position of constraints in terms of the
+ * local numbering of degrees of freedom on a cell. The first number
+ * stores the distance from one constrained degree of freedom to the
+ * next. This allows to identify the position of constrained DoFs as we
+ * loop through the local degrees of freedom of the cell when reading
+ * from or writing to a vector. The second number stores the index of
+ * the constraint weights, stored in the variable constraint_pool_data.
*/
std::vector<std::pair<unsigned short,unsigned short> > constraint_indicator;
/**
- * This stores the parallel partitioning that
- * can be used to set up vectors. The
- * partitioner includes the description of the
- * local range in the vector, and also
- * includes how the ghosts look like. This
- * enables initialization of vectors based on
- * the DoFInfo field.
+ * This stores the parallel partitioning that can be used to set up
+ * vectors. The partitioner includes the description of the local range
+ * in the vector, and also includes how the ghosts look like. This
+ * enables initialization of vectors based on the DoFInfo field.
*/
std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> vector_partitioner;
/**
- * This stores a (sorted) list of all locally
- * owned degrees of freedom that are
- * constrained.
+ * This stores a (sorted) list of all locally owned degrees of freedom
+ * that are constrained.
*/
std::vector<unsigned int> constrained_dofs;
/**
- * Stores the rowstart indices of the
- * compressed row storage in the @p
- * dof_indices_plain
- * fields.
+ * Stores the rowstart indices of the compressed row storage in the @p
+ * dof_indices_plain fields.
*/
std::vector<unsigned int> row_starts_plain_indices;
/**
- * Stores the indices of the degrees of
- * freedom for each cell. This array does not
- * include the indirect contributions from
- * constraints, which are included in @p
- * dof_indices. Because of variable lengths of
- * rows, this would be a vector of a
- * vector. However, we use one contiguous
- * memory region and store the rowstart in the
- * variable @p row_starts_plain_indices.
+ * Stores the indices of the degrees of freedom for each cell. This
+ * array does not include the indirect contributions from constraints,
+ * which are included in @p dof_indices. Because of variable lengths of
+ * rows, this would be a vector of a vector. However, we use one
+ * contiguous memory region and store the rowstart in the variable @p
+ * row_starts_plain_indices.
*/
std::vector<unsigned int> plain_dof_indices;
/**
- * Stores the number of components in the
- * DoFHandler where the indices have been read
- * from.
+ * Stores the number of components in the DoFHandler where the indices
+ * have been read from.
*/
unsigned int n_components;
/**
- * Stores the number of degrees of freedom per
- * cell.
+ * Stores the number of degrees of freedom per cell.
*/
std::vector<unsigned int> dofs_per_cell;
/**
- * Stores the number of degrees of freedom per
- * face.
+ * Stores the number of degrees of freedom per face.
*/
std::vector<unsigned int> dofs_per_face;
/**
- * Informs on whether plain indices are
- * cached.
+ * Informs on whether plain indices are cached.
*/
bool store_plain_indices;
/**
- * Stores the index of the active
- * finite element in the hp case.
+ * Stores the index of the active finite element in the hp case.
*/
std::vector<unsigned int> cell_active_fe_index;
/**
- * Stores the maximum degree of
- * different finite elements for the
- * hp case.
+ * Stores the maximum degree of different finite elements for the hp
+ * case.
*/
unsigned int max_fe_index;
/**
- * This variable stores the dofs per cell and
- * the finite element degree associated for
- * all fe indices in the underlying element
- * for easier access to data in the hp case.
+ * This variable stores the dofs per cell and the finite element degree
+ * associated for all fe indices in the underlying element for easier
+ * access to data in the hp case.
*/
std::vector<std::pair<unsigned int,unsigned int> > fe_index_conversion;
/**
- * Temporarily stores the numbers of ghosts
- * during setup. Cleared when calling @p
- * assign_ghosts. Then, all information is
- * collected by the partitioner.
+ * Temporarily stores the numbers of ghosts during setup. Cleared when
+ * calling @p assign_ghosts. Then, all information is collected by the
+ * partitioner.
*/
std::vector<unsigned int> ghost_dofs;
};
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;
}
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;
}
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;
}
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;
}
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,
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;
}
}
};
/**
- * A struct that takes entries describing
- * a constraint and puts them into a
- * sorted list where duplicates are
- * filtered out
+ * A struct that takes entries describing a constraint and puts them into
+ * a sorted list where duplicates are filtered out
*/
template <typename Number>
struct ConstraintValues
ConstraintValues();
/**
- * This function inserts some
- * constrained entries to the
- * collection of all
- * values. It stores the
- * (reordered) numbering of
- * the dofs (according to the
- * ordering that matches with
- * the function) in
- * new_indices, and returns
- * the storage position the
- * double array for access
- * later on.
+ * This function inserts some constrained entries to the collection of
+ * all values. It stores the (reordered) numbering of the dofs
+ * (according to the ordering that matches with the function) in
+ * new_indices, and returns the storage position the double array for
+ * access later on.
*/
unsigned short
insert_entries (const std::vector<std::pair<unsigned int,double> > &entries);
ConstraintComparator());
for (unsigned int j=0; j<constraint_entries.size(); j++)
{
- // copy the indices of the constraint entries after
- // sorting.
+ // copy the indices of the constraint entries after sorting.
constraint_indices[j] = constraint_entries[j].first;
- // one_constraint takes the weights of the
- // constraint
+ // one_constraint takes the weights of the constraint
next_constraint.first[j] = constraint_entries[j].second;
}
}
next_constraint.second = constraints.size();
- // check whether or not constraint is already
- // in pool. the initial implementation
- // computed a hash value based on the truncated
- // array (to given accuracy around 1e-13) in
- // order to easily detect different arrays and
- // then made a fine-grained check when the
- // hash values were equal. this was quite
- // lenghty and now we use a std::map with a
- // user-defined comparator to compare floating
- // point arrays to a tolerance 1e-13.
+ // check whether or not constraint is already in pool. the initial
+ // implementation computed a hash value based on the truncated array (to
+ // given accuracy around 1e-13) in order to easily detect different
+ // arrays and then made a fine-grained check when the hash values were
+ // equal. this was quite lenghty and now we use a std::map with a
+ // user-defined comparator to compare floating point arrays to a
+ // tolerance 1e-13.
std::pair<typename std::map<std::vector<double>, unsigned int,
FPArrayComparator<double> >::iterator,
bool> it = constraints.insert(next_constraint);
else
insert_position = next_constraint.second;
- // we want to store the result as a short
- // variable, so we have to make sure that the
- // result does not exceed the limits when
- // casting.
+ // we want to store the result as a short variable, so we have to make
+ // sure that the result does not exceed the limits when casting.
Assert(insert_position < (1<<(8*sizeof(unsigned short))),
ExcInternalError());
return static_cast<unsigned short>(insert_position);
// dof is constrained
if (entries_ptr != 0)
{
- // in case we want to access plain indices, we
- // need to know about the location of
- // constrained indices as well (all the other
- // indices are collected by the cases below)
+ // in case we want to access plain indices, we need to know
+ // about the location of constrained indices as well (all the
+ // other indices are collected by the cases below)
if (current_dof < first_owned || current_dof >= last_owned)
{
ghost_dofs.push_back (current_dof);
cell_at_boundary = true;
}
- // check whether this dof is identity
- // constrained to another dof. then we can
- // simply insert that dof and there is no need
+ // check whether this dof is identity constrained to another
+ // dof. then we can simply insert that dof and there is no need
// to actually resolve the constraint entries
const std::vector<std::pair<unsigned int,double> >
&entries = *entries_ptr;
// reset constraint iterator for next round
constraint_iterator.first = 0;
- // add the local_to_global indices computed in
- // the insert_entries function. transform the
- // index to local index space or mark it as
- // ghost if necessary
+ // add the local_to_global indices computed in the
+ // insert_entries function. transform the index to local index
+ // space or mark it as ghost if necessary
if (n_entries > 0)
{
const std::vector<unsigned int> &constraint_indices =
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);
}
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))
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());
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)
IndexSet ghost_indices (vector_partitioner->size());
if (n_ghosts > 0)
{
- // since we need to go back to the
- // local_to_global indices and replace the
- // temporary numbering of ghosts by the real
- // number in the index set, we need to store
- // these values
+ // since we need to go back to the local_to_global indices and
+ // replace the temporary numbering of ghosts by the real number in
+ // the index set, we need to store these values
std::vector<std::pair<unsigned int,unsigned int> > ghost_origin(n_ghosts);
for (unsigned int i=0; i<n_ghosts; ++i)
{
ghost_origin.back().first+1);
ghost_indices.compress();
- // make sure that we got the correct local
- // numbering of the ghost dofs. the ghost
- // index set should store the same number
+ // make sure that we got the correct local numbering of the ghost
+ // dofs. the ghost index set should store the same number
{
AssertDimension (n_unique_ghosts, ghost_indices.n_elements());
for (unsigned int i=0; i<n_ghosts; ++i)
ExcInternalError());
}
- // apply correct numbering for ghost indices:
- // We previously just enumerated them
- // according to their appearance in the
- // local_to_global structure. Above, we
- // derived a relation between this enumeration
- // and the actual number
+ // apply correct numbering for ghost indices: We previously just
+ // enumerated them according to their appearance in the
+ // local_to_global structure. Above, we derived a relation between
+ // this enumeration and the actual number
const unsigned int n_boundary_cells = boundary_cells.size();
for (unsigned int i=0; i<n_boundary_cells; ++i)
{
std::vector<unsigned int> new_ghosts;
ghost_dofs.swap(new_ghosts);
- // set the ghost indices now. need to cast
- // away constness here, but that is uncritical
- // since we reset the Partitioner in the same
- // initialize call as this call here.
+ // set the ghost indices now. need to cast away constness here, but that
+ // is uncritical since we reset the Partitioner in the same initialize
+ // call as this call here.
Utilities::MPI::Partitioner *vec_part =
const_cast<Utilities::MPI::Partitioner *>(vector_partitioner.get());
vec_part->set_ghost_indices (ghost_indices);
for (unsigned int j=0; j<n_comp; ++j)
fe_indices[j]=cell_active_fe_index[renumbering[position_cell+j]];
- // by construction, all cells should
- // have the same fe index.
+ // by construction, all cells should have the same fe index.
for (unsigned int j=1; j<n_comp; ++j)
Assert (fe_indices[j] == fe_indices[0], ExcInternalError());
new_plain_indices.reserve (plain_dof_indices.size());
}
- // copy the indices and the constraint
- // indicators to the new data field: Store the
- // indices in a way so that adjacent data
- // fields in local vectors are adjacent, i.e.,
- // first dof index 0 for all vectors, then dof
- // index 1 for all vectors, and so on. This
- // involves some extra resorting.
+ // copy the indices and the constraint indicators to the new data field:
+ // Store the indices in a way so that adjacent data fields in local
+ // vectors are adjacent, i.e., first dof index 0 for all vectors, then
+ // dof index 1 for all vectors, and so on. This involves some extra
+ // resorting.
std::vector<const unsigned int *> glob_indices (vectorization_length);
std::vector<const unsigned int *> plain_glob_indices (vectorization_length);
std::vector<const std::pair<unsigned short,unsigned short>*>
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)
{
new_rowstart_plain.swap (row_starts_plain_indices);
#ifdef DEBUG
- // sanity check 1: all indices should be
- // smaller than the number of dofs locally
- // owned plus the number of ghosts
+ // sanity check 1: all indices should be smaller than the number of dofs
+ // locally owned plus the number of ghosts
const unsigned int index_range = (vector_partitioner->local_range().second-
vector_partitioner->local_range().first)
+ vector_partitioner->ghost_indices().n_elements();
for (unsigned int i=0; i<dof_indices.size(); ++i)
AssertIndexRange (dof_indices[i], index_range);
- // sanity check 2: for the constraint
- // indicators, the first index should be
- // smaller than the number of indices in the
- // row, and the second index should be smaller
- // than the number of constraints in the
+ // sanity check 2: for the constraint indicators, the first index should
+ // be smaller than the number of indices in the row, and the second
+ // index should be smaller than the number of constraints in the
// constraint pool.
for (unsigned int row=0; row<size_info.n_macro_cells; ++row)
{
}
}
- // sanity check 3: all non-boundary cells
- // should have indices that only refer to the
- // locally owned range
+ // sanity check 3: all non-boundary cells should have indices that only
+ // refer to the locally owned range
const unsigned int local_size = (vector_partitioner->local_range().second-
vector_partitioner->local_range().first);
for (unsigned int row=0; row<size_info.boundary_cells_start; ++row)
void DoFInfo::guess_block_size (const SizeInfo &size_info,
TaskInfo &task_info)
{
- // user did not say a positive number, so we
- // have to guess
+ // user did not say a positive number, so we have to guess
if (task_info.block_size == 0)
{
- // we would like to have enough work to do, so
- // as first guess, try to get 50 times as many
- // chunks as we have threads on the
- // system.
+ // we would like to have enough work to do, so as first guess, try
+ // to get 50 times as many chunks as we have threads on the system.
task_info.block_size =
size_info.n_macro_cells / (multithread_info.n_default_threads * 50);
- // if there are too few degrees of freedom per
- // cell, need to increase the block size
+ // if there are too few degrees of freedom per cell, need to
+ // increase the block size
const unsigned int minimum_parallel_grain_size = 500;
if (dofs_per_cell[0] * task_info.block_size <
minimum_parallel_grain_size)
guess_block_size (size_info, task_info);
- // set up partitions. if we just use coloring
- // without partitions, do nothing here, assume
- // all cells to belong to the zero partition
- // (that we otherwise use for MPI boundary
- // cells)
+ // set up partitions. if we just use coloring without partitions, do
+ // nothing here, assume all cells to belong to the zero partition (that
+ // we otherwise use for MPI boundary cells)
unsigned int partition = 0, start_up = 0, counter = 0;
unsigned int start_nonboundary = numbers::invalid_unsigned_int;
bool work = true;
task_info.block_size_last = task_info.block_size;
}
- // assume that all FEs have the same
- // connectivity graph, so take the zeroth FE
+ // assume that all FEs have the same connectivity graph, so take the
+ // zeroth FE
task_info.n_blocks = (size_info.n_macro_cells+task_info.block_size-1)/
task_info.block_size;
task_info.block_size_last = size_info.n_macro_cells-
(task_info.block_size*(task_info.n_blocks-1));
- // create the connectivity graph with internal
- // blocking
+ // create the connectivity graph with internal blocking
CompressedSimpleSparsityPattern connectivity;
make_connectivity_graph (size_info, task_info, renumbering,irregular_cells,
true, connectivity);
// Create cell-block partitioning.
- // For each block of cells, this variable
- // saves to which partitions the block
- // belongs. Initialize all to n_macro_cells to
- // mark them as not yet assigned a partition.
+ // For each block of cells, this variable saves to which partitions the
+ // block belongs. Initialize all to n_macro_cells to mark them as not
+ // yet assigned a partition.
std::vector<unsigned int> cell_partition(task_info.n_blocks,
size_info.n_macro_cells);
std::vector<unsigned int> neighbor_list;
std::vector<unsigned int> neighbor_neighbor_list;
- // In element j of this variable, one puts the
- // old number of the block that should be the
- // jth block in the new numeration.
+ // In element j of this variable, one puts the old number of the block
+ // that should be the jth block in the new numeration.
std::vector<unsigned int> partition_list (task_info.n_blocks,0);
std::vector<unsigned int> partition_color_list(task_info.n_blocks,0);
- // This vector points to the start of each
- // partition.
+ // This vector points to the start of each partition.
std::vector<unsigned int> partition_blocks (2,0);
std::vector<unsigned int> cell_color(task_info.n_blocks,
size_info.n_macro_cells);
std::vector<bool> color_finder;
- // this performs a classical breath-first
- // search in the connectivity graph of the
- // cell chunks
+ // this performs a classical breath-first search in the connectivity
+ // graph of the cell chunks
while (work)
{
- // put all cells up to begin_inner_cells into
- // first partition. if the numbers do not add
- // up exactly, assign an additional block
+ // put all cells up to begin_inner_cells into first partition. if
+ // the numbers do not add up exactly, assign an additional block
if (start_nonboundary>0 && start_up == start_nonboundary)
{
unsigned int n_blocks = ((start_nonboundary+task_info.block_size-1)
}
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);
neighbor_neighbor_list.resize(0);
}
- // One has to check if the graph is not
- // connected so we have to find another
- // partition.
+ // One has to check if the graph is not connected so we have to find
+ // another partition.
work = false;
for (unsigned int j=start_up; j<task_info.n_blocks; ++j)
if (cell_partition[j] == size_info.n_macro_cells)
unsigned int cell = partition_list[k];
unsigned int no_neighbors = connectivity.row_length(cell);
- // In the worst case, each neighbor has a
- // different color. So we find at least one
- // available color between 0 and no_neighbors.
+ // In the worst case, each neighbor has a different color. So we
+ // find at least one available color between 0 and no_neighbors.
color_finder.resize(no_neighbors+1);
for (unsigned int j=0; j<=no_neighbors; ++j)
color_finder[j]=true;
end = connectivity.row_end(cell);
for (; neighbor!=end ; ++neighbor)
{
- // Mark the color that a neighbor within the
- // partition has as taken
+ // Mark the color that a neighbor within the partition has
+ // as taken
if (cell_partition[*neighbor] == part &&
cell_color[*neighbor] <= no_neighbors)
color_finder[cell_color[*neighbor]] = false;
}
- // Choose the smallest color that is not taken
- // for the block
+ // Choose the smallest color that is not taken for the block
cell_color[cell]=0;
while (color_finder[cell_color[cell]] == false)
cell_color[cell]++;
if (cell_color[cell] > max_color)
max_color = cell_color[cell];
}
- // Reorder within partition: First, all blocks
- // that belong the 0 and then so on until
- // those with color max (Note that the smaller
- // the number the larger the partition)
+ // Reorder within partition: First, all blocks that belong the 0 and
+ // then so on until those with color max (Note that the smaller the
+ // number the larger the partition)
for (unsigned int color=0; color<=max_color; color++)
{
task_info.partition_color_blocks_data.push_back(color_counter);
partition_list = renumbering;
- // in debug mode, check that the partition
- // color list is one-to-one
+ // in debug mode, check that the partition color list is one-to-one
#ifdef DEBUG
{
std::vector<unsigned int> sorted_pc_list (partition_color_list);
}
#endif
- // set the start list for each block and
- // compute the renumbering of cells
+ // set the start list for each block and compute the renumbering of
+ // cells
std::vector<unsigned int> block_start(size_info.n_macro_cells+1);
std::vector<unsigned int> irregular(size_info.n_macro_cells);
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-
task_info.position_short_block = task_info.n_blocks-1;
unsigned int cluster_size = task_info.block_size*vectorization_length;
- // create the connectivity graph without
- // internal blocking
+ // create the connectivity graph without internal blocking
CompressedSimpleSparsityPattern connectivity;
make_connectivity_graph (size_info, task_info, renumbering,irregular_cells,
false, connectivity);
// Create cell-block partitioning.
- // For each block of cells, this variable
- // saves to which partitions the block
- // belongs. Initialize all to n_macro_cells to
- // mark them as not yet assigned a partition.
+ // For each block of cells, this variable saves to which partitions the
+ // block belongs. Initialize all to n_macro_cells to mark them as not
+ // yet assigned a partition.
std::vector<unsigned int> cell_partition (size_info.n_active_cells,
size_info.n_active_cells);
std::vector<unsigned int> neighbor_list;
std::vector<unsigned int> neighbor_neighbor_list;
- // In element j of this variable, one puts the
- // old number of the block that should be the
- // jth block in the new numeration.
+ // In element j of this variable, one puts the old number of the block
+ // that should be the jth block in the new numeration.
std::vector<unsigned int> partition_list(size_info.n_active_cells,0);
std::vector<unsigned int> partition_partition_list(size_info.n_active_cells,0);
- // This vector points to the start of each
- // partition.
+ // This vector points to the start of each partition.
std::vector<unsigned int> partition_size(2,0);
unsigned int partition = 0,start_up=0,counter=0;
bool work = true;
unsigned int remainder = cluster_size;
- // this performs a classical breath-first
- // search in the connectivity graph of the
- // cells under the restriction that the size
- // of the partitions should be a multiple of
- // the given block size
+ // this performs a classical breath-first search in the connectivity
+ // graph of the cells under the restriction that the size of the
+ // partitions should be a multiple of the given block size
while (work)
{
- // put the cells with neighbors on remote MPI
- // processes up front
+ // put the cells with neighbors on remote MPI processes up front
if (start_nonboundary>0)
{
for (unsigned int cell=0; cell<start_nonboundary; ++cell)
if (remainder == cluster_size)
remainder = 0;
- // adjust end of boundary cells to the
- // remainder
+ // adjust end of boundary cells to the remainder
size_info.boundary_cells_end += (remainder+vectorization_length-1)/vectorization_length;
}
else
{
- // To start up, set the start_up cell to
- // partition and list all its neighbors.
+ // To start up, set the start_up cell to partition and list all
+ // its neighbors.
cell_partition[start_up] = partition;
neighbor_list.push_back(start_up);
partition_list[counter++] = start_up;
neighbor_neighbor_list.resize(0);
}
not_connect:
- // One has to check if the graph is not
- // connected so we have to find another
- // partition.
+ // One has to check if the graph is not connected so we have to find
+ // another partition.
work = false;
for (unsigned int j=start_up; j<size_info.n_active_cells; ++j)
if (cell_partition[j] == size_info.n_active_cells)
{
// Create partitioning within partitions.
- // For each block of cells, this variable
- // saves to which partitions the block
- // belongs. Initialize all to n_macro_cells to
- // mark them as not yet assigned a partition.
+ // For each block of cells, this variable saves to which partitions
+ // the block belongs. Initialize all to n_macro_cells to mark them as
+ // not yet assigned a partition.
std::vector<unsigned int> cell_partition_l2(size_info.n_active_cells,
size_info.n_active_cells);
task_info.partition_color_blocks_row_index.resize(partition+1,0);
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;
int index_before = neighbor_neighbor_list.size(),
index = index_before;
{
- // put the cells into separate lists for each
- // FE index within one partition-partition
+ // put the cells into separate lists for each FE index
+ // within one partition-partition
missing_macros = 0;
std::vector<unsigned int> remaining_per_macro_cell
(max_fe_index);
[cell]]].
push_back(partition_partition_list[cell]);
}
- // check how many more cells are needed in the
- // lists
+ // check how many more cells are needed in the lists
for (unsigned int j=0; j<max_fe_index; j++)
{
remaining_per_macro_cell[j] =
missing_macros = task_info.block_size -
(missing_macros%task_info.block_size);
- // now we realized that there are some cells
- // missing.
+ // now we realized that there are some cells missing.
while (missing_macros>0 || filled == false)
{
if (index==0)
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),
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)
{
}
if (hp_bool == true)
{
- // set the renumbering according to their
- // active FE index within one
- // partition-partition which was implicitly
- // assumed above
+ // set the renumbering according to their active FE
+ // index within one partition-partition which was
+ // implicitly assumed above
cell = counter - partition_counter;
for (unsigned int j=0; j<max_fe_index; j++)
{
}
- // internal sparsity pattern structure that is
- // specialized to the situation in the setup
- // of the cells the dofs belong to: allows
- // only insertion of entries at the end
namespace internal
{
- // rudimentary version of a vector that keeps
- // entries always ordered
+ // rudimentary version of a vector that keeps entries always ordered
class ordered_vector : public std::vector<unsigned int>
{
public:
}
- // insert a given entry. dat is a pointer
- // within this vector (make sure that it
- // really stays there, otherwise there will be
- // an error)
+ // insert a given entry. dat is a pointer within this vector (the user
+ // needs to make sure that it really stays there)
void insert (const unsigned int entry,
std::vector<unsigned int>::iterator &dat)
{
if (dat == end())
{
push_back(entry);
+ dat = end();
}
else if (*dat > entry)
{
- this->std::vector<unsigned int>::insert (dat, entry);
+ dat = this->std::vector<unsigned int>::insert (dat, entry);
+ ++dat;
}
- ++dat;
+ else
+ ++dat;
}
};
}
void DoFInfo::renumber_dofs (std::vector<unsigned int> &renumbering)
{
- // first renumber all locally owned degrees of
- // freedom
+ // first renumber all locally owned degrees of freedom
AssertDimension (vector_partitioner->local_size(),
vector_partitioner->size());
const unsigned int local_size = vector_partitioner->local_size();
for (unsigned int i=0; i<constrained_dofs.size(); ++i)
new_constrained_dofs[i] = renumbering[constrained_dofs[i]];
- // the new constrained DoFs should be sorted
- // already as they are not contained in
- // dof_indices and then get contiguous numbers
+ // the new constrained DoFs should be sorted already as they are not
+ // contained in dof_indices and then get contiguous numbers
#ifdef DEBUG
for (unsigned int i=1; i<new_constrained_dofs.size(); ++i)
Assert (new_constrained_dofs[i] > new_constrained_dofs[i-1], ExcInternalError());
// $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
struct MappingInfo
{
/**
- * Determines how many bits of an unsigned int
- * are used to distinguish the cell types
- * (Cartesian, with constant Jacobian, or
- * general)
+ * Determines how many bits of an unsigned int are used to distinguish
+ * the cell types (Cartesian, with constant Jacobian, or general)
*/
static const std::size_t n_cell_type_bits = 2;
/**
- * Determines how many types of different
- * cells can be detected at most. Corresponds
- * to the number of bits we reserved for it.
+ * Determines how many types of different cells can be detected at
+ * most. Corresponds to the number of bits we reserved for it.
*/
static const unsigned int n_cell_types = 1U<<n_cell_type_bits;
MappingInfo();
/**
- * Computes the information in the given
- * cells. The cells are specified by the level
- * and the index within the level (as given by
- * CellIterator::level() and
- * CellIterator::index(), in order to allow
- * for different kinds of iterators,
- * e.g. standard DoFHandler, multigrid, etc.)
- * on a fixed Triangulation. In addition, a
- * mapping and several quadrature formulas are
- * given.
+ * Computes the information in the given cells. The cells are specified
+ * by the level and the index within the level (as given by
+ * CellIterator::level() and CellIterator::index(), in order to allow
+ * for different kinds of iterators, e.g. standard DoFHandler,
+ * multigrid, etc.) on a fixed Triangulation. In addition, a mapping
+ * and several quadrature formulas are given.
*/
void initialize (const dealii::Triangulation<dim> &tria,
const std::vector<std::pair<unsigned int,unsigned int> > &cells,
const UpdateFlags update_flags);
/**
- * Helper function to determine which update
- * flags must be set in the internal functions
- * to initialize all data as requested by the
- * user.
+ * Helper function to determine which update flags must be set in the
+ * internal functions to initialize all data as requested by the user.
*/
UpdateFlags
compute_update_flags (const UpdateFlags update_flags,
const std::vector<dealii::hp::QCollection<1> > &quad) const;
/**
- * Returns the type of a given cell as
- * detected during initialization.
+ * Returns the type of a given cell as detected during initialization.
*/
CellType get_cell_type (const unsigned int cell_chunk_no) const;
/**
- * Returns the type of a given cell as
- * detected during initialization.
+ * Returns the type of a given cell as detected during initialization.
*/
unsigned int get_cell_data_index (const unsigned int cell_chunk_no) const;
void clear ();
/**
- * Returns the memory consumption of this
- * class in bytes.
+ * Returns the memory consumption of this class in bytes.
*/
std::size_t memory_consumption() const;
/**
- * Prints a detailed summary of memory
- * consumption in the different structures of
- * this class to the given output stream.
+ * Prints a detailed summary of memory consumption in the different
+ * structures of this class to the given output stream.
*/
template <typename STREAM>
void print_memory_consumption(STREAM &out,
const SizeInfo &size_info) const;
/**
- * Stores whether a cell is Cartesian, has
- * constant transform data (Jacobians) or is
- * general. cell_type % 4 gives this
- * information (0: Cartesian, 1: constant
- * Jacobian throughout cell, 2: general cell),
- * and cell_type / 4 gives the index in the
- * data field of where to find the information
- * in the fields Jacobian and JxW values
- * (except for quadrature points, for which
- * the index runs as usual).
+ * Stores whether a cell is Cartesian, has constant transform data
+ * (Jacobians) or is general. cell_type % 4 gives this information (0:
+ * Cartesian, 1: constant Jacobian throughout cell, 2: general cell),
+ * and cell_type / 4 gives the index in the data field of where to find
+ * the information in the fields Jacobian and JxW values (except for
+ * quadrature points, for which the index runs as usual).
*/
std::vector<unsigned int> cell_type;
/**
- * The first field stores the inverse Jacobian
- * for Cartesian cells: There, it is a
- * diagonal rank-2 tensor, so we actually just
- * store a rank-1 tensor. It is the same on
- * all cells, therefore we only store it once
- * per cell, and use similarities from one
- * cell to another, too (on structured meshes,
- * there are usually many cells with the same
+ * The first field stores the inverse Jacobian for Cartesian cells:
+ * There, it is a diagonal rank-2 tensor, so we actually just store a
+ * rank-1 tensor. It is the same on all cells, therefore we only store
+ * it once per cell, and use similarities from one cell to another, too
+ * (on structured meshes, there are usually many cells with the same
* Jacobian).
*
- * The second field stores the Jacobian
- * determinant for Cartesian cells (without
- * the quadrature weight, which depends on the
- * quadrature point, whereas the determinant
- * is the same on each quadrature point).
+ * The second field stores the Jacobian determinant for Cartesian cells
+ * (without the quadrature weight, which depends on the quadrature
+ * point, whereas the determinant is the same on each quadrature point).
*/
AlignedVector<std::pair<Tensor<1,dim,VectorizedArray<Number> >,
VectorizedArray<Number> > > cartesian_data;
/**
- * The first field stores the Jacobian for
- * non-Cartesian cells where all the Jacobians
- * on the cell are the same (i.e., constant,
- * which comes from a linear transformation
- * from unit to real cell). Also use
- * similarities from one cell to another (on
- * structured meshes, there are usually many
- * cells with the same Jacobian).
+ * The first field stores the Jacobian for non-Cartesian cells where all
+ * the Jacobians on the cell are the same (i.e., constant, which comes
+ * from a linear transformation from unit to real cell). Also use
+ * similarities from one cell to another (on structured meshes, there
+ * are usually many cells with the same Jacobian).
*
- * The second field stores the Jacobian
- * determinant for non-Cartesian cells with
- * constant Jacobian throughout the cell
- * (without the quadrature weight, which
- * depends on the quadrature point, whereas
- * the determinant is the same on each
- * quadrature point).
+ * The second field stores the Jacobian determinant for non-Cartesian
+ * cells with constant Jacobian throughout the cell (without the
+ * quadrature weight, which depends on the quadrature point, whereas the
+ * determinant is the same on each quadrature point).
*/
AlignedVector<std::pair<Tensor<2,dim,VectorizedArray<Number> >,
VectorizedArray<Number> > > affine_data;
/**
- * Definition of a structure that stores data
- * that depends on the quadrature formula (if
- * we have more than one quadrature formula on
- * a given problem, these fields will be
- * different)
+ * Definition of a structure that stores data that depends on the
+ * quadrature formula (if we have more than one quadrature formula on a
+ * given problem, these fields will be different)
*/
struct MappingInfoDependent
{
/**
- * This field stores the row starts for the
- * inverse Jacobian transformations,
- * quadrature weights and second derivatives.
+ * This field stores the row starts for the inverse Jacobian
+ * transformations, quadrature weights and second derivatives.
*/
std::vector<unsigned int> rowstart_jacobians;
/**
- * This field stores the inverse Jacobian
- * transformation from unit to real cell,
- * which is needed for most gradient
- * transformations (corresponds to
- * FEValues::inverse_jacobian) for general
- * cells.
+ * This field stores the inverse Jacobian transformation from unit to
+ * real cell, which is needed for most gradient transformations
+ * (corresponds to FEValues::inverse_jacobian) for general cells.
*/
AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians;
/**
- * This field stores the Jacobian
- * determinant times the quadrature weights
- * (JxW in deal.II speak) for general cells.
+ * This field stores the Jacobian determinant times the quadrature
+ * weights (JxW in deal.II speak) for general cells.
*/
AlignedVector<VectorizedArray<Number> > JxW_values;
/**
- * Stores the diagonal part of the gradient of
- * the inverse Jacobian transformation. The
- * first index runs over the derivatives
- * $\partial^2/\partial x_i^2$, the second
- * over the space coordinate. Needed for
- * computing the Laplacian of FE functions on
- * the real cell. Uses a separate storage from
- * the off-diagonal part $\partial^2/\partial
- * x_i \partial x_j, i\neq j$ because that is
+ * Stores the diagonal part of the gradient of the inverse Jacobian
+ * transformation. The first index runs over the derivatives
+ * $\partial^2/\partial x_i^2$, the second over the space
+ * coordinate. Needed for computing the Laplacian of FE functions on
+ * the real cell. Uses a separate storage from the off-diagonal part
+ * $\partial^2/\partial x_i \partial x_j, i\neq j$ because that is
* only needed for computing a full Hessian.
*/
AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians_grad_diag;
/**
- * Stores the off-diagonal part of the
- * gradient of the inverse Jacobian
- * transformation. Because of symmetry, only
- * the upper diagonal part is needed. The
- * first index runs through the derivatives
- * row-wise, i.e., $\partial^2/\partial x_1
- * \partial x_2$ first, then
- * $\partial^2/\partial x_1 \partial x_3$, and
- * so on. The second index is the spatial
- * coordinate. Not filled currently.
+ * Stores the off-diagonal part of the gradient of the inverse
+ * Jacobian transformation. Because of symmetry, only the upper
+ * diagonal part is needed. The first index runs through the
+ * derivatives row-wise, i.e., $\partial^2/\partial x_1 \partial x_2$
+ * first, then $\partial^2/\partial x_1 \partial x_3$, and so on. The
+ * second index is the spatial coordinate. Not filled currently.
*/
AlignedVector<Tensor<1,(dim>1?dim*(dim-1)/2:1),
Tensor<1,dim,VectorizedArray<Number> > > > jacobians_grad_upper;
/**
- * Stores the row start for quadrature points
- * in real coordinates for both types of
- * cells. Note that Cartesian cells will have
- * shorter fields (length is @p n_q_points_1d)
- * than non-Cartesian cells (length is @p
- * n_q_points).
+ * Stores the row start for quadrature points in real coordinates for
+ * both types of cells. Note that Cartesian cells will have shorter
+ * fields (length is @p n_q_points_1d) than non-Cartesian cells
+ * (length is @p n_q_points).
*/
std::vector<unsigned int> rowstart_q_points;
/**
- * Stores the quadrature points in real
- * coordinates for Cartesian cells (does not
- * need to store the full data on all points)
+ * Stores the quadrature points in real coordinates for Cartesian
+ * cells (does not need to store the full data on all points)
*/
AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
/**
- * The dim-dimensional quadrature formula
- * underlying the problem (constructed from a
- * 1D tensor product quadrature formula).
+ * The dim-dimensional quadrature formula underlying the problem
+ * (constructed from a 1D tensor product quadrature formula).
*/
dealii::hp::QCollection<dim> quadrature;
/**
- * The (dim-1)-dimensional quadrature formula
- * corresponding to face evaluation
- * (constructed from a 1D tensor product
- * quadrature formula).
+ * The (dim-1)-dimensional quadrature formula corresponding to face
+ * evaluation (constructed from a 1D tensor product quadrature
+ * formula).
*/
dealii::hp::QCollection<dim-1> face_quadrature;
/**
- * The number of quadrature points for the
- * current quadrature formula.
+ * The number of quadrature points for the current quadrature formula.
*/
std::vector<unsigned int> n_q_points;
/**
- * The number of quadrature points for the
- * current quadrature formula when applied to
- * a face. Only set if the quadrature formula
- * is derived from a tensor product, since it
- * is not defined from the full quadrature
- * formula otherwise.
+ * The number of quadrature points for the current quadrature formula
+ * when applied to a face. Only set if the quadrature formula is
+ * derived from a tensor product, since it is not defined from the
+ * full quadrature formula otherwise.
*/
std::vector<unsigned int> n_q_points_face;
/**
- * The quadrature weights (vectorized data
- * format) on the unit cell.
+ * The quadrature weights (vectorized data format) on the unit cell.
*/
std::vector<AlignedVector<VectorizedArray<Number> > > quadrature_weights;
/**
- * This variable stores the number of
- * quadrature points for all quadrature
- * indices in the underlying element for
- * easier access to data in the hp case.
+ * This variable stores the number of quadrature points for all
+ * quadrature indices in the underlying element for easier access to
+ * data in the hp case.
*/
std::vector<unsigned int> quad_index_conversion;
/**
- * Returns the quadrature index for a given
- * number of quadrature points. If not in hp
- * mode or if the index is not found, this
- * function always returns index 0. Hence,
- * this function does not check whether the
- * given degree is actually present.
+ * Returns the quadrature index for a given number of quadrature
+ * points. If not in hp mode or if the index is not found, this
+ * function always returns index 0. Hence, this function does not
+ * check whether the given degree is actually present.
*/
unsigned int
quad_index_from_n_q_points (const unsigned int n_q_points) const;
/**
- * Prints a detailed summary of memory
- * consumption in the different structures of
- * this class to the given output stream.
+ * Prints a detailed summary of memory consumption in the different
+ * structures of this class to the given output stream.
*/
template <typename STREAM>
void print_memory_consumption(STREAM &out,
};
/**
- * Contains all the stuff that depends on the
- * quadrature formula
+ * Contains all the stuff that depends on the quadrature formula
*/
std::vector<MappingInfoDependent> mapping_data_gen;
/**
- * Stores whether JxW values have been
- * initialized
+ * Stores whether JxW values have been initialized
*/
bool JxW_values_initialized;
/**
- * Stores whether we computed second
- * derivatives.
+ * Stores whether we computed second derivatives.
*/
bool second_derivatives_initialized;
/**
- * Stores whether we computed quadrature
- * points.
+ * Stores whether we computed quadrature points.
*/
bool quadrature_points_initialized;
/**
- * Internal temporary data used for the
- * initialization.
+ * Internal temporary data used for the initialization.
*/
struct CellData
{
};
/**
- * Helper function called internally during
- * the initialize function.
+ * Helper function called internally during the initialize function.
*/
void evaluate_on_cell (const dealii::Triangulation<dim> &tria,
const std::pair<unsigned int,unsigned int> *cells,
// $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
compute_update_flags (const UpdateFlags update_flags,
const std::vector<dealii::hp::QCollection<1> > &quad) const
{
- // this class is build around the evaluation
- // this class is build around the evaluation
- // of inverse gradients, so compute them in
- // any case
+ // this class is build around the evaluation this class is build around
+ // the evaluation of inverse gradients, so compute them in any case
UpdateFlags new_flags = update_inverse_jacobians;
- // if the user requested gradients, need
- // inverse Jacobians
+ // if the user requested gradients, need inverse Jacobians
if (update_flags & update_gradients || update_flags & update_inverse_jacobians)
new_flags |= update_inverse_jacobians;
if (update_flags & update_JxW_values)
new_flags |= update_JxW_values;
- // for Hessian information, need inverse
- // Jacobians and the derivative of Jacobians
- // (these two together will give use the
- // gradients of the inverse Jacobians, which
- // is what we need)
+ // for Hessian information, need inverse Jacobians and the derivative of
+ // Jacobians (these two together will give use the gradients of the
+ // inverse Jacobians, which is what we need)
if (update_flags & update_hessians || update_flags & update_jacobian_grads)
new_flags |= update_jacobian_grads;
if (update_flags & update_quadrature_points)
new_flags |= update_quadrature_points;
- // there is one more thing: if we have a
- // quadrature formula with only one quadrature
- // point on the first component, but more
- // points on later components, we need to have
- // Jacobian gradients anyway in order to
- // determine whether the Jacobian is constant
- // throughout a cell
+ // there is one more thing: if we have a quadrature formula with only
+ // one quadrature point on the first component, but more points on later
+ // components, we need to have Jacobian gradients anyway in order to
+ // determine whether the Jacobian is constant throughout a cell
bool formula_with_one_point = false;
for (unsigned int i=0; i<quad[0].size(); ++i)
if (quad[0][i].size() == 1)
mapping_data_gen.resize (n_quads);
cell_type.resize (n_macro_cells);
- // dummy FE that is used to set up an FEValues
- // object. Do not need the actual finite
- // element because we will only evaluate
- // quantities for the mapping that are
- // independent of the FE
+ // dummy FE that is used to set up an FEValues object. Do not need the
+ // actual finite element because we will only evaluate quantities for
+ // the mapping that are independent of the FE
FE_Nothing<dim> dummy_fe;
UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
if (update_flags & update_quadrature_points)
quadrature_points_initialized = true;
- // when we make comparisons about the size of
- // Jacobians we need to know the approximate
- // size of typical entries in Jacobians. We
- // need to fix the Jacobian size once and for
- // all. We choose the diameter of the first
- // cell (on level zero, which is the best
- // accuracy we can hope for, since diameters
- // on finer levels are computed by differences
- // of nearby cells). If the mesh extends over
- // a certain domain, the precision of double
- // values is essentially limited by this
- // precision.
+ // when we make comparisons about the size of Jacobians we need to know
+ // the approximate size of typical entries in Jacobians. We need to fix
+ // the Jacobian size once and for all. We choose the diameter of the
+ // first cell (on level zero, which is the best accuracy we can hope
+ // for, since diameters on finer levels are computed by differences of
+ // nearby cells). If the mesh extends over a certain domain, the
+ // precision of double values is essentially limited by this precision.
const double jacobian_size = internal::get_jacobian_size(tria);
- // objects that hold the data for up to
- // vectorization_length cells while we fill them up. Only
- // after all vectorization_length cells have been
- // processed, we can insert the data into the
- // data structures of this class
+ // objects that hold the data for up to vectorization_length cells while
+ // we fill them up. Only after all vectorization_length cells have been
+ // processed, we can insert the data into the data structures of this
+ // class
CellData data (jacobian_size);
for (unsigned int my_q=0; my_q<n_quads; ++my_q)
if (n_hp_quads > 1)
current_data.quad_index_conversion[q] = n_q_points;
- // To walk on the diagonal for lexicographic
- // ordering, we have to jump one index ahead
- // in each direction. For direction 0, this is
- // just the next point, for direction 1, it
- // means adding n_q_points_1d, and so on.
+ // To walk on the diagonal for lexicographic ordering, we have
+ // to jump one index ahead in each direction. For direction 0,
+ // this is just the next point, for direction 1, it means adding
+ // n_q_points_1d, and so on.
step_size_cartesian[q] = 0;
unsigned int factor = 1;
for (unsigned int d=0; d<dim; ++d)
}
}
- // if there are no cells, there is nothing to
- // do
+ // if there are no cells, there is nothing to do
if (cells.size() == 0)
continue;
Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
Tensor<1,dim,VectorizedArray<Number> > tmp;
- // encodes the cell types of the current
- // cell. Since several cells must be
- // considered together, this variable holds
- // the individual info of the last chunk of
- // cells
+ // encodes the cell types of the current cell. Since several cells
+ // must be considered together, this variable holds the individual
+ // info of the last chunk of cells
CellType cell_t [vectorization_length],
cell_t_prev [vectorization_length];
for (unsigned int j=0; j<vectorization_length; ++j)
cell_t_prev[j] = undefined;
- // fe_values object that is used to compute
- // the mapping data. for the hp case there
- // might be more than one finite
- // element. since we manually select the
- // active FE index and not via a
- // hp::DoFHandler<dim>::active_cell_iterator,
- // we need to manually select the correct
- // finite element, so just hold a vector of
+ // fe_values object that is used to compute the mapping data. for
+ // the hp case there might be more than one finite element. since we
+ // manually select the active FE index and not via a
+ // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
+ // select the correct finite element, so just hold a vector of
// FEValues
std::vector<std_cxx1x::shared_ptr<FEValues<dim> > >
fe_values (current_data.quadrature.size());
(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)
{
// loop over all cells
for (unsigned int cell=0; cell<n_macro_cells; ++cell)
{
- // GENERAL OUTLINE: First generate the data in
- // format "number" for vectorization_length cells, and
- // then find the most general type of cell for
- // appropriate vectorized formats. then fill
- // this data in
+ // GENERAL OUTLINE: First generate the data in format "number"
+ // for vectorization_length cells, and then find the most
+ // general type of cell for appropriate vectorized formats. then
+ // fill this data in
const unsigned int fe_index = active_fe_index.size() > 0 ?
active_fe_index[cell] : 0;
const unsigned int n_q_points = current_data.n_q_points[fe_index];
FEValues<dim> &fe_val = *fe_values[fe_index];
data.resize (n_q_points);
- // if the fe index has changed from the
- // previous cell, set the old cell type to
- // invalid (otherwise, we might detect
+ // if the fe index has changed from the previous cell, set the
+ // old cell type to invalid (otherwise, we might detect
// similarity due to some cells further ahead)
if (cell > 0 && active_fe_index.size() > 0 &&
active_fe_index[cell] != active_fe_index[cell-1])
evaluate_on_cell (tria, &cells[cell*vectorization_length],
cell, my_q, cell_t_prev, cell_t, fe_val, data);
- // now reorder the data into vectorized
- // types. if we are here for the first time,
- // we need to find out whether the Jacobian
- // allows for some simplification (Cartesian,
- // affine) taking vectorization_length cell together and
- // we have to insert that data into the
- // respective fields. Also, we have to
- // compress different cell indicators into one
- // structure.
+ // now reorder the data into vectorized types. if we are here
+ // for the first time, we need to find out whether the Jacobian
+ // allows for some simplification (Cartesian, affine) taking
+ // vectorization_length cell together and we have to insert that
+ // data into the respective fields. Also, we have to compress
+ // different cell indicators into one structure.
if (my_q == 0)
{
- // find the most general cell type (most
- // general type is 2 (general cell))
+ // find the most general cell type (most general type is 2
+ // (general cell))
CellType most_general_type = cartesian;
for (unsigned int j=0; j<vectorization_length; ++j)
if (cell_t[j] > most_general_type)
AssertIndexRange (most_general_type, 3);
unsigned int insert_position = numbers::invalid_unsigned_int;
- // Cartesian cell with diagonal Jacobian: only
- // insert the diagonal of the inverse and the
- // Jacobian determinant. We do this by using
- // an std::map that collects pointers to all
- // Cartesian Jacobians. We need a pointer in
- // the std::map because it cannot store data
- // based on VectorizedArray (alignment
- // issue). We circumvent the problem by
- // temporarily filling the next value into the
- // cartesian_data field and, in case we did an
- // insertion, the data is already in the
- // correct place.
+ // Cartesian cell with diagonal Jacobian: only insert the
+ // diagonal of the inverse and the Jacobian determinant. We
+ // do this by using an std::map that collects pointers to
+ // all Cartesian Jacobians. We need a pointer in the
+ // std::map because it cannot store data based on
+ // VectorizedArray (alignment issue). We circumvent the
+ // problem by temporarily filling the next value into the
+ // cartesian_data field and, in case we did an insertion,
+ // the data is already in the correct place.
if (most_general_type == cartesian)
{
std::pair<Tensor<1,dim,VEC_ARRAY>,unsigned int> new_entry;
insert_position = new_entry.second;
}
- // Constant Jacobian case. same strategy as
- // before, but with other data fields
+ // Constant Jacobian case. same strategy as before, but with
+ // other data fields
else if (most_general_type == affine)
{
std::pair<Tensor<2,dim,VEC_ARRAY>,unsigned int> new_entry;
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());
} // 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 =
current_data.jacobians.push_back (transpose(invert(jac)));
const Tensor<2,dim,VectorizedArray<Number> > &inv_jac = current_data.jacobians.back();
- // TODO: deal.II does not use abs on
- // determinants. Is there an assumption
- // somewhere that negative determinants are
- // not allowed?
+ // TODO: deal.II does not use abs on determinants. Is
+ // there an assumption somewhere that negative
+ // determinants are not allowed?
if (update_flags & update_JxW_values)
current_data.JxW_values.push_back
(std::abs(det) * current_data.quadrature_weights[fe_index][q]);
if (update_flags & update_jacobian_grads)
{
- // for second derivatives on the real cell,
- // need also the gradient of the inverse
- // Jacobian J. This involves some calculus and
- // is done vectorized. This is very cheap
- // compared to what fe_values does (in early
- // 2011). If L is the gradient of the jacobian
- // on the unit cell, the gradient of the
- // inverse is given by (multidimensional
- // calculus) - J * (J * L) * J (the third J is
- // because we need to transform the gradient L
- // from the unit to the real cell, and then
- // apply the inverse Jacobian). Compare this
- // with 1D with j(x) = 1/k(phi(x)), where j =
- // phi' is the inverse of the jacobian and k
- // is the derivative of the jacobian on the
- // unit cell. Then j' = phi' k'/k^2 = j k'
- // j^2.
+ // for second derivatives on the real cell, need
+ // also the gradient of the inverse Jacobian J. This
+ // involves some calculus and is done
+ // vectorized. This is very cheap compared to what
+ // fe_values does (in early 2011). If L is the
+ // gradient of the jacobian on the unit cell, the
+ // gradient of the inverse is given by
+ // (multidimensional calculus) - J * (J * L) * J
+ // (the third J is because we need to transform the
+ // gradient L from the unit to the real cell, and
+ // then apply the inverse Jacobian). Compare this
+ // with 1D with j(x) = 1/k(phi(x)), where j = phi'
+ // is the inverse of the jacobian and k is the
+ // derivative of the jacobian on the unit cell. Then
+ // j' = phi' k'/k^2 = j k' j^2.
// compute: jac_grad = J*grad_unit(J^-1)
for (unsigned int d=0; d<dim; ++d)
if (update_flags & update_quadrature_points)
{
- // eventually we turn to the quadrature points
- // that we can compress in case we have
- // Cartesian cells. we also need to reorder
- // them into arrays of vectorized data types.
- // first go through the cells and find out how
- // much memory we need to allocate for the
- // quadrature points. We store n_q_points for
- // all cells but Cartesian cells. For
- // Cartesian cells, only need to store the
- // values on a diagonal through the cell
- // (n_q_points_1d). This will give (with some
- // little indexing) the location of all
+ // eventually we turn to the quadrature points that we can
+ // compress in case we have Cartesian cells. we also need to
+ // reorder them into arrays of vectorized data types. first
+ // go through the cells and find out how much memory we need
+ // to allocate for the quadrature points. We store
+ // n_q_points for all cells but Cartesian cells. For
+ // Cartesian cells, only need to store the values on a
+ // diagonal through the cell (n_q_points_1d). This will give
+ // (with some little indexing) the location of all
// quadrature points.
const unsigned int old_size =
current_data.quadrature_points.size();
current_data.rowstart_q_points[n_macro_cells] =
current_data.quadrature_points.size();
- // finally, need to invert and
- // transpose the Jacobians in the
- // cartesian_data and affine_data
- // fields and compute the JxW value.
+ // finally, fill the accumulated data for Cartesian and affine cells
+ // into cartesian_data and affine_data, invert and transpose the
+ // Jacobians, and compute the JxW value.
if (my_q == 0)
{
cartesian_data.resize(cartesians.size());
const unsigned int n_q_points = fe_val.n_quadrature_points;
const UpdateFlags update_flags = fe_val.get_update_flags();
- // this should be the same value as
- // used in HashValue::scaling (but we
+ // this should be the same value as used in HashValue::scaling (but we
// not have that field here)
const double zero_tolerance_double = data.jac_size *
std::numeric_limits<double>::epsilon() * 1024.;
fe_val.reinit(cell_it);
cell_t[j] = undefined;
- // extract quadrature points and store them
- // temporarily. if we have Cartesian cells, we
- // can compress the indices
+ // extract quadrature points and store them temporarily. if we have
+ // Cartesian cells, we can compress the indices
if (update_flags & update_quadrature_points)
for (unsigned int q=0; q<n_q_points; ++q)
{
data.quadrature_points[q][d][j] = point[d];
}
- // if this is not the first quadrature formula
- // and we already have determined that this
- // cell is either Cartesian or with constant
+ // if this is not the first quadrature formula and we already have
+ // determined that this cell is either Cartesian or with constant
// Jacobian, we have nothing more to do.
if (my_q > 0 && (get_cell_type(cell) == cartesian
|| get_cell_type(cell) == affine) )
continue;
- // first round: if the transformation is
- // detected to be the same as on the old cell,
- // we only need to copy over the data.
+ // first round: if the transformation is detected to be the same as
+ // on the old cell, we only need to copy over the data.
if (fe_val.get_cell_similarity() == CellSimilarity::translation
&&
my_q == 0)
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;
break;
}
- // check whether the Jacobian is diagonal to
- // machine accuracy
+ // check whether the Jacobian is diagonal to machine
+ // accuracy
bool cell_cartesian = jacobian_constant;
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
// Cartesian cell
if (cell_t[j] == cartesian)
{
- // set Jacobian into diagonal and clear
- // off-diagonal part
+ // set Jacobian into diagonal and clear off-diagonal part
for (unsigned int d=0; d<dim; ++d)
{
data.const_jac[d][d][j] = jac_0[d][d];
// general cell case
- // go through all quadrature points and fill
- // in the data into the temporary data
- // structures with slots for the vectorized
- // data types
+ // go through all quadrature points and fill in the data into the
+ // temporary data structures with slots for the vectorized data
+ // types
for (unsigned int q=0; q<n_q_points; ++q)
{
- // compress out very small numbers which are
- // only noise. Then it is cleaner to use zero
- // straight away (though it does not save any
- // memory)
+ // compress out very small numbers which are only noise. Then it
+ // is cleaner to use zero straight away (though it does not save
+ // any memory)
const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
data.general_jac[q][d][e][j] =
std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
- // need to do some calculus based on the
- // gradient of the Jacobian, in order to find
- // the gradient of the inverse Jacobian which
- // is needed in user code. however, we would
- // like to perform that on vectorized data
- // types instead of doubles or floats. to this
- // end, copy the gradients first
+ // need to do some calculus based on the gradient of the
+ // Jacobian, in order to find the gradient of the inverse
+ // Jacobian which is needed in user code. however, we would like
+ // to perform that on vectorized data types instead of doubles
+ // or floats. to this end, copy the gradients first
if (update_flags & update_jacobian_grads)
{
const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
}
}
- } // end loop over all entries in vectorization (vectorization_length cells)
+ } // end loop over all entries in vectorization (vectorization_length
+ // cells)
// set information for next cell
for (unsigned int j=0; j<vectorization_length; ++j)
(STREAM &out,
const SizeInfo &size_info) const
{
- // print_memory_statistics involves
- // global communication, so we can
- // disable the check here only if no
- // processor has any such data
+ // print_memory_statistics involves global communication, so we can
+ // disable the check here only if no processor has any such data
#if DEAL_II_COMPILER_SUPPORTS_MPI
unsigned int general_size_glob = 0, general_size_loc = jacobians.size();
MPI_Allreduce (&general_size_loc, &general_size_glob, 1, MPI_UNSIGNED,
/**
* This class collects all the data that is stored for the matrix free
- * implementation. The storage scheme is tailored towards several
- * loops performed with the same data, i.e., typically doing many
- * matrix-vector products or residual computations on the same
- * mesh. The class is used in step-37 and step-48.
+ * implementation. The storage scheme is tailored towards several loops
+ * performed with the same data, i.e., typically doing many matrix-vector
+ * products or residual computations on the same mesh. The class is used in
+ * step-37 and step-48.
*
* This class does not implement any operations involving finite element basis
* functions, i.e., regarding the operation performed on the cells. For these
public:
/**
- * Collects the options for initialization of
- * the MatrixFree class. The first
- * parameter specifies the MPI communicator to
- * be used, the second the parallelization
- * options in shared memory (task-based
- * parallelism, where one can choose between
- * no parallelism and three schemes that avoid
- * that cells with access to the same vector
- * entries are accessed simultaneously), the
- * third with the block size for task parallel
- * scheduling, the fourth the update flags
- * that should be stored by this class.
+ * Collects the options for initialization of the MatrixFree class. The
+ * first parameter specifies the MPI communicator to be used, the second the
+ * parallelization options in shared memory (task-based parallelism, where
+ * one can choose between no parallelism and three schemes that avoid that
+ * cells with access to the same vector entries are accessed
+ * simultaneously), the third with the block size for task parallel
+ * scheduling, the fourth the update flags that should be stored by this
+ * class.
*
- * The fifth parameter specifies the level in
- * the triangulation from which the indices
- * are to be used. If the level is set to
- * numbers::invalid_unsigned_int, the active
- * cells are traversed, and otherwise the
- * cells in the given level. This option has
- * no effect in case a DoFHandler or
- * hp::DoFHandler is given.
+ * The fifth parameter specifies the level in the triangulation from which
+ * the indices are to be used. If the level is set to
+ * numbers::invalid_unsigned_int, the active cells are traversed, and
+ * otherwise the cells in the given level. This option has no effect in case
+ * a DoFHandler or hp::DoFHandler is given.
*
- * The parameter @p initialize_plain_indices
- * indicates whether the DoFInfo class should
- * also allow for access to vectors without
- * resolving constraints.
+ * The parameter @p initialize_plain_indices indicates whether the DoFInfo
+ * class should also allow for access to vectors without resolving
+ * constraints.
*
- * The last two parameters allow the user to
- * disable some of the initialization
- * processes. For example, if only the
- * scheduling that avoids touching the same
- * vector/matrix indices simultaneously is to
- * be found, the mapping needs not be
- * initialized. Likewise, if the mapping has
- * changed from one iteration to the next but
- * the topology has not (like when using a
- * deforming mesh with MappingQEulerian), it
- * suffices to initialize the mapping only.
+ * The last two parameters allow the user to disable some of the
+ * initialization processes. For example, if only the scheduling that avoids
+ * touching the same vector/matrix indices simultaneously is to be found,
+ * the mapping needs not be initialized. Likewise, if the mapping has
+ * changed from one iteration to the next but the topology has not (like
+ * when using a deforming mesh with MappingQEulerian), it suffices to
+ * initialize the mapping only.
*/
struct AdditionalData
{
{};
/**
- * Sets the MPI communicator that the parallel
- * layout of the operator should be based
- * upon. Defaults to MPI_COMM_SELF, but should
- * be set to a communicator similar to the one
- * used for a distributed triangulation in
- * order to inform this class over all cells
- * that are present.
+ * Sets the MPI communicator that the parallel layout of the operator
+ * should be based upon. Defaults to MPI_COMM_SELF, but should be set to a
+ * communicator similar to the one used for a distributed triangulation in
+ * order to inform this class over all cells that are present.
*/
MPI_Comm mpi_communicator;
/**
- * Sets the scheme for task parallelism. There
- * are four options available. If set to @p
- * none, the operator application is done in
- * serial without shared memory
- * parallelism. If this class is used together
- * with MPI and MPI is also used for
- * parallelism within the nodes, this flag
- * should be set to @p none. The default value
- * is @p partition_partition, i.e. we actually
- * use multithreading with the first option
- * described below.
+ * Sets the scheme for task parallelism. There are four options
+ * available. If set to @p none, the operator application is done in
+ * serial without shared memory parallelism. If this class is used
+ * together with MPI and MPI is also used for parallelism within the
+ * nodes, this flag should be set to @p none. The default value is @p
+ * partition_partition, i.e. we actually use multithreading with the first
+ * option described below.
*
- * The first option @p partition_partition is
- * to partition the cells on two levels in
- * onion-skin-like partitions and forming
- * chunks of tasks_block_size after the
- * partitioning. The partitioning finds sets
- * of independent cells that enable working in
- * parallel without accessing the same vector
- * entries at the same time.
+ * The first option @p partition_partition is to partition the cells on
+ * two levels in onion-skin-like partitions and forming chunks of
+ * tasks_block_size after the partitioning. The partitioning finds sets of
+ * independent cells that enable working in parallel without accessing the
+ * same vector entries at the same time.
*
- * The second option @p partition_color is to
- * use a partition on the global level and
- * color cells within the partitions (where
- * all chunks within a color are
- * independent). Here, the subdivision into
- * chunks of cells is done before the
- * partitioning, which might give worse
- * partitions but better cache performance if
- * degrees of freedom are not renumbered.
+ * The second option @p partition_color is to use a partition on the
+ * global level and color cells within the partitions (where all chunks
+ * within a color are independent). Here, the subdivision into chunks of
+ * cells is done before the partitioning, which might give worse
+ * partitions but better cache performance if degrees of freedom are not
+ * renumbered.
*
- * The third option @p color is to use a
- * traditional algorithm of coloring on the
- * global level. This scheme is a special case
- * of the second option where only one
- * partition is present. Note that for
- * problems with hanging nodes, there are
- * quite many colors (50 or more in 3D), which
- * might degrade parallel performance (bad
- * cache behavior, many synchronization
- * points).
+ * The third option @p color is to use a traditional algorithm of coloring
+ * on the global level. This scheme is a special case of the second option
+ * where only one partition is present. Note that for problems with
+ * hanging nodes, there are quite many colors (50 or more in 3D), which
+ * might degrade parallel performance (bad cache behavior, many
+ * synchronization points).
*/
TasksParallelScheme tasks_parallel_scheme;
/**
- * Sets the number of so-called macro cells
- * that should form one partition. If zero
- * size is given, the class tries to find a
- * good size for the blocks based on the
- * number of threads stored in
- * multithread_info.n_default_threads and the
- * number of cells present. Otherwise, the
- * given number is used. If the given number
- * is larger than one third of the number of
- * total cells, this means no
- * parallelism. Note that in the case
- * vectorization is used, a macro cell
+ * Sets the number of so-called macro cells that should form one
+ * partition. If zero size is given, the class tries to find a good size
+ * for the blocks based on the number of threads stored in
+ * multithread_info.n_default_threads and the number of cells
+ * present. Otherwise, the given number is used. If the given number is
+ * larger than one third of the number of total cells, this means no
+ * parallelism. Note that in the case vectorization is used, a macro cell
* consists of more than one physical cell.
*/
unsigned int tasks_block_size;
/**
- * This flag is used to determine which
- * quantities should be cached. This class can
- * cache data needed for gradient computations
- * (inverse Jacobians), Jacobian determinants
- * (JxW), quadrature points as well as data
- * for Hessians (derivative of Jacobians). By
- * default, only data for gradients and
- * Jacobian determinants times quadrature
- * weights, JxW, are cached. If quadrature
- * points or second derivatives are needed,
- * they must be specified by this field (even
- * though second derivatives might still be
- * evaluated on Cartesian cells without this
- * option set here, since there the Jacobian
- * describes the mapping completely).
+ * This flag is used to determine which quantities should be cached. This
+ * class can cache data needed for gradient computations (inverse
+ * Jacobians), Jacobian determinants (JxW), quadrature points as well as
+ * data for Hessians (derivative of Jacobians). By default, only data for
+ * gradients and Jacobian determinants times quadrature weights, JxW, are
+ * cached. If quadrature points or second derivatives are needed, they
+ * must be specified by this field (even though second derivatives might
+ * still be evaluated on Cartesian cells without this option set here,
+ * since there the Jacobian describes the mapping completely).
*/
UpdateFlags mapping_update_flags;
/**
- * If working on a MGDoFHandler, this option
- * can be used to define whether we work on a
- * certain level of the mesh, and not the
- * active cells. If set to
- * invalid_unsigned_int (which is the default
- * value), the active cells are gone through,
- * otherwise the level given by this
- * parameter.
+ * If working on a MGDoFHandler, this option can be used to define whether
+ * we work on a certain level of the mesh, and not the active cells. If
+ * set to invalid_unsigned_int (which is the default value), the active
+ * cells are gone through, otherwise the level given by this parameter.
*/
unsigned int level_mg_handler;
/**
- * Controls whether to allow reading from
- * vectors without resolving constraints,
- * i.e., just read the local values of the
- * vector. By default, this option is
- * disabled, so if you want to use
- * FEEvaluationBase::read_dof_values_plain,
- * this flag needs to be set.
+ * Controls whether to allow reading from vectors without resolving
+ * constraints, i.e., just read the local values of the vector. By
+ * default, this option is disabled, so if you want to use
+ * FEEvaluationBase::read_dof_values_plain, this flag needs to be set.
*/
bool store_plain_indices;
/**
- * Option to control whether the indices
- * stored in the DoFHandler should be read and
- * the pattern for task parallelism should be
- * set up in the initialize method of
- * MatrixFree. Defaults to true. Can be
- * disabled in case the mapping should be
- * recomputed (e.g. when using a deforming
- * mesh described through MappingEulerian) but
- * the topology of cells has remained the
- * same.
+ * Option to control whether the indices stored in the DoFHandler should
+ * be read and the pattern for task parallelism should be set up in the
+ * initialize method of MatrixFree. Defaults to true. Can be disabled in
+ * case the mapping should be recomputed (e.g. when using a deforming mesh
+ * described through MappingEulerian) but the topology of cells has
+ * remained the same.
*/
bool initialize_indices;
/**
- * Option to control whether the mapping
- * information should be computed in the
- * initialize method of
- * MatrixFree. Defaults to true. Can be
- * disabled when only some indices should be
- * set up (e.g. when only a set of independent
- * cells should be computed).
+ * Option to control whether the mapping information should be computed in
+ * the initialize method of MatrixFree. Defaults to true. Can be disabled
+ * when only some indices should be set up (e.g. when only a set of
+ * independent cells should be computed).
*/
bool initialize_mapping;
};
*/
//@{
/**
- * Default empty constructor. Does
- * nothing.
+ * Default empty constructor. Does nothing.
*/
MatrixFree ();
~MatrixFree();
/**
- * Extracts the information needed to
- * perform loops over cells. The
- * DoFHandler and ConstraintMatrix
- * describe the layout of degrees of
- * freedom, the DoFHandler and the
- * mapping describe the transformations
- * from unit to real cell, and the finite
- * element underlying the DoFHandler
- * together with the quadrature formula
- * describe the local operations. Note
- * that the finite element underlying the
- * DoFHandler must either be scalar or
- * contain several copies of the same
- * element. Mixing several different
- * elements into one FESystem is not
- * allowed. In that case, use the
- * initialization function with several
- * DoFHandler arguments.
+ * Extracts the information needed to perform loops over cells. The
+ * DoFHandler and ConstraintMatrix describe the layout of degrees of freedom,
+ * the DoFHandler and the mapping describe the transformations from unit to
+ * real cell, and the finite element underlying the DoFHandler together with
+ * the quadrature formula describe the local operations. Note that the finite
+ * element underlying the DoFHandler must either be scalar or contain several
+ * copies of the same element. Mixing several different elements into one
+ * FESystem is not allowed. In that case, use the initialization function
+ * with several DoFHandler arguments.
*
- * The @p IndexSet @p locally_owned_dofs
- * is used to specify the parallel
- * partitioning with MPI. Usually, this
- * needs not be specified, and the other
- * initialization function without and @p
- * IndexSet description can be used,
- * which gets the partitioning
- * information builtin into the
- * DoFHandler.
+ * The @p IndexSet @p locally_owned_dofs is used to specify the parallel
+ * partitioning with MPI. Usually, this needs not be specified, and the other
+ * initialization function without and @p IndexSet description can be used,
+ * which gets the partitioning information builtin into the DoFHandler.
*/
template <typename DH, typename Quadrature>
void reinit (const Mapping<dim> &mapping,
const AdditionalData additional_data = AdditionalData());
/**
- * Initializes the data structures. Same
- * as above, but with index set stored in
- * the DoFHandler for describing the
- * locally owned degrees of freedom.
+ * Initializes the data structures. Same as above, but with index set stored
+ * in the DoFHandler for describing the locally owned degrees of freedom.
*/
template <typename DH, typename Quadrature>
void reinit (const Mapping<dim> &mapping,
const AdditionalData additional_data = AdditionalData());
/**
- * Initializes the data structures. Same
- * as above, but with mapping @p
+ * Initializes the data structures. Same as above, but with mapping @p
* MappingQ1.
*/
template <typename DH, typename Quadrature>
const AdditionalData additional_data = AdditionalData());
/**
- * Extracts the information needed to
- * perform loops over cells. The
- * DoFHandler and ConstraintMatrix
- * describe the layout of degrees of
- * freedom, the DoFHandler and the
- * mapping describe the transformations
- * from unit to real cell, and the finite
- * element underlying the DoFHandler
- * together with the quadrature formula
- * describe the local operations. As
- * opposed to the scalar case treated
- * with the other initialization
- * functions, this function allows for
- * problems with two or more different
- * finite elements. The DoFHandlers to
- * each element must be passed as
- * pointers to the initialization
- * function. Note that the finite element
- * underlying an DoFHandler must either
- * be scalar or contain several copies of
- * the same element. Mixing several
- * different elements into one @p
- * FE_System is not allowed.
+ * Extracts the information needed to perform loops over cells. The
+ * DoFHandler and ConstraintMatrix describe the layout of degrees of freedom,
+ * the DoFHandler and the mapping describe the transformations from unit to
+ * real cell, and the finite element underlying the DoFHandler together with
+ * the quadrature formula describe the local operations. As opposed to the
+ * scalar case treated with the other initialization functions, this function
+ * allows for problems with two or more different finite elements. The
+ * DoFHandlers to each element must be passed as pointers to the
+ * initialization function. Note that the finite element underlying an
+ * DoFHandler must either be scalar or contain several copies of the same
+ * element. Mixing several different elements into one @p FE_System is not
+ * allowed.
*
- * This function also allows for using
- * several quadrature formulas, e.g. when
- * the description contains independent
- * integrations of elements of different
- * degrees. However, the number of
- * different quadrature formulas can be
- * sets independently from the number of
- * DoFHandlers, when several elements are
- * always integrated with the same
- * quadrature formula.
+ * This function also allows for using several quadrature formulas, e.g. when
+ * the description contains independent integrations of elements of different
+ * degrees. However, the number of different quadrature formulas can be sets
+ * independently from the number of DoFHandlers, when several elements are
+ * always integrated with the same quadrature formula.
*
- * The @p IndexSet @p locally_owned_dofs
- * is used to specify the parallel
- * partitioning with MPI. Usually, this
- * needs not be specified, and the other
- * initialization function without and @p
- * IndexSet description can be used,
- * which gets the partitioning
- * information from the DoFHandler. This
- * is the most general initialization
- * function.
+ * The @p IndexSet @p locally_owned_dofs is used to specify the parallel
+ * partitioning with MPI. Usually, this needs not be specified, and the other
+ * initialization function without and @p IndexSet description can be used,
+ * which gets the partitioning information from the DoFHandler. This is the
+ * most general initialization function.
*/
template <typename DH, typename Quadrature>
void reinit (const Mapping<dim> &mapping,
const AdditionalData additional_data = AdditionalData());
/**
- * Initializes the data structures. Same
- * as before, but now the index set
- * description of the locally owned range
- * of degrees of freedom is taken from
+ * Initializes the data structures. Same as before, but now the index set
+ * description of the locally owned range of degrees of freedom is taken from
* the DoFHandler.
*/
template <typename DH, typename Quadrature>
const AdditionalData additional_data = AdditionalData());
/**
- * Initializes the data structures. Same
- * as above, but with mapping @p
+ * Initializes the data structures. Same as above, but with mapping @p
* MappingQ1.
*/
template <typename DH, typename Quadrature>
const AdditionalData additional_data = AdditionalData());
/**
- * Initializes the data structures. Same
- * as before, but now the index set
- * description of the locally owned range
- * of degrees of freedom is taken from
- * the DoFHandler. Moreover, only a
- * single quadrature formula is used, as
- * might be necessary when several
- * components in a vector-valued problem
- * are integrated together based on the
- * same quadrature formula.
+ * Initializes the data structures. Same as before, but now the index set
+ * description of the locally owned range of degrees of freedom is taken from
+ * the DoFHandler. Moreover, only a single quadrature formula is used, as
+ * might be necessary when several components in a vector-valued problem are
+ * integrated together based on the same quadrature formula.
*/
template <typename DH, typename Quadrature>
void reinit (const Mapping<dim> &mapping,
const AdditionalData additional_data = AdditionalData());
/**
- * Initializes the data structures. Same
- * as above, but with mapping @p
+ * Initializes the data structures. Same as above, but with mapping @p
* MappingQ1.
*/
template <typename DH, typename Quadrature>
const AdditionalData additional_data = AdditionalData());
/**
- * Copy function. Creates a deep copy of all
- * data structures. It is usually enough to
- * keep the data for different operations
- * once, so this function should not be needed
- * very often.
+ * Copy function. Creates a deep copy of all data structures. It is usually
+ * enough to keep the data for different operations once, so this function
+ * should not be needed very often.
*/
void copy_from (const MatrixFree<dim,Number> &matrix_free_base);
/**
- * Clears all data fields and brings the class
- * into a condition similar to after having
- * called the default constructor.
+ * Clears all data fields and brings the class into a condition similar to
+ * after having called the default constructor.
*/
void clear();
*/
//@{
/**
- * This method runs the loop over all
- * cells (in parallel) and performs
- * the MPI data exchange on the source
- * vector and destination vector. The
- * first argument indicates a function
- * object that has the following
- * signature: <code>cell_operation
- * (const MatrixFree<dim,Number>
- * &, OutVector &, InVector &,
- * std::pair<unsigned int,unsigned
- * int> &)</code>, where the first
- * argument passes the data of the
- * calling class and the last argument
- * defines the range of cells which
- * should be worked on (typically more
- * than one cell should be worked on
- * in order to reduce overheads). One
- * can pass a pointer to an object in
- * this place if it has an
- * <code>operator()</code> with the
- * correct set of arguments since such
- * a pointer can be converted to the
- * function object.
+ * This method runs the loop over all cells (in parallel) and performs the
+ * MPI data exchange on the source vector and destination vector. The first
+ * argument indicates a function object that has the following signature:
+ * <code>cell_operation (const MatrixFree<dim,Number> &, OutVector &,
+ * InVector &, std::pair<unsigned int,unsigned int> &)</code>, where the
+ * first argument passes the data of the calling class and the last argument
+ * defines the range of cells which should be worked on (typically more than
+ * one cell should be worked on in order to reduce overheads). One can pass
+ * a pointer to an object in this place if it has an <code>operator()</code>
+ * with the correct set of arguments since such a pointer can be converted
+ * to the function object.
*/
template <typename OutVector, typename InVector>
void cell_loop (const std_cxx1x::function<void (const MatrixFree<dim,Number> &,
const InVector &src) const;
/**
- * This is the second variant to run the loop
- * over all cells, now providing a function
- * pointer to a member function of class @p
- * CLASS with the signature
- * <code>cell_operation (const
- * MatrixFree<dim,Number> &, OutVector &,
- * InVector &, std::pair<unsigned int,unsigned
- * int>&)const</code>. This method obviates the need
- * to call std_cxx1x::bind to bind the class
- * into the given function in case the local
- * function needs to access data in the class
- * (i.e., it is a non-static member function).
+ * This is the second variant to run the loop over all cells, now providing
+ * a function pointer to a member function of class @p CLASS with the
+ * signature <code>cell_operation (const MatrixFree<dim,Number> &, OutVector
+ * &, InVector &, std::pair<unsigned int,unsigned int>&)const</code>. This
+ * method obviates the need to call std_cxx1x::bind to bind the class into
+ * the given function in case the local function needs to access data in the
+ * class (i.e., it is a non-static member function).
*/
template <typename CLASS, typename OutVector, typename InVector>
void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
const InVector &src) const;
/**
- * Same as above, but for class member
- * functions which are non-const.
+ * Same as above, but for class member functions which are non-const.
*/
template <typename CLASS, typename OutVector, typename InVector>
void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
const InVector &src) const;
/**
- * In the hp adaptive case, a subrange of
- * cells as computed during the cell loop
- * might contain elements of different
- * degrees. Use this function to compute what
- * the subrange for an individual finite
- * element degree is. The finite element
- * degree is associated to the vector
- * component given in the function call.
+ * In the hp adaptive case, a subrange of cells as computed during the cell
+ * loop might contain elements of different degrees. Use this function to
+ * compute what the subrange for an individual finite element degree is. The
+ * finite element degree is associated to the vector component given in the
+ * function call.
*/
std::pair<unsigned int,unsigned int>
create_cell_subrange_hp (const std::pair<unsigned int,unsigned int> &range,
const unsigned int vector_component = 0) const;
/**
- * In the hp adaptive case, a subrange of
- * cells as computed during the cell loop
- * might contain elements of different
- * degrees. Use this function to compute what
- * the subrange for a given index the hp
- * finite element, as opposed to the finite
- * element degree in the other function.
+ * In the hp adaptive case, a subrange of cells as computed during the cell
+ * loop might contain elements of different degrees. Use this function to
+ * compute what the subrange for a given index the hp finite element, as
+ * opposed to the finite element degree in the other function.
*/
std::pair<unsigned int,unsigned int>
create_cell_subrange_hp_by_index (const std::pair<unsigned int,unsigned int> &range,
*/
//@{
/**
- * Initialize function for a general
- * vector. The length of the vector is equal
- * to the total number of degrees in the
- * DoFHandler. If the vector is of class
- * parallel::distributed::Vector@<Number@>, the ghost
- * entries are set accordingly. For
- * vector-valued problems with several
- * DoFHandlers underlying this class, the
- * parameter @p vector_component defines which
+ * Initialize function for a general vector. The length of the vector is
+ * equal to the total number of degrees in the DoFHandler. If the vector is
+ * of class parallel::distributed::Vector@<Number@>, the ghost entries are
+ * set accordingly. For vector-valued problems with several DoFHandlers
+ * underlying this class, the parameter @p vector_component defines which
* component is to be used.
*/
template <typename VectorType>
const unsigned int vector_component=0) const;
/**
- * Initialize function for a distributed
- * vector. The length of the vector is equal
- * to the total number of degrees in the
- * DoFHandler. If the vector is of class
- * parallel::distributed::Vector@<Number@>, the ghost
- * entries are set accordingly. For
- * vector-valued problems with several
- * DoFHandlers underlying this class, the
- * parameter @p vector_component defines which
+ * Initialize function for a distributed vector. The length of the vector is
+ * equal to the total number of degrees in the DoFHandler. If the vector is
+ * of class parallel::distributed::Vector@<Number@>, the ghost entries are
+ * set accordingly. For vector-valued problems with several DoFHandlers
+ * underlying this class, the parameter @p vector_component defines which
* component is to be used.
*/
template <typename Number2>
const unsigned int vector_component=0) const;
/**
- * Returns the partitioner that represents the
- * locally owned data and the ghost indices
- * where access is needed to for the cell
- * loop. The partitioner is constructed from
- * the locally owned dofs and ghost dofs given
- * by the respective fields. If you want to
- * have specific information about these
- * objects, you can query them with the
- * respective access functions. If you just
- * want to initialize a (parallel) vector, you
- * should usually prefer this data structure
- * as the data exchange information can be
- * reused from one vector to another.
+ * Returns the partitioner that represents the locally owned data and the
+ * ghost indices where access is needed to for the cell loop. The
+ * partitioner is constructed from the locally owned dofs and ghost dofs
+ * given by the respective fields. If you want to have specific information
+ * about these objects, you can query them with the respective access
+ * functions. If you just want to initialize a (parallel) vector, you should
+ * usually prefer this data structure as the data exchange information can
+ * be reused from one vector to another.
*/
const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &
get_vector_partitioner (const unsigned int vector_component=0) const;
/**
- * Returns the set of cells that are
- * oned by the processor.
+ * Returns the set of cells that are oned by the processor.
*/
const IndexSet &
get_locally_owned_set (const unsigned int fe_component = 0) const;
/**
- * Returns the set of ghost cells
- * needed but not owned by the
- * processor.
+ * Returns the set of ghost cells needed but not owned by the processor.
*/
const IndexSet &
get_ghost_set (const unsigned int fe_component = 0) const;
/**
- * Returns a list of all degrees of freedom
- * that are constrained. The list is returned
- * in local index space for the locally owned
- * range of the vector, not in global
- * numbers. In addition, it only returns the
- * indices for degrees of freedom that are
- * owned locally, not for ghosts.
+ * Returns a list of all degrees of freedom that are constrained. The list
+ * is returned in local index space for the locally owned range of the
+ * vector, not in global numbers. In addition, it only returns the indices
+ * for degrees of freedom that are owned locally, not for ghosts.
*/
const std::vector<unsigned int> &
get_constrained_dofs (const unsigned int fe_component = 0) const;
/**
- * Calls renumber_dofs function in dof
- * info which renumbers the the
- * degrees of freedom according to the
- * ordering for parallelization.
+ * Calls renumber_dofs function in dof info which renumbers the the degrees
+ * of freedom according to the ordering for parallelization.
*/
void renumber_dofs (std::vector<unsigned int> &renumbering,
const unsigned int vector_component = 0);
*/
//@{
/**
- * Returns the number of different DoFHandlers
- * specified at initialization.
+ * Returns the number of different DoFHandlers specified at initialization.
*/
unsigned int n_components () const;
/**
- * Returns the number of cells this structure
- * is based on. If you are using a usual
- * DoFHandler, it corresponds to the number of
- * (locally owned) active cells. Note that
- * most data structures in this class do not
- * directly act on this number but rather on
- * n_macro_cells() which gives the number of
- * cells as seen when lumping several cells
- * together with vectorization.
+ * Returns the number of cells this structure is based on. If you are using
+ * a usual DoFHandler, it corresponds to the number of (locally owned)
+ * active cells. Note that most data structures in this class do not
+ * directly act on this number but rather on n_macro_cells() which gives the
+ * number of cells as seen when lumping several cells together with
+ * vectorization.
*/
unsigned int n_physical_cells () const;
/**
- * Returns the number of macro cells
- * that this structure works on, i.e.,
- * the number of cell chunks that are
- * worked on after the application of
- * vectorization which in general
- * works on several cells at once. The
- * cell range in @p cell_loop runs
- * from zero to n_macro_cells()
- * (exclusive), so this is the
- * appropriate size if you want to
- * store arrays of data for all cells
- * to be worked on. This number is
- * approximately
- * n_physical_cells()/VectorizedArray::n_array_elements
- * (depending on how many cell chunks
- * that do not get filled up
- * completely).
+ * Returns the number of macro cells that this structure works on, i.e., the
+ * number of cell chunks that are worked on after the application of
+ * vectorization which in general works on several cells at once. The cell
+ * range in @p cell_loop runs from zero to n_macro_cells() (exclusive), so
+ * this is the appropriate size if you want to store arrays of data for all
+ * cells to be worked on. This number is approximately
+ * n_physical_cells()/VectorizedArray::n_array_elements (depending on how
+ * many cell chunks that do not get filled up completely).
*/
unsigned int n_macro_cells () const;
/**
- * In case this structure was built based on a
- * DoFHandler, this returns the
- * DoFHandler. Note that this function returns
- * an exception in case the structure was
- * based on MGDoFHandler and a level has been
- * specified in InitializationOption.
+ * In case this structure was built based on a DoFHandler, this returns the
+ * DoFHandler. Note that this function returns an exception in case the
+ * structure was based on MGDoFHandler and a level has been specified in
+ * InitializationOption.
*/
const DoFHandler<dim> &
get_dof_handler (const unsigned int fe_component = 0) const;
/**
- * In case this structure was built based on a
- * DoFHandler, this returns the
- * DoFHandler. Note that this function returns
- * an exception in case the structure was
- * based on MGDoFHandler and a level has been
- * specified in InitializationOption.
+ * In case this structure was built based on a DoFHandler, this returns the
+ * DoFHandler. Note that this function returns an exception in case the
+ * structure was based on MGDoFHandler and a level has been specified in
+ * InitializationOption.
*/
const MGDoFHandler<dim> &
get_mg_dof_handler (const unsigned int fe_component = 0) const;
/**
- * This returns the cell iterator in deal.II
- * speak to a given cell in the renumbering of
- * this structure. This function returns an
- * exception in case the structure was
- * constructed based on an MGDoFHandler with
- * level specified, as these cells are in
- * general not active.
+ * This returns the cell iterator in deal.II speak to a given cell in the
+ * renumbering of this structure. This function returns an exception in case
+ * the structure was constructed based on an MGDoFHandler with level
+ * specified, as these cells are in general not active.
*
- * Note that the cell iterators in deal.II go
- * through cells differently to what the cell
- * loop of this class does. This is because
- * several cells are worked on together
- * (vectorization), and since cells with
- * neighbors on different MPI processors need
- * to be accessed at a certain time when
- * accessing remote data and overlapping
- * communication with computation.
+ * Note that the cell iterators in deal.II go through cells differently to
+ * what the cell loop of this class does. This is because several cells are
+ * worked on together (vectorization), and since cells with neighbors on
+ * different MPI processors need to be accessed at a certain time when
+ * accessing remote data and overlapping communication with computation.
*/
typename DoFHandler<dim>::active_cell_iterator
get_cell_iterator (const unsigned int macro_cell_number,
const unsigned int fe_component = 0) const;
/**
- * This returns the cell iterator in deal.II
- * speak to a given cell in the renumbering of
- * this structure. This function returns an
- * exception in case the structure was
- * constructed based on a DoFHandler and not
- * MGDoFHandler.
+ * This returns the cell iterator in deal.II speak to a given cell in the
+ * renumbering of this structure. This function returns an exception in case
+ * the structure was constructed based on a DoFHandler and not MGDoFHandler.
*
- * Note that the cell iterators in deal.II go
- * through cells differently to what the cell
- * loop of this class does. This is because
- * several cells are worked on together
- * (vectorization), and since cells with
- * neighbors on different MPI processors need
- * to be accessed at a certain time when
- * accessing remote data and overlapping
- * communication with computation.
+ * Note that the cell iterators in deal.II go through cells differently to
+ * what the cell loop of this class does. This is because several cells are
+ * worked on together (vectorization), and since cells with neighbors on
+ * different MPI processors need to be accessed at a certain time when
+ * accessing remote data and overlapping communication with computation.
*/
typename MGDoFHandler<dim>::cell_iterator
get_mg_cell_iterator (const unsigned int macro_cell_number,
const unsigned int fe_component = 0) const;
/**
- * This returns the cell iterator in deal.II
- * speak to a given cell in the renumbering of
- * this structure. This function returns an
- * exception in case the structure was
- * not constructed based on an hp::DoFHandler.
+ * This returns the cell iterator in deal.II speak to a given cell in the
+ * renumbering of this structure. This function returns an exception in case
+ * the structure was not constructed based on an hp::DoFHandler.
*
- * Note that the cell iterators in deal.II go
- * through cells differently to what the cell
- * loop of this class does. This is because
- * several cells are worked on together
- * (vectorization), and since cells with
- * neighbors on different MPI processors need
- * to be accessed at a certain time when
- * accessing remote data and overlapping
- * communication with computation.
+ * Note that the cell iterators in deal.II go through cells differently to
+ * what the cell loop of this class does. This is because several cells are
+ * worked on together (vectorization), and since cells with neighbors on
+ * different MPI processors need to be accessed at a certain time when
+ * accessing remote data and overlapping communication with computation.
*/
typename hp::DoFHandler<dim>::active_cell_iterator
get_hp_cell_iterator (const unsigned int macro_cell_number,
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
const unsigned int hp_active_fe_index = 0) const;
/**
- * Returns the quadrature rule for
- * given hp index.
+ * Returns the quadrature rule for given hp index.
*/
const Quadrature<dim> &
get_quadrature (const unsigned int quad_index = 0,
const unsigned int hp_active_fe_index = 0) const;
/**
- * Returns the quadrature rule for
- * given hp index.
+ * Returns the quadrature rule for given hp index.
*/
const Quadrature<dim-1> &
get_face_quadrature (const unsigned int quad_index = 0,
const unsigned int hp_active_fe_index = 0) const;
/**
- * Queries whether or not the
- * indexation has been set.
+ * Queries whether or not the indexation has been set.
*/
bool indices_initialized () const;
/**
- * Queries whether or not the
- * geometry-related information for
- * the cells has been set.
+ * Queries whether or not the geometry-related information for the cells has
+ * been set.
*/
bool mapping_initialized () const;
/**
- * Returns an approximation of the memory
- * consumption of this class in bytes.
+ * Returns an approximation of the memory consumption of this class in
+ * bytes.
*/
std::size_t memory_consumption() const;
/**
- * Prints a detailed summary of memory
- * consumption in the different structures of
- * this class to the given output stream.
+ * Prints a detailed summary of memory consumption in the different
+ * structures of this class to the given output stream.
*/
template <typename STREAM>
void print_memory_consumption(STREAM &out) const;
/**
- * Prints a summary of this class to the given
- * output stream. It is focused on the
- * indices, and does not print all the data
- * stored.
+ * Prints a summary of this class to the given output stream. It is focused
+ * on the indices, and does not print all the data stored.
*/
void print (std::ostream &out) const;
get_size_info () const;
/*
- * Returns geometry-dependent
- * information on the cells.
+ * Returns geometry-dependent information on the cells.
*/
const internal::MatrixFreeFunctions::MappingInfo<dim,Number> &
get_mapping_info () const;
/**
- * Returns information on indexation
- * degrees of freedom.
+ * Returns information on indexation degrees of freedom.
*/
const internal::MatrixFreeFunctions::DoFInfo &
get_dof_info (const unsigned int fe_component = 0) const;
/**
- * Returns the number of weights in the
- * constraint pool.
+ * Returns the number of weights in the constraint pool.
*/
unsigned int n_constraint_pool_entries() const;
/**
- * Returns a pointer to the first
- * number in the constraint pool data
- * with index @p pool_index (to
- * be used together with @p
- * constraint_pool_end()).
+ * Returns a pointer to the first number in the constraint pool data with
+ * index @p pool_index (to be used together with @p constraint_pool_end()).
*/
const Number *
constraint_pool_begin (const unsigned int pool_index) const;
/**
- * Returns a pointer to one past the
- * last number in the constraint pool
- * data with index @p pool_index (to
- * be used together with @p
+ * Returns a pointer to one past the last number in the constraint pool data
+ * with index @p pool_index (to be used together with @p
* constraint_pool_begin()).
*/
const Number *
constraint_pool_end (const unsigned int pool_index) const;
/**
- * Returns the unit cell information
- * for given hp index.
+ * Returns the unit cell information for given hp index.
*/
const internal::MatrixFreeFunctions::ShapeInfo<Number> &
get_shape_info (const unsigned int fe_component = 0,
private:
/**
- * This is the actual reinit function
- * that sets up the indices for the
+ * This is the actual reinit function that sets up the indices for the
* DoFHandler and MGDoFHandler case.
*/
template <typename DoFHandler>
const AdditionalData additional_data);
/**
- * Same as before but for hp::DoFHandler
- * instead of generic DoFHandler type.
+ * Same as before but for hp::DoFHandler instead of generic DoFHandler type.
*/
void internal_reinit (const Mapping<dim> &mapping,
const std::vector<const hp::DoFHandler<dim>*> &dof_handler,
const AdditionalData additional_data);
/**
- * Initializes the fields in DoFInfo
- * together with the constraint pool
- * that holds all different weights in
- * the constraints (not part of
- * DoFInfo because several DoFInfo
- * classes can have the same weights
- * which consequently only need to be
- * stored once).
+ * Initializes the fields in DoFInfo together with the constraint pool that
+ * holds all different weights in the constraints (not part of DoFInfo
+ * because several DoFInfo classes can have the same weights which
+ * consequently only need to be stored once).
*/
void
initialize_indices (const std::vector<const ConstraintMatrix *> &constraint,
const std::vector<IndexSet> &locally_owned_set);
/**
- * Initializes the DoFHandlers based on a
- * DoFHandler<dim> argument.
+ * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
*/
void initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handlers,
const unsigned int level);
/**
- * Initializes the DoFHandlers based on a
- * DoFHandler<dim> argument.
+ * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
*/
void initialize_dof_handlers (const std::vector<const MGDoFHandler<dim>*> &dof_handlers,
const unsigned int level);
/**
- * Initializes the DoFHandlers based on a
- * DoFHandler<dim> argument.
+ * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
*/
void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
const unsigned int level);
/**
- * This struct defines which DoFHandler has
- * actually been given at construction, in
- * order to define the correct behavior when
- * querying the underlying DoFHandler.
+ * This struct defines which DoFHandler has actually been given at
+ * construction, in order to define the correct behavior when querying the
+ * underlying DoFHandler.
*/
struct DoFHandlers
{
};
/**
- * Pointers to the DoFHandlers underlying the
- * current problem.
+ * Pointers to the DoFHandlers underlying the current problem.
*/
DoFHandlers dof_handlers;
/**
- * Contains the information about degrees of
- * freedom on the individual cells and
- * constraints.
+ * Contains the information about degrees of freedom on the individual cells
+ * and constraints.
*/
std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
/**
- * Contains the weights for constraints stored
- * in DoFInfo. Filled into a separate field
- * since several vector components might share
- * similar weights, which reduces memory
- * consumption. Moreover, it obviates template
- * arguments on DoFInfo and keeps it a plain
- * field of indices only.
+ * Contains the weights for constraints stored in DoFInfo. Filled into a
+ * separate field since several vector components might share similar
+ * weights, which reduces memory consumption. Moreover, it obviates template
+ * arguments on DoFInfo and keeps it a plain field of indices only.
*/
std::vector<Number> constraint_pool_data;
/**
- * Contains an indicator to the start
- * of the ith index in the constraint
+ * Contains an indicator to the start of the ith index in the constraint
* pool data.
*/
std::vector<unsigned int> constraint_pool_row_index;
/**
- * Holds information on transformation of
- * cells from reference cell to real cell that
- * is needed for evaluating integrals.
+ * Holds information on transformation of cells from reference cell to real
+ * cell that is needed for evaluating integrals.
*/
internal::MatrixFreeFunctions::MappingInfo<dim,Number> mapping_info;
/**
- * Contains shape value information on the
- * unit cell.
+ * Contains shape value information on the unit cell.
*/
Table<4,internal::MatrixFreeFunctions::ShapeInfo<Number> > shape_info;
/**
- * Describes how the cells are gone
- * through. With the cell level (first index
- * in this field) and the index within the
- * level, one can reconstruct a deal.II
- * cell iterator and use all the traditional
- * things deal.II offers to do with cell
- * iterators.
+ * Describes how the cells are gone through. With the cell level (first
+ * index in this field) and the index within the level, one can reconstruct
+ * a deal.II cell iterator and use all the traditional things deal.II offers
+ * to do with cell iterators.
*/
std::vector<std::pair<unsigned int,unsigned int> > cell_level_index;
/**
- * Stores how many cells we have, how many
- * cells that we see after applying
- * vectorization (i.e., the number of macro
- * cells), and MPI-related stuff.
+ * Stores how many cells we have, how many cells that we see after applying
+ * vectorization (i.e., the number of macro cells), and MPI-related stuff.
*/
internal::MatrixFreeFunctions::SizeInfo size_info;
/**
- * Information regarding the shared memory
- * parallelization.
+ * Information regarding the shared memory parallelization.
*/
internal::MatrixFreeFunctions::TaskInfo task_info;
/**
- * Stores whether indices have been
- * initialized.
+ * Stores whether indices have been initialized.
*/
bool indices_are_initialized;
/**
- * Stores whether indices have been
- * initialized.
+ * Stores whether indices have been initialized.
*/
bool mapping_is_initialized;
};
MatrixFree<dim,Number>::constraint_pool_begin (const unsigned int row) const
{
AssertIndexRange (row, constraint_pool_row_index.size()-1);
- return &constraint_pool_data[0] + constraint_pool_row_index[row];
+ return constraint_pool_data.empty() ? 0 :
+ &constraint_pool_data[0] + constraint_pool_row_index[row];
}
MatrixFree<dim,Number>::constraint_pool_end (const unsigned int row) const
{
AssertIndexRange (row, constraint_pool_row_index.size()-1);
- return &constraint_pool_data[0] + constraint_pool_row_index[row+1];
+ return constraint_pool_data.empty() ? 0 :
+ &constraint_pool_data[0] + constraint_pool_row_index[row+1];
}
return range;
else
{
- // the range over which we are searching must be ordered,
- // otherwise we got a range that spans over too many cells
+ // the range over which we are searching must be ordered, otherwise we
+ // got a range that spans over too many cells
#ifdef DEBUG
for (unsigned int i=range.first+1; i<range.second; ++i)
Assert (fe_indices[i] >= fe_indices[i-1],
#endif
std::pair<unsigned int,unsigned int> return_range;
return_range.first =
- std::lower_bound(&fe_indices[range.first], &fe_indices[range.second],
- fe_index)
+ std::lower_bound(&fe_indices[0] + range.first,
+ &fe_indices[0] + range.second, fe_index)
-&fe_indices[0] ;
return_range.second =
- std::lower_bound(&fe_indices[return_range.first],
- &fe_indices[range.second],
- fe_index+1)-&fe_indices[0];
+ std::lower_bound(&fe_indices[0] + return_range.first,
+ &fe_indices[0] + range.second,
+ fe_index + 1)-&fe_indices[0];
Assert(return_range.first >= range.first &&
return_range.second <= range.second, ExcInternalError());
return return_range;
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];
}
}
const std::vector<Quad> &quad,
const MatrixFree<dim,Number>::AdditionalData additional_data)
{
- // find out whether we use a hp Quadrature or
- // a standard quadrature
+ // find out whether we use a hp Quadrature or a standard quadrature
std::vector<hp::QCollection<1> > quad_hp;
for (unsigned int q=0; q<quad.size(); ++q)
quad_hp.push_back (hp::QCollection<1>(quad[q]));
// ------------------------------ implementation of cell_loop ---------------
-// internal helper functions that define how
-// to call MPI data exchange functions: for
-// generic vectors, do nothing at all. For
-// distributed vectors, can call
-// update_ghost_values_start function and so on. If we
-// have collections of vectors, just do the
-// individual functions of the components
+// internal helper functions that define how to call MPI data exchange
+// functions: for generic vectors, do nothing at all. For distributed vectors,
+// can call update_ghost_values_start function and so on. If we have
+// collections of vectors, just do the individual functions of the components
namespace internal
{
template<typename VectorStruct>
#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
{
#if DEAL_II_USE_MT==1
- // Use multithreading if so requested and if
- // there is enough work to do in parallel (the
- // code might hang if there are less than two
- // chunks!)
+ // Use multithreading if so requested and if there is enough work to do in
+ // parallel (the code might hang if there are less than two chunks!)
if (task_info.use_multithreading == true && task_info.n_blocks > 3)
{
- // to simplify the function calls, bind away
- // all arguments except the cell range
+ // to simplify the function calls, bind away all arguments except the
+ // cell range
typedef
std_cxx1x::function<void (const std::pair<unsigned int,unsigned int> &range)>
Worker;
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;
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());
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;
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;
OutVector &dst,
const InVector &src) const
{
- // here, use std_cxx1x::bind to hand a
- // function handler with the appropriate
+ // here, use std_cxx1x::bind to hand a function handler with the appropriate
// argument to the other loop function
std_cxx1x::function<void (const MatrixFree<dim,Number> &,
OutVector &,
OutVector &dst,
const InVector &src) const
{
- // here, use std_cxx1x::bind to hand a
- // function handler with the appropriate
+ // here, use std_cxx1x::bind to hand a function handler with the appropriate
// argument to the other loop function
std_cxx1x::function<void (const MatrixFree<dim,Number> &,
OutVector &,
for (unsigned int no=0; no<dof_handler.size(); ++no)
dof_info[no].store_plain_indices = additional_data.store_plain_indices;
- // initialize the basic multithreading
- // information that needs to be passed to the
- // DoFInfo structure
+ // initialize the basic multithreading information that needs to be
+ // passed to the DoFInfo structure
#if DEAL_II_USE_MT == 1
if (additional_data.tasks_parallel_scheme != AdditionalData::none)
{
#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));
}
}
- // 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,
for (unsigned int no=0; no<dof_handler.size(); ++no)
dof_info[no].store_plain_indices = additional_data.store_plain_indices;
- // initialize the basic multithreading
- // information that needs to be passed to the
- // DoFInfo structure
+ // initialize the basic multithreading information that needs to be
+ // passed to the DoFInfo structure
#if DEAL_II_USE_MT == 1
if (additional_data.tasks_parallel_scheme != AdditionalData::none)
{
#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;
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,
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;
dof_info.resize (dof_handlers.n_dof_handlers);
- // go through cells on zeroth level and then
- // successively step down into children. This
- // gives a z-ordering of the cells, which is
- // beneficial when setting up neighboring
- // relations between cells for thread
- // parallelization
+ // go through cells on zeroth level and then successively step down into
+ // children. This gives a z-ordering of the cells, which is beneficial when
+ // setting up neighboring relations between cells for thread parallelization
const unsigned int n_mpi_procs = size_info.n_procs;
const unsigned int my_pid = size_info.my_pid;
- // if we have no level given, use the same as
- // for the standard DoFHandler, otherwise we
- // must loop through the respective level
+ // if we have no level given, use the same as for the standard DoFHandler,
+ // otherwise we must loop through the respective level
const Triangulation<dim> &tria = dof_handlers.mg_dof_handler[0]->get_tria();
if (level == numbers::invalid_unsigned_int)
dof_info.resize (dof_handlers.n_dof_handlers);
- // go through cells on zeroth level and then
- // successively step down into children. This
- // gives a z-ordering of the cells, which is
- // beneficial when setting up neighboring
- // relations between cells for thread
- // parallelization
+ // go through cells on zeroth level and then successively step down into
+ // children. This gives a z-ordering of the cells, which is beneficial when
+ // setting up neighboring relations between cells for thread parallelization
const unsigned int n_mpi_procs = size_info.n_procs;
const unsigned int my_pid = size_info.my_pid;
- // if we have no level given, use the same as
- // for the standard DoFHandler, otherwise we
- // must loop through the respective level
+ // if we have no level given, use the same as for the standard DoFHandler,
+ // otherwise we must loop through the respective level
const Triangulation<dim> &tria = dof_handler[0]->get_tria();
if (n_mpi_procs == 1)
ExcMessage ("MatrixFree only works for DoFHandler with one base element"));
const unsigned int n_fe_components = fe.element_multiplicity (0);
- // cache number of finite elements and
- // dofs_per_cell
+ // cache number of finite elements and dofs_per_cell
dof_info[no].dofs_per_cell.push_back (fe.dofs_per_cell);
dof_info[no].dofs_per_face.push_back (fe.dofs_per_face);
dof_info[no].n_components = n_fe_components;
- // get permutation that gives lexicographic
- // renumbering of the cell dofs
- // renumber (this is necessary for FE_Q, for
- // example, since there the vertex DoFs come
- // first, which is incompatible with the
- // lexicographic ordering necessary to apply
- // tensor products efficiently)
+ // get permutation that gives lexicographic renumbering of the cell
+ // dofs renumber (this is necessary for FE_Q, for example, since
+ // there the vertex DoFs come first, which is incompatible with the
+ // lexicographic ordering necessary to apply tensor products
+ // efficiently)
const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>
(&fe.base_element(0));
- // This class currently only works for
- // elements derived from
- // FE_Poly<TensorProductPolynomials<dim>,dim,dim>.
- // For any other element, the dynamic cast
- // above will fail and give fe_poly == 0.
+ // This class currently only works for elements derived from
+ // FE_Poly<TensorProductPolynomials<dim>,dim,dim>. For any other
+ // element, the dynamic cast above will fail and give fe_poly == 0.
Assert (fe_poly != 0, ExcNotImplemented());
if (n_fe_components == 1)
{
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,
ghost_dofs[no].reserve (locally_owned_set[no].n_elements()/10+1);
}
- // extract all the global indices associated
- // with the computation, and form the ghost
- // indices
+ // extract all the global indices associated with the computation, and form
+ // the ghost indices
std::vector<unsigned int> boundary_cells;
for (unsigned int counter = 0 ; counter < n_active_cells ; ++counter)
{
bool cell_at_boundary = false;
for (unsigned int no=0; no<n_fe; ++no)
{
- // OK, read indices from standard DoFHandler
- // or active indices in MGDoFHandler. That is
- // the usual stuff
+ // OK, read indices from standard DoFHandler or active indices in
+ // MGDoFHandler. That is the usual stuff
if (dof_handlers.active_dof_handler == DoFHandlers::usual ||
(dof_handlers.active_dof_handler == DoFHandlers::multigrid &&
dof_handlers.level == numbers::invalid_unsigned_int))
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)
{
}
}
- // 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);
}
for (unsigned int no=0; no<n_fe; ++no)
dof_info[no].assign_ghosts (boundary_cells);
- // reorganize the indices in order to overlap
- // communication in MPI with computations:
- // Place all cells with ghost indices into one
- // chunk. Also reorder cells so that we can
- // parallelize by threads
+ // reorganize the indices in order to overlap communication in MPI with
+ // computations: Place all cells with ghost indices into one chunk. Also
+ // reorder cells so that we can parallelize by threads
std::vector<unsigned int> renumbering;
if (task_info.use_multithreading == true)
{
}
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)
irregular_cells);
}
- // Finally perform the renumbering. We also
- // want to group several cells together to one
- // "macro-cell" for vectorization (where the
- // arithmetic operations will then be done
- // simultaneously).
+ // Finally perform the renumbering. We also want to group several cells
+ // together to one "macro-cell" for vectorization (where the arithmetic
+ // operations will then be done simultaneously).
#ifdef DEBUG
{
std::vector<unsigned int> sorted_renumbering (renumbering);
cell_level_index.push_back
(cell_level_index_old[renumbering[position_cell+j]]);
- // generate a cell and level index
- // also when we have not filled up
- // vectorization_length cells. This is needed for
- // MappingInfo when the transformation
- // data is initialized. We just set
- // the value to the last valid cell in
- // that case.
+ // generate a cell and level index also when we have not filled up
+ // vectorization_length cells. This is needed for MappingInfo when the
+ // transformation data is initialized. We just set the value to the
+ // last valid cell in that case.
for (unsigned int j=n_comp; j<vectorization_length; ++j)
cell_level_index.push_back
(cell_level_index_old[renumbering[position_cell+n_comp-1]]);
AssertDimension (cell_level_index.size(),size_info.n_macro_cells*vectorization_length);
}
- // set constraint pool from the std::map and
- // reorder the indices
+ // set constraint pool from the std::map and reorder the indices
typename std::map<std::vector<double>, unsigned int,
internal::MatrixFreeFunctions::FPArrayComparator<double> >::iterator
it = constraint_values.constraints.begin(),
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);
(vectorization_length - n_boundary_cells%vectorization_length)%vectorization_length;
if (fillup_needed > 0 && n_boundary_cells < n_active_cells)
{
- // fill additional cells into the list of
- // boundary cells to get a balanced number. Go
- // through the indices successively until we
+ // fill additional cells into the list of boundary cells to get a
+ // balanced number. Go through the indices successively until we
// found enough indices
std::vector<unsigned int> new_boundary_cells;
new_boundary_cells.reserve (n_max_boundary_cells);
{
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;
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;
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;
+ /* ------------------------------------------------------------------ */
+
template <typename Number>
FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
:
projects / dealii-svn.git / commitdiff