Assert (sparsity.n_cols() == n_dofs,
ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
- // If we have a distributed::Triangulation only
- // allow locally_owned subdomain. Not setting a
- // subdomain is also okay, because we skip ghost
- // cells in the loop below.
+ // If we have a distributed::Triangulation only allow locally_owned
+ // subdomain. Not setting a subdomain is also okay, because we skip
+ // ghost cells in the loop below.
Assert (
(dof.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id)
||
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
- // In case we work with a distributed
- // sparsity pattern of Trilinos type, we
- // only have to do the work if the
- // current cell is owned by the calling
- // processor. Otherwise, just continue.
+ // In case we work with a distributed sparsity pattern of Trilinos
+ // type, we only have to do the work if the current cell is owned by
+ // the calling processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
if (((subdomain_id == types::invalid_subdomain_id)
||
dofs_on_this_cell.resize (dofs_per_cell);
cell->get_dof_indices (dofs_on_this_cell);
- // make sparsity pattern for this
- // cell. if no constraints pattern was
- // given, then the following call acts
- // as if simply no constraints existed
+ // make sparsity pattern for this cell. if no constraints pattern
+ // was given, then the following call acts as if simply no
+ // constraints existed
constraints.add_entries_local_to_global (dofs_on_this_cell,
sparsity,
keep_constrained_dofs);
Assert (couplings.n_cols() == dof.get_fe().n_components(),
ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components()));
- // If we have a distributed::Triangulation only
- // allow locally_owned subdomain. Not setting a
- // subdomain is also okay, because we skip ghost
- // cells in the loop below.
+ // If we have a distributed::Triangulation only allow locally_owned
+ // subdomain. Not setting a subdomain is also okay, because we skip
+ // ghost cells in the loop below.
Assert (
(dof.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id)
||
const hp::FECollection<DH::dimension,DH::space_dimension> fe_collection (dof.get_fe());
- // first, for each finite element, build a
- // mask for each dof, not like the one
- // given which represents components. make
- // sure we do the right thing also with
- // respect to non-primitive shape
- // functions, which takes some additional
- // thought
+ // first, for each finite element, build a mask for each dof, not like
+ // the one given which represents components. make sure we do the right
+ // thing also with respect to non-primitive shape functions, which
+ // takes some additional thought
std::vector<Table<2,bool> > dof_mask(fe_collection.size());
- // check whether the table of couplings
- // contains only true arguments, i.e., we
- // do not exclude any index. that is the
- // easy case, since we don't have to set
- // up the tables
+ // check whether the table of couplings contains only true arguments,
+ // i.e., we do not exclude any index. that is the easy case, since we
+ // don't have to set up the tables
bool need_dof_mask = false;
for (unsigned int i=0; i<couplings.n_rows(); ++i)
for (unsigned int j=0; j<couplings.n_cols(); ++j)
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
- // In case we work with a distributed
- // sparsity pattern of Trilinos type, we
- // only have to do the work if the
- // current cell is owned by the calling
- // processor. Otherwise, just continue.
+ // In case we work with a distributed sparsity pattern of Trilinos
+ // type, we only have to do the work if the current cell is owned by
+ // the calling processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
if (((subdomain_id == types::invalid_subdomain_id)
||
cell->get_dof_indices (dofs_on_this_cell);
- // make sparsity pattern for this
- // cell. if no constraints pattern was
- // given, then the following call acts
- // as if simply no constraints existed
+ // make sparsity pattern for this cell. if no constraints pattern
+ // was given, then the following call acts as if simply no
+ // constraints existed
constraints.add_entries_local_to_global (dofs_on_this_cell,
sparsity,
keep_constrained_dofs,
{
if (DH::dimension == 1)
{
- // there are only 2 boundary
- // indicators in 1d, so it is no
- // performance problem to call the
- // other function
+ // there are only 2 boundary indicators in 1d, so it is no
+ // performance problem to call the other function
typename DH::FunctionMap boundary_indicators;
boundary_indicators[0] = 0;
boundary_indicators[1] = 0;
std::vector<unsigned int> dofs_on_this_face;
dofs_on_this_face.reserve (max_dofs_per_face(dof));
- // loop over all faces to check
- // whether they are at a
- // boundary. note that we need not
- // take special care of single
- // lines (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every boundary line is
- // also part of a boundary face.
+ // loop over all faces to check whether they are at a boundary. note
+ // that we need not take special care of single lines (using
+ // @p{cell->has_boundary_lines}), since we do not support boundaries of
+ // dimension dim-2, and so every boundary line is also part of a
+ // boundary face.
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
for (; cell!=endc; ++cell)
{
if (DH::dimension == 1)
{
- // first check left, then right
- // boundary point
+ // first check left, then right boundary point
for (unsigned int direction=0; direction<2; ++direction)
{
- // if this boundary is not
- // requested, then go on with next one
+ // if this boundary is not requested, then go on with next one
if (boundary_indicators.find(direction) ==
boundary_indicators.end())
continue;
- // find active cell at that
- // boundary: first go to
- // left/right, then to children
+ // find active cell at that boundary: first go to left/right,
+ // then to children
typename DH::cell_iterator cell = dof.begin(0);
while (!cell->at_boundary(direction))
cell = cell->neighbor(direction);
const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
std::vector<unsigned int> boundary_dof_boundary_indices (dofs_per_vertex);
- // next get boundary mapped dof
- // indices of boundary dofs
+ // next get boundary mapped dof indices of boundary dofs
for (unsigned int i=0; i<dofs_per_vertex; ++i)
boundary_dof_boundary_indices[i]
= dof_to_boundary_mapping[cell->vertex_dof_index(direction,i)];
template <class DH, class SparsityPattern>
void
- make_flux_sparsity_pattern (const DH &dof,
- SparsityPattern &sparsity,
- const ConstraintMatrix &constraints,
- const bool keep_constrained_dofs,
+ make_flux_sparsity_pattern (const DH &dof,
+ SparsityPattern &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
const types::subdomain_id subdomain_id)
+
+ // TODO: QA: reduce the indentation level of this method..., Maier 2012
+
{
const unsigned int n_dofs = dof.n_dofs();
AssertDimension (sparsity.n_rows(), n_dofs);
AssertDimension (sparsity.n_cols(), n_dofs);
- // If we have a distributed::Triangulation only
- // allow locally_owned subdomain. Not setting a
- // subdomain is also okay, because we skip ghost
- // cells in the loop below.
+ // If we have a distributed::Triangulation only allow locally_owned
+ // subdomain. Not setting a subdomain is also okay, because we skip
+ // ghost cells in the loop below.
Assert (
(dof.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id)
||
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
- // TODO: in an old implementation, we used
- // user flags before to tag faces that were
- // already touched. this way, we could reduce
- // the work a little bit. now, we instead add
- // only data from one side. this should be OK,
- // but we need to actually verify it.
-
- // In case we work with a distributed
- // sparsity pattern of Trilinos type, we
- // only have to do the work if the
- // current cell is owned by the calling
- // processor. Otherwise, just continue.
+ // TODO: in an old implementation, we used user flags before to tag
+ // faces that were already touched. this way, we could reduce the work
+ // a little bit. now, we instead add only data from one side. this
+ // should be OK, but we need to actually verify it.
+
+ // In case we work with a distributed sparsity pattern of Trilinos
+ // type, we only have to do the work if the current cell is owned by
+ // the calling processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
if (((subdomain_id == types::invalid_subdomain_id)
||
dofs_on_this_cell.resize (n_dofs_on_this_cell);
cell->get_dof_indices (dofs_on_this_cell);
- // make sparsity pattern for this
- // cell. if no constraints pattern was
- // given, then the following call acts
- // as if simply no constraints existed
+ // make sparsity pattern for this cell. if no constraints pattern
+ // was given, then the following call acts as if simply no
+ // constraints existed
constraints.add_entries_local_to_global (dofs_on_this_cell,
sparsity,
keep_constrained_dofs);
}
else
{
- // Refinement edges are taken care of
- // by coarser cells
+ // Refinement edges are taken care of by coarser
+ // cells
- // TODO: in the distributed case, we miss out
- // the constraints when the neighbor cell is
- // coarser, but only the current cell is owned
- // locally!
+ // TODO: in the distributed case, we miss out the
+ // constraints when the neighbor cell is coarser, but
+ // only the current cell is owned locally!
if (cell->neighbor_is_coarser(face))
continue;
(dofs_on_this_cell, dofs_on_other_cell,
sparsity, keep_constrained_dofs);
- // only need to add these in case the neighbor
- // cell is not locally owned - otherwise, we
- // touch each face twice and hence put the
- // indices the other way around
+ // only need to add these in case the neighbor cell
+ // is not locally owned - otherwise, we touch each
+ // face twice and hence put the indices the other way
+ // around
if (cell->neighbor(face)->subdomain_id() !=
cell->subdomain_id())
constraints.add_entries_local_to_global
}
+
template <class DH, class SparsityPattern>
void
- make_flux_sparsity_pattern (
- const DH &dof,
- SparsityPattern &sparsity)
+ make_flux_sparsity_pattern (const DH &dof,
+ SparsityPattern &sparsity)
{
ConstraintMatrix constraints;
make_flux_sparsity_pattern (dof, sparsity, constraints);
template <int dim, int spacedim>
Table<2,Coupling>
- dof_couplings_from_component_couplings
- (const FiniteElement<dim,spacedim> &fe,
- const Table<2,Coupling> &component_couplings)
+ dof_couplings_from_component_couplings (const FiniteElement<dim,spacedim> &fe,
+ const Table<2,Coupling> &component_couplings)
{
Assert(component_couplings.n_rows() == fe.n_components(),
ExcDimensionMismatch(component_couplings.n_rows(),
}
+
template <int dim, int spacedim>
std::vector<Table<2,Coupling> >
dof_couplings_from_component_couplings
-
namespace internal
{
namespace
{
- // implementation of the same function in
- // namespace DoFTools for non-hp
- // DoFHandlers
+ // implementation of the same function in namespace DoFTools for
+ // non-hp DoFHandlers
template <class DH, class SparsityPattern>
void
make_flux_sparsity_pattern (const DH &dof,
{
typename DH::cell_iterator
neighbor = cell->neighbor(face);
- // Refinement edges are taken care of
- // by coarser cells
+ // Refinement edges are taken care of by coarser
+ // cells
if (cell->neighbor_is_coarser(face))
continue;
}
- // implementation of the same function in
- // namespace DoFTools for non-hp
- // DoFHandlers
+ // implementation of the same function in namespace DoFTools for
+ // non-hp DoFHandlers
template <int dim, int spacedim, class SparsityPattern>
void
make_flux_sparsity_pattern (const dealii::hp::DoFHandler<dim,spacedim> &dof,
const Table<2,Coupling> &int_mask,
const Table<2,Coupling> &flux_mask)
{
- // while the implementation above is
- // quite optimized and caches a lot of
- // data (see e.g. the int/flux_dof_mask
- // tables), this is no longer practical
- // for the hp version since we would
- // have to have it for all combinations
- // of elements in the
- // hp::FECollection. consequently, the
- // implementation here is simpler and
- // probably less efficient but at least
- // readable...
+ // while the implementation above is quite optimized and caches a
+ // lot of data (see e.g. the int/flux_dof_mask tables), this is no
+ // longer practical for the hp version since we would have to have
+ // it for all combinations of elements in the hp::FECollection.
+ // consequently, the implementation here is simpler and probably
+ // less efficient but at least readable...
const dealii::hp::FECollection<dim,spacedim> &fe = dof.get_fe();
{
dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell);
cell->get_dof_indices (dofs_on_this_cell);
+
// make sparsity pattern for this cell
for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
for (unsigned int j=0; j<cell->get_fe().dofs_per_cell; ++j)
{
typename dealii::hp::DoFHandler<dim,spacedim>::cell_iterator
neighbor = cell->neighbor(face);
- // Refinement edges are taken care of
- // by coarser cells
+
+ // Refinement edges are taken care of by coarser cells
if (cell->neighbor_is_coarser(face))
continue;
const Table<2,Coupling> &int_mask,
const Table<2,Coupling> &flux_mask)
{
- // do the error checking and frame code
- // here, and then pass on to more
- // specialized functions in the internal
- // namespace
+ // do the error checking and frame code here, and then pass on to more
+ // specialized functions in the internal namespace
const unsigned int n_dofs = dof.n_dofs();
const unsigned int n_comp = dof.get_fe().n_components();
Assert (flux_mask.n_cols() == n_comp,
ExcDimensionMismatch (flux_mask.n_cols(), n_comp));
-
- // Clear user flags because we will
- // need them. But first we save
- // them and make sure that we
- // restore them later such that at
- // the end of this function the
- // Triangulation will be in the
- // same state as it was at the
- // beginning of this function.
+ // Clear user flags because we will need them. But first we save them
+ // and make sure that we restore them later such that at the end of
+ // this function the Triangulation will be in the same state as it was
+ // at the beginning of this function.
std::vector<bool> user_flags;
dof.get_tria().save_user_flags(user_flags);
const_cast<Triangulation<DH::dimension,DH::space_dimension> &>(dof.get_tria()).clear_user_flags ();
-
-
-
namespace internal
{
namespace
{
- inline
- bool
+ inline bool
check_master_dof_list (const FullMatrix<double> &face_interpolation_matrix,
const std::vector<unsigned int> &master_dof_list)
{
for (unsigned int j=0; j<N; ++j)
tmp(i,j) = face_interpolation_matrix (master_dof_list[i], j);
- // then use the algorithm
- // from
- // FullMatrix::gauss_jordan
- // on this matrix to find out
- // whether it is
- // singular. the algorithm
- // there does piviting and at
- // the end swaps rows back
- // into their proper order --
- // we omit this step here,
- // since we don't care about
- // the inverse matrix, all we
- // care about is whether the
- // matrix is regular or
- // singular
-
- // first get an estimate of the
- // size of the elements of this
- // matrix, for later checks whether
- // the pivot element is large
- // enough, or whether we have to
- // fear that the matrix is not
+ // then use the algorithm from FullMatrix::gauss_jordan on this
+ // matrix to find out whether it is singular. the algorithm there
+ // does piviting and at the end swaps rows back into their proper
+ // order -- we omit this step here, since we don't care about the
+ // inverse matrix, all we care about is whether the matrix is
+ // regular or singular
+
+ // first get an estimate of the size of the elements of this
+ // matrix, for later checks whether the pivot element is large
+ // enough, or whether we have to fear that the matrix is not
// regular
double diagonal_sum = 0;
for (unsigned int i=0; i<N; ++i)
diagonal_sum += std::fabs(tmp(i,i));
const double typical_diagonal_element = diagonal_sum/N;
- // initialize the array that holds
- // the permutations that we find
+ // initialize the array that holds the permutations that we find
// during pivot search
std::vector<unsigned int> p(N);
for (unsigned int i=0; i<N; ++i)
for (unsigned int j=0; j<N; ++j)
{
- // pivot search: search that
- // part of the line on and
- // right of the diagonal for
- // the largest element
+ // pivot search: search that part of the line on and right of
+ // the diagonal for the largest element
double max = std::fabs(tmp(j,j));
unsigned int r = j;
for (unsigned int i=j+1; i<N; ++i)
r = i;
}
}
- // check whether the
- // pivot is too small. if
- // that is the case, then
- // the matrix is singular
- // and we shouldn't use
- // this set of master
- // dofs
+ // check whether the pivot is too small. if that is the case,
+ // then the matrix is singular and we shouldn't use this set of
+ // master dofs
if (max < 1.e-12*typical_diagonal_element)
return false;
tmp(j,j) = hr;
}
- // everything went fine, so
- // we can accept this set of
- // master dofs (at least as
- // far as they have already
- // been collected)
+ // everything went fine, so we can accept this set of master dofs
+ // (at least as far as they have already been collected)
return true;
}
+
/**
- * When restricting, on a face, the
- * degrees of freedom of fe1 to the
- * space described by fe2 (for example
- * for the complex case described in
- * the @ref hp_paper "hp paper"), we have to select
- * fe2.dofs_per_face out of the
- * fe1.dofs_per_face face DoFs as the
- * master DoFs, and the rest become
- * slave dofs. This function selects
- * which ones will be masters, and
- * which ones will be slaves.
+ * When restricting, on a face, the degrees of freedom of fe1 to the
+ * space described by fe2 (for example for the complex case described
+ * in the @ref hp_paper "hp paper"), we have to select
+ * fe2.dofs_per_face out of the fe1.dofs_per_face face DoFs as the
+ * master DoFs, and the rest become slave dofs. This function selects
+ * which ones will be masters, and which ones will be slaves.
*
- * The function assumes that
- * master_dofs already has size
- * fe1.dofs_per_face. After the
- * function, exactly fe2.dofs_per_face
+ * The function assumes that master_dofs already has size
+ * fe1.dofs_per_face. After the function, exactly fe2.dofs_per_face
* entries will be true.
*
- * The function is a bit
- * complicated since it has to
- * figure out a set a DoFs so
- * that the corresponding rows
- * in the face interpolation
- * matrix are all linearly
- * independent. we have a good
- * heuristic (see the function
- * body) for selecting these
- * rows, but there are cases
- * where this fails and we have
- * to pick them
- * differently. what we do is
- * to run the heuristic and
- * then go back to determine
- * whether we have a set of
- * rows with full row rank. if
- * this isn't the case, go back
- * and select dofs differently
+ * The function is a bit complicated since it has to figure out a set
+ * a DoFs so that the corresponding rows in the face interpolation
+ * matrix are all linearly independent. we have a good heuristic (see
+ * the function body) for selecting these rows, but there are cases
+ * where this fails and we have to pick them differently. what we do
+ * is to run the heuristic and then go back to determine whether we
+ * have a set of rows with full row rank. if this isn't the case, go
+ * back and select dofs differently
*/
template <int dim, int spacedim>
void
(fe2.dofs_per_quad <= fe1.dofs_per_quad),
ExcInternalError());
- // the idea here is to designate as
- // many DoFs in fe1 per object
- // (vertex, line, quad) as master as
- // there are such dofs in fe2
- // (indices are int, because we want
- // to avoid the 'unsigned int < 0 is
- // always false warning for the cases
- // at the bottom in 1d and 2d)
+ // the idea here is to designate as many DoFs in fe1 per object
+ // (vertex, line, quad) as master as there are such dofs in fe2
+ // (indices are int, because we want to avoid the 'unsigned int < 0
+ // is always false warning for the cases at the bottom in 1d and
+ // 2d)
//
- // as mentioned in the paper, it is
- // not always easy to find a set of
- // master dofs that produces an
- // invertible matrix. to this end, we
- // check in each step whether the
- // matrix is still invertible and
- // simply discard this dof if the
- // matrix is not invertible anymore.
+ // as mentioned in the paper, it is not always easy to find a set
+ // of master dofs that produces an invertible matrix. to this end,
+ // we check in each step whether the matrix is still invertible and
+ // simply discard this dof if the matrix is not invertible anymore.
//
- // the cases where we did have
- // trouble in the past were with
- // adding more quad dofs when Q3 and
- // Q4 elements meet at a refined face
- // in 3d (see the hp/crash_12 test
- // that tests that we can do exactly
- // this, and failed before we had
- // code to compensate for this
- // case). the other case are system
- // elements: if we have say a Q1Q2 vs
- // a Q2Q3 element, then we can't just
- // take all master dofs on a line
- // from a single base element, since
- // the shape functions of that base
- // element are independent of that of
- // the other one. this latter case
- // shows up when running
- // hp/hp_constraints_q_system_06
+ // the cases where we did have trouble in the past were with adding
+ // more quad dofs when Q3 and Q4 elements meet at a refined face in
+ // 3d (see the hp/crash_12 test that tests that we can do exactly
+ // this, and failed before we had code to compensate for this
+ // case). the other case are system elements: if we have say a Q1Q2
+ // vs a Q2Q3 element, then we can't just take all master dofs on a
+ // line from a single base element, since the shape functions of
+ // that base element are independent of that of the other one. this
+ // latter case shows up when running hp/hp_constraints_q_system_06
+
std::vector<unsigned int> master_dof_list;
unsigned int index = 0;
for (int v=0;
unsigned int i = 0;
while (dofs_added < fe2.dofs_per_vertex)
{
- // make sure that we
- // were able to find
- // a set of master
- // dofs and that the
- // code down below
- // didn't just reject
- // all our efforts
+ // make sure that we were able to find a set of master dofs
+ // and that the code down below didn't just reject all our
+ // efforts
Assert (i < fe1.dofs_per_vertex,
ExcInternalError());
- // tentatively push
- // this vertex dof
+
+ // tentatively push this vertex dof
master_dof_list.push_back (index+i);
- // then see what
- // happens. if it
- // succeeds, fine
+ // then see what happens. if it succeeds, fine
if (check_master_dof_list (face_interpolation_matrix,
master_dof_list)
== true)
++dofs_added;
else
- // well, it
- // didn't. simply
- // pop that dof
- // from the list
- // again and try
- // with the next
- // dof
+ // well, it didn't. simply pop that dof from the list
+ // again and try with the next dof
master_dof_list.pop_back ();
- // forward counter by
- // one
+ // forward counter by one
++i;
}
index += fe1.dofs_per_vertex;
AssertDimension (index, fe1.dofs_per_face);
AssertDimension (master_dof_list.size(), fe2.dofs_per_face);
- // finally copy the list into the
- // mask
+ // finally copy the list into the mask
std::fill (master_dof_mask.begin(), master_dof_mask.end(), false);
for (std::vector<unsigned int>::const_iterator i=master_dof_list.begin();
i!=master_dof_list.end(); ++i)
/**
- * Make sure that the mask exists that
- * determines which dofs will be the
- * masters on refined faces where an
- * fe1 and a fe2 meet.
+ * Make sure that the mask exists that determines which dofs will be
+ * the masters on refined faces where an fe1 and a fe2 meet.
*/
template <int dim, int spacedim>
void
}
+
/**
- * Make sure that the given @p
- * face_interpolation_matrix
- * pointer points to a valid
- * matrix. If the pointer is zero
- * beforehand, create an entry
- * with the correct data. If it
- * is nonzero, don't touch it.
+ * Make sure that the given @p face_interpolation_matrix pointer
+ * points to a valid matrix. If the pointer is zero beforehand,
+ * create an entry with the correct data. If it is nonzero, don't
+ * touch it.
*/
template <int dim, int spacedim>
void
/**
- * Same, but for subface
- * interpolation matrices.
+ * Same, but for subface interpolation matrices.
*/
template <int dim, int spacedim>
void
}
+
/**
- * Given the face interpolation
- * matrix between two elements,
- * split it into its master and
- * slave parts and invert the
- * master part as explained in
- * the @ref hp_paper "hp paper".
+ * Given the face interpolation matrix between two elements, split it
+ * into its master and slave parts and invert the master part as
+ * explained in the @ref hp_paper "hp paper".
*/
void
ensure_existence_of_split_face_matrix (const FullMatrix<double> &face_interpolation_matrix,
AssertDimension (nth_master_dof, n_master_dofs);
AssertDimension (nth_slave_dof, n_dofs-n_master_dofs);
-//TODO[WB]: We should make sure very small entries are removed after inversion
+ //TODO[WB]: We should make sure very small entries are removed after inversion
split_matrix->first.gauss_jordan ();
}
}
- // a template that can
- // determine statically whether
- // a given DoFHandler class
- // supports different finite
- // element elements
+ // a template that can determine statically whether a given
+ // DoFHandler class supports different finite element elements
template <typename>
struct DoFHandlerSupportsDifferentFEs
{
/**
- * A function that returns how
- * many different finite
- * elements a dof handler
- * uses. This is one for non-hp
- * DoFHandlers and
- * dof_handler.get_fe().size()
- * for the hp-versions.
+ * A function that returns how many different finite elements a dof
+ * handler uses. This is one for non-hp DoFHandlers and
+ * dof_handler.get_fe().size() for the hp-versions.
*/
template <int dim, int spacedim>
unsigned int
/**
- * For a given face belonging
- * to an active cell that
- * borders to a more refined
- * cell, return the fe_index of
- * the most dominating finite
- * element used on any of the
- * face's subfaces.
+ * For a given face belonging to an active cell that borders to a
+ * more refined cell, return the fe_index of the most dominating
+ * finite element used on any of the face's subfaces.
*/
template <typename face_iterator>
unsigned int
for (; dominating_subface_no<face->n_children();
++dominating_subface_no)
{
- // each of the subfaces
- // can have only a single
- // fe_index associated
- // with them, since there
- // is no cell on the
- // other side
+ // each of the subfaces can have only a single fe_index
+ // associated with them, since there is no cell on the other
+ // side
Assert (face->child(dominating_subface_no)
->n_active_fe_indices()
== 1,
this_subface_fe.compare_for_face_domination(that_subface_fe);
}
- // see if the element
- // on this subface is
- // able to dominate
- // the ones on all
- // other subfaces,
- // and if so take it
+ // see if the element on this subface is able to dominate the
+ // ones on all other subfaces, and if so take it
if ((domination == FiniteElementDomination::this_element_dominates)
||
(domination == FiniteElementDomination::either_element_can_dominate))
break;
}
- // check that we have
- // found one such subface
+ // check that we have found one such subface
Assert (dominating_subface_no < face->n_children(),
ExcNotImplemented());
- // return the finite element
- // index used on it. note
- // that only a single fe can
- // be active on such subfaces
+ // return the finite element index used on it. note that only a
+ // single fe can be active on such subfaces
return face->child (dominating_subface_no)->nth_active_fe_index(0);
}
/**
- * Copy constraints into a constraint
- * matrix object.
+ * Copy constraints into a constraint matrix object.
*
- * This function removes zero
- * constraints and those, which
- * constrain a DoF which was
- * already eliminated in one of
- * the previous steps of the hp
- * hanging node procedure.
+ * This function removes zero constraints and those, which constrain
+ * a DoF which was already eliminated in one of the previous steps of
+ * the hp hanging node procedure.
*
- * It also suppresses very small
- * entries in the constraint matrix to
- * avoid making the sparsity pattern
- * fuller than necessary.
+ * It also suppresses very small entries in the constraint matrix to
+ * avoid making the sparsity pattern fuller than necessary.
*/
void
filter_constraints (const std::vector<unsigned int> &master_dofs,
const unsigned int n_master_dofs = master_dofs.size ();
const unsigned int n_slave_dofs = slave_dofs.size ();
- // check for a couple
- // conditions that happened
- // in parallel distributed
- // mode
+ // check for a couple conditions that happened in parallel
+ // distributed mode
for (unsigned int row=0; row!=n_slave_dofs; ++row)
Assert (slave_dofs[row] != numbers::invalid_unsigned_int,
ExcInternalError());
{
bool constraint_already_satisfied = false;
- // Check if we have an identity
- // constraint, which is already
- // satisfied by unification of
- // the corresponding global dof
+ // Check if we have an identity constraint, which is already
+ // satisfied by unification of the corresponding global dof
// indices
for (unsigned int i=0; i<n_master_dofs; ++i)
if (face_constraints (row,i) == 1.0)
if (constraint_already_satisfied == false)
{
- // add up the absolute
- // values of all
- // constraints in this line
- // to get a measure of
- // their absolute size
+ // add up the absolute values of all constraints in this
+ // line to get a measure of their absolute size
double abs_sum = 0;
for (unsigned int i=0; i<n_master_dofs; ++i)
abs_sum += std::abs (face_constraints(row,i));
- // then enter those
- // constraints that are
- // larger than
- // 1e-14*abs_sum. everything
- // else probably originated
- // from inexact inversion
- // of matrices and similar
- // effects. having those
- // constraints in here will
- // only lead to problems
- // because it makes
- // sparsity patterns fuller
- // than necessary without
- // producing any
+ // then enter those constraints that are larger than
+ // 1e-14*abs_sum. everything else probably originated
+ // from inexact inversion of matrices and similar
+ // effects. having those constraints in here will only
+ // lead to problems because it makes sparsity patterns
+ // fuller than necessary without producing any
// significant effect
constraints.add_line (slave_dofs[row]);
for (unsigned int i=0; i<n_master_dofs; ++i)
make_hp_hanging_node_constraints (const dealii::DoFHandler<1> &,
ConstraintMatrix &)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
ConstraintMatrix &,
dealii::internal::int2type<1>)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &,
ConstraintMatrix &)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
ConstraintMatrix &,
dealii::internal::int2type<1>)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/,
ConstraintMatrix & /*constraints*/)
{
- // we may have to compute
- // constraints for
- // vertices. gotta think about
+ // we may have to compute constraints for vertices. gotta think about
// that a bit more
-//TODO[WB]: think about what to do here...
+
+ //TODO[WB]: think about what to do here...
}
ConstraintMatrix & /*constraints*/,
dealii::internal::int2type<1>)
{
- // we may have to compute
- // constraints for
- // vertices. gotta think about
+ // we may have to compute constraints for vertices. gotta think about
// that a bit more
-//TODO[WB]: think about what to do here...
+
+ //TODO[WB]: think about what to do here...
}
make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &,
ConstraintMatrix &)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
ConstraintMatrix &,
dealii::internal::int2type<1>)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
make_hp_hanging_node_constraints (const dealii::DoFHandler<1,3> &,
ConstraintMatrix &)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
void
ConstraintMatrix &,
dealii::internal::int2type<1>)
{
- // nothing to do for regular
- // dof handlers in 1d
+ // nothing to do for regular dof handlers in 1d
}
// // constraints for
// // vertices. gotta think about
// // that a bit more
+//
// //TODO[WB]: think about what to do here...
// }
//#endif
std::vector<unsigned int> dofs_on_mother;
std::vector<unsigned int> dofs_on_children;
- // loop over all lines; only on
- // lines there can be constraints.
- // We do so by looping over all
- // active cells and checking
- // whether any of the faces are
- // refined which can only be from
- // the neighboring cell because
- // this one is active. In that
- // case, the face is subject to
+ // loop over all lines; only on lines there can be constraints. We do
+ // so by looping over all active cells and checking whether any of
+ // the faces are refined which can only be from the neighboring cell
+ // because this one is active. In that case, the face is subject to
// constraints
//
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
+ // note that even though we may visit a face twice if the neighboring
+ // cells are equally refined, we can only visit each face with
+ // hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
- // artificial cells can at
- // best neighbor ghost cells,
- // but we're not interested
- // in these interfaces
+ // artificial cells can at best neighbor ghost cells, but we're not
+ // interested in these interfaces
if (!cell->is_artificial ())
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->has_children())
{
- // in any case, faces
- // can have at most two
- // active fe indices,
- // but here the face
- // can have only one
- // (namely the same as
- // that from the cell
- // we're sitting on),
- // and each of the
- // children can have
- // only one as
- // well. check this
+ // in any case, faces can have at most two active fe
+ // indices, but here the face can have only one (namely the
+ // same as that from the cell we're sitting on), and each
+ // of the children can have only one as well. check this
Assert (cell->face(face)->n_active_fe_indices() == 1,
ExcInternalError());
Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
ExcInternalError());
- // right now, all that
- // is implemented is
- // the case that both
- // sides use the same
- // fe
+ // right now, all that is implemented is the case that both
+ // sides use the same fe
for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
Assert (cell->face(face)->child(c)
const typename DH::line_iterator this_face = cell->face(face);
- // fill the dofs indices. Use same
- // enumeration scheme as in
+ // fill the dofs indices. Use same enumeration scheme as in
// @p{FiniteElement::constraints()}
unsigned int next_index = 0;
for (unsigned int vertex=0; vertex<2; ++vertex)
= this_face->child(child)->dof_index(dof, fe_index);
AssertDimension (next_index, dofs_on_children.size());
- // for each row in the constraint
- // matrix for this line:
+ // for each row in the constraint matrix for this line:
for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
{
constraints.add_line (dofs_on_children[row]);
}
else
{
- // this face has no
- // children, but it
- // could still be
- // that it is shared
- // by two cells that
- // use a different fe
- // index. check a
- // couple of things,
- // but ignore the
- // case that the
- // neighbor is an
- // artificial cell
+ // this face has no children, but it could still be that it
+ // is shared by two cells that use a different fe index.
+ // check a couple of things, but ignore the case that the
+ // neighbor is an artificial cell
if (!cell->at_boundary(face) &&
!cell->neighbor(face)->is_artificial())
{
std::vector<unsigned int> dofs_on_mother;
std::vector<unsigned int> dofs_on_children;
- // loop over all quads; only on
- // quads there can be constraints.
- // We do so by looping over all
- // active cells and checking
- // whether any of the faces are
- // refined which can only be from
- // the neighboring cell because
- // this one is active. In that
- // case, the face is subject to
+ // loop over all quads; only on quads there can be constraints. We do
+ // so by looping over all active cells and checking whether any of
+ // the faces are refined which can only be from the neighboring cell
+ // because this one is active. In that case, the face is subject to
// constraints
//
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
+ // note that even though we may visit a face twice if the neighboring
+ // cells are equally refined, we can only visit each face with
+ // hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
- // artificial cells can at
- // best neighbor ghost cells,
- // but we're not interested
- // in these interfaces
+ // artificial cells can at best neighbor ghost cells, but we're not
+ // interested in these interfaces
if (!cell->is_artificial ())
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->has_children())
{
- // first of all, make sure that
- // we treat a case which is
- // possible, i.e. either no dofs
- // on the face at all or no
+ // first of all, make sure that we treat a case which is
+ // possible, i.e. either no dofs on the face at all or no
// anisotropic refinement
if (cell->get_fe().dofs_per_face == 0)
continue;
Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
ExcNotImplemented());
- // in any case, faces
- // can have at most two
- // active fe indices,
- // but here the face
- // can have only one
- // (namely the same as
- // that from the cell
- // we're sitting on),
- // and each of the
- // children can have
- // only one as
- // well. check this
+ // in any case, faces can have at most two active fe
+ // indices, but here the face can have only one (namely the
+ // same as that from the cell we're sitting on), and each
+ // of the children can have only one as well. check this
AssertDimension (cell->face(face)->n_active_fe_indices(), 1);
Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
== true,
for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
AssertDimension (cell->face(face)->child(c)->n_active_fe_indices(), 1);
- // right now, all that
- // is implemented is
- // the case that both
- // sides use the same
- // fe, and not only
- // that but also that
- // all lines bounding
- // this face and the
- // children have the
+ // right now, all that is implemented is the case that both
+ // sides use the same fe, and not only that but also that
+ // all lines bounding this face and the children have the
// same fe
for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+
12*fe.dofs_per_line+
4*fe.dofs_per_quad);
-//TODO[TL]: think about this and the following in case of anisotropic refinement
+
+ //TODO[TL]: think about this and the following in case of anisotropic refinement
dofs_on_mother.resize (n_dofs_on_mother);
dofs_on_children.resize (n_dofs_on_children);
const typename DH::face_iterator this_face = cell->face(face);
- // fill the dofs indices. Use same
- // enumeration scheme as in
+ // fill the dofs indices. Use same enumeration scheme as in
// @p{FiniteElement::constraints()}
unsigned int next_index = 0;
for (unsigned int vertex=0; vertex<4; ++vertex)
next_index = 0;
- // assert some consistency
- // assumptions
- //TODO[TL]: think about this in case of anisotropic refinement
+ // assert some consistency assumptions
+
+ //TODO[TL]: think about this in case of anisotropic
+ //refinement
+
Assert (dof_handler.get_tria().get_anisotropic_refinement_flag() ||
((this_face->child(0)->vertex_index(3) ==
this_face->child(1)->vertex_index(2)) &&
dofs_on_children[next_index++]
= this_face->child(0)->vertex_dof_index(3,dof);
- // dof numbers on the centers of
- // the lines bounding this face
+ // dof numbers on the centers of the lines bounding this
+ // face
for (unsigned int line=0; line<4; ++line)
for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
dofs_on_children[next_index++]
= this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index);
- // next the dofs on the lines interior
- // to the face; the order of these
- // lines is laid down in the
- // FiniteElement class documentation
+ // next the dofs on the lines interior to the face; the
+ // order of these lines is laid down in the FiniteElement
+ // class documentation
for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
dofs_on_children[next_index++]
= this_face->child(0)->line(1)->dof_index(dof, fe_index);
dofs_on_children[next_index++]
= this_face->line(line)->child(child)->dof_index(dof, fe_index);
- // finally, for the dofs interior
- // to the four child faces
+ // finally, for the dofs interior to the four child faces
for (unsigned int child=0; child<4; ++child)
for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
dofs_on_children[next_index++]
= this_face->child(child)->dof_index(dof, fe_index);
AssertDimension (next_index, dofs_on_children.size());
- // for each row in the constraint
- // matrix for this line:
+ // for each row in the constraint matrix for this line:
for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
{
constraints.add_line (dofs_on_children[row]);
}
else
{
- // this face has no
- // children, but it
- // could still be
- // that it is shared
- // by two cells that
- // use a different fe
- // index. check a
- // couple of things,
- // but ignore the
- // case that the
- // neighbor is an
- // artificial cell
+ // this face has no children, but it could still be that it
+ // is shared by two cells that use a different fe index.
+ // check a couple of things, but ignore the case that the
+ // neighbor is an artificial cell
if (!cell->at_boundary(face) &&
!cell->neighbor(face)->is_artificial())
{
make_hp_hanging_node_constraints (const DH &dof_handler,
ConstraintMatrix &constraints)
{
- // note: this function is going
- // to be hard to understand if
- // you haven't read the hp
- // paper. however, we try to
- // follow the notation laid out
- // there, so go read the paper
- // before you try to understand
+ // note: this function is going to be hard to understand if you
+ // haven't read the hp paper. however, we try to follow the notation
+ // laid out there, so go read the paper before you try to understand
// what is going on here
+
const unsigned int dim = DH::dimension;
const unsigned int spacedim = DH::space_dimension;
- // a matrix to be used for
- // constraints below. declared
- // here and simply resized down
- // below to avoid permanent
- // re-allocation of memory
+ // a matrix to be used for constraints below. declared here and
+ // simply resized down below to avoid permanent re-allocation of
+ // memory
FullMatrix<double> constraint_matrix;
- // similarly have arrays that
- // will hold master and slave
- // dof numbers, as well as a
- // scratch array needed for the
- // complicated case below
+ // similarly have arrays that will hold master and slave dof numbers,
+ // as well as a scratch array needed for the complicated case below
std::vector<unsigned int> master_dofs;
std::vector<unsigned int> slave_dofs;
std::vector<unsigned int> scratch_dofs;
- // caches for the face and
- // subface interpolation
- // matrices between different
- // (or the same) finite
- // elements. we compute them
- // only once, namely the first
- // time they are needed, and
- // then just reuse them
+ // caches for the face and subface interpolation matrices between
+ // different (or the same) finite elements. we compute them only
+ // once, namely the first time they are needed, and then just reuse
+ // them
Table<2,std_cxx1x::shared_ptr<FullMatrix<double> > >
face_interpolation_matrices (n_finite_elements (dof_handler),
n_finite_elements (dof_handler));
n_finite_elements (dof_handler),
GeometryInfo<dim>::max_children_per_face);
- // similarly have a cache for
- // the matrices that are split
- // into their master and slave
- // parts, and for which the
- // master part is
- // inverted. these two matrices
- // are derived from the face
- // interpolation matrix as
- // described in the @ref hp_paper "hp paper"
+ // similarly have a cache for the matrices that are split into their
+ // master and slave parts, and for which the master part is inverted.
+ // these two matrices are derived from the face interpolation matrix
+ // as described in the @ref hp_paper "hp paper"
Table<2,std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > > >
split_face_interpolation_matrices (n_finite_elements (dof_handler),
n_finite_elements (dof_handler));
- // finally, for each pair of finite
- // elements, have a mask that states
- // which of the degrees of freedom on
- // the coarse side of a refined face
- // will act as master dofs.
+ // finally, for each pair of finite elements, have a mask that states
+ // which of the degrees of freedom on the coarse side of a refined
+ // face will act as master dofs.
Table<2,std_cxx1x::shared_ptr<std::vector<bool> > >
master_dof_masks (n_finite_elements (dof_handler),
n_finite_elements (dof_handler));
// loop over all faces
//
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
+ // note that even though we may visit a face twice if the neighboring
+ // cells are equally refined, we can only visit each face with
+ // hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
- // artificial cells can at
- // best neighbor ghost cells,
- // but we're not interested
- // in these interfaces
+ // artificial cells can at best neighbor ghost cells, but we're not
+ // interested in these interfaces
if (!cell->is_artificial ())
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->has_children())
{
- // first of all, make sure that
- // we treat a case which is
- // possible, i.e. either no dofs
- // on the face at all or no
+ // first of all, make sure that we treat a case which is
+ // possible, i.e. either no dofs on the face at all or no
// anisotropic refinement
if (cell->get_fe().dofs_per_face == 0)
continue;
Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
ExcNotImplemented());
- // so now we've found a
- // face of an active
- // cell that has
- // children. that means
- // that there are
- // hanging nodes here.
-
- // in any case, faces
- // can have at most two
- // sets of active fe
- // indices, but here
- // the face can have
- // only one (namely the
- // same as that from
- // the cell we're
- // sitting on), and
- // each of the children
- // can have only one as
- // well. check this
+ // so now we've found a face of an active cell that has
+ // children. that means that there are hanging nodes here.
+
+ // in any case, faces can have at most two sets of active
+ // fe indices, but here the face can have only one (namely
+ // the same as that from the cell we're sitting on), and
+ // each of the children can have only one as well. check
+ // this
Assert (cell->face(face)->n_active_fe_indices() == 1,
ExcInternalError());
Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
ExcInternalError());
- // first find out
- // whether we can
- // constrain each of
- // the subfaces to the
- // mother face. in the
- // lingo of the hp
- // paper, this would be
- // the simple
- // case. note that we
- // can short-circuit
- // this decision if the
- // dof_handler doesn't
+ // first find out whether we can constrain each of the
+ // subfaces to the mother face. in the lingo of the hp
+ // paper, this would be the simple case. note that we can
+ // short-circuit this decision if the dof_handler doesn't
// support hp at all
//
- // ignore all
- // interfaces with
- // artificial cells
+ // ignore all interfaces with artificial cells
FiniteElementDomination::Domination
mother_face_dominates = FiniteElementDomination::either_element_can_dominate;
case FiniteElementDomination::this_element_dominates:
case FiniteElementDomination::either_element_can_dominate:
{
- // Case 1 (the
- // simple case
- // and the only
- // case that can
- // happen for
- // non-hp
- // DoFHandlers):
- // The coarse
- // element
- // dominates the
- // elements on
- // the subfaces
- // (or they are
+ // Case 1 (the simple case and the only case that can
+ // happen for non-hp DoFHandlers): The coarse element
+ // dominates the elements on the subfaces (or they are
// all the same)
//
- // so we are
- // going to
- // constrain
- // the DoFs on
- // the face
- // children
- // against the
- // DoFs on the
- // face itself
+ // so we are going to constrain the DoFs on the face
+ // children against the DoFs on the face itself
master_dofs.resize (cell->get_fe().dofs_per_face);
cell->face(face)->get_dof_indices (master_dofs,
cell->active_fe_index ());
- // Now create
- // constraint matrix
- // for the subfaces and
- // assemble it. ignore
- // all interfaces with
- // artificial cells
- // because we can only
- // get to such
- // interfaces if the
- // current cell is a
- // ghost cell. also
- // ignore the interface
- // if the neighboring
- // cell is a ghost cell
- // in 2d: what we would
- // compute here are the
- // constraints on the
- // ghost cell's DoFs,
- // but we are not
- // interested in those:
- // we only want
- // constraints on
- // *locally active*
- // DoFs, not on
- // *locally relevant*
- // DoFs. However, in 3d
- // we must still
- // compute those
- // constraints because
- // it might happen that
- // a constraint is
- // related to an edge
- // where the hanging
- // node is only
- // detected if we also
- // look between ghosts
+ // Now create constraint matrix for the subfaces and
+ // assemble it. ignore all interfaces with artificial
+ // cells because we can only get to such interfaces if
+ // the current cell is a ghost cell. also ignore the
+ // interface if the neighboring cell is a ghost cell in
+ // 2d: what we would compute here are the constraints
+ // on the ghost cell's DoFs, but we are not interested
+ // in those: we only want constraints on *locally
+ // active* DoFs, not on *locally relevant* DoFs.
+ // However, in 3d we must still compute those
+ // constraints because it might happen that a
+ // constraint is related to an edge where the hanging
+ // node is only detected if we also look between ghosts
for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
{
if (cell->neighbor_child_on_subface (face, c)->is_artificial()
const unsigned int
subface_fe_index = subface->nth_active_fe_index(0);
- // we sometime run
- // into the
- // situation where
- // for example on
- // one big cell we
- // have a FE_Q(1)
- // and on the
- // subfaces we have
- // a mixture of
- // FE_Q(1) and
- // FE_Nothing. In
- // that case, the
- // face domination
- // is
- // either_element_can_dominate
- // for the whole
- // collection of
- // subfaces, but on
- // the particular
- // subface between
- // FE_Q(1) and
- // FE_Nothing,
- // there are no
- // constraints that
- // we need to take
- // care of. in that
- // case, just
- // continue
+ // we sometime run into the situation where for
+ // example on one big cell we have a FE_Q(1) and on
+ // the subfaces we have a mixture of FE_Q(1) and
+ // FE_Nothing. In that case, the face domination is
+ // either_element_can_dominate for the whole
+ // collection of subfaces, but on the particular
+ // subface between FE_Q(1) and FE_Nothing, there
+ // are no constraints that we need to take care of.
+ // in that case, just continue
if (cell->get_fe().compare_for_face_domination
(subface->get_fe(subface_fe_index))
==
FiniteElementDomination::no_requirements)
continue;
- // Same procedure as for the
- // mother cell. Extract the face
- // DoFs from the cell DoFs.
+ // Same procedure as for the mother cell. Extract
+ // the face DoFs from the cell DoFs.
slave_dofs.resize (subface->get_fe(subface_fe_index)
.dofs_per_face);
subface->get_dof_indices (slave_dofs, subface_fe_index);
Assert (slave_dofs[i] != numbers::invalid_unsigned_int,
ExcInternalError());
- // Now create the
- // element constraint
- // for this subface.
+ // Now create the element constraint for this
+ // subface.
//
- // As a side remark,
- // one may wonder the
- // following:
- // neighbor_child is
- // clearly computed
- // correctly,
- // i.e. taking into
- // account
- // face_orientation
- // (just look at the
- // implementation of
- // that
- // function). however,
- // we don't care about
- // this here, when we
- // ask for
- // subface_interpolation
- // on subface c. the
- // question rather is:
- // do we have to
- // translate 'c' here
- // as well?
+ // As a side remark, one may wonder the following:
+ // neighbor_child is clearly computed correctly,
+ // i.e. taking into account face_orientation (just
+ // look at the implementation of that function).
+ // however, we don't care about this here, when we
+ // ask for subface_interpolation on subface c. the
+ // question rather is: do we have to translate 'c'
+ // here as well?
//
- // the answer is in
- // fact 'no'. if one
- // does that, results
- // are wrong:
- // constraints are
- // added twice for the
- // same pair of nodes
- // but with differing
- // weights. in
- // addition, one can
- // look at the
- // deal.II/project_*_03
- // tests that look at
- // exactly this case:
- // there, we have a
- // mesh with at least
- // one
- // face_orientation==false
- // and hanging nodes,
- // and the results of
- // those tests show
- // that the result of
- // projection verifies
- // the approximation
- // properties of a
- // finite element onto
- // that mesh
+ // the answer is in fact 'no'. if one does that,
+ // results are wrong: constraints are added twice
+ // for the same pair of nodes but with differing
+ // weights. in addition, one can look at the
+ // deal.II/project_*_03 tests that look at exactly
+ // this case: there, we have a mesh with at least
+ // one face_orientation==false and hanging nodes,
+ // and the results of those tests show that the
+ // result of projection verifies the approximation
+ // properties of a finite element onto that mesh
ensure_existence_of_subface_matrix
(cell->get_fe(),
subface->get_fe(subface_fe_index),
subface_interpolation_matrices
[cell->active_fe_index()][subface_fe_index][c]);
- // Add constraints to global constraint
- // matrix.
+ // Add constraints to global constraint matrix.
filter_constraints (master_dofs,
slave_dofs,
*(subface_interpolation_matrices
case FiniteElementDomination::other_element_dominates:
case FiniteElementDomination::neither_element_dominates:
{
- // Case 2 (the "complex"
- // case): at least one
- // (the neither_... case)
- // of the finer elements
- // or all of them (the
- // other_... case) is
- // dominating. See the hp
- // paper for a way how to
- // deal with this
- // situation
+ // Case 2 (the "complex" case): at least one (the
+ // neither_... case) of the finer elements or all of
+ // them (the other_... case) is dominating. See the hp
+ // paper for a way how to deal with this situation
//
- // since this is
- // something that
- // can only
- // happen for hp
- // dof handlers,
- // add a check
- // here...
+ // since this is something that can only happen for hp
+ // dof handlers, add a check here...
Assert (DoFHandlerSupportsDifferentFEs<DH>::value == true,
ExcInternalError());
- // we first have
- // to find the
- // finite element
- // that is able
- // to generate a
- // space that all
- // the other ones
- // can be
- // constrained to
+ // we first have to find the finite element that is
+ // able to generate a space that all the other ones can
+ // be constrained to
const unsigned int dominating_fe_index
= get_most_dominating_subface_fe_index (cell->face(face));
const FiniteElement<dim,spacedim> &dominating_fe
= dof_handler.get_fe()[dominating_fe_index];
- // check also
- // that it is
- // able to
- // constrain the
- // mother
- // face. it
- // should be, or
- // we wouldn't
- // have gotten
- // into the
- // branch for the
- // 'complex' case
+ // check also that it is able to constrain the mother
+ // face. it should be, or we wouldn't have gotten into
+ // the branch for the 'complex' case
Assert ((dominating_fe.compare_for_face_domination
(cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0)))
== FiniteElementDomination::this_element_dominates)
ExcInternalError());
- // first get the
- // interpolation matrix
- // from the mother to the
- // virtual dofs
+ // first get the interpolation matrix from the mother
+ // to the virtual dofs
Assert (dominating_fe.dofs_per_face <=
cell->get_fe().dofs_per_face,
ExcInternalError());
face_interpolation_matrices
[dominating_fe_index][cell->active_fe_index()]);
- // split this matrix into
- // master and slave
- // components. invert the
- // master component
+ // split this matrix into master and slave components.
+ // invert the master component
ensure_existence_of_master_dof_mask
(cell->get_fe(),
dominating_fe,
= (split_face_interpolation_matrices
[dominating_fe_index][cell->active_fe_index()]->second);
- // now compute
- // the constraint
- // matrix as the
- // product
- // between the
- // inverse matrix
- // and the slave
- // part
+ // now compute the constraint matrix as the product
+ // between the inverse matrix and the slave part
constraint_matrix.reinit (cell->get_fe().dofs_per_face -
dominating_fe.dofs_per_face,
dominating_fe.dofs_per_face);
.mmult (constraint_matrix,
restrict_mother_to_virtual_master_inv);
- // then figure
- // out the global
- // numbers of
- // master and
- // slave dofs and
- // apply
- // constraints
+ // then figure out the global numbers of master and
+ // slave dofs and apply constraints
scratch_dofs.resize (cell->get_fe().dofs_per_face);
cell->face(face)->get_dof_indices (scratch_dofs,
cell->active_fe_index ());
- // split dofs into master
- // and slave components
+ // split dofs into master and slave components
master_dofs.clear ();
slave_dofs.clear ();
for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
- // next we have
- // to deal with
- // the
- // subfaces. do
- // as discussed
- // in the hp
- // paper
+ // next we have to deal with the subfaces. do as
+ // discussed in the hp paper
for (unsigned int sf=0;
sf<cell->face(face)->n_children(); ++sf)
{
- // ignore
- // interfaces
- // with
- // artificial
- // cells as
- // well as
- // interfaces
- // between
- // ghost
- // cells in 2d
+ // ignore interfaces with artificial cells as well
+ // as interfaces between ghost cells in 2d
if (cell->neighbor_child_on_subface (face, sf)->is_artificial()
||
(dim==2 && cell->is_ghost()
const FiniteElement<dim,spacedim> &subface_fe
= dof_handler.get_fe()[subface_fe_index];
- // first get the
- // interpolation
- // matrix from the
- // subface to the
- // virtual dofs
+ // first get the interpolation matrix from the
+ // subface to the virtual dofs
Assert (dominating_fe.dofs_per_face <=
subface_fe.dofs_per_face,
ExcInternalError());
}
case FiniteElementDomination::no_requirements:
- // there
- // are no
- // continuity
- // requirements
- // between
- // the two
- // elements. record
- // no
- // constraints
+ // there are no continuity requirements between the two
+ // elements. record no constraints
break;
default:
}
else
{
- // this face has no
- // children, but it
- // could still be that
- // it is shared by two
- // cells that use a
- // different fe index
+ // this face has no children, but it could still be that it
+ // is shared by two cells that use a different fe index
Assert (cell->face(face)
->fe_index_is_active(cell->active_fe_index()) == true,
ExcInternalError());
- // see if there is a
- // neighbor that is
- // an artificial
- // cell. in that
- // case, we're not
- // interested in this
- // interface. we test
- // this case first
- // since artificial
- // cells may not have
- // an active_fe_index
- // set, etc
+ // see if there is a neighbor that is an artificial cell.
+ // in that case, we're not interested in this interface. we
+ // test this case first since artificial cells may not have
+ // an active_fe_index set, etc
if (!cell->at_boundary(face)
&&
cell->neighbor(face)->is_artificial())
continue;
- // Only if there is
- // a neighbor with
- // a different
- // active_fe_index
- // and the same h-level,
- // some action has
- // to be taken.
+ // Only if there is a neighbor with a different
+ // active_fe_index and the same h-level, some action has to
+ // be taken.
if ((DoFHandlerSupportsDifferentFEs<DH>::value == true)
&&
!cell->face(face)->at_boundary ()
{
const typename DH::cell_iterator neighbor = cell->neighbor (face);
- // see which side of the
- // face we have to
- // constrain
+ // see which side of the face we have to constrain
switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ()))
{
case FiniteElementDomination::this_element_dominates:
{
- // Get DoFs on
- // dominating and
- // dominated side of
- // the face
+ // Get DoFs on dominating and dominated side of the
+ // face
master_dofs.resize (cell->get_fe().dofs_per_face);
cell->face(face)->get_dof_indices (master_dofs,
cell->active_fe_index ());
cell->face(face)->get_dof_indices (slave_dofs,
neighbor->active_fe_index ());
- // break if the n_master_dofs == 0,
- // because we are attempting to
- // constrain to an element that has
+ // break if the n_master_dofs == 0, because we are
+ // attempting to constrain to an element that has
// has no face dofs
if (master_dofs.size() == 0) break;
- // make sure
- // the element
- // constraints
- // for this
- // face are
- // available
+ // make sure the element constraints for this face
+ // are available
ensure_existence_of_face_matrix
(cell->get_fe(),
neighbor->get_fe(),
face_interpolation_matrices
[cell->active_fe_index()][neighbor->active_fe_index()]);
- // Add constraints to global constraint
- // matrix.
+ // Add constraints to global constraint matrix.
filter_constraints (master_dofs,
slave_dofs,
*(face_interpolation_matrices
case FiniteElementDomination::other_element_dominates:
{
- // we don't do anything
- // here since we will
- // come back to this
- // face from the other
- // cell, at which time
- // we will fall into
- // the first case
- // clause above
+ // we don't do anything here since we will come
+ // back to this face from the other cell, at which
+ // time we will fall into the first case clause
+ // above
break;
}
case FiniteElementDomination::either_element_can_dominate:
{
- // it appears as if
- // neither element has
- // any constraints on
- // its neighbor. this
- // may be because
- // neither element has
- // any DoFs on faces at
- // all. or that the two
- // elements are
- // actually the same,
- // although they happen
- // to run under
- // different fe_indices
- // (this is what
- // happens in
- // hp/hp_hanging_nodes_01
- // for example).
+ // it appears as if neither element has any
+ // constraints on its neighbor. this may be because
+ // neither element has any DoFs on faces at all. or
+ // that the two elements are actually the same,
+ // although they happen to run under different
+ // fe_indices (this is what happens in
+ // hp/hp_hanging_nodes_01 for example).
//
- // another possibility
- // is what happens in
- // crash_13. there, we
- // have
- // FESystem(FE_Q(1),FE_DGQ(0))
- // vs. FESystem(FE_Q(1),FE_DGQ(1)).
- // neither of them
- // dominates the
- // other. the point is
- // that it doesn't
- // matter, since
- // hopefully for this
- // case, both sides'
- // dofs should have
- // been unified.
+ // another possibility is what happens in crash_13.
+ // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs.
+ // FESystem(FE_Q(1),FE_DGQ(1)). neither of them
+ // dominates the other. the point is that it
+ // doesn't matter, since hopefully for this case,
+ // both sides' dofs should have been unified.
//
- // make sure this is
- // actually true. this
- // actually only
- // matters, of course,
- // if either of the two
- // finite elements
- // actually do have
- // dofs on the face
+ // make sure this is actually true. this actually
+ // only matters, of course, if either of the two
+ // finite elements actually do have dofs on the
+ // face
if ((cell->get_fe().dofs_per_face != 0)
||
(cell->neighbor(face)->get_fe().dofs_per_face != 0))
cell->neighbor(face)->get_fe().dofs_per_face,
ExcNotImplemented());
- // (ab)use the master
- // and slave dofs
- // arrays for a
- // moment here
+ // (ab)use the master and slave dofs arrays for
+ // a moment here
master_dofs.resize (cell->get_fe().dofs_per_face);
cell->face(face)
->get_dof_indices (master_dofs,
case FiniteElementDomination::neither_element_dominates:
{
- // we don't presently
- // know what exactly to
- // do here. it isn't quite
- // clear what exactly we
- // would have to do
- // here. sit tight until
- // someone trips over the
- // following statement
- // and see what exactly
- // is going on
+ // we don't presently know what exactly to do here.
+ // it isn't quite clear what exactly we would have
+ // to do here. sit tight until someone trips over
+ // the following statement and see what exactly is
+ // going on
Assert (false, ExcNotImplemented());
break;
}
}
default:
- // we shouldn't get
- // here
+ // we shouldn't get here
Assert (false, ExcInternalError());
}
}
make_hanging_node_constraints (const DH &dof_handler,
ConstraintMatrix &constraints)
{
- // Decide whether to use the
- // new or old make_hanging_node_constraints
- // function. If all the FiniteElement
- // or all elements in a FECollection support
- // the new face constraint matrix, the
- // new code will be used.
- // Otherwise, the old implementation is used
- // for the moment.
+ // Decide whether to use the new or old make_hanging_node_constraints
+ // function. If all the FiniteElement or all elements in a FECollection
+ // support the new face constraint matrix, the new code will be used.
+ // Otherwise, the old implementation is used for the moment.
if (dof_handler.get_fe().hp_constraints_are_implemented ())
internal::
make_hp_hanging_node_constraints (dof_handler,
// A lookup table on how to go through the child faces depending on the
// orientation:
- static const int lookup_table_2d[2][2]
- // flip
- = { {0, 1}, // false
- {1, 0}, // true
+ static const int lookup_table_2d[2][2] =
+ {
+ // flip:
+ {0, 1}, // false
+ {1, 0}, // true
};
- static const int lookup_table_3d[2][2][2][4]
- // orientation flip rotation
- = { { { {0, 2, 1, 3}, // false false false
- {2, 3, 0, 1},
- }, // false false true
- { {3, 1, 2, 0}, // false true false
- {1, 0, 3, 2},
+ static const int lookup_table_3d[2][2][2][4] =
+ {
+ // orientation flip rotation
+ { { {0, 2, 1, 3}, // false false false
+ {2, 3, 0, 1}, // false false true
+ },
+ { {3, 1, 2, 0}, // false true false
+ {1, 0, 3, 2}, // false true true
+ },
+ },
+ { { {0, 1, 2, 3}, // true false false
+ {1, 3, 0, 2}, // true false true
},
- }, // false true true
- { { {0, 1, 2, 3}, // true false false
- {1, 3, 0, 2},
- }, // true false true
- { {3, 2, 1, 0}, // true true false
- {2, 0, 3, 1},
+ { {3, 2, 1, 0}, // true true false
+ {2, 0, 3, 1}, // true true true
},
- }, // true true true
+ },
};
// In the case that both faces have children, we loop over all
// face_flip,
// face_rotation)
// But this will give us an index on a cell - something we cannot work
- // with. But luckily we can match them back :-]
+ // with directly. But luckily we can match them back :-]
std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
// Build up a cell to face index for face_2:
for (unsigned int i = 0; i < dofs_per_face; ++i)
{
- const unsigned int cell_index = fe.
- face_to_cell_index(i,
- 0, // It doesn't really matter, just assume
- // we're on the first face...
- true, false, false // default orientation
- );
+ const unsigned int cell_index =
+ fe.face_to_cell_index(i,
+ 0, /* It doesn't really matter, just assume
+ * we're on the first face...
+ */
+ true, false, false // default orientation
+ );
cell_to_rotated_face_index[cell_index] = i;
}
// Query the correct face_index on face_2 respecting the given
// orientation:
const unsigned int j =
- cell_to_rotated_face_index[fe.
- face_to_cell_index(i,
- 0, // It doesn't really matter, just assume
- // we're on the first face...
- face_orientation,
- face_flip,
- face_rotation)];
+ cell_to_rotated_face_index[fe.face_to_cell_index(i,
+ 0, /* It doesn't really matter, just assume
+ * we're on the first face...
+ */
+ face_orientation,
+ face_flip,
+ face_rotation)];
// And finally constrain the two DoFs respecting component_mask:
- if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components())
- ||
- (component_mask[fe.face_system_to_component_index(i).first] == true))
+ if (component_mask.n_selected_components(fe.n_components()) == fe.n_components()
+ || component_mask[fe.face_system_to_component_index(i).first])
{
if (!constraint_matrix.is_constrained(dofs_1[i]))
{
Assert (0<=direction && direction<space_dim,
ExcIndexRange (direction, 0, space_dim));
Assert ((dynamic_cast<const parallel::distributed::Triangulation<DH::dimension,DH::space_dimension>*>
- (&dof_handler.get_tria())
- ==
- 0),
+ (&dof_handler.get_tria()) == 0),
ExcMessage ("This function can not be used with distributed triangulations."
"See the documentation for more information."));
ExcNotImplemented());
Assert ((dynamic_cast<const parallel::distributed::Triangulation<DH::dimension,DH::space_dimension>*>
- (&dof_handler.get_tria())
- ==
- 0),
+ (&dof_handler.get_tria()) == 0),
ExcMessage ("This function can not be used with distributed triangulations."
"See the documentation for more information."));
}
- // this internal function assigns to each dof
- // the respective component of the vector
- // system.
+ // this internal function assigns to each dof the respective component
+ // of the vector system.
//
- // the output array dofs_by_component
- // lists for each dof the corresponding
- // vector component. if the DoFHandler is
- // based on a parallel distributed
- // triangulation then the output array is
- // index by
+ // the output array dofs_by_component lists for each dof the
+ // corresponding vector component. if the DoFHandler is based on a
+ // parallel distributed triangulation then the output array is index by
// dof.locally_owned_dofs().index_within_set(indices[i])
//
- // if an element is non-primitive then we assign to each degree of freedom
- // the following component:
+ // if an element is non-primitive then we assign to each degree of
+ // freedom the following component:
// - if the nonzero components that belong to a shape function are not
// selected in the component_mask, then the shape function is assigned
// to the first nonzero vector component that corresponds to this
ExcDimensionMismatch(dofs_by_component.size(),
dof.n_locally_owned_dofs()));
- // next set up a table for the degrees
- // of freedom on each of the cells
- // (regardless of the fact whether it
- // is listed in the component_select
- // argument or not)
+ // next set up a table for the degrees of freedom on each of the
+ // cells (regardless of the fact whether it is listed in the
+ // component_select argument or not)
//
- // for each element 'f' of the
- // FECollection,
- // local_component_association[f][d]
- // then returns the vector component
- // that degree of freedom 'd' belongs
- // to
+ // for each element 'f' of the FECollection,
+ // local_component_association[f][d] then returns the vector
+ // component that degree of freedom 'd' belongs to
std::vector<std::vector<unsigned char> >
local_component_association (fe_collection.size());
for (unsigned int f=0; f<fe_collection.size(); ++f)
= get_local_component_association (fe, component_mask);
}
- // then loop over all cells and do
- // the work
+ // then loop over all cells and do the work
std::vector<unsigned int> indices;
for (typename DH::active_cell_iterator c=dof.begin_active();
c!=dof.end(); ++ c)
}
- // this is the function corresponding to the one above but working
- // on blocks instead of components.
+ // this is the function corresponding to the one above but working on
+ // blocks instead of components.
//
- // the output array dofs_by_block
- // lists for each dof the corresponding
- // vector block. if the DoFHandler is
- // based on a parallel distributed
- // triangulation then the output array is
- // index by
+ // the output array dofs_by_block lists for each dof the corresponding
+ // vector block. if the DoFHandler is based on a parallel distributed
+ // triangulation then the output array is index by
// dof.locally_owned_dofs().index_within_set(indices[i])
template <class DH>
inline
ExcDimensionMismatch(dofs_by_block.size(),
dof.n_locally_owned_dofs()));
- // next set up a table for the degrees
- // of freedom on each of the cells
- // (regardless of the fact whether it
- // is listed in the component_select
- // argument or not)
+ // next set up a table for the degrees of freedom on each of the
+ // cells (regardless of the fact whether it is listed in the
+ // component_select argument or not)
//
- // for each element 'f' of the
- // FECollection,
- // local_block_association[f][d]
- // then returns the vector block
- // that degree of freedom 'd' belongs
- // to
+ // for each element 'f' of the FECollection,
+ // local_block_association[f][d] then returns the vector block that
+ // degree of freedom 'd' belongs to
std::vector<std::vector<unsigned char> > local_block_association
(fe_collection.size());
for (unsigned int f=0; f<fe_collection.size(); ++f)
ExcInternalError());
}
- // then loop over all cells and do
- // the work
+ // then loop over all cells and do the work
std::vector<unsigned int> indices;
for (typename DH::active_cell_iterator c=dof.begin_active();
c!=dof.end(); ++ c)
Assert (fe_is_primitive(dof_handler) == true,
typename FiniteElement<dim>::ExcFENotPrimitive());
- // store a flag whether we should care
- // about different components. this is
- // just a simplification, we could ask
- // for this at every single place
- // equally well
+ // store a flag whether we should care about different components. this
+ // is just a simplification, we could ask for this at every single
+ // place equally well
const bool consider_components = (n_components(dof_handler) != 1);
- // zero out the components that we
- // will touch
+ // zero out the components that we will touch
if (consider_components == false)
dof_data = 0;
else
dof_data(i) = 0;
}
- // count how often we have added a value
- // in the sum for each dof
+ // count how often we have added a value in the sum for each dof
std::vector<unsigned char> touch_count (dof_handler.n_dofs(), 0);
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
cell->get_dof_indices (dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
- // consider this dof only if it
- // is the right component. if there
- // is only one component, short cut
- // the test
+ // consider this dof only if it is the right component. if there
+ // is only one component, short cut the test
if (!consider_components ||
(cell->get_fe().system_to_component_index(i).first == component))
{
- // sum up contribution of the
- // present_cell to this dof
+ // sum up contribution of the present_cell to this dof
dof_data(dof_indices[i]) += cell_data(present_cell);
- // note that we added another
- // summand
+
+ // note that we added another summand
++touch_count[dof_indices[i]];
}
}
- // compute the mean value on all the
- // dofs by dividing with the number
+ // compute the mean value on all the dofs by dividing with the number
// of summands.
for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
{
- // assert that each dof was used
- // at least once. this needs not be
- // the case if the vector has more than
- // one component
+ // assert that each dof was used at least once. this needs not be
+ // the case if the vector has more than one component
Assert (consider_components || (touch_count[i]!=0),
ExcInternalError());
if (touch_count[i] != 0)
Assert(selected_dofs.size() == dof.n_locally_owned_dofs(),
ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs()));
- // two special cases: no component
- // is selected, and all components
- // are selected; both rather
- // stupid, but easy to catch
+ // two special cases: no component is selected, and all components are
+ // selected; both rather stupid, but easy to catch
if (component_mask.n_selected_components(n_components(dof)) == 0)
{
std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false);
// preset all values by false
std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false);
- // get the component association of each DoF
- // and then select the ones that match the
- // given set of blocks
+ // get the component association of each DoF and then select the ones
+ // that match the given set of blocks
std::vector<unsigned char> dofs_by_component (dof.n_locally_owned_dofs());
internal::get_component_association (dof, component_mask,
dofs_by_component);
}
-//TODO: Unify the following two functions with the non-hp case
+ // TODO: Unify the following two functions with the non-hp case
template <int dim, int spacedim>
void
Assert(selected_dofs.size() == dof.n_dofs(),
ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs()));
- // two special cases: no component
- // is selected, and all components
- // are selected; both rather
- // stupid, but easy to catch
+ // two special cases: no component is selected, and all components are
+ // selected; both rather stupid, but easy to catch
if (component_mask.n_selected_components(n_components(dof)) == 0)
{
std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
// preset all values by false
std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
- // get the component association of each DoF
- // and then select the ones that match the
- // given set of components
+ // get the component association of each DoF and then select the ones
+ // that match the given set of components
std::vector<unsigned char> dofs_by_component (dof.n_dofs());
internal::get_component_association (dof, component_mask,
dofs_by_component);
Assert(selected_dofs.size() == dof.n_dofs(level),
ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
- // two special cases: no component
- // is selected, and all components
- // are selected, both rather
- // stupid, but easy to catch
+ // two special cases: no component is selected, and all components are
+ // selected, both rather stupid, but easy to catch
if (component_mask.n_selected_components(n_components(dof)) == 0)
{
std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
// preset all values by false
std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
- // next set up a table for the
- // degrees of freedom on each of
- // the cells whether it is
- // something interesting or not
+ // next set up a table for the degrees of freedom on each of the cells
+ // whether it is something interesting or not
std::vector<unsigned char> local_component_asssociation
= internal::get_local_component_association (fe, component_mask);
std::vector<bool> local_selected_dofs (fe.dofs_per_cell);
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
local_selected_dofs[i] = component_mask[local_component_asssociation[i]];
- // then loop over all cells and do
- // work
+ // then loop over all cells and do work
std::vector<unsigned int> indices(fe.dofs_per_cell);
typename MGDoFHandler<dim,spacedim>::cell_iterator c;
for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
selected_dofs.clear ();
selected_dofs.set_size(dof_handler.n_dofs());
- // let's see whether we have to
- // check for certain boundary
- // indicators or whether we can
- // accept all
+ // let's see whether we have to check for certain boundary indicators
+ // or whether we can accept all
const bool check_boundary_indicator = (boundary_indicators.size() != 0);
- // also see whether we have to
- // check whether a certain vector
- // component is selected, or all
+ // also see whether we have to check whether a certain vector component
+ // is selected, or all
const bool check_vector_component
= ((component_mask.represents_the_all_selected_mask() == false)
||
std::vector<unsigned int> face_dof_indices;
face_dof_indices.reserve (max_dofs_per_face(dof_handler));
- // now loop over all cells and
- // check whether their faces are at
- // the boundary. note that we need
- // not take special care of single
- // lines being at the boundary
- // (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every isolated boundary
- // line is also part of a boundary
- // face which we will be visiting
- // sooner or later
+ // now loop over all cells and check whether their faces are at the
+ // boundary. note that we need not take special care of single lines
+ // being at the boundary (using @p{cell->has_boundary_lines}), since we
+ // do not support boundaries of dimension dim-2, and so every isolated
+ // boundary line is also part of a boundary face which we will be
+ // visiting sooner or later
for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell)
- // only work on cells that are either locally owned or at
- // least ghost cells
+
+ // only work on cells that are either locally owned or at least ghost
+ // cells
if (cell->is_artificial() == false)
for (unsigned int face=0;
face<GeometryInfo<DH::dimension>::faces_per_cell; ++face)
if (!check_vector_component)
selected_dofs.add_index (face_dof_indices[i]);
else
- // check for
- // component is
- // required. somewhat
- // tricky as usual
- // for the case that
- // the shape function
- // is non-primitive,
- // but use usual
- // convention (see
- // docs)
+ // check for component is required. somewhat tricky as
+ // usual for the case that the shape function is
+ // non-primitive, but use usual convention (see docs)
{
- // first get at the
- // cell-global
- // number of a face
- // dof, to ask the
- // fe certain
- // questions
+ // first get at the cell-global number of a face dof,
+ // to ask the fe certain questions
const unsigned int cell_index
= (dim == 1 ?
i
Assert (boundary_indicators.find (numbers::internal_face_boundary_id) == boundary_indicators.end(),
ExcInvalidBoundaryIndicator());
- // let's see whether we have to
- // check for certain boundary
- // indicators or whether we can
- // accept all
+ // let's see whether we have to check for certain boundary indicators
+ // or whether we can accept all
const bool check_boundary_indicator = (boundary_indicators.size() != 0);
- // also see whether we have to
- // check whether a certain vector
- // component is selected, or all
+ // also see whether we have to check whether a certain vector component
+ // is selected, or all
const bool check_vector_component
= (component_mask.represents_the_all_selected_mask() == false);
- // clear and reset array by default
- // values
+ // clear and reset array by default values
selected_dofs.clear ();
selected_dofs.resize (dof_handler.n_dofs(), false);
std::vector<unsigned int> cell_dof_indices;
cell_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // now loop over all cells and
- // check whether their faces are at
- // the boundary. note that we need
- // not take special care of single
- // lines being at the boundary
- // (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every isolated boundary
- // line is also part of a boundary
- // face which we will be visiting
- // sooner or later
+ // now loop over all cells and check whether their faces are at the
+ // boundary. note that we need not take special care of single lines
+ // being at the boundary (using @p{cell->has_boundary_lines}), since we
+ // do not support boundaries of dimension dim-2, and so every isolated
+ // boundary line is also part of a boundary face which we will be
+ // visiting sooner or later
for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell)
for (unsigned int face=0;
if (!check_vector_component)
selected_dofs[cell_dof_indices[i]] = true;
else
- // check for
- // component is
- // required. somewhat
- // tricky as usual
- // for the case that
- // the shape function
- // is non-primitive,
- // but use usual
- // convention (see
- // docs)
+ // check for component is required. somewhat tricky
+ // as usual for the case that the shape function is
+ // non-primitive, but use usual convention (see docs)
{
if (fe.is_primitive (i))
selected_dofs[cell_dof_indices[i]]
const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
- // this function is similar to the
- // make_sparsity_pattern function,
+ // this function is similar to the make_sparsity_pattern function,
// see there for more information
typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
cell = dof_handler.begin_active(),
const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
- // this function is similar to the
- // make_sparsity_pattern function,
+ // this function is similar to the make_sparsity_pattern function,
// see there for more information
typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true;
- // dof numbers on the centers of
- // the lines bounding this face
+ // dof numbers on the centers of the lines bounding this
+ // face
for (unsigned int line=0; line<4; ++line)
for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true;
- // next the dofs on the lines interior
- // to the face; the order of these
- // lines is laid down in the
- // FiniteElement class documentation
+ // next the dofs on the lines interior to the face; the
+ // order of these lines is laid down in the FiniteElement
+ // class documentation
for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
selected_dofs[face->child(0)->line(1)->dof_index(dof)] = true;
for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true;
- // finally, for the dofs interior
- // to the four child faces
+ // finally, for the dofs interior to the four child faces
for (unsigned int child=0; child<4; ++child)
for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
selected_dofs[face->child(child)->dof_index(dof)] = true;
std::vector<unsigned int> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // this function is similar to the
- // make_sparsity_pattern function,
- // see there for more information
+ // this function is similar to the make_sparsity_pattern function, see
+ // there for more information
typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
- // add the DoF on the adjacent ghost cells
- // to the IndexSet, cache them in a
- // set. need to check each dof manually
- // because we can't be sure that the dof
- // range of locally_owned_dofs is really
- // contiguous.
+ // add the DoF on the adjacent ghost cells to the IndexSet, cache them
+ // in a set. need to check each dof manually because we can't be sure
+ // that the dof range of locally_owned_dofs is really contiguous.
std::vector<unsigned int> dof_indices;
std::set<unsigned int> global_dof_indices;
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
- // add the DoF on the adjacent ghost cells
- // to the IndexSet, cache them in a
- // set. need to check each dof manually
- // because we can't be sure that the dof
- // range of locally_owned_dofs is really
- // contiguous.
+ // add the DoF on the adjacent ghost cells to the IndexSet, cache them
+ // in a set. need to check each dof manually because we can't be sure
+ // that the dof range of locally_owned_dofs is really contiguous.
std::vector<unsigned int> dof_indices;
std::set<unsigned int> global_dof_indices;
n_selected_dofs += std::count (dofs_by_component.begin(),
dofs_by_component.end(), i);
- // First count the number of dofs
- // in the current component.
+ // First count the number of dofs in the current component.
constant_modes.resize (n_components_selected, std::vector<bool>(n_selected_dofs,
false));
std::vector<unsigned int> component_list (n_components, 0);
get_subdomain_association (const DH &dof_handler,
std::vector<types::subdomain_id> &subdomain_association)
{
- // if the Triangulation is distributed, the
- // only thing we can usefully ask is for
- // its locally owned subdomain
+ // if the Triangulation is distributed, the only thing we can usefully
+ // ask is for its locally owned subdomain
Assert ((dynamic_cast<const parallel::distributed::
Triangulation<DH::dimension,DH::space_dimension> *>
(&dof_handler.get_tria()) == 0),
std::vector<unsigned int> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // pseudo-randomly assign variables
- // which lie on the interface
- // between subdomains to each of
- // the two or more
+ // pseudo-randomly assign variables which lie on the interface between
+ // subdomains to each of the two or more
bool coin_flip = true;
- // loop over all cells and record
- // which subdomain a DoF belongs
- // to. toss a coin in case it is on
- // an interface
+ // loop over all cells and record which subdomain a DoF belongs to.
+ // toss a coin in case it is on an interface
typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
local_dof_indices.resize (dofs_per_cell);
cell->get_dof_indices (local_dof_indices);
- // set subdomain ids. if dofs
- // already have their values
- // set then they must be on
- // partition interfaces. in
- // that case randomly assign
- // them to either the previous
- // association or the current
- // one, where we take "random"
- // to be "once this way once
- // that way"
+ // set subdomain ids. if dofs already have their values set then
+ // they must be on partition interfaces. in that case randomly
+ // assign them to either the previous association or the current
+ // one, where we take "random" to be "once this way once that way"
for (unsigned int i=0; i<dofs_per_cell; ++i)
if (subdomain_association[local_dof_indices[i]] ==
numbers::invalid_unsigned_int)
const types::subdomain_id subdomain)
{
- // If we have a distributed::Triangulation only
- // allow locally_owned subdomain.
+ // If we have a distributed::Triangulation only allow locally_owned
+ // subdomain.
Assert (
(dof_handler.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id)
||
std::vector<unsigned int> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // first generate an unsorted list of all
- // indices which we fill from the back. could
- // also insert them directly into the
- // IndexSet, but that inserts indices in the
- // middle, which is an O(n^2) algorithm and
- // hence too expensive. Could also use
- // std::set, but that is in general more
- // expensive than a vector
+ // first generate an unsorted list of all indices which we fill from
+ // the back. could also insert them directly into the IndexSet, but
+ // that inserts indices in the middle, which is an O(n^2) algorithm and
+ // hence too expensive. Could also use std::set, but that is in general
+ // more expensive than a vector
std::vector<unsigned int> subdomain_indices;
typename DH::active_cell_iterator
dof_handler.get_fe().n_components()));
std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
- // in debug mode, make sure that there are
- // some cells at least with this subdomain
- // id
+ // in debug mode, make sure that there are some cells at least with
+ // this subdomain id
#ifdef DEBUG
{
bool found = false;
namespace internal
{
-//TODO: why is this function so complicated? It would be nice to have
-// comments that explain why we can't just loop over all components
-// and count the entries in dofs_by_component that have this
-// component's index
+ // TODO: why is this function so complicated? It would be nice to have
+ // comments that explain why we can't just loop over all components and
+ // count the entries in dofs_by_component that have this component's
+ // index
template <int dim, int spacedim>
void
resolve_components (const FiniteElement<dim,spacedim> &fe,
dofs_by_component.end(),
component);
- // if we have non-primitive FEs and want all
- // components to show the number of dofs, need
- // to copy the result to those components
+ // if we have non-primitive FEs and want all components
+ // to show the number of dofs, need to copy the result to
+ // those components
if (!base.is_primitive() && !only_once)
for (unsigned int dd=1; dd<d; ++dd)
dofs_per_component[target_component[component-d+dd]] =
}
}
+
template <int dim, int spacedim>
void
resolve_components (const hp::FECollection<dim,spacedim> &fe_collection,
std::vector<unsigned int> &dofs_per_component,
unsigned int &component)
{
- // assert that all elements in the collection have the same
- // structure (base elements and multiplicity, components per base
- // element) and then simply call the function above
+ // assert that all elements in the collection have the same structure
+ // (base elements and multiplicity, components per base element) and
+ // then simply call the function above
for (unsigned int fe=1; fe<fe_collection.size(); ++fe)
{
Assert (fe_collection[fe].n_components() == fe_collection[0].n_components(),
}
+
namespace internal
{
namespace
std::fill (dofs_per_component.begin(), dofs_per_component.end(), 0U);
- // If the empty vector was given as
- // default argument, set up this
+ // If the empty vector was given as default argument, set up this
// vector as identity.
if (target_component.size()==0)
{
AssertDimension (dofs_per_component.size(), n_target_components);
- // special case for only one
- // component. treat this first
- // since it does not require any
- // computations
+ // special case for only one component. treat this first since it does
+ // not require any computations
if (n_components == 1)
{
dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
}
- // otherwise determine the number
- // of dofs in each component
- // separately. do so in parallel
+ // otherwise determine the number of dofs in each component separately.
+ // do so in parallel
std::vector<unsigned char> dofs_by_component (dof_handler.n_locally_owned_dofs());
internal::get_component_association (dof_handler, ComponentMask(),
dofs_by_component);
only_once, dofs_per_component, component);
Assert (n_components == component, ExcInternalError());
- // finally sanity check. this is
- // only valid if the finite element
- // is actually primitive, so
- // exclude other elements from this
+ // finally sanity check. this is only valid if the finite element is
+ // actually primitive, so exclude other elements from this
Assert ((internal::all_elements_are_primitive(dof_handler.get_fe()) == false)
||
(std::accumulate (dofs_per_component.begin(),
const FiniteElement<DH::dimension,DH::space_dimension> &fe = fe_collection[this_fe];
std::fill (dofs_per_block.begin(), dofs_per_block.end(), 0U);
- // If the empty vector was given as
- // default argument, set up this
+ // If the empty vector was given as default argument, set up this
// vector as identity.
if (target_block.size()==0)
{
AssertDimension (dofs_per_block.size(), n_target_blocks);
- // special case for only one
- // block. treat this first
- // since it does not require any
- // computations
+ // special case for only one block. treat this first since it does
+ // not require any computations
if (n_blocks == 1)
{
dofs_per_block[0] = dof_handler.n_dofs();
return;
}
- // otherwise determine the number
- // of dofs in each block
- // separately.
+ // otherwise determine the number of dofs in each block separately.
std::vector<unsigned char> dofs_by_block (dof_handler.n_locally_owned_dofs());
internal::get_block_association (dof_handler, dofs_by_block);
#ifdef DEAL_II_USE_P4EST
#if DEAL_II_COMPILER_SUPPORTS_MPI
- // if we are working on a parallel
- // mesh, we now need to collect
- // this information from all
- // processors
+ // if we are working on a parallel mesh, we now need to collect
+ // this information from all processors
if (const parallel::distributed::Triangulation<DH::dimension,DH::space_dimension> *tria
= (dynamic_cast<const parallel::distributed::Triangulation<DH::dimension,DH::space_dimension>*>
(&dof_handler.get_tria())))
namespace
{
/**
- * This is a function that is
- * called by the _2 function and
- * that operates on a range of
- * cells only. It is used to
- * split up the whole range of
- * cells into chunks which are
- * then worked on in parallel, if
- * multithreading is available.
+ * This is a function that is called by the _2 function and that
+ * operates on a range of cells only. It is used to split up the
+ * whole range of cells into chunks which are then worked on in
+ * parallel, if multithreading is available.
*/
template <int dim, int spacedim>
void
const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &begin,
const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &end)
{
- // aliases to the finite elements
- // used by the dof handlers:
+ // aliases to the finite elements used by the dof handlers:
const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe();
- // for each cell on the parameter grid:
- // find out which degrees of freedom on the
- // fine grid correspond in which way to
- // the degrees of freedom on the parameter
- // grid
+ // for each cell on the parameter grid: find out which degrees of
+ // freedom on the fine grid correspond in which way to the degrees
+ // of freedom on the parameter grid
//
- // since for continuous FEs some
- // dofs exist on more than one
- // cell, we have to track which
- // ones were already visited. the
- // problem is that if we visit a
- // dof first on one cell and
- // compute its weight with respect
- // to some global dofs to be
- // non-zero, and later visit the
- // dof again on another cell and
- // (since we are on another cell)
- // recompute the weights with
- // respect to the same dofs as
- // above to be zero now, we have to
- // preserve them. we therefore
- // overwrite all weights if they
- // are nonzero and do not enforce
- // zero weights since that might be
- // only due to the fact that we are
- // on another cell.
+ // since for continuous FEs some dofs exist on more than one cell,
+ // we have to track which ones were already visited. the problem is
+ // that if we visit a dof first on one cell and compute its weight
+ // with respect to some global dofs to be non-zero, and later visit
+ // the dof again on another cell and (since we are on another cell)
+ // recompute the weights with respect to the same dofs as above to
+ // be zero now, we have to preserve them. we therefore overwrite
+ // all weights if they are nonzero and do not enforce zero weights
+ // since that might be only due to the fact that we are on another
+ // cell.
//
// example:
// coarse grid
// | | | | |
// *--x--y--*--*
//
- // when on cell 1, we compute the
- // weights of dof 'x' to be 1/2
- // from parameter dofs 0 and 1,
- // respectively. however, when
- // later we are on cell 2, we again
- // compute the prolongation of
- // shape function 1 restricted to
- // cell 2 to the globla grid and
- // find that the weight of global
- // dof 'x' now is zero. however, we
- // should not overwrite the old
- // value.
+ // when on cell 1, we compute the weights of dof 'x' to be 1/2 from
+ // parameter dofs 0 and 1, respectively. however, when later we are
+ // on cell 2, we again compute the prolongation of shape function 1
+ // restricted to cell 2 to the globla grid and find that the weight
+ // of global dof 'x' now is zero. however, we should not overwrite
+ // the old value.
//
- // we therefore always only set
- // nonzero values. why adding up is
- // not useful: dof 'y' would get
- // weight 1 from parameter dof 1 on
- // both cells 1 and 2, but the
- // correct weight is nevertheless
- // only 1.
-
- // vector to hold the representation of
- // a single degree of freedom on the
- // coarse grid (for the selected fe)
- // on the fine grid
+ // we therefore always only set nonzero values. why adding up is
+ // not useful: dof 'y' would get weight 1 from parameter dof 1 on
+ // both cells 1 and 2, but the correct weight is nevertheless only
+ // 1.
+
+ // vector to hold the representation of a single degree of freedom
+ // on the coarse grid (for the selected fe) on the fine grid
const unsigned int n_fine_dofs = weight_mapping.size();
dealii::Vector<double> global_parameter_representation (n_fine_dofs);
for (cell=begin; cell!=end; ++cell)
{
- // get the global indices of the
- // parameter dofs on this parameter
- // grid cell
+ // get the global indices of the parameter dofs on this
+ // parameter grid cell
cell->get_dof_indices (parameter_dof_indices);
- // loop over all dofs on this
- // cell and check whether they
- // are interesting for us
+ // loop over all dofs on this cell and check whether they are
+ // interesting for us
for (unsigned int local_dof=0;
local_dof<coarse_fe.dofs_per_cell;
++local_dof)
==
coarse_component)
{
- // the how-many-th
- // parameter is this on
- // this cell?
+ // the how-many-th parameter is this on this cell?
const unsigned int local_parameter_dof
= coarse_fe.system_to_component_index(local_dof).second;
global_parameter_representation = 0;
// distribute the representation of
- // @p{local_parameter_dof} on the
- // parameter grid cell @p{cell} to
- // the global data space
+ // @p{local_parameter_dof} on the parameter grid cell
+ // @p{cell} to the global data space
coarse_to_fine_grid_map[cell]->
set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
global_parameter_representation);
- // now that we've got the global
- // representation of each parameter
- // dof, we've only got to clobber the
- // non-zero entries in that vector and
- // store the result
+ // now that we've got the global representation of each
+ // parameter dof, we've only got to clobber the non-zero
+ // entries in that vector and store the result
//
- // what we have learned: if entry @p{i}
- // of the global vector holds the value
- // @p{v[i]}, then this is the weight with
- // which the present dof contributes
- // to @p{i}. there may be several such
- // @p{i}s and their weights' sum should
- // be one. Then, @p{v[i]} should be
- // equal to @p{\sum_j w_{ij} p[j]} with
- // @p{p[j]} be the values of the degrees
- // of freedom on the coarse grid. we
- // can thus compute constraints which
- // link the degrees of freedom @p{v[i]}
- // on the fine grid to those on the
- // coarse grid, @p{p[j]}. Now to use
- // these as real constraints, rather
- // than as additional equations, we
- // have to identify representants
- // among the @p{i} for each @p{j}. this will
- // be done by simply taking the first
- // @p{i} for which @p{w_{ij}==1}.
+ // what we have learned: if entry @p{i} of the global
+ // vector holds the value @p{v[i]}, then this is the
+ // weight with which the present dof contributes to
+ // @p{i}. there may be several such @p{i}s and their
+ // weights' sum should be one. Then, @p{v[i]} should be
+ // equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the
+ // values of the degrees of freedom on the coarse grid.
+ // we can thus compute constraints which link the degrees
+ // of freedom @p{v[i]} on the fine grid to those on the
+ // coarse grid, @p{p[j]}. Now to use these as real
+ // constraints, rather than as additional equations, we
+ // have to identify representants among the @p{i} for
+ // each @p{j}. this will be done by simply taking the
+ // first @p{i} for which @p{w_{ij}==1}.
//
- // guard modification of
- // the weights array by a
- // Mutex. since it should
- // happen rather rarely
- // that there are several
- // threads operating on
- // different intergrid
- // weights, have only one
- // mutex for all of them
+ // guard modification of the weights array by a Mutex.
+ // since it should happen rather rarely that there are
+ // several threads operating on different intergrid
+ // weights, have only one mutex for all of them
static Threads::ThreadMutex mutex;
Threads::ThreadMutex::ScopedLock lock (mutex);
for (unsigned int i=0; i<global_parameter_representation.size(); ++i)
- // set this weight if it belongs
- // to a parameter dof.
+ // set this weight if it belongs to a parameter dof.
if (weight_mapping[i] != -1)
{
- // only overwrite old
- // value if not by
- // zero
+ // only overwrite old value if not by zero
if (global_parameter_representation(i) != 0)
{
const unsigned int wi = parameter_dof_indices[local_dof],
/**
- * This is a helper function that
- * is used in the computation of
- * integrid constraints. See the
- * function for a thorough
- * description of how it works.
+ * This is a helper function that is used in the computation of
+ * integrid constraints. See the function for a thorough description
+ * of how it works.
*/
template <int dim, int spacedim>
void
const std::vector<int> &weight_mapping,
std::vector<std::map<unsigned int,float> > &weights)
{
- // simply distribute the range of
- // cells to different threads
+ // simply distribute the range of cells to different threads
typedef typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator active_cell_iterator;
std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
cell_intervals = Threads::split_range<active_cell_iterator> (coarse_grid.begin_active(),
coarse_grid.end(),
multithread_info.n_default_threads);
-//TODO: use WorkStream here
+ // TODO: use WorkStream here
+
Threads::TaskGroup<> tasks;
void (*fun_ptr) (const dealii::DoFHandler<dim,spacedim> &,
const unsigned int ,
/**
- * This is a helper function that
- * is used in the computation of
- * integrid constraints. See the
- * function for a thorough
- * description of how it works.
+ * This is a helper function that is used in the computation of
+ * integrid constraints. See the function for a thorough description
+ * of how it works.
*/
template <int dim, int spacedim>
unsigned int
std::vector<std::map<unsigned int, float> > &weights,
std::vector<int> &weight_mapping)
{
- // aliases to the finite elements
- // used by the dof handlers:
+ // aliases to the finite elements used by the dof handlers:
const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe(),
&fine_fe = fine_grid.get_fe();
// local numbers of dofs
const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell;
- // alias the number of dofs per
- // cell belonging to the
- // coarse_component which is to be
- // the restriction of the fine
+ // alias the number of dofs per cell belonging to the
+ // coarse_component which is to be the restriction of the fine
// grid:
const unsigned int coarse_dofs_per_cell_component
= coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell;
- // Try to find out whether the
- // grids stem from the same coarse
- // grid. This is a rather crude
- // test, but better than nothing
+ // Try to find out whether the grids stem from the same coarse
+ // grid. This is a rather crude test, but better than nothing
Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0),
ExcGridsDontMatch());
- // check whether the map correlates
- // the right objects
+ // check whether the map correlates the right objects
Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
ExcGridsDontMatch ());
Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
ExcGridsDontMatch ());
- // check whether component numbers
- // are valid
+ // check whether component numbers are valid
AssertIndexRange (coarse_component,coarse_fe.n_components());
AssertIndexRange (fine_component, fine_fe.n_components());
- // check whether respective finite
- // elements are equal
+
+ // check whether respective finite elements are equal
Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first)
==
fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first),
ExcFiniteElementsDontMatch());
#ifdef DEBUG
- // if in debug mode, check whether
- // the coarse grid is indeed
- // coarser everywhere than the fine
- // grid
+ // if in debug mode, check whether the coarse grid is indeed
+ // coarser everywhere than the fine grid
for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
cell=coarse_grid.begin_active();
cell != coarse_grid.end(); ++cell)
ExcGridNotCoarser());
#endif
-
-
/*
- * From here on: the term `parameter' refers to the selected component
- * on the coarse grid and its analogon on the fine grid. The naming of
- * variables containing this term is due to the fact that
- * `selected_component' is longer, but also due to the fact that the
- * code of this function was initially written for a program where the
- * component which we wanted to match between grids was actually the
- * `parameter' variable.
+ * From here on: the term `parameter' refers to the selected
+ * component on the coarse grid and its analogon on the fine grid.
+ * The naming of variables containing this term is due to the fact
+ * that `selected_component' is longer, but also due to the fact
+ * that the code of this function was initially written for a
+ * program where the component which we wanted to match between
+ * grids was actually the `parameter' variable.
*
- * Likewise, the terms `parameter grid' and `state grid' refer to the
- * coarse and fine grids, respectively.
+ * Likewise, the terms `parameter grid' and `state grid' refer to
+ * the coarse and fine grids, respectively.
*
- * Changing the names of variables would in principle be a good idea,
- * but would not make things simpler and would be another source of
- * errors. If anyone feels like doing so: patches would be welcome!
+ * Changing the names of variables would in principle be a good
+ * idea, but would not make things simpler and would be another
+ * source of errors. If anyone feels like doing so: patches would
+ * be welcome!
*/
- // set up vectors of cell-local
- // data; each vector represents one
- // degree of freedom of the
- // coarse-grid variable in the
- // fine-grid element
+ // set up vectors of cell-local data; each vector represents one
+ // degree of freedom of the coarse-grid variable in the fine-grid
+ // element
std::vector<dealii::Vector<double> >
parameter_dofs (coarse_dofs_per_cell_component,
dealii::Vector<double>(fine_dofs_per_cell));
- // for each coarse dof: find its
- // position within the fine element
- // and set this value to one in the
- // respective vector (all other values
- // are zero by construction)
+ // for each coarse dof: find its position within the fine element
+ // and set this value to one in the respective vector (all other
+ // values are zero by construction)
for (unsigned int local_coarse_dof=0;
local_coarse_dof<coarse_dofs_per_cell_component;
++local_coarse_dof)
};
- // find out how many DoFs there are
- // on the grids belonging to the
+ // find out how many DoFs there are on the grids belonging to the
// components we want to match
unsigned int n_parameters_on_fine_grid=0;
if (true)
{
- // have a flag for each dof on
- // the fine grid and set it
- // to true if this is an
- // interesting dof. finally count
- // how many true's there
+ // have a flag for each dof on the fine grid and set it to true
+ // if this is an interesting dof. finally count how many true's
+ // there
std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
std::vector<unsigned int> local_dof_indices (fine_fe.dofs_per_cell);
{
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
- // if this DoF is a
- // parameter dof and has
- // not yet been numbered,
- // then do so
+ // if this DoF is a parameter dof and has not yet been
+ // numbered, then do so
if ((fine_fe.system_to_component_index(i).first == fine_component) &&
(weight_mapping[local_dof_indices[i]] == -1))
{
};
- // for each cell on the parameter grid:
- // find out which degrees of freedom on the
- // fine grid correspond in which way to
- // the degrees of freedom on the parameter
- // grid
+ // for each cell on the parameter grid: find out which degrees of
+ // freedom on the fine grid correspond in which way to the degrees
+ // of freedom on the parameter grid
//
- // do this in a separate function
- // to allow for multithreading
- // there. see this function also if
- // you want to read more
- // information on the algorithm
- // used.
+ // do this in a separate function to allow for multithreading
+ // there. see this function also if you want to read more
+ // information on the algorithm used.
compute_intergrid_weights_2 (coarse_grid, coarse_component,
coarse_to_fine_grid_map, parameter_dofs,
weight_mapping, weights);
- // ok, now we have all weights for each
- // dof on the fine grid. if in debug
- // mode lets see if everything went smooth,
- // i.e. each dof has sum of weights one
+ // ok, now we have all weights for each dof on the fine grid. if in
+ // debug mode lets see if everything went smooth, i.e. each dof has
+ // sum of weights one
//
- // in other words this means that
- // if the sum of all shape
- // functions on the parameter grid
- // is one (which is always the
- // case), then the representation
- // on the state grid should be as
- // well (division of unity)
+ // in other words this means that if the sum of all shape functions
+ // on the parameter grid is one (which is always the case), then
+ // the representation on the state grid should be as well (division
+ // of unity)
//
- // if the parameter grid has more
- // than one component, then the
- // respective dofs of the other
- // components have sum of weights
- // zero, of course. we do not
- // explicitly ask which component
- // a dof belongs to, but this at
- // least tests some errors
+ // if the parameter grid has more than one component, then the
+ // respective dofs of the other components have sum of weights
+ // zero, of course. we do not explicitly ask which component a dof
+ // belongs to, but this at least tests some errors
#ifdef DEBUG
for (unsigned int col=0; col<n_parameters_on_fine_grid; ++col)
{
const InterGridMap<DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
ConstraintMatrix &constraints)
{
- // store the weights with which a dof
- // on the parameter grid contributes
- // to a dof on the fine grid. see the
- // long doc below for more info
+ // store the weights with which a dof on the parameter grid contributes
+ // to a dof on the fine grid. see the long doc below for more info
//
- // allocate as many rows as there are
- // parameter dofs on the coarse grid
- // and as many columns as there are
- // parameter dofs on the fine grid.
+ // allocate as many rows as there are parameter dofs on the coarse grid
+ // and as many columns as there are parameter dofs on the fine grid.
//
- // weight_mapping is used to map the
- // global (fine grid) parameter dof
+ // weight_mapping is used to map the global (fine grid) parameter dof
// indices to the columns
//
- // in the original implementation,
- // the weights array was actually
- // of FullMatrix<double> type. this
- // wasted huge amounts of memory,
- // but was fast. nonetheless, since
- // the memory consumption was
- // quadratic in the number of
- // degrees of freedom, this was not
- // very practical, so we now use a
- // vector of rows of the matrix,
- // and in each row a vector of
- // pairs (colnum,value). this seems
- // like the best tradeoff between
- // memory and speed, as it is now
- // linear in memory and still fast
+ // in the original implementation, the weights array was actually of
+ // FullMatrix<double> type. this wasted huge amounts of memory, but was
+ // fast. nonetheless, since the memory consumption was quadratic in the
+ // number of degrees of freedom, this was not very practical, so we now
+ // use a vector of rows of the matrix, and in each row a vector of
+ // pairs (colnum,value). this seems like the best tradeoff between
+ // memory and speed, as it is now linear in memory and still fast
// enough.
//
- // to save some memory and since
- // the weights are usually
- // (negative) powers of 2, we
- // choose the value type of the
- // matrix to be @p{float} rather
- // than @p{double}.
+ // to save some memory and since the weights are usually (negative)
+ // powers of 2, we choose the value type of the matrix to be @p{float}
+ // rather than @p{double}.
std::vector<std::map<unsigned int, float> > weights;
- // this is this mapping. there is one
- // entry for each dof on the fine grid;
- // if it is a parameter dof, then its
- // value is the column in weights for
- // that parameter dof, if it is any
- // other dof, then its value is -1,
- // indicating an error
+ // this is this mapping. there is one entry for each dof on the fine
+ // grid; if it is a parameter dof, then its value is the column in
+ // weights for that parameter dof, if it is any other dof, then its
+ // value is -1, indicating an error
std::vector<int> weight_mapping;
const unsigned int n_parameters_on_fine_grid
n_fine_dofs = fine_grid.n_dofs();
- // get an array in which we store
- // which dof on the coarse grid is
- // a parameter and which is not
+ // get an array in which we store which dof on the coarse grid is a
+ // parameter and which is not
std::vector<bool> coarse_dof_is_parameter (coarse_grid.n_dofs());
if (true)
{
extract_dofs (coarse_grid, ComponentMask(mask), coarse_dof_is_parameter);
}
- // now we know that the weights in
- // each row constitute a
- // constraint. enter this into the
- // constraints object
+ // now we know that the weights in each row constitute a constraint.
+ // enter this into the constraints object
//
- // first task: for each parameter
- // dof on the parameter grid, find
- // a representant on the fine,
- // global grid. this is possible
- // since we use conforming finite
- // element. we take this
- // representant to be the first
- // element in this row with weight
- // identical to one. the
- // representant will become an
- // unconstrained degree of freedom,
- // while all others will be
- // constrained to this dof (and
- // possibly others)
+ // first task: for each parameter dof on the parameter grid, find a
+ // representant on the fine, global grid. this is possible since we use
+ // conforming finite element. we take this representant to be the first
+ // element in this row with weight identical to one. the representant
+ // will become an unconstrained degree of freedom, while all others
+ // will be constrained to this dof (and possibly others)
std::vector<int> representants(n_coarse_dofs, -1);
for (unsigned int parameter_dof=0; parameter_dof<n_coarse_dofs;
++parameter_dof)
if (coarse_dof_is_parameter[parameter_dof] == true)
{
- // if this is the line of a
- // parameter dof on the
- // coarse grid, then it
- // should have at least one
- // dependent node on the fine
+ // if this is the line of a parameter dof on the coarse grid,
+ // then it should have at least one dependent node on the fine
// grid
Assert (weights[parameter_dof].size() > 0, ExcInternalError());
- // find the column where the
- // representant is mentioned
+ // find the column where the representant is mentioned
std::map<unsigned int,float>::const_iterator i = weights[parameter_dof].begin();
for (; i!=weights[parameter_dof].end(); ++i)
if (i->second == 1)
Assert (i!=weights[parameter_dof].end(), ExcInternalError());
const unsigned int column = i->first;
- // now we know in which column of
- // weights the representant is, but
- // we don't know its global index. get
- // it using the inverse operation of
- // the weight_mapping
+ // now we know in which column of weights the representant is,
+ // but we don't know its global index. get it using the inverse
+ // operation of the weight_mapping
unsigned int global_dof=0;
for (; global_dof<weight_mapping.size(); ++global_dof)
if (weight_mapping[global_dof] == static_cast<int>(column))
break;
Assert (global_dof < weight_mapping.size(), ExcInternalError());
- // now enter the representants global
- // index into our list
+ // now enter the representants global index into our list
representants[parameter_dof] = global_dof;
}
else
{
- // consistency check: if this
- // is no parameter dof on the
- // coarse grid, then the
- // respective row must be
- // empty!
+ // consistency check: if this is no parameter dof on the coarse
+ // grid, then the respective row must be empty!
Assert (weights[parameter_dof].size() == 0, ExcInternalError());
};
- // note for people that want to
- // optimize this function: the
- // largest part of the computing
- // time is spent in the following,
- // rather innocent block of
- // code. basically, it must be the
- // ConstraintMatrix::add_entry call
- // which takes the bulk of the
- // time, but it is not known to the
- // author how to make it faster...
+ // note for people that want to optimize this function: the largest
+ // part of the computing time is spent in the following, rather
+ // innocent block of code. basically, it must be the
+ // ConstraintMatrix::add_entry call which takes the bulk of the time,
+ // but it is not known to the author how to make it faster...
std::vector<std::pair<unsigned int,double> > constraint_line;
for (unsigned int global_dof=0; global_dof<n_fine_dofs; ++global_dof)
if (weight_mapping[global_dof] != -1)
- // this global dof is a parameter
- // dof, so it may carry a constraint
- // note that for each global dof,
- // the sum of weights shall be one,
- // so we can find out whether this
- // dof is constrained in the following
- // way: if the only weight in this row
- // is a one, and the representant for
- // the parameter dof of the line in
- // which this one is is the present
- // dof, then we consider this dof
- // to be unconstrained. otherwise,
- // all other dofs are constrained
+ // this global dof is a parameter dof, so it may carry a constraint
+ // note that for each global dof, the sum of weights shall be one,
+ // so we can find out whether this dof is constrained in the
+ // following way: if the only weight in this row is a one, and the
+ // representant for the parameter dof of the line in which this one
+ // is is the present dof, then we consider this dof to be
+ // unconstrained. otherwise, all other dofs are constrained
{
const unsigned int col = weight_mapping[global_dof];
Assert (col < n_parameters_on_fine_grid, ExcInternalError());
if ((col_entry->second == 1) &&
(representants[first_used_row] == static_cast<int>(global_dof)))
- // dof unconstrained or
- // constrained to itself
- // (in case this cell is
- // mapped to itself, rather
- // than to children of
+ // dof unconstrained or constrained to itself (in case this
+ // cell is mapped to itself, rather than to children of
// itself)
continue;
}
const InterGridMap<DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
std::vector<std::map<unsigned int, float> > &transfer_representation)
{
- // store the weights with which a dof
- // on the parameter grid contributes
- // to a dof on the fine grid. see the
- // long doc below for more info
+ // store the weights with which a dof on the parameter grid contributes
+ // to a dof on the fine grid. see the long doc below for more info
//
- // allocate as many rows as there are
- // parameter dofs on the coarse grid
- // and as many columns as there are
- // parameter dofs on the fine grid.
+ // allocate as many rows as there are parameter dofs on the coarse grid
+ // and as many columns as there are parameter dofs on the fine grid.
//
- // weight_mapping is used to map the
- // global (fine grid) parameter dof
+ // weight_mapping is used to map the global (fine grid) parameter dof
// indices to the columns
//
- // in the original implementation,
- // the weights array was actually
- // of FullMatrix<double> type. this
- // wasted huge amounts of memory,
- // but was fast. nonetheless, since
- // the memory consumption was
- // quadratic in the number of
- // degrees of freedom, this was not
- // very practical, so we now use a
- // vector of rows of the matrix,
- // and in each row a vector of
- // pairs (colnum,value). this seems
- // like the best tradeoff between
- // memory and speed, as it is now
- // linear in memory and still fast
+ // in the original implementation, the weights array was actually of
+ // FullMatrix<double> type. this wasted huge amounts of memory, but was
+ // fast. nonetheless, since the memory consumption was quadratic in the
+ // number of degrees of freedom, this was not very practical, so we now
+ // use a vector of rows of the matrix, and in each row a vector of
+ // pairs (colnum,value). this seems like the best tradeoff between
+ // memory and speed, as it is now linear in memory and still fast
// enough.
//
- // to save some memory and since
- // the weights are usually
- // (negative) powers of 2, we
- // choose the value type of the
- // matrix to be @p{float} rather
- // than @p{double}.
+ // to save some memory and since the weights are usually (negative)
+ // powers of 2, we choose the value type of the matrix to be @p{float}
+ // rather than @p{double}.
std::vector<std::map<unsigned int, float> > weights;
- // this is this mapping. there is one
- // entry for each dof on the fine grid;
- // if it is a parameter dof, then its
- // value is the column in weights for
- // that parameter dof, if it is any
- // other dof, then its value is -1,
- // indicating an error
+ // this is this mapping. there is one entry for each dof on the fine
+ // grid; if it is a parameter dof, then its value is the column in
+ // weights for that parameter dof, if it is any other dof, then its
+ // value is -1, indicating an error
std::vector<int> weight_mapping;
internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
coarse_to_fine_grid_map,
weights, weight_mapping);
- // now compute the requested
- // representation
+ // now compute the requested representation
const unsigned int n_global_parm_dofs
= std::count_if (weight_mapping.begin(), weight_mapping.end(),
std::bind2nd (std::not_equal_to<int> (), -1));
- // first construct the inverse
- // mapping of weight_mapping
+ // first construct the inverse mapping of weight_mapping
std::vector<unsigned int> inverse_weight_mapping (n_global_parm_dofs,
DoFHandler<dim,spacedim>::invalid_dof_index);
for (unsigned int i=0; i<weight_mapping.size(); ++i)
{
const unsigned int parameter_dof = weight_mapping[i];
- // if this global dof is a
- // parameter
+ // if this global dof is a parameter
if (parameter_dof != numbers::invalid_unsigned_int)
{
Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
};
};
- // next copy over weights array
- // and replace respective
- // numbers
+ // next copy over weights array and replace respective numbers
const unsigned int n_rows = weight_mapping.size();
transfer_representation.clear ();
dofs_on_face.reserve (max_dofs_per_face(dof_handler));
unsigned int next_boundary_index = 0;
- // now loop over all cells and
- // check whether their faces are at
- // the boundary. note that we need
- // not take special care of single
- // lines being at the boundary
- // (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every isolated boundary
- // line is also part of a boundary
- // face which we will be visiting
- // sooner or later
+ // now loop over all cells and check whether their faces are at the
+ // boundary. note that we need not take special care of single lines
+ // being at the boundary (using @p{cell->has_boundary_lines}), since we
+ // do not support boundaries of dimension dim-2, and so every isolated
+ // boundary line is also part of a boundary face which we will be
+ // visiting sooner or later
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index)
{
- // check whether every fe in the collection
- // has support points
+ // check whether every fe in the collection has support points
Assert(fe_collection[fe_index].has_support_points(),
typename FiniteElement<dim>::ExcFEHasNoSupportPoints());
q_coll_dummy.push_back(
fe_collection[fe_index].get_unit_support_points()));
}
- // now loop over all cells and
- // enquire the support points on
- // each of these. we use dummy
- // quadrature formulas where the
- // quadrature points are located at
- // the unit support points to
- // enquire the location of the
- // support points in real space
+ // Now loop over all cells and enquire the support points on each
+ // of these. we use dummy quadrature formulas where the quadrature
+ // points are located at the unit support points to enquire the
+ // location of the support points in real space.
//
- // the weights of the quadrature
- // rule have been set to invalid values
- // by the used constructor.
+ // The weights of the quadrature rule have been set to invalid
+ // values by the used constructor.
hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe_collection,
q_coll_dummy, update_quadrature_points);
typename DH::active_cell_iterator cell =
ExcMessage ("This function can not be used with distributed triangulations."
"See the documentation for more information."));
- //Let the internal function do all the work,
- //just make sure that it gets a MappingCollection
+ // Let the internal function do all the work, just make sure that it
+ // gets a MappingCollection
const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
internal::map_dofs_to_support_points (mapping_collection,
ExcMessage ("This function can not be used with distributed triangulations."
"See the documentation for more information."));
- //Let the internal function do all the work,
- //just make sure that it gets a MappingCollection
+ // Let the internal function do all the work, just make sure that it
+ // gets a MappingCollection
internal::map_dofs_to_support_points (mapping,
dof_handler,
support_points);
{
support_points.clear();
- //Let the internal function do all the work,
- //just make sure that it gets a MappingCollection
+ // Let the internal function do all the work, just make sure that it
+ // gets a MappingCollection
const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
internal::map_dofs_to_support_points (mapping_collection,
{
support_points.clear();
- //Let the internal function do all the work,
- //just make sure that it gets a MappingCollection
+ // Let the internal function do all the work, just make sure that it
+ // gets a MappingCollection
internal::map_dofs_to_support_points (mapping,
dof_handler,
support_points);
||
(face->boundary_indicator() == boundary_indicator)))
{
- // get indices and physical
- // location on this face
+ // get indices and physical location on this face
face_dofs.resize (fe.dofs_per_face);
face->get_dof_indices (face_dofs, cell->active_fe_index());
- // enter those dofs into the list
- // that match the component
+ // enter those dofs into the list that match the component
// signature.
for (unsigned int i=0; i<face_dofs.size(); ++i)
{
- // Find out if a dof
- // has a contribution
- // in this component,
- // and if so, add it
- // to the list
+ // Find out if a dof has a contribution in this
+ // component, and if so, add it to the list
const ComponentMask &nonzero_component_array
= cell->get_fe().get_nonzero_components (i);
bool nonzero = false;
if (interior_only)
{
- // Exclude degrees of
- // freedom on faces
- // opposite to the vertex
+ // Exclude degrees of freedom on faces opposite to the vertex
exclude.resize(fe.dofs_per_cell);
std::fill(exclude.begin(), exclude.end(), false);
const unsigned int dpf = fe.dofs_per_face;
for (unsigned int child=0; child<pcell->n_children(); ++child)
{
cell = pcell->child(child);
- // For hp, only this line
- // here would have to be
- // replaced.
+
+ // For hp, only this line here would have to be replaced.
const FiniteElement<DH::dimension> &fe = dof_handler.get_fe();
const unsigned int n_dofs = fe.dofs_per_cell;
indices.resize(n_dofs);
if (interior_dofs_only)
{
- // Eliminate dofs on
- // faces of the child
- // which are on faces
+ // Eliminate dofs on faces of the child which are on faces
// of the parent
const unsigned int dpf = fe.dofs_per_face;
for (unsigned int i=0; i<dpf; ++i)
exclude[fe.face_to_cell_index(i,face)] = true;
}
- // Now remove all
- // degrees of freedom
- // on the domain
- // boundary from the
- // exclusion list
+
+ // Now remove all degrees of freedom on the domain boundary
+ // from the exclusion list
if (boundary_dofs)
for (unsigned int face=0; face< GeometryInfo<DH::dimension>::faces_per_cell; ++face)
if (cell->at_boundary(face))
typename DH::cell_iterator cell;
typename DH::cell_iterator endc = dof_handler.end(level);
- // Vector mapping from vertex
- // index in the triangulation to
- // consecutive block indices on
- // this level
- // The number of cells at a vertex
+ // Vector mapping from vertex index in the triangulation to consecutive
+ // block indices on this level The number of cells at a vertex
std::vector<unsigned int> vertex_cell_count(dof_handler.get_tria().n_vertices(), 0);
+
// Is a vertex at the boundary?
std::vector<bool> vertex_boundary(dof_handler.get_tria().n_vertices(), false);
+
std::vector<unsigned int> vertex_mapping(dof_handler.get_tria().n_vertices(),
numbers::invalid_unsigned_int);
- // Estimate for the number of
- // dofs at this point
+
+ // Estimate for the number of dofs at this point
std::vector<unsigned int> vertex_dof_count(dof_handler.get_tria().n_vertices(), 0);
- // Identify all vertices active
- // on this level and remember
- // some data about them
+ // Identify all vertices active on this level and remember some data
+ // about them
for (cell=dof_handler.begin(level); cell != endc; ++cell)
for (unsigned int v=0; v<GeometryInfo<DH::dimension>::vertices_per_cell; ++v)
{
vertex_boundary[vg] = true;
}
}
- // From now on, only vertices
- // with positive dof count are
- // "in".
+ // From now on, only vertices with positive dof count are "in".
- // Remove vertices at boundaries
- // or in corners
+ // Remove vertices at boundaries or in corners
for (unsigned int vg=0; vg<vertex_dof_count.size(); ++vg)
if ((!single_cell_patches && vertex_cell_count[vg] < 2)
||
(!boundary_patches && vertex_boundary[vg]))
vertex_dof_count[vg] = 0;
- // Create a mapping from all
- // vertices to the ones used here
+ // Create a mapping from all vertices to the ones used here
unsigned int i=0;
for (unsigned int vg=0; vg<vertex_mapping.size(); ++vg)
if (vertex_dof_count[vg] != 0)
if (vertex_dof_count[vg] != 0)
vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
- // Now that we have all the data,
- // we reduce it to the part we
- // actually want
+ // Now that we have all the data, we reduce it to the part we actually
+ // want
vertex_dof_count.resize(i);
- // At this point, the list of
- // patches is ready. Now we enter
- // the dofs into the sparsity
- // pattern.
+ // At this point, the list of patches is ready. Now we enter the dofs
+ // into the sparsity pattern.
block_list.reinit(vertex_dof_count.size(), dof_handler.n_dofs(level), vertex_dof_count);
std::vector<unsigned int> indices;
if (interior_only)
{
- // Exclude degrees of
- // freedom on faces
- // opposite to the vertex
+ // Exclude degrees of freedom on faces opposite to the
+ // vertex
exclude.resize(fe.dofs_per_cell);
std::fill(exclude.begin(), exclude.end(), false);
const unsigned int dpf = fe.dofs_per_face;
projects / dealii-svn.git / commitdiff