face_1->get_dof_indices(dofs_1, face_1_index);
face_2->get_dof_indices(dofs_2, face_2_index);
- // .. and constrain them (respecting
- // component_mask):
+ // .. and constrain them (respecting
+ // component_mask):
for (unsigned int i = 0; i < dofs_per_face; ++i) {
- if (component_mask.size() == 0 ||
- component_mask[fe.face_system_to_component_index(i).first]) {
- if(!constraint_matrix.is_constrained(dofs_1[i])) {
- constraint_matrix.add_line(dofs_1[i]);
- constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0);
+ if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components())
+ ||
+ (component_mask[fe.face_system_to_component_index(i).first] == true))
+ {
+ if(!constraint_matrix.is_constrained(dofs_1[i]))
+ {
+ constraint_matrix.add_line(dofs_1[i]);
+ constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0);
+ }
}
- }
}
}
const types::boundary_id boundary_component,
const int direction,
dealii::Tensor<1,DH::space_dimension>
- &offset,
+ &offset,
dealii::ConstraintMatrix &constraint_matrix,
- const std::vector<bool> &component_mask)
+ const ComponentMask &component_mask)
{
static const int space_dim = DH::space_dimension;
Assert (0<=direction && direction<space_dim,
- template<int dim, int spacedim>
+ template <int dim, int spacedim>
+ void
+ extract_dofs (const DoFHandler<dim,spacedim> &dof,
+ const BlockMask &block_mask,
+ std::vector<bool> &selected_dofs)
+ {
+ // simply forward to the function that works based on a component mask
+ extract_dofs (dof, dof.get_fe().component_mask (block_mask),
+ selected_dofs);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ extract_dofs (const hp::DoFHandler<dim,spacedim> &dof,
+ const BlockMask &block_mask,
+ std::vector<bool> &selected_dofs)
+ {
+ // simply forward to the function that works based on a component mask
+ extract_dofs (dof, dof.get_fe().component_mask (block_mask),
+ selected_dofs);
+ }
+
+
+
+ template<class DH>
void
- extract_level_dofs(
- const unsigned int level,
- const DH &dof,
- const std::vector<bool> &component_select,
- std::vector<bool> &selected_dofs,
- const bool count_by_blocks)
+ extract_level_dofs(const unsigned int level,
- const MGDoFHandler<dim,spacedim> &dof,
++ const DH &dof,
+ const ComponentMask &component_mask,
+ std::vector<bool> &selected_dofs)
{
- const FiniteElement<dim,spacedim>& fe = dof.get_fe();
+ const FiniteElement<DH::dimension,DH::space_dimension>& fe = dof.get_fe();
- if (count_by_blocks == true)
- {
- Assert(component_select.size() == fe.n_blocks(),
- ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
- }
- else
- {
- Assert(component_select.size() == fe.n_components(),
- ExcDimensionMismatch(component_select.size(), fe.n_components()));
- }
-
+ Assert(component_mask.represents_n_components(n_components(dof)),
+ ExcMessage ("The given component mask is not sized correctly to represent the "
+ "components of the given finite element."));
Assert(selected_dofs.size() == dof.n_dofs(level),
ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
// 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)
+ typename DH::cell_iterator c;
- for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
++ for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
{
c->get_mg_dof_indices(indices);
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- template<int dim, int spacedim>
++ template<class DH>
+ void
+ extract_level_dofs(const unsigned int level,
- const MGDoFHandler<dim,spacedim> &dof,
++ const DH &dof,
+ const BlockMask &block_mask,
+ std::vector<bool> &selected_dofs)
+ {
+ // simply defer to the other extract_level_dofs() function
+ extract_level_dofs (level, dof, dof.get_fe().component_mask(block_mask),
+ selected_dofs);
+ }
+
+
+
template <class DH>
void
extract_boundary_dofs (const DH &dof_handler,
"See the documentation for more information."));
IndexSet indices;
- extract_boundary_dofs (dof_handler, component_select,
+ extract_boundary_dofs (dof_handler, component_mask,
indices, boundary_indicators);
- // 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);
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
- :
- (dim == 2 ?
- (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ = (dim == 1 ?
+ i
+ :
+ (dim == 2 ?
+ (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ :
+ (dim == 3 ?
+ (i<4*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)));
+ (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]
+ if (component_mask[fe.face_system_to_component_index(i).first]
- == true)
+ == true)
selected_dofs.add_index (face_dof_indices[i]);
}
else // not primitive
template void DoFTools::extract_dofs<deal_II_dimension>
(const hp::DoFHandler<deal_II_dimension>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+ const ComponentMask &,
+ std::vector<bool>&);
+
+ template void DoFTools::extract_dofs<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension>&,
+ const BlockMask &,
+ std::vector<bool>&);
+
-template void DoFTools::extract_level_dofs<deal_II_dimension>
++template void DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension> >
++(const unsigned int level,
++ const DoFHandler<deal_II_dimension>&,
++ const ComponentMask &,
++ std::vector<bool>&);
+
+template void DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension> >
- (const unsigned int, const DoFHandler<deal_II_dimension>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
++(const unsigned int level,
++ const DoFHandler<deal_II_dimension>&,
++ const BlockMask &,
++ std::vector<bool>&);
+
+template void DoFTools::extract_level_dofs<MGDoFHandler<deal_II_dimension> >
- (const unsigned int, const MGDoFHandler<deal_II_dimension>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
+ const MGDoFHandler<deal_II_dimension>&,
+ const ComponentMask &,
+ std::vector<bool>&);
+
-template void DoFTools::extract_level_dofs<deal_II_dimension>
++template void DoFTools::extract_level_dofs<MGDoFHandler<deal_II_dimension> >
+ (const unsigned int level,
+ const MGDoFHandler<deal_II_dimension>&,
+ const BlockMask &,
+ std::vector<bool>&);
template
void
template
void
-DoFTools::extract_level_dofs<deal_II_dimension, deal_II_dimension+1>
+DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension, deal_II_dimension+1> >
- (const unsigned int, const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
- const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
++ const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
+ const ComponentMask&,
+ std::vector<bool>&);
template
void
-DoFTools::extract_level_dofs<deal_II_dimension, deal_II_dimension+1>
+DoFTools::extract_level_dofs<MGDoFHandler<deal_II_dimension, deal_II_dimension+1> >
- (const unsigned int, const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
+ const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
+ const BlockMask&,
+ std::vector<bool>&);
#endif
template
void
-DoFTools::extract_level_dofs<1,3>
+DoFTools::extract_level_dofs<DoFHandler<1, 3> >
- (const unsigned int, const DoFHandler<1, 3>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
- const MGDoFHandler<1,3>&,
++ const DoFHandler<1,3>&,
+ const ComponentMask &,
+ std::vector<bool>&);
template
void
-DoFTools::extract_level_dofs<1,3>
+DoFTools::extract_level_dofs<MGDoFHandler<1, 3> >
- (const unsigned int, const MGDoFHandler<1, 3>&,
- const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
+ const MGDoFHandler<1,3>&,
+ const BlockMask &,
+ std::vector<bool>&);
#endif
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/base/mg_level_object.h>
+#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe.h>
+ #include <deal.II/fe/component_mask.h>
#include <vector>
#include <algorithm>
// otherwise determine the number
// of dofs in each component
// separately. do so in parallel
- std::vector<std::vector<bool> >
- dofs_in_component (n_components,
- std::vector<bool>(dof_handler.n_dofs(l),
- false));
- std::vector<std::vector<bool> >
- component_select (n_components,
- std::vector<bool>(n_components, false));
- Threads::TaskGroup<> tasks;
- for (unsigned int i=0; i<n_components; ++i)
- {
- void (*fun_ptr) (const unsigned int level,
- const MGDoFHandler<dim,spacedim> &,
- const std::vector<bool> &,
- std::vector<bool> &,
- bool)
- = &DoFTools::template extract_level_dofs<MGDoFHandler<dim, spacedim> >;
- component_select[i][i] = true;
- tasks += Threads::new_task (fun_ptr,
- l, dof_handler,
- component_select[i],
- dofs_in_component[i], false);
- }
- tasks.join_all();
+ std::vector<std::vector<bool> >
+ dofs_in_component (n_components,
+ std::vector<bool>(dof_handler.n_dofs(l),
+ false));
+ std::vector<ComponentMask> component_select (n_components);
+ Threads::TaskGroup<> tasks;
+ for (unsigned int i=0; i<n_components; ++i)
+ {
+ void (*fun_ptr) (const unsigned int level,
+ const MGDoFHandler<dim,spacedim> &,
+ const ComponentMask &,
+ std::vector<bool> &)
- = &DoFTools::template extract_level_dofs<dim>;
++ = &DoFTools::template extract_level_dofs<MGDoFHandler<dim,spacedim> >;
+
+ std::vector<bool> tmp(n_components, false);
+ tmp[i] = true;
+ component_select[i] = ComponentMask(tmp);
+
+ tasks += Threads::new_task (fun_ptr,
+ l, dof_handler,
+ component_select[i],
+ dofs_in_component[i]);
+ }
+ tasks.join_all();
// next count what we got
- unsigned int component = 0;
- for (unsigned int b=0;b<fe.n_base_elements();++b)
- {
- const FiniteElement<dim>& base = fe.base_element(b);
+ unsigned int component = 0;
+ for (unsigned int b=0;b<fe.n_base_elements();++b)
+ {
+ const FiniteElement<dim>& base = fe.base_element(b);
// Dimension of base element
- unsigned int d = base.n_components();
-
- for (unsigned int m=0;m<fe.element_multiplicity(b);++m)
- {
- for (unsigned int dd=0;dd<d;++dd)
- {
- if (base.is_primitive() || (!only_once || dd==0))
- result[l][target_component[component]]
- += std::count(dofs_in_component[component].begin(),
- dofs_in_component[component].end(),
- true);
- ++component;
- }
- }
- }
+ unsigned int d = base.n_components();
+
+ for (unsigned int m=0;m<fe.element_multiplicity(b);++m)
+ {
+ for (unsigned int dd=0;dd<d;++dd)
+ {
+ if (base.is_primitive() || (!only_once || dd==0))
+ result[l][target_component[component]]
+ += std::count(dofs_in_component[component].begin(),
+ dofs_in_component[component].end(),
+ true);
+ ++component;
+ }
+ }
+ }
// finally sanity check
- Assert (!dof_handler.get_fe().is_primitive()
- ||
- std::accumulate (result[l].begin(),
- result[l].end(), 0U)
- ==
- dof_handler.n_dofs(l),
- ExcInternalError());
- }
+ Assert (!dof_handler.get_fe().is_primitive()
+ ||
+ std::accumulate (result[l].begin(),
+ result[l].end(), 0U)
+ ==
+ dof_handler.n_dofs(l),
+ ExcInternalError());
+ }
}
}
for (unsigned int l=0;l<n_levels;++l)
{
std::vector<std::vector<bool> >
- dofs_in_block (n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
+ dofs_in_block (n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
- std::vector<std::vector<bool> >
- block_select (n_blocks, std::vector<bool>(n_blocks, false));
+ std::vector<BlockMask> block_select (n_blocks);
Threads::TaskGroup<> tasks;
for (unsigned int i=0; i<n_blocks; ++i)
- {
- void (*fun_ptr) (const unsigned int level,
- const DH&,
- const std::vector<bool>&,
- std::vector<bool>&,
- bool)
- = &DoFTools::template extract_level_dofs<DH>;
- block_select[i][i] = true;
- tasks += Threads::new_task (fun_ptr,
- l, dof_handler, block_select[i],
- dofs_in_block[i], true);
- };
+ {
+ void (*fun_ptr) (const unsigned int level,
- const MGDoFHandler<dim,spacedim>&,
++ const DH&,
+ const BlockMask &,
+ std::vector<bool>&)
- = &DoFTools::template extract_level_dofs<dim>;
++ = &DoFTools::template extract_level_dofs<DH>;
+
+ std::vector<bool> tmp(n_blocks, false);
+ tmp[i] = true;
+ block_select[i] = tmp;
+
+ tasks += Threads::new_task (fun_ptr,
+ l, dof_handler, block_select[i],
+ dofs_in_block[i]);
+ }
tasks.join_all ();
// next count what we got
std::vector<unsigned int> local_dofs;
local_dofs.reserve (DoFTools::max_dofs_per_face(dof));
- std::fill (local_dofs.begin (), local_dofs.end (),
+ std::fill (local_dofs.begin (),
+ local_dofs.end (),
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ DoFHandler<dim,spacedim>::invalid_dof_index);
// First, deal with the simpler
// case when we have to identify
// all boundary dofs
- if (component_mask.size() == 0)
+ if (component_mask.n_selected_components(n_components) == n_components)
{
typename MGDoFHandler<dim,spacedim>::cell_iterator
- cell = dof.begin(),
- endc = dof.end();
+ cell = dof.begin(),
+ endc = dof.end();
for (; cell!=endc; ++cell)
- {
- const FiniteElement<dim> &fe = cell->get_fe();
- const unsigned int level = cell->level();
- local_dofs.resize(fe.dofs_per_face);
-
- for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- if (cell->at_boundary(face_no) == true)
- {
- const typename MGDoFHandler<dim,spacedim>::face_iterator
- face = cell->face(face_no);
- const types::boundary_id bi = face->boundary_indicator();
+ {
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int level = cell->level();
+ local_dofs.resize(fe.dofs_per_face);
+
+ for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ if (cell->at_boundary(face_no) == true)
+ {
+ const typename MGDoFHandler<dim,spacedim>::face_iterator
+ face = cell->face(face_no);
+ const types::boundary_id bi = face->boundary_indicator();
// Face is listed in
// boundary map
- if (function_map.find(bi) != function_map.end())
- {
- face->get_mg_dof_indices(level, local_dofs);
- for (unsigned int i=0;i<fe.dofs_per_face;++i)
- boundary_indices[level].insert(local_dofs[i]);
- }
- }
- }
+ if (function_map.find(bi) != function_map.end())
+ {
+ face->get_mg_dof_indices(level, local_dofs);
+ for (unsigned int i=0;i<fe.dofs_per_face;++i)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
+ }
}
else
{
- Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
+ Assert (component_mask.n_selected_components(n_components) > 0,
- ExcMessage("It's probably worthwhile to select at least one component."));
+ ExcMessage("It's probably worthwhile to select at least one component."));
typename MGDoFHandler<dim,spacedim>::cell_iterator
- cell = dof.begin(),
- endc = dof.end();
+ cell = dof.begin(),
+ endc = dof.end();
for (; cell!=endc; ++cell)
- for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- {
- if (cell->at_boundary(face_no) == false)
- continue;
+ for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ if (cell->at_boundary(face_no) == false)
+ continue;
- const FiniteElement<dim> &fe = cell->get_fe();
- const unsigned int level = cell->level();
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int level = cell->level();
// we can presently deal only with
// primitive elements for boundary
// that are non-zero for the
// components we are interested in,
// are in fact primitive
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- {
- const std::vector<bool> &nonzero_component_array
- = cell->get_fe().get_nonzero_components (i);
- for (unsigned int c=0; c<n_components; ++c)
- if ((nonzero_component_array[c] == true)
- &&
- (component_mask[c] == true))
- Assert (cell->get_fe().is_primitive (i),
- ExcMessage ("This function can only deal with requested boundary "
- "values that correspond to primitive (scalar) base "
- "elements"));
- }
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ {
+ const ComponentMask &nonzero_component_array
+ = cell->get_fe().get_nonzero_components (i);
+ for (unsigned int c=0; c<n_components; ++c)
+ if ((nonzero_component_array[c] == true)
+ &&
+ (component_mask[c] == true))
+ Assert (cell->get_fe().is_primitive (i),
+ ExcMessage ("This function can only deal with requested boundary "
+ "values that correspond to primitive (scalar) base "
+ "elements"));
+ }
- typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
- const types::boundary_id boundary_component = face->boundary_indicator();
- if (function_map.find(boundary_component) != function_map.end())
+ typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
+ const types::boundary_id boundary_component = face->boundary_indicator();
+ if (function_map.find(boundary_component) != function_map.end())
// face is of the right component
- {
+ {
// get indices, physical location and
// boundary values of dofs on this
// face
// components. if shape function is non-primitive, then we will ignore
// the result in the following anyway, otherwise there's only one
// non-zero component which we will use
- component = (std::find (fe.get_nonzero_components(cell_i).begin(),
- fe.get_nonzero_components(cell_i).end(),
- true)
- -
- fe.get_nonzero_components(cell_i).begin());
- }
+ component = fe.get_nonzero_components(cell_i).first_selected_component();
+ }
- if (component_mask[component] == true)
- boundary_indices[level].insert(local_dofs[i]);
- }
- }
- else
- for (unsigned int i=0; i<local_dofs.size(); ++i)
- boundary_indices[level].insert(local_dofs[i]);
- }
- }
+ if (component_mask[component] == true)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
+ else
+ for (unsigned int i=0; i<local_dofs.size(); ++i)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
}
}
template <int dim, int spacedim>
void
+ make_boundary_list(
+ const DoFHandler<dim,spacedim>& dof,
+ const typename FunctionMap<dim>::type& function_map,
+ std::vector<std::set<unsigned int> >& boundary_indices,
+ const std::vector<bool>& component_mask)
+ {
+ // if for whatever reason we were
+ // passed an empty map, return
+ // immediately
+ if (function_map.size() == 0)
+ return;
+
+ const unsigned int n_levels = dof.get_tria().n_levels();
+
+
- const std::vector<bool> &nonzero_component_array
- = cell->get_fe().get_nonzero_components (i);
+ const unsigned int n_components = DoFTools::n_components(dof);
+ const bool fe_is_system = (n_components != 1);
+
+ AssertDimension (boundary_indices.size(), n_levels);
+
+ std::vector<unsigned int> local_dofs;
+ local_dofs.reserve (DoFTools::max_dofs_per_face(dof));
+ std::fill (local_dofs.begin (), local_dofs.end (),
+ DoFHandler<dim,spacedim>::invalid_dof_index);
+
+ // First, deal with the simpler
+ // case when we have to identify
+ // all boundary dofs
+ if (component_mask.size() == 0)
+ {
+ typename DoFHandler<dim,spacedim>::cell_iterator
+ cell = dof.begin(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ {
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int level = cell->level();
+ local_dofs.resize(fe.dofs_per_face);
+
+ for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ if (cell->at_boundary(face_no))
+ {
+ const typename DoFHandler<dim,spacedim>::face_iterator
+ face = cell->face(face_no);
+ const unsigned char bi = face->boundary_indicator();
+ // Face is listed in
+ // boundary map
+ if (function_map.find(bi) != function_map.end())
+ {
+ face->get_mg_dof_indices(level, local_dofs);
+ for (unsigned int i=0;i<fe.dofs_per_face;++i)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
+ }
+ }
+ else
+ {
+ Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
+ ExcMessage("It's probably worthwhile to select at least one component."));
+
+ typename DoFHandler<dim,spacedim>::cell_iterator
+ cell = dof.begin(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ if (!(cell->at_boundary(face_no)))
+ continue;
+
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int level = cell->level();
+
+ // we can presently deal only with
+ // primitive elements for boundary
+ // values. this does not preclude
+ // us using non-primitive elements
+ // in components that we aren't
+ // interested in, however. make
+ // sure that all shape functions
+ // that are non-zero for the
+ // components we are interested in,
+ // are in fact primitive
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ {
- component = (std::find (fe.get_nonzero_components(cell_i).begin(),
- fe.get_nonzero_components(cell_i).end(),
- true)
- -
- fe.get_nonzero_components(cell_i).begin());
++ const ComponentMask &nonzero_component_array
++ = cell->get_fe().get_nonzero_components (i);
+ for (unsigned int c=0; c<n_components; ++c)
+ if ((nonzero_component_array[c] == true)
+ &&
+ (component_mask[c] == true))
+ Assert (cell->get_fe().is_primitive (i),
+ ExcMessage ("This function can only deal with requested boundary "
+ "values that correspond to primitive (scalar) base "
+ "elements"));
+ }
+
+ typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
+ const unsigned char boundary_component = face->boundary_indicator();
+ if (function_map.find(boundary_component) != function_map.end())
+ // face is of the right component
+ {
+ // get indices, physical location and
+ // boundary values of dofs on this
+ // face
+ local_dofs.resize (fe.dofs_per_face);
+ face->get_mg_dof_indices (level, local_dofs);
+ if (fe_is_system)
+ {
+ // enter those dofs
+ // into the list that
+ // match the
+ // component
+ // signature. avoid
+ // the usual
+ // complication that
+ // we can't just use
+ // *_system_to_component_index
+ // for non-primitive
+ // FEs
+ for (unsigned int i=0; i<local_dofs.size(); ++i)
+ {
+ unsigned int component;
+ if (fe.is_primitive())
+ component = fe.face_system_to_component_index(i).first;
+ else
+ {
+ // non-primitive
+ // case. make
+ // sure that
+ // this
+ // particular
+ // shape
+ // function
+ // _is_
+ // primitive,
+ // and get at
+ // it's
+ // component. use
+ // usual
+ // trick to
+ // transfer
+ // face dof
+ // index to
+ // cell dof
+
+ // index
+ const unsigned int cell_i
+ = (dim == 1 ?
+ i
+ :
+ (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)));
+ Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
+
+ // make sure
+ // that if
+ // this is
+ // not a
+ // primitive
+ // shape function,
+ // then all
+ // the
+ // corresponding
+ // components
+ // in the
+ // mask are
+ // not set
+// if (!fe.is_primitive(cell_i))
+// for (unsigned int c=0; c<n_components; ++c)
+// if (fe.get_nonzero_components(cell_i)[c])
+// Assert (component_mask[c] == false,
+// ExcFENotPrimitive());
+
+// let's pick the first of possibly more than one non-zero
+// components. if shape function is non-primitive, then we will ignore
+// the result in the following anyway, otherwise there's only one
+// non-zero component which we will use
++ component = fe.get_nonzero_components(cell_i).first_selected_component();
+ }
+
+ if (component_mask[component] == true)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
+ else
+ for (unsigned int i=0; i<local_dofs.size(); ++i)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
make_boundary_list(const MGDoFHandler<dim,spacedim>& dof,
- const typename FunctionMap<dim>::type& function_map,
- std::vector<IndexSet>& boundary_indices,
- const std::vector<bool>& component_mask)
+ const typename FunctionMap<dim>::type& function_map,
+ std::vector<IndexSet>& boundary_indices,
+ const ComponentMask & component_mask)
{
Assert (boundary_indices.size() == dof.get_tria().n_levels(),
- ExcDimensionMismatch (boundary_indices.size(),
- dof.get_tria().n_levels()));
+ ExcDimensionMismatch (boundary_indices.size(),
+ dof.get_tria().n_levels()));
std::vector<std::set<unsigned int> >
my_boundary_indices (dof.get_tria().n_levels());
projects / dealii-svn.git / commitdiff