||
(subdomain_id == dof.get_tria().locally_owned_subdomain()),
ExcMessage ("For parallel::distributed::Triangulation objects and "
- "associated DoF handler objects, asking for any subdomain other "
- "than the locally owned one does not make sense."));
+ "associated DoF handler objects, asking for any subdomain other "
+ "than the locally owned one does not make sense."));
- std::vector<unsigned int> dofs_on_this_cell;
+ std::vector<types::global_dof_index> dofs_on_this_cell;
dofs_on_this_cell.reserve (max_dofs_per_cell(dof));
typename DH::active_cell_iterator cell = dof.begin_active(),
- endc = dof.end();
+ endc = dof.end();
// In case we work with a distributed
// sparsity pattern of Trilinos type, we
if (need_dof_mask == true)
for (unsigned int f=0; f<fe_collection.size(); ++f)
- {
- const unsigned int dofs_per_cell = fe_collection[f].dofs_per_cell;
-
- dof_mask[f].reinit (dofs_per_cell, dofs_per_cell);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (fe_collection[f].is_primitive(i) &&
- fe_collection[f].is_primitive(j))
- dof_mask[f](i,j)
- = (couplings(fe_collection[f].system_to_component_index(i).first,
- fe_collection[f].system_to_component_index(j).first) != none);
- else
- {
- const unsigned int first_nonzero_comp_i
- = (std::find (fe_collection[f].get_nonzero_components(i).begin(),
- fe_collection[f].get_nonzero_components(i).end(),
- true)
- -
- fe_collection[f].get_nonzero_components(i).begin());
- const unsigned int first_nonzero_comp_j
- = (std::find (fe_collection[f].get_nonzero_components(j).begin(),
- fe_collection[f].get_nonzero_components(j).end(),
- true)
- -
- fe_collection[f].get_nonzero_components(j).begin());
- Assert (first_nonzero_comp_i < fe_collection[f].n_components(),
- ExcInternalError());
- Assert (first_nonzero_comp_j < fe_collection[f].n_components(),
- ExcInternalError());
-
- dof_mask[f](i,j)
- = (couplings(first_nonzero_comp_i,first_nonzero_comp_j) != none);
- }
- }
+ {
+ const unsigned int dofs_per_cell = fe_collection[f].dofs_per_cell;
+
+ dof_mask[f].reinit (dofs_per_cell, dofs_per_cell);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (fe_collection[f].is_primitive(i) &&
+ fe_collection[f].is_primitive(j))
+ dof_mask[f](i,j)
+ = (couplings(fe_collection[f].system_to_component_index(i).first,
+ fe_collection[f].system_to_component_index(j).first) != none);
+ else
+ {
+ const unsigned int first_nonzero_comp_i
+ = fe_collection[f].get_nonzero_components(i).first_selected_component();
+ const unsigned int first_nonzero_comp_j
+ = fe_collection[f].get_nonzero_components(j).first_selected_component();
+ Assert (first_nonzero_comp_i < fe_collection[f].n_components(),
+ ExcInternalError());
+ Assert (first_nonzero_comp_j < fe_collection[f].n_components(),
+ ExcInternalError());
+
+ dof_mask[f](i,j)
+ = (couplings(first_nonzero_comp_i,first_nonzero_comp_j) != none);
+ }
+ }
- std::vector<unsigned int> dofs_on_this_cell(fe_collection.max_dofs_per_cell());
+ std::vector<types::global_dof_index> dofs_on_this_cell(fe_collection.max_dofs_per_cell());
typename DH::active_cell_iterator cell = dof.begin_active(),
- endc = dof.end();
+ endc = dof.end();
// In case we work with a distributed
// sparsity pattern of Trilinos type, we
const std::list<std::pair<typename DH::cell_iterator,
-- typename DH::cell_iterator> >
++ typename DH::cell_iterator> >
cell_list
= GridTools::get_finest_common_cells (dof_row, dof_col);
for (; cell_iter!=cell_list.end(); ++cell_iter)
{
- const typename DH::cell_iterator cell_row = cell_iter->first;
- const typename DH::cell_iterator cell_col = cell_iter->second;
-
- if (!cell_row->has_children() && !cell_col->has_children())
+ const typename DH::cell_iterator cell_row = cell_iter->first;
+ const typename DH::cell_iterator cell_col = cell_iter->second;
+
- if (!cell_row->has_children() && !cell_col->has_children())
- {
- const unsigned int dofs_per_cell_row =
- cell_row->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col->get_fe().dofs_per_cell;
- std::vector<unsigned int>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<unsigned int>
- local_dof_indices_col(dofs_per_cell_col);
- cell_row->get_dof_indices (local_dof_indices_row);
- cell_col->get_dof_indices (local_dof_indices_col);
- for (unsigned int i=0; i<dofs_per_cell_row; ++i)
- sparsity.add_entries (local_dof_indices_row[i],
- local_dof_indices_col.begin(),
- local_dof_indices_col.end());
- }
- else if (cell_row->has_children())
- {
- const std::vector<typename DH::active_cell_iterator >
- child_cells = GridTools::get_active_child_cells<DH> (cell_row);
- for (unsigned int i=0; i<child_cells.size(); i++)
- {
- const typename DH::active_cell_iterator
- cell_row_child = child_cells[i];
- const unsigned int dofs_per_cell_row =
- cell_row_child->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col->get_fe().dofs_per_cell;
- std::vector<unsigned int>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<unsigned int>
- local_dof_indices_col(dofs_per_cell_col);
- cell_row_child->get_dof_indices (local_dof_indices_row);
- cell_col->get_dof_indices (local_dof_indices_col);
- for (unsigned int i=0; i<dofs_per_cell_row; ++i)
- sparsity.add_entries (local_dof_indices_row[i],
- local_dof_indices_col.begin(),
- local_dof_indices_col.end());
- }
- }
- else
- {
- std::vector<typename DH::active_cell_iterator>
- child_cells = GridTools::get_active_child_cells<DH> (cell_col);
- for (unsigned int i=0; i<child_cells.size(); i++)
- {
- const typename DH::active_cell_iterator
- cell_col_child = child_cells[i];
- const unsigned int dofs_per_cell_row =
- cell_row->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col_child->get_fe().dofs_per_cell;
- std::vector<unsigned int>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<unsigned int>
- local_dof_indices_col(dofs_per_cell_col);
- cell_row->get_dof_indices (local_dof_indices_row);
- cell_col_child->get_dof_indices (local_dof_indices_col);
- for (unsigned int i=0; i<dofs_per_cell_row; ++i)
- sparsity.add_entries (local_dof_indices_row[i],
- local_dof_indices_col.begin(),
- local_dof_indices_col.end());
- }
- }
++ if (!cell_row->has_children() && !cell_col->has_children())
++ {
++ const unsigned int dofs_per_cell_row =
++ cell_row->get_fe().dofs_per_cell;
++ const unsigned int dofs_per_cell_col =
++ cell_col->get_fe().dofs_per_cell;
++ std::vector<types::global_dof_index>
++ local_dof_indices_row(dofs_per_cell_row);
++ std::vector<types::global_dof_index>
++ local_dof_indices_col(dofs_per_cell_col);
++ cell_row->get_dof_indices (local_dof_indices_row);
++ cell_col->get_dof_indices (local_dof_indices_col);
++ for (unsigned int i=0; i<dofs_per_cell_row; ++i)
++ sparsity.add_entries (local_dof_indices_row[i],
++ local_dof_indices_col.begin(),
++ local_dof_indices_col.end());
++ }
++ else if (cell_row->has_children())
++ {
++ const std::vector<typename DH::active_cell_iterator >
++ child_cells = GridTools::get_active_child_cells<DH> (cell_row);
++ for (unsigned int i=0; i<child_cells.size(); i++)
+ {
- const unsigned int dofs_per_cell_row =
- cell_row->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col->get_fe().dofs_per_cell;
- std::vector<types::global_dof_index>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<types::global_dof_index>
- local_dof_indices_col(dofs_per_cell_col);
- cell_row->get_dof_indices (local_dof_indices_row);
- cell_col->get_dof_indices (local_dof_indices_col);
- for (unsigned int i=0; i<dofs_per_cell_row; ++i)
- sparsity.add_entries (local_dof_indices_row[i],
- local_dof_indices_col.begin(),
- local_dof_indices_col.end());
++ const typename DH::active_cell_iterator
++ cell_row_child = child_cells[i];
++ const unsigned int dofs_per_cell_row =
++ cell_row_child->get_fe().dofs_per_cell;
++ const unsigned int dofs_per_cell_col =
++ cell_col->get_fe().dofs_per_cell;
++ std::vector<types::global_dof_index>
++ local_dof_indices_row(dofs_per_cell_row);
++ std::vector<types::global_dof_index>
++ local_dof_indices_col(dofs_per_cell_col);
++ cell_row_child->get_dof_indices (local_dof_indices_row);
++ cell_col->get_dof_indices (local_dof_indices_col);
++ for (unsigned int i=0; i<dofs_per_cell_row; ++i)
++ sparsity.add_entries (local_dof_indices_row[i],
++ local_dof_indices_col.begin(),
++ local_dof_indices_col.end());
+ }
- else if (cell_row->has_children())
++ }
++ else
++ {
++ std::vector<typename DH::active_cell_iterator>
++ child_cells = GridTools::get_active_child_cells<DH> (cell_col);
++ for (unsigned int i=0; i<child_cells.size(); i++)
+ {
- const std::vector<typename DH::active_cell_iterator >
- child_cells = GridTools::get_active_child_cells<DH> (cell_row);
- for (unsigned int i=0; i<child_cells.size(); i++)
- {
- const typename DH::active_cell_iterator
- cell_row_child = child_cells[i];
- const unsigned int dofs_per_cell_row =
- cell_row_child->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col->get_fe().dofs_per_cell;
- std::vector<types::global_dof_index>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<types::global_dof_index>
- local_dof_indices_col(dofs_per_cell_col);
- cell_row_child->get_dof_indices (local_dof_indices_row);
- cell_col->get_dof_indices (local_dof_indices_col);
- for (unsigned int i=0; i<dofs_per_cell_row; ++i)
- sparsity.add_entries (local_dof_indices_row[i],
- local_dof_indices_col.begin(),
- local_dof_indices_col.end());
- }
- }
- else
- {
- std::vector<typename DH::active_cell_iterator>
- child_cells = GridTools::get_active_child_cells<DH> (cell_col);
- for (unsigned int i=0; i<child_cells.size(); i++)
- {
- const typename DH::active_cell_iterator
- cell_col_child = child_cells[i];
- const unsigned int dofs_per_cell_row =
- cell_row->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col_child->get_fe().dofs_per_cell;
- std::vector<types::global_dof_index>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<types::global_dof_index>
- local_dof_indices_col(dofs_per_cell_col);
- cell_row->get_dof_indices (local_dof_indices_row);
- cell_col_child->get_dof_indices (local_dof_indices_col);
- for (unsigned int i=0; i<dofs_per_cell_row; ++i)
- sparsity.add_entries (local_dof_indices_row[i],
- local_dof_indices_col.begin(),
- local_dof_indices_col.end());
- }
++ const typename DH::active_cell_iterator
++ cell_col_child = child_cells[i];
++ const unsigned int dofs_per_cell_row =
++ cell_row->get_fe().dofs_per_cell;
++ const unsigned int dofs_per_cell_col =
++ cell_col_child->get_fe().dofs_per_cell;
++ std::vector<types::global_dof_index>
++ local_dof_indices_row(dofs_per_cell_row);
++ std::vector<types::global_dof_index>
++ local_dof_indices_col(dofs_per_cell_col);
++ cell_row->get_dof_indices (local_dof_indices_row);
++ cell_col_child->get_dof_indices (local_dof_indices_col);
++ for (unsigned int i=0; i<dofs_per_cell_row; ++i)
++ sparsity.add_entries (local_dof_indices_row[i],
++ local_dof_indices_col.begin(),
++ local_dof_indices_col.end());
+ }
++ }
}
}
#ifdef DEBUG
if (sparsity.n_rows() != 0)
{
- unsigned int max_element = 0;
- for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != DH::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- AssertDimension (max_element, sparsity.n_rows()-1);
- unsigned int max_element = 0;
- for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != DH::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- AssertDimension (max_element, sparsity.n_rows()-1);
++ types::global_dof_index max_element = 0;
++ for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
++ i!=dof_to_boundary_mapping.end(); ++i)
++ if ((*i != DH::invalid_dof_index) &&
++ (*i > max_element))
++ max_element = *i;
++ AssertDimension (max_element, sparsity.n_rows()-1);
};
#endif
// 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);
- while (!cell->active())
- cell = cell->child(direction);
+ typename DH::cell_iterator cell = dof.begin(0);
+ while (!cell->at_boundary(direction))
+ cell = cell->neighbor(direction);
+ while (!cell->active())
+ cell = cell->child(direction);
- const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
- std::vector<types::global_dof_index> boundary_dof_boundary_indices (dofs_per_vertex);
- const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
- std::vector<unsigned int> boundary_dof_boundary_indices (dofs_per_vertex);
++ const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
++ std::vector<types::global_dof_index> boundary_dof_boundary_indices (dofs_per_vertex);
// next get boundary mapped dof
// indices of boundary dofs
#ifdef DEBUG
if (sparsity.n_rows() != 0)
{
- unsigned int max_element = 0;
- for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != DH::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- AssertDimension (max_element, sparsity.n_rows()-1);
- unsigned int max_element = 0;
- for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != DH::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- AssertDimension (max_element, sparsity.n_rows()-1);
++ types::global_dof_index max_element = 0;
++ for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
++ i!=dof_to_boundary_mapping.end(); ++i)
++ if ((*i != DH::invalid_dof_index) &&
++ (*i > max_element))
++ max_element = *i;
++ AssertDimension (max_element, sparsity.n_rows()-1);
};
#endif
- std::vector<unsigned int> dofs_on_this_face;
+ std::vector<types::global_dof_index> dofs_on_this_face;
dofs_on_this_face.reserve (max_dofs_per_face(dof));
typename DH::active_cell_iterator cell = dof.begin_active(),
- endc = dof.end();
+ endc = dof.end();
for (; cell!=endc; ++cell)
for (unsigned int f=0; f<GeometryInfo<DH::dimension>::faces_per_cell; ++f)
- if (boundary_indicators.find(cell->face(f)->boundary_indicator()) !=
- boundary_indicators.end())
- {
- const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
- dofs_on_this_face.resize (dofs_per_face);
- cell->face(f)->get_dof_indices (dofs_on_this_face,
- cell->active_fe_index());
+ if (boundary_indicators.find(cell->face(f)->boundary_indicator()) !=
+ boundary_indicators.end())
+ {
+ const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
+ dofs_on_this_face.resize (dofs_per_face);
+ cell->face(f)->get_dof_indices (dofs_on_this_face,
+ cell->active_fe_index());
// make sparsity pattern for this cell
- for (unsigned int i=0; i<dofs_per_face; ++i)
- for (unsigned int j=0; j<dofs_per_face; ++j)
- sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
- dof_to_boundary_mapping[dofs_on_this_face[j]]);
- }
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ sparsity.add (dof_to_boundary_mapping[dofs_on_this_face[i]],
+ dof_to_boundary_mapping[dofs_on_this_face[j]]);
+ }
}
||
(subdomain_id == dof.get_tria().locally_owned_subdomain()),
ExcMessage ("For parallel::distributed::Triangulation objects and "
- "associated DoF handler objects, asking for any subdomain other "
- "than the locally owned one does not make sense."));
+ "associated DoF handler objects, asking for any subdomain other "
+ "than the locally owned one does not make sense."));
- std::vector<unsigned int> dofs_on_this_cell;
- std::vector<unsigned int> dofs_on_other_cell;
+ std::vector<types::global_dof_index> dofs_on_this_cell;
+ std::vector<types::global_dof_index> dofs_on_other_cell;
dofs_on_this_cell.reserve (max_dofs_per_cell(dof));
dofs_on_other_cell.reserve (max_dofs_per_cell(dof));
typename DH::active_cell_iterator cell = dof.begin_active(),
template <class DH, class SparsityPattern>
void
make_flux_sparsity_pattern (const DH &dof,
- SparsityPattern &sparsity,
- const Table<2,Coupling> &int_mask,
- const Table<2,Coupling> &flux_mask)
+ SparsityPattern &sparsity,
+ const Table<2,Coupling> &int_mask,
+ const Table<2,Coupling> &flux_mask)
{
- const FiniteElement<DH::dimension,DH::space_dimension> &fe = dof.get_fe();
+ const FiniteElement<DH::dimension,DH::space_dimension> &fe = dof.get_fe();
- std::vector<types::global_dof_index> dofs_on_this_cell(fe.dofs_per_cell);
- std::vector<types::global_dof_index> dofs_on_other_cell(fe.dofs_per_cell);
- std::vector<unsigned int> dofs_on_this_cell(fe.dofs_per_cell);
- std::vector<unsigned int> dofs_on_other_cell(fe.dofs_per_cell);
++ std::vector<types::global_dof_index> dofs_on_this_cell(fe.dofs_per_cell);
++ std::vector<types::global_dof_index> dofs_on_other_cell(fe.dofs_per_cell);
- const Table<2,Coupling>
- int_dof_mask = dof_couplings_from_component_couplings(fe, int_mask),
- flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask);
+ const Table<2,Coupling>
+ int_dof_mask = dof_couplings_from_component_couplings(fe, int_mask),
+ flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask);
- Table<2,bool> support_on_face(fe.dofs_per_cell,
- GeometryInfo<DH::dimension>::faces_per_cell);
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int f=0; f<GeometryInfo<DH::dimension>::faces_per_cell;++f)
- support_on_face(i,f) = fe.has_support_on_face(i,f);
+ Table<2,bool> support_on_face(fe.dofs_per_cell,
+ GeometryInfo<DH::dimension>::faces_per_cell);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int f=0; f<GeometryInfo<DH::dimension>::faces_per_cell;++f)
+ support_on_face(i,f) = fe.has_support_on_face(i,f);
- typename DH::active_cell_iterator cell = dof.begin_active(),
- endc = dof.end();
- for (; cell!=endc; ++cell)
+ typename DH::active_cell_iterator cell = dof.begin_active(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
- {
- cell->get_dof_indices (dofs_on_this_cell);
+ {
+ cell->get_dof_indices (dofs_on_this_cell);
// make sparsity pattern for this cell
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- if (int_dof_mask(i,j) != none)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ if (int_dof_mask(i,j) != none)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
// Loop over all interior neighbors
- for (unsigned int face = 0;
- face < GeometryInfo<DH::dimension>::faces_per_cell;
- ++face)
- {
- const typename DH::face_iterator
- cell_face = cell->face(face);
- if (cell_face->user_flag_set ())
- continue;
-
- if (cell->at_boundary (face) )
- {
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- {
- const bool i_non_zero_i = support_on_face (i, face);
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- {
- const bool j_non_zero_i = support_on_face (j, face);
-
- if ((flux_dof_mask(i,j) == always)
- ||
- (flux_dof_mask(i,j) == nonzero
- &&
- i_non_zero_i
- &&
- j_non_zero_i))
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- }
- }
- }
- else
- {
- typename DH::cell_iterator
- neighbor = cell->neighbor(face);
+ for (unsigned int face = 0;
+ face < GeometryInfo<DH::dimension>::faces_per_cell;
+ ++face)
+ {
+ const typename DH::face_iterator
+ cell_face = cell->face(face);
+ if (cell_face->user_flag_set ())
+ continue;
+
+ if (cell->at_boundary (face) )
+ {
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ const bool i_non_zero_i = support_on_face (i, face);
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ {
+ const bool j_non_zero_i = support_on_face (j, face);
+
+ if ((flux_dof_mask(i,j) == always)
+ ||
+ (flux_dof_mask(i,j) == nonzero
+ &&
+ i_non_zero_i
+ &&
+ j_non_zero_i))
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ }
+ }
+ }
+ else
+ {
+ typename DH::cell_iterator
+ neighbor = cell->neighbor(face);
// Refinement edges are taken care of
// by coarser cells
- if (cell->neighbor_is_coarser(face))
- continue;
-
- typename DH::face_iterator cell_face = cell->face(face);
- const unsigned int
- neighbor_face = cell->neighbor_of_neighbor(face);
-
- if (cell_face->has_children())
- {
- for (unsigned int sub_nr = 0;
- sub_nr != cell_face->n_children();
- ++sub_nr)
- {
- const typename DH::cell_iterator
- sub_neighbor
- = cell->neighbor_child_on_subface (face, sub_nr);
-
- sub_neighbor->get_dof_indices (dofs_on_other_cell);
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- {
- const bool i_non_zero_i = support_on_face (i, face);
- const bool i_non_zero_e = support_on_face (i, neighbor_face);
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- {
- const bool j_non_zero_i = support_on_face (j, face);
- const bool j_non_zero_e = support_on_face (j, neighbor_face);
-
- if (flux_dof_mask(i,j) == always)
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
- else if (flux_dof_mask(i,j) == nonzero)
- {
- if (i_non_zero_i && j_non_zero_e)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- if (i_non_zero_e && j_non_zero_i)
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- if (i_non_zero_i && j_non_zero_i)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- if (i_non_zero_e && j_non_zero_e)
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
-
- if (flux_dof_mask(j,i) == always)
- {
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- else if (flux_dof_mask(j,i) == nonzero)
- {
- if (j_non_zero_i && i_non_zero_e)
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- if (j_non_zero_e && i_non_zero_i)
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- if (j_non_zero_i && i_non_zero_i)
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- if (j_non_zero_e && i_non_zero_e)
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- }
- }
- sub_neighbor->face(neighbor_face)->set_user_flag ();
- }
- }
- else
- {
- neighbor->get_dof_indices (dofs_on_other_cell);
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- {
- const bool i_non_zero_i = support_on_face (i, face);
- const bool i_non_zero_e = support_on_face (i, neighbor_face);
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- {
- const bool j_non_zero_i = support_on_face (j, face);
- const bool j_non_zero_e = support_on_face (j, neighbor_face);
- if (flux_dof_mask(i,j) == always)
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
- if (flux_dof_mask(i,j) == nonzero)
- {
- if (i_non_zero_i && j_non_zero_e)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- if (i_non_zero_e && j_non_zero_i)
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- if (i_non_zero_i && j_non_zero_i)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- if (i_non_zero_e && j_non_zero_e)
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
-
- if (flux_dof_mask(j,i) == always)
- {
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- if (flux_dof_mask(j,i) == nonzero)
- {
- if (j_non_zero_i && i_non_zero_e)
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- if (j_non_zero_e && i_non_zero_i)
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- if (j_non_zero_i && i_non_zero_i)
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- if (j_non_zero_e && i_non_zero_e)
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- }
- }
- neighbor->face(neighbor_face)->set_user_flag ();
- }
- }
- }
+ if (cell->neighbor_is_coarser(face))
+ continue;
+
+ typename DH::face_iterator cell_face = cell->face(face);
+ const unsigned int
+ neighbor_face = cell->neighbor_of_neighbor(face);
+
+ if (cell_face->has_children())
+ {
+ for (unsigned int sub_nr = 0;
+ sub_nr != cell_face->n_children();
+ ++sub_nr)
+ {
+ const typename DH::cell_iterator
+ sub_neighbor
+ = cell->neighbor_child_on_subface (face, sub_nr);
+
+ sub_neighbor->get_dof_indices (dofs_on_other_cell);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ const bool i_non_zero_i = support_on_face (i, face);
+ const bool i_non_zero_e = support_on_face (i, neighbor_face);
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ {
+ const bool j_non_zero_i = support_on_face (j, face);
+ const bool j_non_zero_e = support_on_face (j, neighbor_face);
+
+ if (flux_dof_mask(i,j) == always)
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+ else if (flux_dof_mask(i,j) == nonzero)
+ {
+ if (i_non_zero_i && j_non_zero_e)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ if (i_non_zero_e && j_non_zero_i)
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ if (i_non_zero_i && j_non_zero_i)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ if (i_non_zero_e && j_non_zero_e)
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+
+ if (flux_dof_mask(j,i) == always)
+ {
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ else if (flux_dof_mask(j,i) == nonzero)
+ {
+ if (j_non_zero_i && i_non_zero_e)
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ if (j_non_zero_e && i_non_zero_i)
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ if (j_non_zero_i && i_non_zero_i)
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ if (j_non_zero_e && i_non_zero_e)
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ }
+ }
+ sub_neighbor->face(neighbor_face)->set_user_flag ();
+ }
+ }
+ else
+ {
+ neighbor->get_dof_indices (dofs_on_other_cell);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ const bool i_non_zero_i = support_on_face (i, face);
+ const bool i_non_zero_e = support_on_face (i, neighbor_face);
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ {
+ const bool j_non_zero_i = support_on_face (j, face);
+ const bool j_non_zero_e = support_on_face (j, neighbor_face);
+ if (flux_dof_mask(i,j) == always)
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+ if (flux_dof_mask(i,j) == nonzero)
+ {
+ if (i_non_zero_i && j_non_zero_e)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ if (i_non_zero_e && j_non_zero_i)
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ if (i_non_zero_i && j_non_zero_i)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ if (i_non_zero_e && j_non_zero_e)
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+
+ if (flux_dof_mask(j,i) == always)
+ {
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ if (flux_dof_mask(j,i) == nonzero)
+ {
+ if (j_non_zero_i && i_non_zero_e)
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ if (j_non_zero_e && i_non_zero_i)
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ if (j_non_zero_i && i_non_zero_i)
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ if (j_non_zero_e && i_non_zero_e)
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ }
+ }
+ neighbor->face(neighbor_face)->set_user_flag ();
+ }
+ }
+ }
}
}
// probably less efficient but at least
// readable...
- const dealii::hp::FECollection<dim,spacedim> &fe = dof.get_fe();
+ const dealii::hp::FECollection<dim,spacedim> &fe = dof.get_fe();
- std::vector<types::global_dof_index> dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof));
- std::vector<types::global_dof_index> dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof));
- std::vector<unsigned int> dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof));
- std::vector<unsigned int> dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof));
++ std::vector<types::global_dof_index> dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof));
++ std::vector<types::global_dof_index> dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof));
-
const std::vector<Table<2,Coupling> >
- int_dof_mask
- = dof_couplings_from_component_couplings(fe, int_mask);
+ const std::vector<Table<2,Coupling> >
+ int_dof_mask
+ = dof_couplings_from_component_couplings(fe, int_mask);
- typename dealii::hp::DoFHandler<dim,spacedim>::active_cell_iterator
- cell = dof.begin_active(),
- endc = dof.end();
- for (; cell!=endc; ++cell)
- {
- dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell);
- cell->get_dof_indices (dofs_on_this_cell);
+ typename dealii::hp::DoFHandler<dim,spacedim>::active_cell_iterator
+ cell = dof.begin_active(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ {
+ 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)
- if (int_dof_mask[cell->active_fe_index()](i,j) != none)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
+ 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)
+ if (int_dof_mask[cell->active_fe_index()](i,j) != none)
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
// Loop over all interior neighbors
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell;
- ++face)
- {
- const typename dealii::hp::DoFHandler<dim,spacedim>::face_iterator
- cell_face = cell->face(face);
- if (cell_face->user_flag_set ())
- continue;
-
- if (cell->at_boundary (face) )
- {
- 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)
- if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
- cell->get_fe().system_to_component_index(j).first)
- == always)
- ||
- (flux_mask(cell->get_fe().system_to_component_index(i).first,
- cell->get_fe().system_to_component_index(j).first)
- == nonzero))
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- }
- else
- {
- typename dealii::hp::DoFHandler<dim,spacedim>::cell_iterator
- neighbor = cell->neighbor(face);
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ {
+ const typename dealii::hp::DoFHandler<dim,spacedim>::face_iterator
+ cell_face = cell->face(face);
+ if (cell_face->user_flag_set ())
+ continue;
+
+ if (cell->at_boundary (face) )
+ {
+ 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)
+ if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
+ cell->get_fe().system_to_component_index(j).first)
+ == always)
+ ||
+ (flux_mask(cell->get_fe().system_to_component_index(i).first,
+ cell->get_fe().system_to_component_index(j).first)
+ == nonzero))
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ }
+ else
+ {
+ typename dealii::hp::DoFHandler<dim,spacedim>::cell_iterator
+ neighbor = cell->neighbor(face);
// Refinement edges are taken care of
// by coarser cells
- if (cell->neighbor_is_coarser(face))
- continue;
-
- typename dealii::hp::DoFHandler<dim,spacedim>::face_iterator
- cell_face = cell->face(face);
- const unsigned int
- neighbor_face = cell->neighbor_of_neighbor(face);
-
- if (cell_face->has_children())
- {
- for (unsigned int sub_nr = 0;
- sub_nr != cell_face->n_children();
- ++sub_nr)
- {
- const typename dealii::hp::DoFHandler<dim,spacedim>::cell_iterator
- sub_neighbor
- = cell->neighbor_child_on_subface (face, sub_nr);
-
- dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell);
- sub_neighbor->get_dof_indices (dofs_on_other_cell);
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<sub_neighbor->get_fe().dofs_per_cell;
- ++j)
- {
- if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
- sub_neighbor->get_fe().system_to_component_index(j).first)
- == always)
- ||
- (flux_mask(cell->get_fe().system_to_component_index(i).first,
- sub_neighbor->get_fe().system_to_component_index(j).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
-
- if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == always)
- ||
- (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- }
- }
- sub_neighbor->face(neighbor_face)->set_user_flag ();
- }
- }
- else
- {
- dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell);
- neighbor->get_dof_indices (dofs_on_other_cell);
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<neighbor->get_fe().dofs_per_cell; ++j)
- {
- if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
- neighbor->get_fe().system_to_component_index(j).first)
- == always)
- ||
- (flux_mask(cell->get_fe().system_to_component_index(i).first,
- neighbor->get_fe().system_to_component_index(j).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
-
- if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == always)
- ||
- (flux_mask(neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- }
- }
- neighbor->face(neighbor_face)->set_user_flag ();
- }
- }
- }
- }
+ if (cell->neighbor_is_coarser(face))
+ continue;
+
+ typename dealii::hp::DoFHandler<dim,spacedim>::face_iterator
+ cell_face = cell->face(face);
+ const unsigned int
+ neighbor_face = cell->neighbor_of_neighbor(face);
+
+ if (cell_face->has_children())
+ {
+ for (unsigned int sub_nr = 0;
+ sub_nr != cell_face->n_children();
+ ++sub_nr)
+ {
+ const typename dealii::hp::DoFHandler<dim,spacedim>::cell_iterator
+ sub_neighbor
+ = cell->neighbor_child_on_subface (face, sub_nr);
+
+ dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell);
+ sub_neighbor->get_dof_indices (dofs_on_other_cell);
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<sub_neighbor->get_fe().dofs_per_cell;
+ ++j)
+ {
+ if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
+ sub_neighbor->get_fe().system_to_component_index(j).first)
+ == always)
+ ||
+ (flux_mask(cell->get_fe().system_to_component_index(i).first,
+ sub_neighbor->get_fe().system_to_component_index(j).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+
+ if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == always)
+ ||
+ (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ }
+ }
+ sub_neighbor->face(neighbor_face)->set_user_flag ();
+ }
+ }
+ else
+ {
+ dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell);
+ neighbor->get_dof_indices (dofs_on_other_cell);
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<neighbor->get_fe().dofs_per_cell; ++j)
+ {
+ if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
+ neighbor->get_fe().system_to_component_index(j).first)
+ == always)
+ ||
+ (flux_mask(cell->get_fe().system_to_component_index(i).first,
+ neighbor->get_fe().system_to_component_index(j).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+
+ if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == always)
+ ||
+ (flux_mask(neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ }
+ }
+ neighbor->face(neighbor_face)->set_user_flag ();
+ }
+ }
+ }
+ }
}
}
inline
bool
check_master_dof_list (const FullMatrix<double> &face_interpolation_matrix,
- const std::vector<types::global_dof_index> &master_dof_list)
- const std::vector<unsigned int> &master_dof_list)
++ const std::vector<types::global_dof_index> &master_dof_list)
{
- const unsigned int N = master_dof_list.size();
+ const unsigned int N = master_dof_list.size();
- FullMatrix<double> tmp (N,N);
- for (unsigned int i=0; i<N; ++i)
- for (unsigned int j=0; j<N; ++j)
- tmp(i,j) = face_interpolation_matrix (master_dof_list[i], j);
+ FullMatrix<double> tmp (N,N);
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ tmp(i,j) = face_interpolation_matrix (master_dof_list[i], j);
// then use the algorithm
// from
// the other one. this latter case
// shows up when running
// hp/hp_constraints_q_system_06
- std::vector<types::global_dof_index> master_dof_list;
- unsigned int index = 0;
- for (int v=0;
- v<static_cast<signed int>(GeometryInfo<dim>::vertices_per_face);
- ++v)
- std::vector<unsigned int> master_dof_list;
- unsigned int index = 0;
- for (int v=0;
- v<static_cast<signed int>(GeometryInfo<dim>::vertices_per_face);
- ++v)
- {
- unsigned int dofs_added = 0;
- unsigned int i = 0;
- while (dofs_added < fe2.dofs_per_vertex)
- {
++ std::vector<types::global_dof_index> master_dof_list;
++ unsigned int index = 0;
++ for (int v=0;
++ v<static_cast<signed int>(GeometryInfo<dim>::vertices_per_face);
++ ++v)
++ {
++ unsigned int dofs_added = 0;
++ unsigned int i = 0;
++ while (dofs_added < fe2.dofs_per_vertex)
+ {
- unsigned int dofs_added = 0;
- unsigned int i = 0;
- while (dofs_added < fe2.dofs_per_vertex)
- {
// make sure that we
// were able to find
// a set of master
// finally copy the list into the
// mask
- std::fill (master_dof_mask.begin(), master_dof_mask.end(), false);
- for (std::vector<types::global_dof_index>::const_iterator i=master_dof_list.begin();
- i!=master_dof_list.end(); ++i)
- master_dof_mask[*i] = true;
- 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)
- master_dof_mask[*i] = true;
++ std::fill (master_dof_mask.begin(), master_dof_mask.end(), false);
++ for (std::vector<types::global_dof_index>::const_iterator i=master_dof_list.begin();
++ i!=master_dof_list.end(); ++i)
++ master_dof_mask[*i] = true;
}
* fuller than necessary.
*/
void
- filter_constraints (const std::vector<unsigned int> &master_dofs,
- const std::vector<unsigned int> &slave_dofs,
- const FullMatrix<double> &face_constraints,
- ConstraintMatrix &constraints)
+ filter_constraints (const std::vector<types::global_dof_index> &master_dofs,
- const std::vector<types::global_dof_index> &slave_dofs,
- const FullMatrix<double> &face_constraints,
- ConstraintMatrix &constraints)
++ const std::vector<types::global_dof_index> &slave_dofs,
++ const FullMatrix<double> &face_constraints,
++ ConstraintMatrix &constraints)
{
- Assert (face_constraints.n () == master_dofs.size (),
- ExcDimensionMismatch(master_dofs.size (),
- face_constraints.n()));
- Assert (face_constraints.m () == slave_dofs.size (),
- ExcDimensionMismatch(slave_dofs.size (),
- face_constraints.m()));
+ Assert (face_constraints.n () == master_dofs.size (),
+ ExcDimensionMismatch(master_dofs.size (),
+ face_constraints.n()));
+ Assert (face_constraints.m () == slave_dofs.size (),
+ ExcDimensionMismatch(slave_dofs.size (),
+ face_constraints.m()));
- const unsigned int n_master_dofs = master_dofs.size ();
- const unsigned int n_slave_dofs = slave_dofs.size ();
+ 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
}
return;
}
-- // .. otherwise we should be in the case
-- // were both faces are active and have
-- // no children ..
++ // .. otherwise we should be in the case
++ // were both faces are active and have
++ // no children ..
Assert (!face_1->has_children() && !face_2->has_children(),
ExcNotImplemented());
Assert (face_1->n_active_fe_indices() == 1 && face_2->n_active_fe_indices() == 1,
ExcInternalError());
-- // .. then we match the
-- // corresponding DoFs of both faces ..
++ // .. then we match the
++ // corresponding DoFs of both faces ..
const unsigned int face_1_index = face_1->nth_active_fe_index(0);
const unsigned int face_2_index = face_2->nth_active_fe_index(0);
Assert ( face_1->get_fe(face_1_index)
-- == face_2->get_fe(face_1_index),
-- ExcMessage ("Matching periodic cells need to use the same finite element"));
++ == face_2->get_fe(face_1_index),
++ ExcMessage ("Matching periodic cells need to use the same finite element"));
const dealii::FiniteElement<dim> &fe = face_1->get_fe(face_1_index);
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,
typedef typename DH::cell_iterator CellIterator;
-- // We collect matching periodic cells on
-- // the coarsest level:
++ // We collect matching periodic cells on
++ // the coarsest level:
std::map<CellIterator, CellIterator>
matched_cells =
-- GridTools::collect_periodic_cell_pairs(dof_handler.begin(0),
-- dof_handler.end(0),
-- boundary_component,
-- direction,
-- offset);
--
-- // And apply the low level
-- // make_periodicity_constraints function
-- // to every matching pair:
++ GridTools::collect_periodic_cell_pairs(dof_handler.begin(0),
++ dof_handler.end(0),
++ boundary_component,
++ direction,
++ offset);
++
++ // And apply the low level
++ // make_periodicity_constraints function
++ // to every matching pair:
for (typename std::map<CellIterator, CellIterator>::iterator it = matched_cells.begin();
- it != matched_cells.end(); ++it)
+ it != matched_cells.end(); ++it)
{
- typedef typename DH::face_iterator FaceIterator;
- FaceIterator face_1 = it->first->face(2*direction);
- FaceIterator face_2 = it->second->face(2*direction+1);
+ typedef typename DH::face_iterator FaceIterator;
+ FaceIterator face_1 = it->first->face(2*direction);
+ FaceIterator face_2 = it->second->face(2*direction+1);
- Assert(face_1->at_boundary() && face_2->at_boundary(),
- ExcInternalError());
+ Assert(face_1->at_boundary() && face_2->at_boundary(),
+ ExcInternalError());
- Assert (face_1->boundary_indicator() == boundary_component &&
- face_2->boundary_indicator() == boundary_component,
- ExcInternalError());
+ Assert (face_1->boundary_indicator() == boundary_component &&
+ face_2->boundary_indicator() == boundary_component,
+ ExcInternalError());
- make_periodicity_constraints(face_1,
- face_2,
- constraint_matrix,
- component_mask);
+ make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask);
}
}
for (unsigned int f=0; f<fe_collection.size(); ++f)
{
const FiniteElement<DH::dimension,DH::space_dimension> &fe =
- fe_collection[f];
- local_component_association[f].resize(fe.dofs_per_cell);
- if (sort_by_blocks == true)
- // compute the block each
- // local dof belongs to
- {
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- local_component_association[f][i]
- = fe.system_to_block_index(i).first;
- }
- else
- // compute the component each
- // local dof belongs to
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- if (fe.is_primitive(i))
- local_component_association[f][i] =
- fe.system_to_component_index(i).first;
- else
- // if this shape function is
- // not primitive, then we have
- // to work harder. we have to
- // find out whether _any_ of
- // the vector components of
- // this element is selected or
- // not
- //
- // to do so, get the a list of
- // nonzero elements and see which are
- // actually active
- {
- const unsigned int first_comp =
- (std::find(fe.get_nonzero_components(i).begin(),
- fe.get_nonzero_components(i).end(),
- true) -
- fe.get_nonzero_components(i).begin());
- const unsigned int end_comp =
- (std::find(fe.get_nonzero_components(i).begin()+first_comp,
- fe.get_nonzero_components(i).end(),
- false)-
- fe.get_nonzero_components(i).begin());
-
- // now check whether any of
- // the components in between
- // is set
- if (component_select.size() == 0 ||
- (component_select[first_comp] == true ||
- std::count(component_select.begin()+first_comp,
- component_select.begin()+end_comp, true) == 0))
- local_component_association[f][i] = first_comp;
- else
- for (unsigned int c=first_comp; c<end_comp; ++c)
- if (component_select[c] == true)
- {
- local_component_association[f][i] = c;
- break;
- }
- }
- }
-
- // then loop over all cells and do
- // the work
+ fe_collection[f];
+ local_component_association[f]
+ = get_local_component_association (fe, component_mask);
+ }
+
+ // 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)
+ if (c->is_locally_owned())
+ {
+ const unsigned int fe_index = c->active_fe_index();
+ const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell;
+ indices.resize(dofs_per_cell);
+ c->get_dof_indices(indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (dof.locally_owned_dofs().is_element(indices[i]))
+ dofs_by_component[dof.locally_owned_dofs().index_within_set(indices[i])]
+ = local_component_association[fe_index][i];
+ }
+ }
+
+
+ // 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
+ // dof.locally_owned_dofs().index_within_set(indices[i])
+ template <class DH>
+ inline
+ void
+ get_block_association (const DH &dof,
+ std::vector<unsigned char> &dofs_by_block)
+ {
+ const dealii::hp::FECollection<DH::dimension,DH::space_dimension>
+ fe_collection (dof.get_fe());
+ Assert (fe_collection.n_components() < 256, ExcNotImplemented());
+ Assert (dofs_by_block.size() == dof.n_locally_owned_dofs(),
+ 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)
+ //
+ // 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)
+ {
+ const FiniteElement<DH::dimension,DH::space_dimension> &fe =
+ fe_collection[f];
+ local_block_association[f].resize(fe.dofs_per_cell,
+ (unsigned char)(-1));
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ local_block_association[f][i] = fe.system_to_block_index(i).first;
+
+ Assert (std::find (local_block_association[f].begin(),
+ local_block_association[f].end(),
+ (unsigned char)(-1))
+ ==
+ local_block_association[f].end(),
+ ExcInternalError());
+ }
+
+ // then loop over all cells and do
+ // the work
- std::vector<unsigned int> indices;
+ std::vector<types::global_dof_index> indices;
for (typename DH::active_cell_iterator c=dof.begin_active();
- c!=dof.end(); ++ c)
- if (c->is_locally_owned())
- {
- const unsigned int fe_index = c->active_fe_index();
- const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell;
- indices.resize(dofs_per_cell);
- c->get_dof_indices(indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- if (dof.locally_owned_dofs().is_element(indices[i]))
- dofs_by_component[dof.locally_owned_dofs().index_within_set(indices[i])]
- = local_component_association[fe_index][i];
- }
+ c!=dof.end(); ++ c)
+ if (c->is_locally_owned())
+ {
+ const unsigned int fe_index = c->active_fe_index();
+ const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell;
+ indices.resize(dofs_per_cell);
+ c->get_dof_indices(indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (dof.locally_owned_dofs().is_element(indices[i]))
+ dofs_by_block[dof.locally_owned_dofs().index_within_set(indices[i])]
+ = local_block_association[fe_index][i];
+ }
}
}
std::vector<unsigned char> touch_count (dof_handler.n_dofs(), 0);
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
- std::vector<unsigned int> dof_indices;
+ std::vector<types::global_dof_index> dof_indices;
dof_indices.reserve (max_dofs_per_cell(dof_handler));
for (unsigned int present_cell = 0; cell!=endc; ++cell, ++present_cell)
- template<int dim, int spacedim>
+ 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)));
+ ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
// two special cases: no component
// is selected, and all components
// then loop over all cells and do
// work
- std::vector<unsigned int> indices(fe.dofs_per_cell);
- typename MGDoFHandler<dim,spacedim>::cell_iterator c;
+ std::vector<types::global_dof_index> indices(fe.dofs_per_cell);
+ typename DH::cell_iterator 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)
- selected_dofs[indices[i]] = local_selected_dofs[i];
+ c->get_mg_dof_indices(indices);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ selected_dofs[indices[i]] = local_selected_dofs[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);
-- // then convert the values computed above to the binary vector
++ // then convert the values computed above to the binary vector
indices.fill_binary_vector(selected_dofs);
}
ExcInvalidBoundaryIndicator());
const unsigned int dim=DH::dimension;
-- // first reset output argument
++ // first reset output argument
selected_dofs.clear ();
selected_dofs.set_size(dof_handler.n_dofs());
// check whether a certain vector
// component is selected, or all
const bool check_vector_component
- = (component_select != std::vector<bool>(component_select.size(),
- true));
+ = ((component_mask.represents_the_all_selected_mask() == false)
+ ||
+ (component_mask.n_selected_components(n_components(dof_handler)) !=
+ n_components(dof_handler)));
- std::vector<unsigned int> face_dof_indices;
+ std::vector<types::global_dof_index> face_dof_indices;
face_dof_indices.reserve (max_dofs_per_face(dof_handler));
// now loop over all cells and
// 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)
++ 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()))
++ != boundary_indicators.end()))
{
const FiniteElement<DH::dimension, DH::space_dimension> &fe = cell->get_fe();
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]
- == true)
+ if (component_mask[fe.face_system_to_component_index(i).first]
- == true)
++ == true)
selected_dofs.add_index (face_dof_indices[i]);
}
else // not primitive
// 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;
+ std::vector<types::global_dof_index> dof_indices;
+ std::set<types::global_dof_index> global_dof_indices;
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
- {
- dof_indices.resize(cell->get_fe().dofs_per_cell);
- cell->get_dof_indices(dof_indices);
+ {
+ dof_indices.resize(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(dof_indices);
- for (std::vector<types::global_dof_index>::iterator it=dof_indices.begin();
- it!=dof_indices.end();
- ++it)
- if (!dof_set.is_element(*it))
- global_dof_indices.insert(*it);
- }
- for (std::vector<unsigned int>::iterator it=dof_indices.begin();
- it!=dof_indices.end();
- ++it)
- if (!dof_set.is_element(*it))
- global_dof_indices.insert(*it);
- }
++ for (std::vector<types::global_dof_index>::iterator it=dof_indices.begin();
++ it!=dof_indices.end();
++ ++it)
++ if (!dof_set.is_element(*it))
++ global_dof_indices.insert(*it);
++ }
dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
// 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;
+ std::vector<types::global_dof_index> dof_indices;
+ std::set<types::global_dof_index> global_dof_indices;
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->is_ghost())
- {
- dof_indices.resize(cell->get_fe().dofs_per_cell);
- cell->get_dof_indices(dof_indices);
+ {
+ dof_indices.resize(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(dof_indices);
- for (std::vector<types::global_dof_index>::iterator it=dof_indices.begin();
- it!=dof_indices.end();
- ++it)
- if (!dof_set.is_element(*it))
- global_dof_indices.insert(*it);
- }
- for (std::vector<unsigned int>::iterator it=dof_indices.begin();
- it!=dof_indices.end();
- ++it)
- if (!dof_set.is_element(*it))
- global_dof_indices.insert(*it);
- }
++ for (std::vector<types::global_dof_index>::iterator it=dof_indices.begin();
++ it!=dof_indices.end();
++ ++it)
++ if (!dof_set.is_element(*it))
++ global_dof_indices.insert(*it);
++ }
dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
// preset all values by an invalid value
std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(),
- types::invalid_subdomain_id);
+ types::invalid_subdomain_id);
- std::vector<unsigned int> local_dof_indices;
+ std::vector<types::global_dof_index> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
// pseudo-randomly assign variables
IndexSet index_set (dof_handler.n_dofs());
- std::vector<unsigned int> local_dof_indices;
+ std::vector<types::global_dof_index> 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
- std::vector<unsigned int> subdomain_indices;
++ // 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<types::global_dof_index> subdomain_indices;
typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
local_dof_indices.begin(),
local_dof_indices.end());
}
-- // sort indices and remove duplicates
++ // sort indices and remove duplicates
std::sort (subdomain_indices.begin(), subdomain_indices.end());
subdomain_indices.erase (std::unique(subdomain_indices.begin(),
subdomain_indices.end()),
subdomain_indices.end());
-- // insert into IndexSet
++ // insert into IndexSet
index_set.add_indices (subdomain_indices.begin(), subdomain_indices.end());
index_set.compress ();
template <int dim, int spacedim>
void
resolve_components (const FiniteElement<dim,spacedim>&fe,
- const std::vector<unsigned char> &dofs_by_component,
- const std::vector<unsigned int> &target_component,
- const bool only_once,
- std::vector<types::global_dof_index> &dofs_per_component,
- unsigned int &component)
- const std::vector<unsigned char> &dofs_by_component,
- const std::vector<unsigned int> &target_component,
- const bool only_once,
- std::vector<unsigned int> &dofs_per_component,
- unsigned int &component)
++ const std::vector<unsigned char> &dofs_by_component,
++ const std::vector<unsigned int> &target_component,
++ const bool only_once,
++ std::vector<types::global_dof_index> &dofs_per_component,
++ unsigned int &component)
{
for (unsigned int b=0;b<fe.n_base_elements();++b)
- {
- const FiniteElement<dim,spacedim>& base = fe.base_element(b);
+ {
+ const FiniteElement<dim,spacedim>& 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)
- {
- if (base.n_base_elements() > 1)
- resolve_components(base, dofs_by_component, target_component,
- only_once, dofs_per_component, component);
- else
- {
- for (unsigned int dd=0;dd<d;++dd,++component)
- dofs_per_component[target_component[component]]
- += std::count(dofs_by_component.begin(),
- dofs_by_component.end(),
- component);
+ unsigned int d = base.n_components();
+
+ for (unsigned int m=0;m<fe.element_multiplicity(b);++m)
+ {
+ if (base.n_base_elements() > 1)
+ resolve_components(base, dofs_by_component, target_component,
+ only_once, dofs_per_component, component);
+ else
+ {
+ for (unsigned int dd=0;dd<d;++dd,++component)
+ dofs_per_component[target_component[component]]
+ += std::count(dofs_by_component.begin(),
+ dofs_by_component.end(),
+ component);
// if we have non-primitive FEs and want all
// components to show the number of dofs, need
const std::vector<unsigned char> &dofs_by_component,
const std::vector<unsigned int> &target_component,
const bool only_once,
- std::vector<unsigned int> &dofs_per_component,
+ std::vector<types::global_dof_index> &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(),
-- ExcNotImplemented());
-- Assert (fe_collection[fe].n_base_elements() == fe_collection[0].n_base_elements(),
-- ExcNotImplemented());
-- for (unsigned int b=0;b<fe_collection[0].n_base_elements();++b)
{
-- Assert (fe_collection[fe].base_element(b).n_components() == fe_collection[0].base_element(b).n_components(),
++ Assert (fe_collection[fe].n_components() == fe_collection[0].n_components(),
ExcNotImplemented());
-- Assert (fe_collection[fe].base_element(b).n_base_elements() == fe_collection[0].base_element(b).n_base_elements(),
++ Assert (fe_collection[fe].n_base_elements() == fe_collection[0].n_base_elements(),
ExcNotImplemented());
++ for (unsigned int b=0;b<fe_collection[0].n_base_elements();++b)
++ {
++ Assert (fe_collection[fe].base_element(b).n_components() == fe_collection[0].base_element(b).n_components(),
++ ExcNotImplemented());
++ Assert (fe_collection[fe].base_element(b).n_base_elements() == fe_collection[0].base_element(b).n_base_elements(),
++ ExcNotImplemented());
++ }
}
-- }
resolve_components (fe_collection[0], dofs_by_component,
target_component, only_once, dofs_per_component,
namespace internal
{
namespace {
-- /**
-- * Return true if the given element is primitive.
-- */
++ /**
++ * Return true if the given element is primitive.
++ */
template <int dim, int spacedim>
bool all_elements_are_primitive (const FiniteElement<dim,spacedim> &fe)
{
}
-- /**
-- * Return true if each element of the given element collection is primitive.
-- */
++ /**
++ * Return true if each element of the given element collection is primitive.
++ */
template <int dim, int spacedim>
bool all_elements_are_primitive (const dealii::hp::FECollection<dim,spacedim> &fe_collection)
{
template <int dim, int spacedim>
void
count_dofs_per_component (const DoFHandler<dim,spacedim> &dof_handler,
- std::vector<types::global_dof_index> &dofs_per_component,
- std::vector<unsigned int> target_component)
- std::vector<unsigned int> &dofs_per_component,
- std::vector<unsigned int> target_component)
++ std::vector<types::global_dof_index> &dofs_per_component,
++ std::vector<unsigned int> target_component)
{
count_dofs_per_component (dof_handler, dofs_per_component,
- false, target_component);
+ false, target_component);
}
// 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);
+ const unsigned int n_fine_dofs = weight_mapping.size();
+ dealii::Vector<double> global_parameter_representation (n_fine_dofs);
- typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator cell;
- std::vector<types::global_dof_index> parameter_dof_indices (coarse_fe.dofs_per_cell);
- typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator cell;
- std::vector<unsigned int> parameter_dof_indices (coarse_fe.dofs_per_cell);
++ typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator cell;
++ std::vector<types::global_dof_index> parameter_dof_indices (coarse_fe.dofs_per_cell);
- for (cell=begin; cell!=end; ++cell)
- {
+ for (cell=begin; cell!=end; ++cell)
+ {
// get the global indices of the
// parameter dofs on this parameter
// grid cell
// 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<types::global_dof_index> local_dof_indices (fine_fe.dofs_per_cell);
- std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
- std::vector<unsigned int> local_dof_indices (fine_fe.dofs_per_cell);
++ std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
++ std::vector<types::global_dof_index> local_dof_indices (fine_fe.dofs_per_cell);
- for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
- cell=fine_grid.begin_active();
- cell!=fine_grid.end(); ++cell)
- {
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
- if (fine_fe.system_to_component_index(i).first == fine_component)
- dof_is_interesting[local_dof_indices[i]] = true;
- };
+ for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+ cell=fine_grid.begin_active();
+ cell!=fine_grid.end(); ++cell)
+ {
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+ if (fine_fe.system_to_component_index(i).first == fine_component)
+ dof_is_interesting[local_dof_indices[i]] = true;
+ };
- n_parameters_on_fine_grid = std::count (dof_is_interesting.begin(),
- dof_is_interesting.end(),
- true);
- };
+ n_parameters_on_fine_grid = std::count (dof_is_interesting.begin(),
+ dof_is_interesting.end(),
+ true);
+ };
// set up the weights mapping
- weights.clear ();
- weights.resize (n_coarse_dofs);
+ weights.clear ();
+ weights.resize (n_coarse_dofs);
- weight_mapping.clear ();
- weight_mapping.resize (n_fine_dofs, -1);
+ weight_mapping.clear ();
+ weight_mapping.resize (n_fine_dofs, -1);
- if (true)
- if (true)
- {
- std::vector<unsigned int> local_dof_indices(fine_fe.dofs_per_cell);
- unsigned int next_free_index=0;
- for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
- cell=fine_grid.begin_active();
- cell != fine_grid.end(); ++cell)
- {
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
++ if (true)
++ {
++ std::vector<types::global_dof_index> local_dof_indices(fine_fe.dofs_per_cell);
++ unsigned int next_free_index=0;
++ for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
++ cell=fine_grid.begin_active();
++ cell != fine_grid.end(); ++cell)
+ {
- std::vector<types::global_dof_index> local_dof_indices(fine_fe.dofs_per_cell);
- unsigned int next_free_index=0;
- for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
- cell=fine_grid.begin_active();
- cell != fine_grid.end(); ++cell)
- {
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
++ 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,
// otherwise enter all constraints
- constraints.add_line (global_dof);
-
- constraint_line.clear ();
- for (unsigned int row=first_used_row; row<n_coarse_dofs; ++row)
- {
- const std::map<unsigned int,float>::const_iterator
- j = weights[row].find(col);
- if ((j != weights[row].end()) && (j->second != 0))
- constraint_line.push_back (std::pair<unsigned int,double>(representants[row],
- j->second));
- };
-
- constraints.add_entries (global_dof, constraint_line);
- };
+ constraints.add_line (global_dof);
+
+ constraint_line.clear ();
+ for (unsigned int row=first_used_row; row<n_coarse_dofs; ++row)
+ {
+ const std::map<unsigned int,float>::const_iterator
+ j = weights[row].find(col);
+ if ((j != weights[row].end()) && (j->second != 0))
+ constraint_line.push_back (std::pair<unsigned int,double>(representants[row],
- j->second));
++ j->second));
+ };
+
+ constraints.add_entries (global_dof, constraint_line);
+ };
}
template <class DH>
void
map_dof_to_boundary_indices (const DH &dof_handler,
- std::vector<types::global_dof_index> &mapping)
- std::vector<unsigned int> &mapping)
++ std::vector<types::global_dof_index> &mapping)
{
Assert (&dof_handler.get_fe() != 0, ExcNoFESelected());
mapping.clear ();
mapping.insert (mapping.end(), dof_handler.n_dofs(),
- DH::invalid_dof_index);
+ DH::invalid_dof_index);
- std::vector<unsigned int> dofs_on_face;
+ std::vector<types::global_dof_index> dofs_on_face;
dofs_on_face.reserve (max_dofs_per_face(dof_handler));
unsigned int next_boundary_index = 0;
void
map_dofs_to_support_points(const hp::MappingCollection<DH::dimension, DH::space_dimension> & mapping,
const DH &dof_handler,
- std::map<unsigned int,Point<DH::space_dimension> > &support_points)
+ std::map<types::global_dof_index,Point<DH::space_dimension> > &support_points)
{
-- const unsigned int dim = DH::dimension;
-- const unsigned int spacedim = DH::space_dimension;
-
- hp::FECollection<dim, spacedim> fe_collection(dof_handler.get_fe());
- hp::QCollection<dim> q_coll_dummy;
-
- for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index)
- {
- // 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(
- Quadrature<dim> (
- fe_collection[fe_index].get_unit_support_points()));
- }
++ const unsigned int dim = DH::dimension;
++ const unsigned int spacedim = DH::space_dimension;
- hp::FECollection<dim, spacedim> fe_collection(dof_handler.get_fe());
- hp::QCollection<dim> q_coll_dummy;
- // 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.
- hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe_collection,
- q_coll_dummy, update_quadrature_points);
- typename DH::active_cell_iterator cell =
- dof_handler.begin_active(), endc = dof_handler.end();
-
- std::vector<unsigned int> local_dof_indices;
- for (; cell != endc; ++cell)
- // only work on locally relevant cells
- if (cell->is_artificial() == false)
- {
- hp_fe_values.reinit(cell);
- const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
++ hp::FECollection<dim, spacedim> fe_collection(dof_handler.get_fe());
++ hp::QCollection<dim> q_coll_dummy;
- for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index)
- {
- // 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(
- Quadrature<dim> (
- fe_collection[fe_index].get_unit_support_points()));
- }
- local_dof_indices.resize(cell->get_fe().dofs_per_cell);
- cell->get_dof_indices(local_dof_indices);
++ for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index)
++ {
++ // 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(
++ Quadrature<dim> (
++ 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
- //
- // 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 =
- dof_handler.begin_active(), endc = dof_handler.end();
-
- std::vector<types::global_dof_index> local_dof_indices;
- for (; cell != endc; ++cell)
- // only work on locally relevant cells
- if (cell->is_artificial() == false)
- {
- hp_fe_values.reinit(cell);
- const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
- const std::vector<Point<spacedim> > & points =
- fe_values.get_quadrature_points();
- for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
- // insert the values into the map
- support_points[local_dof_indices[i]] = points[i];
- }
++ // 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.
++ hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe_collection,
++ q_coll_dummy, update_quadrature_points);
++ typename DH::active_cell_iterator cell =
++ dof_handler.begin_active(), endc = dof_handler.end();
++
++ std::vector<types::global_dof_index> local_dof_indices;
++ for (; cell != endc; ++cell)
++ // only work on locally relevant cells
++ if (cell->is_artificial() == false)
++ {
++ hp_fe_values.reinit(cell);
++ const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
+
- local_dof_indices.resize(cell->get_fe().dofs_per_cell);
- cell->get_dof_indices(local_dof_indices);
++ local_dof_indices.resize(cell->get_fe().dofs_per_cell);
++ cell->get_dof_indices(local_dof_indices);
+
- const std::vector<Point<spacedim> > & points =
- fe_values.get_quadrature_points();
- for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
- // insert the values into the map
- support_points[local_dof_indices[i]] = points[i];
- }
++ const std::vector<Point<spacedim> > & points =
++ fe_values.get_quadrature_points();
++ for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
++ // insert the values into the map
++ support_points[local_dof_indices[i]] = points[i];
++ }
}
map_dofs_to_support_points(const hp::MappingCollection<DH::dimension, DH::space_dimension> & mapping,
const DH &dof_handler,
std::vector<Point<DH::space_dimension> > &support_points)
-- {
-- // get the data in the form of the map as above
- std::map<types::global_dof_index,Point<DH::space_dimension> > x_support_points;
- std::map<unsigned int,Point<DH::space_dimension> > x_support_points;
-- map_dofs_to_support_points(mapping, dof_handler, x_support_points);
++ {
++ // get the data in the form of the map as above
++ std::map<types::global_dof_index,Point<DH::space_dimension> > x_support_points;
++ map_dofs_to_support_points(mapping, dof_handler, x_support_points);
-- // now convert from the map to the linear vector. make sure every
-- // entry really appeared in the map
-- for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-- {
-- Assert (x_support_points.find(i) != x_support_points.end(),
-- ExcInternalError());
-- support_points[i] = x_support_points[i];
-- }
-- }
++ // now convert from the map to the linear vector. make sure every
++ // entry really appeared in the map
++ for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
++ {
++ Assert (x_support_points.find(i) != x_support_points.end(),
++ ExcInternalError());
++ support_points[i] = x_support_points[i];
++ }
++ }
}
}
template <int dim, int spacedim>
void
map_dofs_to_support_points (const Mapping<dim,spacedim> &mapping,
-- const DoFHandler<dim,spacedim> &dof_handler,
-- std::vector<Point<spacedim> > &support_points)
++ const DoFHandler<dim,spacedim> &dof_handler,
++
std::vector<Point<spacedim> > &support_points)
{
AssertDimension(support_points.size(), dof_handler.n_dofs());
Assert ((dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
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);
const unsigned int n_components = DoFTools::n_components (dof);
- // set the component mask to either
- // the original value or a vector
- // of trues
- const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
- std::vector<bool> (n_components, true) :
- component_mask_);
- Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
- VectorTools::ExcNoComponentSelected());
+ Assert (component_mask.n_selected_components(n_components) > 0,
+ VectorTools::ExcNoComponentSelected());
// a field to store the indices
- std::vector<unsigned int> face_dofs;
+ std::vector<types::global_dof_index> face_dofs;
face_dofs.reserve (max_dofs_per_face(dof));
typename DH<dim,spacedim>::active_cell_iterator
template <class DH, class Sparsity>
void make_cell_patches(
-- Sparsity& block_list,
-- const DH& dof_handler,
-- const unsigned int level,
-- const std::vector<bool>& selected_dofs,
- types::global_dof_index offset)
- unsigned int offset)
++ Sparsity& block_list,
++ const DH& dof_handler,
++ const unsigned int level,
++ const std::vector<bool>& selected_dofs,
++ types::global_dof_index offset)
{
typename DH::cell_iterator cell;
typename DH::cell_iterator endc = dof_handler.end(level);
AssertDimension(indices.size(), selected_dofs.size());
for (unsigned int j=0;j<indices.size();++j)
-- {
-- if(selected_dofs.size() == 0)
-- block_list.add(i,indices[j]-offset);
-- else
{
-- if(selected_dofs[j])
++ if(selected_dofs.size() == 0)
block_list.add(i,indices[j]-offset);
++ else
++ {
++ if(selected_dofs[j])
++ block_list.add(i,indices[j]-offset);
++ }
}
-- }
}
}
const DH& dof_handler,
const unsigned int level,
const bool interior_only)
-- {
-- const FiniteElement<DH::dimension>& fe = dof_handler.get_fe();
-- block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level));
-- typename DH::cell_iterator cell;
-- typename DH::cell_iterator endc = dof_handler.end(level);
++ {
++ const FiniteElement<DH::dimension>& fe = dof_handler.get_fe();
++ block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level));
++ typename DH::cell_iterator cell;
++ typename DH::cell_iterator endc = dof_handler.end(level);
-- std::vector<unsigned int> indices;
-- std::vector<bool> exclude;
++ std::vector<unsigned int> indices;
++ std::vector<bool> exclude;
-- for (cell=dof_handler.begin(level); cell != endc;++cell)
++ for (cell=dof_handler.begin(level); cell != endc;++cell)
{
indices.resize(cell->get_fe().dofs_per_cell);
cell->get_mg_dof_indices(indices);
block_list.add(0, indices[j]);
}
}
-- }
++ }
template <class DH>
projects / dealii-svn.git / commitdiff