std::vector<bool> &selected_dofs,
const bool blocks = false);
- /**
- * Extract all degrees of freedom
- * which are at the boundary and
- * belong to specified components
- * of the solution. The function
- * returns its results in the
- * last non-default-valued
- * parameter which contains
- * @p true if a degree of
- * freedom is at the boundary and
- * belongs to one of the selected
- * components, and @p false
- * otherwise.
- *
- * By specifying the
- * @p boundary_indicator
- * variable, you can select which
- * boundary indicators the faces
- * have to have on which the
- * degrees of freedom are located
- * that shall be extracted. If it
- * is an empty list, then all
- * boundary indicators are
- * accepted.
- *
- * The size of @p component_select
- * shall equal the number of
- * components in the finite
- * element used by @p dof. The
- * size of @p selected_dofs shall
- * equal
- * <tt>dof_handler.n_dofs()</tt>. Previous
- * contents of this array or
- * overwritten.
- *
- * Using the usual convention, if
- * a shape function is non-zero
- * in more than one component
- * (i.e. it is non-primitive),
- * then the element in the
- * component mask is used that
- * corresponds to the first
- * non-zero components. Elements
- * in the mask corresponding to
- * later components are ignored.
- *
- * @note This function will not work
- * for DoFHandler objects that are
- * built on a
- * parallel::distributed::Triangulation
- * object. The reasons is that the
- * output argument @p selected_dofs
- * has to have a length equal to <i>all</i>
- * global degrees of freedom.
- * Consequently, this does not scale
- * to very large problems. If you
- * need the functionality of this
- * function for parallel triangulations,
- * then you need to use the other
- * DoFTools::extract_boundary_dofs
- * function.
- */
+ /**
+ * Extract all degrees of freedom
+ * which are at the boundary and
+ * belong to specified components
+ * of the solution. The function
+ * returns its results in the
+ * last non-default-valued
+ * parameter which contains
+ * @p true if a degree of
+ * freedom is at the boundary and
+ * belongs to one of the selected
+ * components, and @p false
+ * otherwise.
+ *
+ * By specifying the
+ * @p boundary_indicator
+ * variable, you can select which
+ * boundary indicators the faces
+ * have to have on which the
+ * degrees of freedom are located
+ * that shall be extracted. If it
+ * is an empty list, then all
+ * boundary indicators are
+ * accepted.
+ *
+ * The size of @p component_select
+ * shall equal the number of
+ * components in the finite
+ * element used by @p dof. The
+ * size of @p selected_dofs shall
+ * equal
+ * <tt>dof_handler.n_dofs()</tt>. Previous
+ * contents of this array or
+ * overwritten.
+ *
+ * Using the usual convention, if
+ * a shape function is non-zero
+ * in more than one component
+ * (i.e. it is non-primitive),
+ * then the element in the
+ * component mask is used that
+ * corresponds to the first
+ * non-zero components. Elements
+ * in the mask corresponding to
+ * later components are ignored.
+ *
+ * @note This function will not work
+ * for DoFHandler objects that are
+ * built on a
+ * parallel::distributed::Triangulation
+ * object. The reasons is that the
+ * output argument @p selected_dofs
+ * has to have a length equal to <i>all</i>
+ * global degrees of freedom.
+ * Consequently, this does not scale
+ * to very large problems. If you
+ * need the functionality of this
+ * function for parallel triangulations,
+ * then you need to use the other
+ * DoFTools::extract_boundary_dofs
+ * function.
+ *
+ * @param dof_handler The object that describes which degrees of freedom
+ * live on which cell
+ * @param component_select A mask denoting the vector components of the
+ * finite element that should be considered.
+ * @param selected_dofs The IndexSet object that is returned and that
+ * will contain the indices of degrees of freedom that are
+ * located on the boundary (and correspond to the selected
+ * vector components and boundary indicators, depending on
+ * the values of the @p component_select and @p boundary_indicators
+ * arguments).
+ * @param boundary_indicators If empty, this function extracts the
+ * indices of the degrees of freedom for all parts of the boundary.
+ * If it is a non-empty list, then the function only considers
+ * boundary faces with the boundary indicators listed in this
+ * argument.
+ */
template <class DH>
void
extract_boundary_dofs (const DH &dof_handler,
std::vector<bool> &selected_dofs,
const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
+ /**
+ * This function does the same as the previous one but it
+ * returns its result as an IndexSet rather than a std::vector@<bool@>.
+ * Thus, it can also be called for DoFHandler objects that are
+ * defined on parallel::distributed::Triangulation objects.
+ *
+ * @note If the DoFHandler object is indeed defined on a
+ * parallel::distributed::Triangulation, then the @p selected_dofs
+ * index set will contain only those degrees of freedom on the
+ * boundary that belong to the locally relevant set (see
+ * @ref GlossLocallyRelevantDof "locally relevant DoFs").
+ *
+ * @param dof_handler The object that describes which degrees of freedom
+ * live on which cell
+ * @param component_select A mask denoting the vector components of the
+ * finite element that should be considered.
+ * @param selected_dofs The IndexSet object that is returned and that
+ * will contain the indices of degrees of freedom that are
+ * located on the boundary (and correspond to the selected
+ * vector components and boundary indicators, depending on
+ * the values of the @p component_select and @p boundary_indicators
+ * arguments).
+ * @param boundary_indicators If empty, this function extracts the
+ * indices of the degrees of freedom for all parts of the boundary.
+ * If it is a non-empty list, then the function only considers
+ * boundary faces with the boundary indicators listed in this
+ * argument.
+ */
+ template <class DH>
+ void
+ extract_boundary_dofs (const DH &dof_handler,
+ const std::vector<bool> &component_select,
+ IndexSet &selected_dofs,
+ const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
+
/**
* This function is similar to
* the extract_boundary_dofs()
std::vector<bool> &selected_dofs,
const std::set<types::boundary_id_t> &boundary_indicators)
{
- AssertDimension (component_select.size(), n_components(dof_handler));
- Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
- ExcInvalidBoundaryIndicator());
Assert ((dynamic_cast<const parallel::distributed::Triangulation<DH::dimension,DH::space_dimension>*>
(&dof_handler.get_tria())
==
ExcMessage ("This function can not be used with distributed triangulations."
"See the documentation for more information."));
+ IndexSet indices;
+ extract_boundary_dofs (dof_handler, component_select,
+ indices, boundary_indicators);
+
+ // clear and reset array by default values
+ selected_dofs.clear ();
+ selected_dofs.resize (dof_handler.n_dofs(), false);
+
+ // then convert the values computed above to the binary vector
+ indices.fill_binary_vector(selected_dofs);
+ }
+
+
+ template <class DH>
+ void
+ extract_boundary_dofs (const DH &dof_handler,
+ const std::vector<bool> &component_select,
+ IndexSet &selected_dofs,
+ const std::set<types::boundary_id_t> &boundary_indicators)
+ {
+ AssertDimension (component_select.size(), n_components(dof_handler));
+ Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
+ ExcInvalidBoundaryIndicator());
const unsigned int dim=DH::dimension;
+ // first reset output argument
+ 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
= (component_select != std::vector<bool>(component_select.size(),
true));
- // clear and reset array by default
- // values
- selected_dofs.clear ();
- selected_dofs.resize (dof_handler.n_dofs(), false);
std::vector<unsigned int> face_dof_indices;
face_dof_indices.reserve (max_dofs_per_face(dof_handler));
// sooner or later
for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell)
- for (unsigned int face=0;
- face<GeometryInfo<DH::dimension>::faces_per_cell; ++face)
- if (cell->at_boundary(face))
- if (! check_boundary_indicator ||
- (boundary_indicators.find (cell->face(face)->boundary_indicator())
- != boundary_indicators.end()))
- {
- const FiniteElement<DH::dimension, DH::space_dimension> &fe = cell->get_fe();
+ // 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 (cell->at_boundary(face))
+ if (! check_boundary_indicator ||
+ (boundary_indicators.find (cell->face(face)->boundary_indicator())
+ != boundary_indicators.end()))
+ {
+ const FiniteElement<DH::dimension, DH::space_dimension> &fe = cell->get_fe();
- const unsigned int dofs_per_face = fe.dofs_per_face;
- face_dof_indices.resize (dofs_per_face);
- cell->face(face)->get_dof_indices (face_dof_indices,
- cell->active_fe_index());
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+ face_dof_indices.resize (dofs_per_face);
+ cell->face(face)->get_dof_indices (face_dof_indices,
+ cell->active_fe_index());
- for (unsigned int i=0; i<fe.dofs_per_face; ++i)
- if (!check_vector_component)
- selected_dofs[face_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)
- {
- // first get at the
- // cell-global
- // number of a face
- // dof, to ask the
- // fe certain
- // questions
- const unsigned int cell_index
+ for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+ 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)
+ {
+ // 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
- :
- (dim == 2 ?
- (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ i
:
- (dim == 3 ?
- (i<4*fe.dofs_per_vertex ?
- i
- :
- (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
- i+4*fe.dofs_per_vertex
- :
- i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
- :
- numbers::invalid_unsigned_int)));
- if (fe.is_primitive (cell_index))
- selected_dofs[face_dof_indices[i]]
- = (component_select[fe.face_system_to_component_index(i).first]
- == true);
- else // not primitive
- {
- const unsigned int first_nonzero_comp
+ (dim == 2 ?
+ (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ :
+ (dim == 3 ?
+ (i<4*fe.dofs_per_vertex ?
+ i
+ :
+ (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+ i+4*fe.dofs_per_vertex
+ :
+ i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ :
+ numbers::invalid_unsigned_int)));
+ if (fe.is_primitive (cell_index))
+ {
+ if (component_select[fe.face_system_to_component_index(i).first]
+ == true)
+ selected_dofs.add_index (face_dof_indices[i]);
+ }
+ else // not primitive
+ {
+ const unsigned int first_nonzero_comp
= (std::find (fe.get_nonzero_components(cell_index).begin(),
fe.get_nonzero_components(cell_index).end(),
true)
- -
- fe.get_nonzero_components(cell_index).begin());
- Assert (first_nonzero_comp < fe.n_components(),
- ExcInternalError());
-
- selected_dofs[face_dof_indices[i]]
- = (component_select[first_nonzero_comp]
- == true);
- }
- }
- }
+ -
+ fe.get_nonzero_components(cell_index).begin());
+ Assert (first_nonzero_comp < fe.n_components(),
+ ExcInternalError());
+
+ if (component_select[first_nonzero_comp] == true)
+ selected_dofs.add_index (face_dof_indices[i]);
+ }
+ }
+ }
}