template <class DH, class SparsityPattern>
void
make_sparsity_pattern (const DH &dof,
- SparsityPattern &sparsity,
- const ConstraintMatrix &constraints,
- const bool keep_constrained_dofs,
- const types::subdomain_id subdomain_id)
+ SparsityPattern &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
+ const types::subdomain_id subdomain_id)
{
const unsigned int n_dofs = dof.n_dofs();
Assert (sparsity.n_rows() == n_dofs,
- ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
+ ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
Assert (sparsity.n_cols() == n_dofs,
- ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
+ ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
// If we have a distributed::Triangulation only
// allow locally_owned subdomain. Not setting a
||
(subdomain_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;
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
// processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
if (((subdomain_id == types::invalid_subdomain_id)
- ||
- (subdomain_id == cell->subdomain_id()))
- &&
+ ||
+ (subdomain_id == cell->subdomain_id()))
+ &&
cell->is_locally_owned())
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- dofs_on_this_cell.resize (dofs_per_cell);
- cell->get_dof_indices (dofs_on_this_cell);
+ {
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ dofs_on_this_cell.resize (dofs_per_cell);
+ cell->get_dof_indices (dofs_on_this_cell);
// make sparsity pattern for this
// cell. if no constraints pattern was
// given, then the following call acts
// as if simply no constraints existed
- constraints.add_entries_local_to_global (dofs_on_this_cell,
- sparsity,
- keep_constrained_dofs);
- }
+ constraints.add_entries_local_to_global (dofs_on_this_cell,
+ sparsity,
+ keep_constrained_dofs);
+ }
}
template <class DH, class SparsityPattern>
void
make_sparsity_pattern (const DH &dof,
- const Table<2,Coupling> &couplings,
- SparsityPattern &sparsity,
- const ConstraintMatrix &constraints,
- const bool keep_constrained_dofs,
- const types::subdomain_id subdomain_id)
+ const Table<2,Coupling> &couplings,
+ SparsityPattern &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
+ const types::subdomain_id subdomain_id)
{
const unsigned int n_dofs = dof.n_dofs();
Assert (sparsity.n_rows() == n_dofs,
- ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
+ ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
Assert (sparsity.n_cols() == n_dofs,
- ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
+ ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
Assert (couplings.n_rows() == dof.get_fe().n_components(),
- ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components()));
+ ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components()));
Assert (couplings.n_cols() == dof.get_fe().n_components(),
- ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components()));
+ ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components()));
// If we have a distributed::Triangulation only
// allow locally_owned subdomain. Not setting a
||
(subdomain_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."));
const hp::FECollection<DH::dimension,DH::space_dimension> fe_collection (dof.get_fe());
bool need_dof_mask = false;
for (unsigned int i=0; i<couplings.n_rows(); ++i)
for (unsigned int j=0; j<couplings.n_cols(); ++j)
- if (couplings(i,j) == none)
- need_dof_mask = true;
+ if (couplings(i,j) == none)
+ need_dof_mask = true;
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
+ = (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);
+ }
+ }
std::vector<unsigned int> 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
// processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
if (((subdomain_id == types::invalid_subdomain_id)
- ||
- (subdomain_id == cell->subdomain_id()))
- &&
+ ||
+ (subdomain_id == cell->subdomain_id()))
+ &&
cell->is_locally_owned())
- {
- const unsigned int fe_index = cell->active_fe_index();
- const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
+ {
+ const unsigned int fe_index = cell->active_fe_index();
+ const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
- dofs_on_this_cell.resize (dofs_per_cell);
- cell->get_dof_indices (dofs_on_this_cell);
+ dofs_on_this_cell.resize (dofs_per_cell);
+ cell->get_dof_indices (dofs_on_this_cell);
// make sparsity pattern for this
// cell. if no constraints pattern was
// given, then the following call acts
// as if simply no constraints existed
- constraints.add_entries_local_to_global (dofs_on_this_cell,
- sparsity,
- keep_constrained_dofs,
- dof_mask[fe_index]);
- }
+ constraints.add_entries_local_to_global (dofs_on_this_cell,
+ sparsity,
+ keep_constrained_dofs,
+ dof_mask[fe_index]);
+ }
}
const unsigned int n_dofs_col = dof_col.n_dofs();
Assert (sparsity.n_rows() == n_dofs_row,
- ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row));
+ ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row));
Assert (sparsity.n_cols() == n_dofs_col,
- ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col));
+ ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col));
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 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());
- }
- }
+ 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());
+ }
+ }
}
}
// indicators in 1d, so it is no
// performance problem to call the
// other function
- typename DH::FunctionMap boundary_indicators;
- boundary_indicators[0] = 0;
- boundary_indicators[1] = 0;
- make_boundary_sparsity_pattern<DH, SparsityPattern> (dof,
- boundary_indicators,
- dof_to_boundary_mapping,
- sparsity);
- return;
+ typename DH::FunctionMap boundary_indicators;
+ boundary_indicators[0] = 0;
+ boundary_indicators[1] = 0;
+ make_boundary_sparsity_pattern<DH, SparsityPattern> (dof,
+ boundary_indicators,
+ dof_to_boundary_mapping,
+ sparsity);
+ return;
}
const unsigned int n_dofs = dof.n_dofs();
#ifdef DEBUG
if (sparsity.n_rows() != 0)
{
- 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);
+ 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);
};
#endif
// and so every boundary line is
// also part of a boundary face.
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 (cell->at_boundary(f))
- {
- 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 (cell->at_boundary(f))
+ {
+ 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]]);
+ }
}
{
// first check left, then right
// boundary point
- for (unsigned int direction=0; direction<2; ++direction)
- {
+ for (unsigned int direction=0; direction<2; ++direction)
+ {
// if this boundary is not
// requested, then go on with next one
- if (boundary_indicators.find(direction) ==
- boundary_indicators.end())
- continue;
+ if (boundary_indicators.find(direction) ==
+ boundary_indicators.end())
+ continue;
// 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<unsigned int> 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);
// next get boundary mapped dof
// indices of boundary dofs
- for (unsigned int i=0; i<dofs_per_vertex; ++i)
- boundary_dof_boundary_indices[i]
- = dof_to_boundary_mapping[cell->vertex_dof_index(direction,i)];
-
- for (unsigned int i=0; i<dofs_per_vertex; ++i)
- sparsity.add_entries (boundary_dof_boundary_indices[i],
- boundary_dof_boundary_indices.begin(),
- boundary_dof_boundary_indices.end());
- };
- return;
+ for (unsigned int i=0; i<dofs_per_vertex; ++i)
+ boundary_dof_boundary_indices[i]
+ = dof_to_boundary_mapping[cell->vertex_dof_index(direction,i)];
+
+ for (unsigned int i=0; i<dofs_per_vertex; ++i)
+ sparsity.add_entries (boundary_dof_boundary_indices[i],
+ boundary_dof_boundary_indices.begin(),
+ boundary_dof_boundary_indices.end());
+ };
+ return;
}
const unsigned int n_dofs = dof.n_dofs();
AssertDimension (dof_to_boundary_mapping.size(), n_dofs);
Assert (boundary_indicators.find(numbers::internal_face_boundary_id) == boundary_indicators.end(),
- typename DH::ExcInvalidBoundaryIndicator());
+ typename DH::ExcInvalidBoundaryIndicator());
Assert (sparsity.n_rows() == dof.n_boundary_dofs (boundary_indicators),
- ExcDimensionMismatch (sparsity.n_rows(), dof.n_boundary_dofs (boundary_indicators)));
+ ExcDimensionMismatch (sparsity.n_rows(), dof.n_boundary_dofs (boundary_indicators)));
Assert (sparsity.n_cols() == dof.n_boundary_dofs (boundary_indicators),
- ExcDimensionMismatch (sparsity.n_cols(), dof.n_boundary_dofs (boundary_indicators)));
+ ExcDimensionMismatch (sparsity.n_cols(), dof.n_boundary_dofs (boundary_indicators)));
#ifdef DEBUG
if (sparsity.n_rows() != 0)
{
- 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);
+ 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);
};
#endif
std::vector<unsigned int> 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]]);
+ }
}
template <class DH, class SparsityPattern>
void
make_flux_sparsity_pattern (const DH &dof,
- SparsityPattern &sparsity,
- const ConstraintMatrix &constraints,
- const bool keep_constrained_dofs,
+ SparsityPattern &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
const types::subdomain_id subdomain_id)
{
const unsigned int n_dofs = dof.n_dofs();
||
(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;
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(),
- endc = dof.end();
+ endc = dof.end();
// TODO: in an old implementation, we used
// user flags before to tag faces that were
// processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
if (((subdomain_id == types::invalid_subdomain_id)
- ||
- (subdomain_id == cell->subdomain_id()))
+ ||
+ (subdomain_id == cell->subdomain_id()))
&&
cell->is_locally_owned())
- {
- const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell;
- dofs_on_this_cell.resize (n_dofs_on_this_cell);
- cell->get_dof_indices (dofs_on_this_cell);
+ {
+ const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell;
+ dofs_on_this_cell.resize (n_dofs_on_this_cell);
+ cell->get_dof_indices (dofs_on_this_cell);
// make sparsity pattern for this
// cell. if no constraints pattern was
// given, then the following call acts
// as if simply no constraints existed
- constraints.add_entries_local_to_global (dofs_on_this_cell,
- sparsity,
- keep_constrained_dofs);
-
- for (unsigned int face = 0;
- face < GeometryInfo<DH::dimension>::faces_per_cell;
- ++face)
- {
- typename DH::face_iterator cell_face = cell->face(face);
- if (! cell->at_boundary(face) )
- {
- typename DH::cell_iterator neighbor = cell->neighbor(face);
-
- if (cell_face->has_children())
- {
- for (unsigned int sub_nr = 0;
- sub_nr != cell_face->number_of_children();
- ++sub_nr)
- {
- const typename DH::cell_iterator
- sub_neighbor
- = cell->neighbor_child_on_subface (face, sub_nr);
-
- const unsigned int n_dofs_on_neighbor
- = sub_neighbor->get_fe().dofs_per_cell;
- dofs_on_other_cell.resize (n_dofs_on_neighbor);
- sub_neighbor->get_dof_indices (dofs_on_other_cell);
-
- constraints.add_entries_local_to_global
- (dofs_on_this_cell, dofs_on_other_cell,
- sparsity, keep_constrained_dofs);
- constraints.add_entries_local_to_global
- (dofs_on_other_cell, dofs_on_this_cell,
- sparsity, keep_constrained_dofs);
- }
- }
- else
- {
+ constraints.add_entries_local_to_global (dofs_on_this_cell,
+ sparsity,
+ keep_constrained_dofs);
+
+ for (unsigned int face = 0;
+ face < GeometryInfo<DH::dimension>::faces_per_cell;
+ ++face)
+ {
+ typename DH::face_iterator cell_face = cell->face(face);
+ if (! cell->at_boundary(face) )
+ {
+ typename DH::cell_iterator neighbor = cell->neighbor(face);
+
+ if (cell_face->has_children())
+ {
+ for (unsigned int sub_nr = 0;
+ sub_nr != cell_face->number_of_children();
+ ++sub_nr)
+ {
+ const typename DH::cell_iterator
+ sub_neighbor
+ = cell->neighbor_child_on_subface (face, sub_nr);
+
+ const unsigned int n_dofs_on_neighbor
+ = sub_neighbor->get_fe().dofs_per_cell;
+ dofs_on_other_cell.resize (n_dofs_on_neighbor);
+ sub_neighbor->get_dof_indices (dofs_on_other_cell);
+
+ constraints.add_entries_local_to_global
+ (dofs_on_this_cell, dofs_on_other_cell,
+ sparsity, keep_constrained_dofs);
+ constraints.add_entries_local_to_global
+ (dofs_on_other_cell, dofs_on_this_cell,
+ sparsity, keep_constrained_dofs);
+ }
+ }
+ else
+ {
// Refinement edges are taken care of
// by coarser cells
// the constraints when the neighbor cell is
// coarser, but only the current cell is owned
// locally!
- if (cell->neighbor_is_coarser(face))
- continue;
+ if (cell->neighbor_is_coarser(face))
+ continue;
- const unsigned int n_dofs_on_neighbor
- = neighbor->get_fe().dofs_per_cell;
- dofs_on_other_cell.resize (n_dofs_on_neighbor);
+ const unsigned int n_dofs_on_neighbor
+ = neighbor->get_fe().dofs_per_cell;
+ dofs_on_other_cell.resize (n_dofs_on_neighbor);
- neighbor->get_dof_indices (dofs_on_other_cell);
+ neighbor->get_dof_indices (dofs_on_other_cell);
- constraints.add_entries_local_to_global
- (dofs_on_this_cell, dofs_on_other_cell,
- sparsity, keep_constrained_dofs);
+ constraints.add_entries_local_to_global
+ (dofs_on_this_cell, dofs_on_other_cell,
+ sparsity, keep_constrained_dofs);
// only need to add these in case the neighbor
// cell is not locally owned - otherwise, we
// touch each face twice and hence put the
// indices the other way around
- if (cell->neighbor(face)->subdomain_id() !=
- cell->subdomain_id())
- constraints.add_entries_local_to_global
- (dofs_on_other_cell, dofs_on_this_cell,
- sparsity, keep_constrained_dofs);
- }
- }
- }
- }
+ if (cell->neighbor(face)->subdomain_id() !=
+ cell->subdomain_id())
+ constraints.add_entries_local_to_global
+ (dofs_on_other_cell, dofs_on_this_cell,
+ sparsity, keep_constrained_dofs);
+ }
+ }
+ }
+ }
}
const Table<2,Coupling> &component_couplings)
{
Assert(component_couplings.n_rows() == fe.n_components(),
- ExcDimensionMismatch(component_couplings.n_rows(),
- fe.n_components()));
+ ExcDimensionMismatch(component_couplings.n_rows(),
+ fe.n_components()));
Assert(component_couplings.n_cols() == fe.n_components(),
- ExcDimensionMismatch(component_couplings.n_cols(),
- fe.n_components()));
+ ExcDimensionMismatch(component_couplings.n_cols(),
+ fe.n_components()));
const unsigned int n_dofs = fe.dofs_per_cell;
for (unsigned int i=0; i<n_dofs; ++i)
{
- const unsigned int ii
- = (fe.is_primitive(i) ?
- fe.system_to_component_index(i).first
- :
- (std::find (fe.get_nonzero_components(i).begin(),
- fe.get_nonzero_components(i).end(),
- true)
- -
- fe.get_nonzero_components(i).begin())
- );
- Assert (ii < fe.n_components(), ExcInternalError());
-
- for (unsigned int j=0; j<n_dofs; ++j)
- {
- const unsigned int jj
- = (fe.is_primitive(j) ?
- fe.system_to_component_index(j).first
- :
- (std::find (fe.get_nonzero_components(j).begin(),
- fe.get_nonzero_components(j).end(),
- true)
- -
- fe.get_nonzero_components(j).begin())
- );
- Assert (jj < fe.n_components(), ExcInternalError());
-
- dof_couplings(i,j) = component_couplings(ii,jj);
- }
+ const unsigned int ii
+ = (fe.is_primitive(i) ?
+ fe.system_to_component_index(i).first
+ :
+ (std::find (fe.get_nonzero_components(i).begin(),
+ fe.get_nonzero_components(i).end(),
+ true)
+ -
+ fe.get_nonzero_components(i).begin())
+ );
+ Assert (ii < fe.n_components(), ExcInternalError());
+
+ for (unsigned int j=0; j<n_dofs; ++j)
+ {
+ const unsigned int jj
+ = (fe.is_primitive(j) ?
+ fe.system_to_component_index(j).first
+ :
+ (std::find (fe.get_nonzero_components(j).begin(),
+ fe.get_nonzero_components(j).end(),
+ true)
+ -
+ fe.get_nonzero_components(j).begin())
+ );
+ Assert (jj < fe.n_components(), ExcInternalError());
+
+ dof_couplings(i,j) = component_couplings(ii,jj);
+ }
}
return dof_couplings;
}
std::vector<Table<2,Coupling> > return_value (fe.size());
for (unsigned int i=0; i<fe.size(); ++i)
return_value[i]
- = dof_couplings_from_component_couplings(fe[i], component_couplings);
+ = dof_couplings_from_component_couplings(fe[i], component_couplings);
return return_value;
}
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<unsigned int> dofs_on_this_cell(fe.dofs_per_cell);
- std::vector<unsigned int> 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);
- 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 ();
+ }
+ }
+ }
}
}
template <int dim, int spacedim, class SparsityPattern>
void
make_flux_sparsity_pattern (const dealii::hp::DoFHandler<dim,spacedim> &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)
{
// while the implementation above is
// quite optimized and caches a lot of
// 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<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<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));
- 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 ();
+ }
+ }
+ }
+ }
}
}
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)
{
// do the error checking and frame code
// here, and then pass on to more
const unsigned int n_comp = dof.get_fe().n_components();
Assert (sparsity.n_rows() == n_dofs,
- ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
+ ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
Assert (sparsity.n_cols() == n_dofs,
- ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
+ ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
Assert (int_mask.n_rows() == n_comp,
- ExcDimensionMismatch (int_mask.n_rows(), n_comp));
+ ExcDimensionMismatch (int_mask.n_rows(), n_comp));
Assert (int_mask.n_cols() == n_comp,
- ExcDimensionMismatch (int_mask.n_cols(), n_comp));
+ ExcDimensionMismatch (int_mask.n_cols(), n_comp));
Assert (flux_mask.n_rows() == n_comp,
- ExcDimensionMismatch (flux_mask.n_rows(), n_comp));
+ ExcDimensionMismatch (flux_mask.n_rows(), n_comp));
Assert (flux_mask.n_cols() == n_comp,
- ExcDimensionMismatch (flux_mask.n_cols(), n_comp));
+ ExcDimensionMismatch (flux_mask.n_cols(), n_comp));
// Clear user flags because we will
const_cast<Triangulation<DH::dimension,DH::space_dimension> &>(dof.get_tria()).clear_user_flags ();
internal::make_flux_sparsity_pattern (dof, sparsity,
- int_mask, flux_mask);
+ int_mask, flux_mask);
// finally restore the user flags
const_cast<Triangulation<DH::dimension,DH::space_dimension> &>(dof.get_tria()).load_user_flags(user_flags);
inline
bool
check_master_dof_list (const FullMatrix<double> &face_interpolation_matrix,
- const std::vector<unsigned int> &master_dof_list)
+ const std::vector<unsigned int> &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
// enough, or whether we have to
// fear that the matrix is not
// regular
- double diagonal_sum = 0;
- for (unsigned int i=0; i<N; ++i)
- diagonal_sum += std::fabs(tmp(i,i));
- const double typical_diagonal_element = diagonal_sum/N;
+ double diagonal_sum = 0;
+ for (unsigned int i=0; i<N; ++i)
+ diagonal_sum += std::fabs(tmp(i,i));
+ const double typical_diagonal_element = diagonal_sum/N;
// initialize the array that holds
// the permutations that we find
// during pivot search
- std::vector<unsigned int> p(N);
- for (unsigned int i=0; i<N; ++i)
- p[i] = i;
+ std::vector<unsigned int> p(N);
+ for (unsigned int i=0; i<N; ++i)
+ p[i] = i;
- for (unsigned int j=0; j<N; ++j)
- {
+ for (unsigned int j=0; j<N; ++j)
+ {
// pivot search: search that
// part of the line on and
// right of the diagonal for
// the largest element
- double max = std::fabs(tmp(j,j));
- unsigned int r = j;
- for (unsigned int i=j+1; i<N; ++i)
- {
- if (std::fabs(tmp(i,j)) > max)
- {
- max = std::fabs(tmp(i,j));
- r = i;
- }
- }
+ double max = std::fabs(tmp(j,j));
+ unsigned int r = j;
+ for (unsigned int i=j+1; i<N; ++i)
+ {
+ if (std::fabs(tmp(i,j)) > max)
+ {
+ max = std::fabs(tmp(i,j));
+ r = i;
+ }
+ }
// check whether the
// pivot is too small. if
// that is the case, then
// and we shouldn't use
// this set of master
// dofs
- if (max < 1.e-12*typical_diagonal_element)
- return false;
+ if (max < 1.e-12*typical_diagonal_element)
+ return false;
// row interchange
- if (r>j)
- {
- for (unsigned int k=0; k<N; ++k)
- std::swap (tmp(j,k), tmp(r,k));
+ if (r>j)
+ {
+ for (unsigned int k=0; k<N; ++k)
+ std::swap (tmp(j,k), tmp(r,k));
- std::swap (p[j], p[r]);
- }
+ std::swap (p[j], p[r]);
+ }
// transformation
- const double hr = 1./tmp(j,j);
- tmp(j,j) = hr;
- for (unsigned int k=0; k<N; ++k)
- {
- if (k==j) continue;
- for (unsigned int i=0; i<N; ++i)
- {
- if (i==j) continue;
- tmp(i,k) -= tmp(i,j)*tmp(j,k)*hr;
- }
- }
- for (unsigned int i=0; i<N; ++i)
- {
- tmp(i,j) *= hr;
- tmp(j,i) *= -hr;
- }
- tmp(j,j) = hr;
- }
+ const double hr = 1./tmp(j,j);
+ tmp(j,j) = hr;
+ for (unsigned int k=0; k<N; ++k)
+ {
+ if (k==j) continue;
+ for (unsigned int i=0; i<N; ++i)
+ {
+ if (i==j) continue;
+ tmp(i,k) -= tmp(i,j)*tmp(j,k)*hr;
+ }
+ }
+ for (unsigned int i=0; i<N; ++i)
+ {
+ tmp(i,j) *= hr;
+ tmp(j,i) *= -hr;
+ }
+ tmp(j,j) = hr;
+ }
// everything went fine, so
// we can accept this set of
// master dofs (at least as
// far as they have already
// been collected)
- return true;
+ return true;
}
template <int dim, int spacedim>
void
select_master_dofs_for_face_restriction (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- const FullMatrix<double> &face_interpolation_matrix,
- std::vector<bool> &master_dof_mask)
+ const FiniteElement<dim,spacedim> &fe2,
+ const FullMatrix<double> &face_interpolation_matrix,
+ std::vector<bool> &master_dof_mask)
{
- Assert (fe1.dofs_per_face >= fe2.dofs_per_face,
- ExcInternalError());
- AssertDimension (master_dof_mask.size(), fe1.dofs_per_face);
-
- Assert (fe2.dofs_per_vertex <= fe1.dofs_per_vertex,
- ExcInternalError());
- Assert (fe2.dofs_per_line <= fe1.dofs_per_line,
- ExcInternalError());
- Assert ((dim < 3)
- ||
- (fe2.dofs_per_quad <= fe1.dofs_per_quad),
- ExcInternalError());
+ Assert (fe1.dofs_per_face >= fe2.dofs_per_face,
+ ExcInternalError());
+ AssertDimension (master_dof_mask.size(), fe1.dofs_per_face);
+
+ Assert (fe2.dofs_per_vertex <= fe1.dofs_per_vertex,
+ ExcInternalError());
+ Assert (fe2.dofs_per_line <= fe1.dofs_per_line,
+ ExcInternalError());
+ Assert ((dim < 3)
+ ||
+ (fe2.dofs_per_quad <= fe1.dofs_per_quad),
+ ExcInternalError());
// the idea here is to designate as
// many DoFs in fe1 per object
// the other one. this latter case
// shows up when running
// hp/hp_constraints_q_system_06
- std::vector<unsigned int> master_dof_list;
- unsigned int index = 0;
- for (int v=0;
- 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<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)
+ {
// make sure that we
// were able to find
// a set of master
// code down below
// didn't just reject
// all our efforts
- Assert (i < fe1.dofs_per_vertex,
- ExcInternalError());
+ Assert (i < fe1.dofs_per_vertex,
+ ExcInternalError());
// tentatively push
// this vertex dof
- master_dof_list.push_back (index+i);
+ master_dof_list.push_back (index+i);
// then see what
// happens. if it
// succeeds, fine
- if (check_master_dof_list (face_interpolation_matrix,
- master_dof_list)
- == true)
- ++dofs_added;
- else
+ if (check_master_dof_list (face_interpolation_matrix,
+ master_dof_list)
+ == true)
+ ++dofs_added;
+ else
// well, it
// didn't. simply
// pop that dof
// again and try
// with the next
// dof
- master_dof_list.pop_back ();
+ master_dof_list.pop_back ();
// forward counter by
// one
- ++i;
- }
- index += fe1.dofs_per_vertex;
- }
-
- for (int l=0;
- l<static_cast<signed int>(GeometryInfo<dim>::lines_per_face);
- ++l)
- {
+ ++i;
+ }
+ index += fe1.dofs_per_vertex;
+ }
+
+ for (int l=0;
+ l<static_cast<signed int>(GeometryInfo<dim>::lines_per_face);
+ ++l)
+ {
// same algorithm as above
- unsigned int dofs_added = 0;
- unsigned int i = 0;
- while (dofs_added < fe2.dofs_per_line)
- {
- Assert (i < fe1.dofs_per_line,
- ExcInternalError());
-
- master_dof_list.push_back (index+i);
- if (check_master_dof_list (face_interpolation_matrix,
- master_dof_list)
- == true)
- ++dofs_added;
- else
- master_dof_list.pop_back ();
-
- ++i;
- }
- index += fe1.dofs_per_line;
- }
-
- for (int q=0;
- q<static_cast<signed int>(GeometryInfo<dim>::quads_per_face);
- ++q)
- {
+ unsigned int dofs_added = 0;
+ unsigned int i = 0;
+ while (dofs_added < fe2.dofs_per_line)
+ {
+ Assert (i < fe1.dofs_per_line,
+ ExcInternalError());
+
+ master_dof_list.push_back (index+i);
+ if (check_master_dof_list (face_interpolation_matrix,
+ master_dof_list)
+ == true)
+ ++dofs_added;
+ else
+ master_dof_list.pop_back ();
+
+ ++i;
+ }
+ index += fe1.dofs_per_line;
+ }
+
+ for (int q=0;
+ q<static_cast<signed int>(GeometryInfo<dim>::quads_per_face);
+ ++q)
+ {
// same algorithm as above
- unsigned int dofs_added = 0;
- unsigned int i = 0;
- while (dofs_added < fe2.dofs_per_quad)
- {
- Assert (i < fe1.dofs_per_quad,
- ExcInternalError());
-
- master_dof_list.push_back (index+i);
- if (check_master_dof_list (face_interpolation_matrix,
- master_dof_list)
- == true)
- ++dofs_added;
- else
- master_dof_list.pop_back ();
-
- ++i;
- }
- index += fe1.dofs_per_quad;
- }
-
- AssertDimension (index, fe1.dofs_per_face);
- AssertDimension (master_dof_list.size(), fe2.dofs_per_face);
+ unsigned int dofs_added = 0;
+ unsigned int i = 0;
+ while (dofs_added < fe2.dofs_per_quad)
+ {
+ Assert (i < fe1.dofs_per_quad,
+ ExcInternalError());
+
+ master_dof_list.push_back (index+i);
+ if (check_master_dof_list (face_interpolation_matrix,
+ master_dof_list)
+ == true)
+ ++dofs_added;
+ else
+ master_dof_list.pop_back ();
+
+ ++i;
+ }
+ index += fe1.dofs_per_quad;
+ }
+
+ AssertDimension (index, fe1.dofs_per_face);
+ AssertDimension (master_dof_list.size(), fe2.dofs_per_face);
// finally copy the list into the
// mask
- std::fill (master_dof_mask.begin(), master_dof_mask.end(), false);
- for (std::vector<unsigned int>::const_iterator i=master_dof_list.begin();
- i!=master_dof_list.end(); ++i)
- 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;
}
template <int dim, int spacedim>
void
ensure_existence_of_master_dof_mask (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- const FullMatrix<double> &face_interpolation_matrix,
- std_cxx1x::shared_ptr<std::vector<bool> > &master_dof_mask)
+ const FiniteElement<dim,spacedim> &fe2,
+ const FullMatrix<double> &face_interpolation_matrix,
+ std_cxx1x::shared_ptr<std::vector<bool> > &master_dof_mask)
{
- if (master_dof_mask == std_cxx1x::shared_ptr<std::vector<bool> >())
- {
- master_dof_mask = std_cxx1x::shared_ptr<std::vector<bool> >
- (new std::vector<bool> (fe1.dofs_per_face));
- select_master_dofs_for_face_restriction (fe1,
- fe2,
- face_interpolation_matrix,
- *master_dof_mask);
- }
+ if (master_dof_mask == std_cxx1x::shared_ptr<std::vector<bool> >())
+ {
+ master_dof_mask = std_cxx1x::shared_ptr<std::vector<bool> >
+ (new std::vector<bool> (fe1.dofs_per_face));
+ select_master_dofs_for_face_restriction (fe1,
+ fe2,
+ face_interpolation_matrix,
+ *master_dof_mask);
+ }
}
template <int dim, int spacedim>
void
ensure_existence_of_face_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- std_cxx1x::shared_ptr<FullMatrix<double> > &matrix)
+ const FiniteElement<dim,spacedim> &fe2,
+ std_cxx1x::shared_ptr<FullMatrix<double> > &matrix)
{
- if (matrix == std_cxx1x::shared_ptr<FullMatrix<double> >())
- {
- matrix = std_cxx1x::shared_ptr<FullMatrix<double> >
- (new FullMatrix<double> (fe2.dofs_per_face,
- fe1.dofs_per_face));
- fe1.get_face_interpolation_matrix (fe2,
- *matrix);
- }
+ if (matrix == std_cxx1x::shared_ptr<FullMatrix<double> >())
+ {
+ matrix = std_cxx1x::shared_ptr<FullMatrix<double> >
+ (new FullMatrix<double> (fe2.dofs_per_face,
+ fe1.dofs_per_face));
+ fe1.get_face_interpolation_matrix (fe2,
+ *matrix);
+ }
}
template <int dim, int spacedim>
void
ensure_existence_of_subface_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- const unsigned int subface,
- std_cxx1x::shared_ptr<FullMatrix<double> > &matrix)
+ const FiniteElement<dim,spacedim> &fe2,
+ const unsigned int subface,
+ std_cxx1x::shared_ptr<FullMatrix<double> > &matrix)
{
- if (matrix == std_cxx1x::shared_ptr<FullMatrix<double> >())
- {
- matrix = std_cxx1x::shared_ptr<FullMatrix<double> >
- (new FullMatrix<double> (fe2.dofs_per_face,
- fe1.dofs_per_face));
- fe1.get_subface_interpolation_matrix (fe2,
- subface,
- *matrix);
- }
+ if (matrix == std_cxx1x::shared_ptr<FullMatrix<double> >())
+ {
+ matrix = std_cxx1x::shared_ptr<FullMatrix<double> >
+ (new FullMatrix<double> (fe2.dofs_per_face,
+ fe1.dofs_per_face));
+ fe1.get_subface_interpolation_matrix (fe2,
+ subface,
+ *matrix);
+ }
}
*/
void
ensure_existence_of_split_face_matrix (const FullMatrix<double> &face_interpolation_matrix,
- const std::vector<bool> &master_dof_mask,
- std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > > &split_matrix)
+ const std::vector<bool> &master_dof_mask,
+ std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > > &split_matrix)
{
- AssertDimension (master_dof_mask.size(), face_interpolation_matrix.m());
- Assert (std::count (master_dof_mask.begin(), master_dof_mask.end(), true) ==
- static_cast<signed int>(face_interpolation_matrix.n()),
- ExcInternalError());
+ AssertDimension (master_dof_mask.size(), face_interpolation_matrix.m());
+ Assert (std::count (master_dof_mask.begin(), master_dof_mask.end(), true) ==
+ static_cast<signed int>(face_interpolation_matrix.n()),
+ ExcInternalError());
- if (split_matrix ==
- std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > >())
- {
- split_matrix
- = std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > >
- (new std::pair<FullMatrix<double>,FullMatrix<double> >());
+ if (split_matrix ==
+ std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > >())
+ {
+ split_matrix
+ = std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > >
+ (new std::pair<FullMatrix<double>,FullMatrix<double> >());
- const unsigned int n_master_dofs = face_interpolation_matrix.n();
- const unsigned int n_dofs = face_interpolation_matrix.m();
+ const unsigned int n_master_dofs = face_interpolation_matrix.n();
+ const unsigned int n_dofs = face_interpolation_matrix.m();
- Assert (n_master_dofs <= n_dofs, ExcInternalError());
+ Assert (n_master_dofs <= n_dofs, ExcInternalError());
// copy and invert the master
// component, copy the slave
// component
- split_matrix->first.reinit (n_master_dofs, n_master_dofs);
- split_matrix->second.reinit (n_dofs-n_master_dofs, n_master_dofs);
-
- unsigned int nth_master_dof = 0,
- nth_slave_dof = 0;
-
- for (unsigned int i=0; i<n_dofs; ++i)
- if (master_dof_mask[i] == true)
- {
- for (unsigned int j=0; j<n_master_dofs; ++j)
- split_matrix->first(nth_master_dof,j)
- = face_interpolation_matrix(i,j);
- ++nth_master_dof;
- }
- else
- {
- for (unsigned int j=0; j<n_master_dofs; ++j)
- split_matrix->second(nth_slave_dof,j)
- = face_interpolation_matrix(i,j);
- ++nth_slave_dof;
- }
-
- AssertDimension (nth_master_dof, n_master_dofs);
- AssertDimension (nth_slave_dof, n_dofs-n_master_dofs);
+ split_matrix->first.reinit (n_master_dofs, n_master_dofs);
+ split_matrix->second.reinit (n_dofs-n_master_dofs, n_master_dofs);
+
+ unsigned int nth_master_dof = 0,
+ nth_slave_dof = 0;
+
+ for (unsigned int i=0; i<n_dofs; ++i)
+ if (master_dof_mask[i] == true)
+ {
+ for (unsigned int j=0; j<n_master_dofs; ++j)
+ split_matrix->first(nth_master_dof,j)
+ = face_interpolation_matrix(i,j);
+ ++nth_master_dof;
+ }
+ else
+ {
+ for (unsigned int j=0; j<n_master_dofs; ++j)
+ split_matrix->second(nth_slave_dof,j)
+ = face_interpolation_matrix(i,j);
+ ++nth_slave_dof;
+ }
+
+ AssertDimension (nth_master_dof, n_master_dofs);
+ AssertDimension (nth_slave_dof, n_dofs-n_master_dofs);
//TODO[WB]: We should make sure very small entries are removed after inversion
- split_matrix->first.gauss_jordan ();
- }
+ split_matrix->first.gauss_jordan ();
+ }
}
template <typename>
struct DoFHandlerSupportsDifferentFEs
{
- static const bool value = true;
+ static const bool value = true;
};
template <int dim, int spacedim>
struct DoFHandlerSupportsDifferentFEs< dealii::DoFHandler<dim,spacedim> >
{
- static const bool value = false;
+ static const bool value = false;
};
template <int dim, int spacedim>
struct DoFHandlerSupportsDifferentFEs< dealii::MGDoFHandler<dim,spacedim> >
{
- static const bool value = false;
+ static const bool value = false;
};
unsigned int
n_finite_elements (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler)
{
- return dof_handler.get_fe().size();
+ return dof_handler.get_fe().size();
}
unsigned int
n_finite_elements (const DH &)
{
- return 1;
+ return 1;
}
unsigned int
get_most_dominating_subface_fe_index (const face_iterator &face)
{
- const unsigned int dim
- = face_iterator::AccessorType::dimension;
- const unsigned int spacedim
- = face_iterator::AccessorType::space_dimension;
-
- unsigned int dominating_subface_no = 0;
- for (; dominating_subface_no<face->n_children();
- ++dominating_subface_no)
- {
+ const unsigned int dim
+ = face_iterator::AccessorType::dimension;
+ const unsigned int spacedim
+ = face_iterator::AccessorType::space_dimension;
+
+ unsigned int dominating_subface_no = 0;
+ for (; dominating_subface_no<face->n_children();
+ ++dominating_subface_no)
+ {
// each of the subfaces
// can have only a single
// fe_index associated
// with them, since there
// is no cell on the
// other side
- Assert (face->child(dominating_subface_no)
- ->n_active_fe_indices()
- == 1,
- ExcInternalError());
-
- const FiniteElement<dim,spacedim> &
- this_subface_fe = (face->child(dominating_subface_no)
- ->get_fe (face->child(dominating_subface_no)
- ->nth_active_fe_index(0)));
-
- FiniteElementDomination::Domination
- domination = FiniteElementDomination::either_element_can_dominate;
- for (unsigned int sf=0; sf<face->n_children(); ++sf)
- if (sf != dominating_subface_no)
- {
- const FiniteElement<dim,spacedim> &
- that_subface_fe = (face->child(sf)
- ->get_fe (face->child(sf)
- ->nth_active_fe_index(0)));
-
- domination = domination &
- this_subface_fe.compare_for_face_domination(that_subface_fe);
- }
+ Assert (face->child(dominating_subface_no)
+ ->n_active_fe_indices()
+ == 1,
+ ExcInternalError());
+
+ const FiniteElement<dim,spacedim> &
+ this_subface_fe = (face->child(dominating_subface_no)
+ ->get_fe (face->child(dominating_subface_no)
+ ->nth_active_fe_index(0)));
+
+ FiniteElementDomination::Domination
+ domination = FiniteElementDomination::either_element_can_dominate;
+ for (unsigned int sf=0; sf<face->n_children(); ++sf)
+ if (sf != dominating_subface_no)
+ {
+ const FiniteElement<dim,spacedim> &
+ that_subface_fe = (face->child(sf)
+ ->get_fe (face->child(sf)
+ ->nth_active_fe_index(0)));
+
+ domination = domination &
+ this_subface_fe.compare_for_face_domination(that_subface_fe);
+ }
// see if the element
// on this subface is
// the ones on all
// other subfaces,
// and if so take it
- if ((domination == FiniteElementDomination::this_element_dominates)
- ||
- (domination == FiniteElementDomination::either_element_can_dominate))
- break;
- }
+ if ((domination == FiniteElementDomination::this_element_dominates)
+ ||
+ (domination == FiniteElementDomination::either_element_can_dominate))
+ break;
+ }
// check that we have
// found one such subface
- Assert (dominating_subface_no < face->n_children(),
- ExcNotImplemented());
+ Assert (dominating_subface_no < face->n_children(),
+ ExcNotImplemented());
// return the finite element
// index used on it. note
// that only a single fe can
// be active on such subfaces
- return face->child (dominating_subface_no)->nth_active_fe_index(0);
+ return face->child (dominating_subface_no)->nth_active_fe_index(0);
}
*/
void
filter_constraints (const std::vector<unsigned int> &master_dofs,
- const std::vector<unsigned int> &slave_dofs,
- const FullMatrix<double> &face_constraints,
- ConstraintMatrix &constraints)
+ const std::vector<unsigned int> &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
// in parallel distributed
// mode
- for (unsigned int row=0; row!=n_slave_dofs; ++row)
- Assert (slave_dofs[row] != numbers::invalid_unsigned_int,
- ExcInternalError());
- for (unsigned int col=0; col!=n_master_dofs; ++col)
- Assert (master_dofs[col] != numbers::invalid_unsigned_int,
- ExcInternalError());
+ for (unsigned int row=0; row!=n_slave_dofs; ++row)
+ Assert (slave_dofs[row] != numbers::invalid_unsigned_int,
+ ExcInternalError());
+ for (unsigned int col=0; col!=n_master_dofs; ++col)
+ Assert (master_dofs[col] != numbers::invalid_unsigned_int,
+ ExcInternalError());
- for (unsigned int row=0; row!=n_slave_dofs; ++row)
- if (constraints.is_constrained (slave_dofs[row]) == false)
- {
- bool constraint_already_satisfied = false;
+ for (unsigned int row=0; row!=n_slave_dofs; ++row)
+ if (constraints.is_constrained (slave_dofs[row]) == false)
+ {
+ bool constraint_already_satisfied = false;
// Check if we have an identity
// constraint, which is already
// satisfied by unification of
// the corresponding global dof
// indices
- for (unsigned int i=0; i<n_master_dofs; ++i)
- if (face_constraints (row,i) == 1.0)
- if (master_dofs[i] == slave_dofs[row])
- {
- constraint_already_satisfied = true;
- break;
- }
-
- if (constraint_already_satisfied == false)
- {
+ for (unsigned int i=0; i<n_master_dofs; ++i)
+ if (face_constraints (row,i) == 1.0)
+ if (master_dofs[i] == slave_dofs[row])
+ {
+ constraint_already_satisfied = true;
+ break;
+ }
+
+ if (constraint_already_satisfied == false)
+ {
// add up the absolute
// values of all
// constraints in this line
// to get a measure of
// their absolute size
- double abs_sum = 0;
- for (unsigned int i=0; i<n_master_dofs; ++i)
- abs_sum += std::abs (face_constraints(row,i));
+ double abs_sum = 0;
+ for (unsigned int i=0; i<n_master_dofs; ++i)
+ abs_sum += std::abs (face_constraints(row,i));
// then enter those
// constraints that are
// than necessary without
// producing any
// significant effect
- constraints.add_line (slave_dofs[row]);
- for (unsigned int i=0; i<n_master_dofs; ++i)
- if ((face_constraints(row,i) != 0)
- &&
- (std::fabs(face_constraints(row,i)) >= 1e-14*abs_sum))
- constraints.add_entry (slave_dofs[row],
- master_dofs[i],
- face_constraints (row,i));
- constraints.set_inhomogeneity (slave_dofs[row], 0.);
- }
- }
+ constraints.add_line (slave_dofs[row]);
+ for (unsigned int i=0; i<n_master_dofs; ++i)
+ if ((face_constraints(row,i) != 0)
+ &&
+ (std::fabs(face_constraints(row,i)) >= 1e-14*abs_sum))
+ constraints.add_entry (slave_dofs[row],
+ master_dofs[i],
+ face_constraints (row,i));
+ constraints.set_inhomogeneity (slave_dofs[row], 0.);
+ }
+ }
}
}
void
make_hp_hanging_node_constraints (const dealii::DoFHandler<1> &,
- ConstraintMatrix &)
+ ConstraintMatrix &)
{
// nothing to do for regular
// dof handlers in 1d
void
make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1> &,
- ConstraintMatrix &,
- dealii::internal::int2type<1>)
+ ConstraintMatrix &,
+ dealii::internal::int2type<1>)
{
// nothing to do for regular
// dof handlers in 1d
void
make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &,
- ConstraintMatrix &)
+ ConstraintMatrix &)
{
// nothing to do for regular
// dof handlers in 1d
void
make_oldstyle_hanging_node_constraints (const dealii::MGDoFHandler<1> &,
- ConstraintMatrix &,
- dealii::internal::int2type<1>)
+ ConstraintMatrix &,
+ dealii::internal::int2type<1>)
{
// nothing to do for regular
// dof handlers in 1d
void
make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/,
- ConstraintMatrix &/*constraints*/)
+ ConstraintMatrix &/*constraints*/)
{
// we may have to compute
// constraints for
void
make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/,
- ConstraintMatrix &/*constraints*/,
- dealii::internal::int2type<1>)
+ ConstraintMatrix &/*constraints*/,
+ dealii::internal::int2type<1>)
{
// we may have to compute
// constraints for
void
make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &,
- ConstraintMatrix &)
+ ConstraintMatrix &)
{
// nothing to do for regular
// dof handlers in 1d
void
make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &,
- ConstraintMatrix &,
- dealii::internal::int2type<1>)
+ ConstraintMatrix &,
+ dealii::internal::int2type<1>)
{
// nothing to do for regular
// dof handlers in 1d
template <class DH>
void
make_oldstyle_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints,
- dealii::internal::int2type<2>)
+ ConstraintMatrix &constraints,
+ dealii::internal::int2type<2>)
{
const unsigned int dim = 2;
// refined, we can only visit each
// face with hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
// artificial cells can at
// best neighbor ghost cells,
// but we're not interested
// in these interfaces
- if (!cell->is_artificial ())
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->has_children())
- {
+ if (!cell->is_artificial ())
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->has_children())
+ {
// in any case, faces
// can have at most two
// active fe indices,
// children can have
// only one as
// well. check this
- Assert (cell->face(face)->n_active_fe_indices() == 1,
- ExcInternalError());
- Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
- == true,
- ExcInternalError());
- for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
- if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
- Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
- ExcInternalError());
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcInternalError());
+ Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+ == true,
+ ExcInternalError());
+ for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+ if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
+ Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
+ ExcInternalError());
// right now, all that
// is implemented is
// the case that both
// sides use the same
// fe
- for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
- if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
- Assert (cell->face(face)->child(c)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcNotImplemented());
+ for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+ if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
+ Assert (cell->face(face)->child(c)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
// ok, start up the work
- const FiniteElement<dim,spacedim> &fe = cell->get_fe();
- const unsigned int fe_index = cell->active_fe_index();
+ const FiniteElement<dim,spacedim> &fe = cell->get_fe();
+ const unsigned int fe_index = cell->active_fe_index();
- const unsigned int
- n_dofs_on_mother = 2*fe.dofs_per_vertex + fe.dofs_per_line,
- n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line;
+ const unsigned int
+ n_dofs_on_mother = 2*fe.dofs_per_vertex + fe.dofs_per_line,
+ n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line;
- dofs_on_mother.resize (n_dofs_on_mother);
- dofs_on_children.resize (n_dofs_on_children);
+ dofs_on_mother.resize (n_dofs_on_mother);
+ dofs_on_children.resize (n_dofs_on_children);
- Assert(n_dofs_on_mother == fe.constraints().n(),
- ExcDimensionMismatch(n_dofs_on_mother,
- fe.constraints().n()));
- Assert(n_dofs_on_children == fe.constraints().m(),
- ExcDimensionMismatch(n_dofs_on_children,
- fe.constraints().m()));
+ Assert(n_dofs_on_mother == fe.constraints().n(),
+ ExcDimensionMismatch(n_dofs_on_mother,
+ fe.constraints().n()));
+ Assert(n_dofs_on_children == fe.constraints().m(),
+ ExcDimensionMismatch(n_dofs_on_children,
+ fe.constraints().m()));
- const typename DH::line_iterator this_face = cell->face(face);
+ const typename DH::line_iterator this_face = cell->face(face);
// fill the dofs indices. Use same
// enumeration scheme as in
// @p{FiniteElement::constraints()}
- unsigned int next_index = 0;
- for (unsigned int vertex=0; vertex<2; ++vertex)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
- fe_index);
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
- AssertDimension (next_index, dofs_on_mother.size());
-
- next_index = 0;
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(0)->vertex_dof_index(1,dof,fe_index);
- for (unsigned int child=0; child<2; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(child)->dof_index(dof, fe_index);
- AssertDimension (next_index, dofs_on_children.size());
+ unsigned int next_index = 0;
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
+ fe_index);
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
+ AssertDimension (next_index, dofs_on_mother.size());
+
+ next_index = 0;
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->vertex_dof_index(1,dof,fe_index);
+ for (unsigned int child=0; child<2; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(child)->dof_index(dof, fe_index);
+ AssertDimension (next_index, dofs_on_children.size());
// for each row in the constraint
// matrix for this line:
- for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
- {
- constraints.add_line (dofs_on_children[row]);
- for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
- constraints.add_entry (dofs_on_children[row],
- dofs_on_mother[i],
- fe.constraints()(row,i));
-
- constraints.set_inhomogeneity (dofs_on_children[row], 0.);
- }
- }
- else
- {
+ for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
+ {
+ constraints.add_line (dofs_on_children[row]);
+ for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
+ constraints.add_entry (dofs_on_children[row],
+ dofs_on_mother[i],
+ fe.constraints()(row,i));
+
+ constraints.set_inhomogeneity (dofs_on_children[row], 0.);
+ }
+ }
+ else
+ {
// this face has no
// children, but it
// could still be
// case that the
// neighbor is an
// artificial cell
- if (!cell->at_boundary(face) &&
- !cell->neighbor(face)->is_artificial())
- {
- Assert (cell->face(face)->n_active_fe_indices() == 1,
- ExcNotImplemented());
- Assert (cell->face(face)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcInternalError());
- }
- }
+ if (!cell->at_boundary(face) &&
+ !cell->neighbor(face)->is_artificial())
+ {
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcInternalError());
+ }
+ }
}
template <class DH>
void
make_oldstyle_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints,
- dealii::internal::int2type<3>)
+ ConstraintMatrix &constraints,
+ dealii::internal::int2type<3>)
{
const unsigned int dim = 3;
// refined, we can only visit each
// face with hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
// artificial cells can at
// best neighbor ghost cells,
// but we're not interested
// in these interfaces
- if (!cell->is_artificial ())
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->has_children())
- {
+ if (!cell->is_artificial ())
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->has_children())
+ {
// first of all, make sure that
// we treat a case which is
// possible, i.e. either no dofs
// on the face at all or no
// anisotropic refinement
- if (cell->get_fe().dofs_per_face == 0)
- continue;
+ if (cell->get_fe().dofs_per_face == 0)
+ continue;
- Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
- ExcNotImplemented());
+ Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
+ ExcNotImplemented());
// in any case, faces
// can have at most two
// children can have
// only one as
// well. check this
- AssertDimension (cell->face(face)->n_active_fe_indices(), 1);
- Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
- == true,
- ExcInternalError());
- for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
- AssertDimension (cell->face(face)->child(c)->n_active_fe_indices(), 1);
+ AssertDimension (cell->face(face)->n_active_fe_indices(), 1);
+ Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+ == true,
+ ExcInternalError());
+ for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+ AssertDimension (cell->face(face)->child(c)->n_active_fe_indices(), 1);
// right now, all that
// is implemented is
// this face and the
// children have the
// same fe
- for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
- if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
- {
- Assert (cell->face(face)->child(c)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcNotImplemented());
- for (unsigned int e=0; e<4; ++e)
- {
- Assert (cell->face(face)->child(c)->line(e)
- ->n_active_fe_indices() == 1,
- ExcNotImplemented());
- Assert (cell->face(face)->child(c)->line(e)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcNotImplemented());
- }
- }
- for (unsigned int e=0; e<4; ++e)
- {
- Assert (cell->face(face)->line(e)
- ->n_active_fe_indices() == 1,
- ExcNotImplemented());
- Assert (cell->face(face)->line(e)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcNotImplemented());
- }
+ for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+ if (!cell->neighbor_child_on_subface(face,c)->is_artificial())
+ {
+ Assert (cell->face(face)->child(c)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+ for (unsigned int e=0; e<4; ++e)
+ {
+ Assert (cell->face(face)->child(c)->line(e)
+ ->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)->child(c)->line(e)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+ }
+ }
+ for (unsigned int e=0; e<4; ++e)
+ {
+ Assert (cell->face(face)->line(e)
+ ->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)->line(e)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+ }
// ok, start up the work
- const FiniteElement<dim> &fe = cell->get_fe();
- const unsigned int fe_index = cell->active_fe_index();
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int fe_index = cell->active_fe_index();
- const unsigned int n_dofs_on_mother = fe.dofs_per_face;
- const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+
- 12*fe.dofs_per_line+
- 4*fe.dofs_per_quad);
+ const unsigned int n_dofs_on_mother = fe.dofs_per_face;
+ const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+
+ 12*fe.dofs_per_line+
+ 4*fe.dofs_per_quad);
//TODO[TL]: think about this and the following in case of anisotropic refinement
- dofs_on_mother.resize (n_dofs_on_mother);
- dofs_on_children.resize (n_dofs_on_children);
+ dofs_on_mother.resize (n_dofs_on_mother);
+ dofs_on_children.resize (n_dofs_on_children);
- Assert(n_dofs_on_mother == fe.constraints().n(),
- ExcDimensionMismatch(n_dofs_on_mother,
- fe.constraints().n()));
- Assert(n_dofs_on_children == fe.constraints().m(),
- ExcDimensionMismatch(n_dofs_on_children,
- fe.constraints().m()));
+ Assert(n_dofs_on_mother == fe.constraints().n(),
+ ExcDimensionMismatch(n_dofs_on_mother,
+ fe.constraints().n()));
+ Assert(n_dofs_on_children == fe.constraints().m(),
+ ExcDimensionMismatch(n_dofs_on_children,
+ fe.constraints().m()));
- const typename DH::face_iterator this_face = cell->face(face);
+ const typename DH::face_iterator this_face = cell->face(face);
// fill the dofs indices. Use same
// enumeration scheme as in
// @p{FiniteElement::constraints()}
- unsigned int next_index = 0;
- for (unsigned int vertex=0; vertex<4; ++vertex)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
- fe_index);
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_mother[next_index++]
- = this_face->line(line)->dof_index(dof, fe_index);
- for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
- dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
- AssertDimension (next_index, dofs_on_mother.size());
-
- next_index = 0;
+ unsigned int next_index = 0;
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
+ fe_index);
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_mother[next_index++]
+ = this_face->line(line)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+ dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
+ AssertDimension (next_index, dofs_on_mother.size());
+
+ next_index = 0;
// assert some consistency
// assumptions
//TODO[TL]: think about this in case of anisotropic refinement
- Assert (dof_handler.get_tria().get_anisotropic_refinement_flag() ||
- ((this_face->child(0)->vertex_index(3) ==
- this_face->child(1)->vertex_index(2)) &&
- (this_face->child(0)->vertex_index(3) ==
- this_face->child(2)->vertex_index(1)) &&
- (this_face->child(0)->vertex_index(3) ==
- this_face->child(3)->vertex_index(0))),
- ExcInternalError());
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(0)->vertex_dof_index(3,dof);
+ Assert (dof_handler.get_tria().get_anisotropic_refinement_flag() ||
+ ((this_face->child(0)->vertex_index(3) ==
+ this_face->child(1)->vertex_index(2)) &&
+ (this_face->child(0)->vertex_index(3) ==
+ this_face->child(2)->vertex_index(1)) &&
+ (this_face->child(0)->vertex_index(3) ==
+ this_face->child(3)->vertex_index(0))),
+ ExcInternalError());
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->vertex_dof_index(3,dof);
// dof numbers on the centers of
// the lines bounding this face
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_children[next_index++]
- = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index);
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index);
// next the dofs on the lines interior
// to the face; the order of these
// lines is laid down in the
// FiniteElement class documentation
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(0)->line(1)->dof_index(dof, fe_index);
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(2)->line(1)->dof_index(dof, fe_index);
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(0)->line(3)->dof_index(dof, fe_index);
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(1)->line(3)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->line(1)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(2)->line(1)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->line(3)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(1)->line(3)->dof_index(dof, fe_index);
// dofs on the bordering lines
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int child=0; child<2; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = this_face->line(line)->child(child)->dof_index(dof, fe_index);
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int child=0; child<2; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->line(line)->child(child)->dof_index(dof, fe_index);
// finally, for the dofs interior
// to the four child faces
- for (unsigned int child=0; child<4; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
- dofs_on_children[next_index++]
- = this_face->child(child)->dof_index(dof, fe_index);
- AssertDimension (next_index, dofs_on_children.size());
+ for (unsigned int child=0; child<4; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(child)->dof_index(dof, fe_index);
+ AssertDimension (next_index, dofs_on_children.size());
// for each row in the constraint
// matrix for this line:
- for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
- {
- constraints.add_line (dofs_on_children[row]);
- for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
- constraints.add_entry (dofs_on_children[row],
- dofs_on_mother[i],
- fe.constraints()(row,i));
-
- constraints.set_inhomogeneity(dofs_on_children[row], 0.);
- }
- }
- else
- {
+ for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
+ {
+ constraints.add_line (dofs_on_children[row]);
+ for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
+ constraints.add_entry (dofs_on_children[row],
+ dofs_on_mother[i],
+ fe.constraints()(row,i));
+
+ constraints.set_inhomogeneity(dofs_on_children[row], 0.);
+ }
+ }
+ else
+ {
// this face has no
// children, but it
// could still be
// case that the
// neighbor is an
// artificial cell
- if (!cell->at_boundary(face) &&
- !cell->neighbor(face)->is_artificial())
- {
- Assert (cell->face(face)->n_active_fe_indices() == 1,
- ExcNotImplemented());
- Assert (cell->face(face)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcInternalError());
- }
- }
+ if (!cell->at_boundary(face) &&
+ !cell->neighbor(face)->is_artificial())
+ {
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcInternalError());
+ }
+ }
}
template <class DH>
void
make_hp_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints)
+ ConstraintMatrix &constraints)
{
// note: this function is going
// to be hard to understand if
// time they are needed, and
// then just reuse them
Table<2,std_cxx1x::shared_ptr<FullMatrix<double> > >
- face_interpolation_matrices (n_finite_elements (dof_handler),
- n_finite_elements (dof_handler));
+ face_interpolation_matrices (n_finite_elements (dof_handler),
+ n_finite_elements (dof_handler));
Table<3,std_cxx1x::shared_ptr<FullMatrix<double> > >
- subface_interpolation_matrices (n_finite_elements (dof_handler),
- n_finite_elements (dof_handler),
- GeometryInfo<dim>::max_children_per_face);
+ subface_interpolation_matrices (n_finite_elements (dof_handler),
+ n_finite_elements (dof_handler),
+ GeometryInfo<dim>::max_children_per_face);
// similarly have a cache for
// the matrices that are split
// interpolation matrix as
// described in the @ref hp_paper "hp paper"
Table<2,std_cxx1x::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > > >
- split_face_interpolation_matrices (n_finite_elements (dof_handler),
- n_finite_elements (dof_handler));
+ split_face_interpolation_matrices (n_finite_elements (dof_handler),
+ n_finite_elements (dof_handler));
// finally, for each pair of finite
// elements, have a mask that states
// the coarse side of a refined face
// will act as master dofs.
Table<2,std_cxx1x::shared_ptr<std::vector<bool> > >
- master_dof_masks (n_finite_elements (dof_handler),
- n_finite_elements (dof_handler));
+ master_dof_masks (n_finite_elements (dof_handler),
+ n_finite_elements (dof_handler));
// loop over all faces
//
// refined, we can only visit each
// face with hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
// artificial cells can at
// best neighbor ghost cells,
// but we're not interested
// in these interfaces
- if (!cell->is_artificial ())
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->has_children())
- {
+ if (!cell->is_artificial ())
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->has_children())
+ {
// first of all, make sure that
// we treat a case which is
// possible, i.e. either no dofs
// on the face at all or no
// anisotropic refinement
- if (cell->get_fe().dofs_per_face == 0)
- continue;
+ if (cell->get_fe().dofs_per_face == 0)
+ continue;
- Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
- ExcNotImplemented());
+ Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
+ ExcNotImplemented());
// so now we've found a
// face of an active
// each of the children
// can have only one as
// well. check this
- Assert (cell->face(face)->n_active_fe_indices() == 1,
- ExcInternalError());
- Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
- == true,
- ExcInternalError());
- for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
- Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
- ExcInternalError());
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcInternalError());
+ Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+ == true,
+ ExcInternalError());
+ for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+ Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
+ ExcInternalError());
// first find out
// whether we can
// ignore all
// interfaces with
// artificial cells
- FiniteElementDomination::Domination
- mother_face_dominates = FiniteElementDomination::either_element_can_dominate;
-
- if (DoFHandlerSupportsDifferentFEs<DH>::value == true)
- for (unsigned int c=0; c<cell->face(face)->number_of_children(); ++c)
- if (!cell->neighbor_child_on_subface (face, c)->is_artificial())
- mother_face_dominates = mother_face_dominates &
- (cell->get_fe().compare_for_face_domination
- (cell->neighbor_child_on_subface (face, c)->get_fe()));
-
- switch (mother_face_dominates)
- {
- case FiniteElementDomination::this_element_dominates:
- case FiniteElementDomination::either_element_can_dominate:
- {
+ FiniteElementDomination::Domination
+ mother_face_dominates = FiniteElementDomination::either_element_can_dominate;
+
+ if (DoFHandlerSupportsDifferentFEs<DH>::value == true)
+ for (unsigned int c=0; c<cell->face(face)->number_of_children(); ++c)
+ if (!cell->neighbor_child_on_subface (face, c)->is_artificial())
+ mother_face_dominates = mother_face_dominates &
+ (cell->get_fe().compare_for_face_domination
+ (cell->neighbor_child_on_subface (face, c)->get_fe()));
+
+ switch (mother_face_dominates)
+ {
+ case FiniteElementDomination::this_element_dominates:
+ case FiniteElementDomination::either_element_can_dominate:
+ {
// Case 1 (the
// simple case
// and the only
// against the
// DoFs on the
// face itself
- master_dofs.resize (cell->get_fe().dofs_per_face);
+ master_dofs.resize (cell->get_fe().dofs_per_face);
- cell->face(face)->get_dof_indices (master_dofs,
- cell->active_fe_index ());
+ cell->face(face)->get_dof_indices (master_dofs,
+ cell->active_fe_index ());
// Now create
// constraint matrix
// node is only
// detected if we also
// look between ghosts
- for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
- {
- if (cell->neighbor_child_on_subface (face, c)->is_artificial()
- ||
- (dim == 2 && cell->neighbor_child_on_subface (face, c)->is_ghost()))
- continue;
+ for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+ {
+ if (cell->neighbor_child_on_subface (face, c)->is_artificial()
+ ||
+ (dim == 2 && cell->neighbor_child_on_subface (face, c)->is_ghost()))
+ continue;
- const typename DH::active_face_iterator
- subface = cell->face(face)->child(c);
+ const typename DH::active_face_iterator
+ subface = cell->face(face)->child(c);
- Assert (subface->n_active_fe_indices() == 1,
- ExcInternalError());
+ Assert (subface->n_active_fe_indices() == 1,
+ ExcInternalError());
- const unsigned int
- subface_fe_index = subface->nth_active_fe_index(0);
+ const unsigned int
+ subface_fe_index = subface->nth_active_fe_index(0);
// we sometime run
// into the
// care of. in that
// case, just
// continue
- if (cell->get_fe().compare_for_face_domination
- (subface->get_fe(subface_fe_index))
- ==
- FiniteElementDomination::no_requirements)
- continue;
+ if (cell->get_fe().compare_for_face_domination
+ (subface->get_fe(subface_fe_index))
+ ==
+ FiniteElementDomination::no_requirements)
+ continue;
// Same procedure as for the
// mother cell. Extract the face
// DoFs from the cell DoFs.
- slave_dofs.resize (subface->get_fe(subface_fe_index)
- .dofs_per_face);
- subface->get_dof_indices (slave_dofs, subface_fe_index);
+ slave_dofs.resize (subface->get_fe(subface_fe_index)
+ .dofs_per_face);
+ subface->get_dof_indices (slave_dofs, subface_fe_index);
- for (unsigned int i=0; i<slave_dofs.size(); ++i)
- Assert (slave_dofs[i] != numbers::invalid_unsigned_int,
- ExcInternalError());
+ for (unsigned int i=0; i<slave_dofs.size(); ++i)
+ Assert (slave_dofs[i] != numbers::invalid_unsigned_int,
+ ExcInternalError());
// Now create the
// element constraint
// properties of a
// finite element onto
// that mesh
- ensure_existence_of_subface_matrix
- (cell->get_fe(),
- subface->get_fe(subface_fe_index),
- c,
- subface_interpolation_matrices
- [cell->active_fe_index()][subface_fe_index][c]);
+ ensure_existence_of_subface_matrix
+ (cell->get_fe(),
+ subface->get_fe(subface_fe_index),
+ c,
+ subface_interpolation_matrices
+ [cell->active_fe_index()][subface_fe_index][c]);
// Add constraints to global constraint
// matrix.
- filter_constraints (master_dofs,
- slave_dofs,
- *(subface_interpolation_matrices
- [cell->active_fe_index()][subface_fe_index][c]),
- constraints);
- }
-
- break;
- }
-
- case FiniteElementDomination::other_element_dominates:
- case FiniteElementDomination::neither_element_dominates:
- {
+ filter_constraints (master_dofs,
+ slave_dofs,
+ *(subface_interpolation_matrices
+ [cell->active_fe_index()][subface_fe_index][c]),
+ constraints);
+ }
+
+ break;
+ }
+
+ case FiniteElementDomination::other_element_dominates:
+ case FiniteElementDomination::neither_element_dominates:
+ {
// Case 2 (the "complex"
// case): at least one
// (the neither_... case)
// dof handlers,
// add a check
// here...
- Assert (DoFHandlerSupportsDifferentFEs<DH>::value == true,
- ExcInternalError());
+ Assert (DoFHandlerSupportsDifferentFEs<DH>::value == true,
+ ExcInternalError());
// we first have
// to find the
// the other ones
// can be
// constrained to
- const unsigned int dominating_fe_index
- = get_most_dominating_subface_fe_index (cell->face(face));
+ const unsigned int dominating_fe_index
+ = get_most_dominating_subface_fe_index (cell->face(face));
- const FiniteElement<dim,spacedim> &dominating_fe
- = dof_handler.get_fe()[dominating_fe_index];
+ const FiniteElement<dim,spacedim> &dominating_fe
+ = dof_handler.get_fe()[dominating_fe_index];
// check also
// that it is
// into the
// branch for the
// 'complex' case
- Assert ((dominating_fe.compare_for_face_domination
- (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0)))
- == FiniteElementDomination::this_element_dominates)
- ||
- (dominating_fe.compare_for_face_domination
- (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0)))
- == FiniteElementDomination::either_element_can_dominate),
- ExcInternalError());
+ Assert ((dominating_fe.compare_for_face_domination
+ (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0)))
+ == FiniteElementDomination::this_element_dominates)
+ ||
+ (dominating_fe.compare_for_face_domination
+ (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0)))
+ == FiniteElementDomination::either_element_can_dominate),
+ ExcInternalError());
// first get the
// interpolation matrix
// from the mother to the
// virtual dofs
- Assert (dominating_fe.dofs_per_face <=
- cell->get_fe().dofs_per_face,
- ExcInternalError());
+ Assert (dominating_fe.dofs_per_face <=
+ cell->get_fe().dofs_per_face,
+ ExcInternalError());
- ensure_existence_of_face_matrix
- (dominating_fe,
- cell->get_fe(),
- face_interpolation_matrices
- [dominating_fe_index][cell->active_fe_index()]);
+ ensure_existence_of_face_matrix
+ (dominating_fe,
+ cell->get_fe(),
+ face_interpolation_matrices
+ [dominating_fe_index][cell->active_fe_index()]);
// split this matrix into
// master and slave
// components. invert the
// master component
- ensure_existence_of_master_dof_mask
- (cell->get_fe(),
- dominating_fe,
- (*face_interpolation_matrices
- [dominating_fe_index]
- [cell->active_fe_index()]),
- master_dof_masks
- [dominating_fe_index]
- [cell->active_fe_index()]);
-
- ensure_existence_of_split_face_matrix
- (*face_interpolation_matrices
- [dominating_fe_index][cell->active_fe_index()],
- (*master_dof_masks
- [dominating_fe_index][cell->active_fe_index()]),
- split_face_interpolation_matrices
- [dominating_fe_index][cell->active_fe_index()]);
-
- const FullMatrix<double> &restrict_mother_to_virtual_master_inv
- = (split_face_interpolation_matrices
- [dominating_fe_index][cell->active_fe_index()]->first);
-
- const FullMatrix<double> &restrict_mother_to_virtual_slave
- = (split_face_interpolation_matrices
- [dominating_fe_index][cell->active_fe_index()]->second);
+ ensure_existence_of_master_dof_mask
+ (cell->get_fe(),
+ dominating_fe,
+ (*face_interpolation_matrices
+ [dominating_fe_index]
+ [cell->active_fe_index()]),
+ master_dof_masks
+ [dominating_fe_index]
+ [cell->active_fe_index()]);
+
+ ensure_existence_of_split_face_matrix
+ (*face_interpolation_matrices
+ [dominating_fe_index][cell->active_fe_index()],
+ (*master_dof_masks
+ [dominating_fe_index][cell->active_fe_index()]),
+ split_face_interpolation_matrices
+ [dominating_fe_index][cell->active_fe_index()]);
+
+ const FullMatrix<double> &restrict_mother_to_virtual_master_inv
+ = (split_face_interpolation_matrices
+ [dominating_fe_index][cell->active_fe_index()]->first);
+
+ const FullMatrix<double> &restrict_mother_to_virtual_slave
+ = (split_face_interpolation_matrices
+ [dominating_fe_index][cell->active_fe_index()]->second);
// now compute
// the constraint
// inverse matrix
// and the slave
// part
- constraint_matrix.reinit (cell->get_fe().dofs_per_face -
- dominating_fe.dofs_per_face,
- dominating_fe.dofs_per_face);
- restrict_mother_to_virtual_slave
- .mmult (constraint_matrix,
- restrict_mother_to_virtual_master_inv);
+ constraint_matrix.reinit (cell->get_fe().dofs_per_face -
+ dominating_fe.dofs_per_face,
+ dominating_fe.dofs_per_face);
+ restrict_mother_to_virtual_slave
+ .mmult (constraint_matrix,
+ restrict_mother_to_virtual_master_inv);
// then figure
// out the global
// slave dofs and
// apply
// constraints
- scratch_dofs.resize (cell->get_fe().dofs_per_face);
- cell->face(face)->get_dof_indices (scratch_dofs,
- cell->active_fe_index ());
+ scratch_dofs.resize (cell->get_fe().dofs_per_face);
+ cell->face(face)->get_dof_indices (scratch_dofs,
+ cell->active_fe_index ());
// split dofs into master
// and slave components
- master_dofs.clear ();
- slave_dofs.clear ();
- for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
- if ((*master_dof_masks
- [dominating_fe_index][cell->active_fe_index()])[i] == true)
- master_dofs.push_back (scratch_dofs[i]);
- else
- slave_dofs.push_back (scratch_dofs[i]);
+ master_dofs.clear ();
+ slave_dofs.clear ();
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
+ if ((*master_dof_masks
+ [dominating_fe_index][cell->active_fe_index()])[i] == true)
+ master_dofs.push_back (scratch_dofs[i]);
+ else
+ slave_dofs.push_back (scratch_dofs[i]);
- AssertDimension (master_dofs.size(), dominating_fe.dofs_per_face);
- AssertDimension (slave_dofs.size(),
- cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face);
+ AssertDimension (master_dofs.size(), dominating_fe.dofs_per_face);
+ AssertDimension (slave_dofs.size(),
+ cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face);
- filter_constraints (master_dofs,
- slave_dofs,
- constraint_matrix,
- constraints);
+ filter_constraints (master_dofs,
+ slave_dofs,
+ constraint_matrix,
+ constraints);
// as discussed
// in the hp
// paper
- for (unsigned int sf=0;
- sf<cell->face(face)->n_children(); ++sf)
- {
+ for (unsigned int sf=0;
+ sf<cell->face(face)->n_children(); ++sf)
+ {
// ignore
// interfaces
// with
// between
// ghost
// cells in 2d
- if (cell->neighbor_child_on_subface (face, sf)->is_artificial()
- ||
- (dim==2 && cell->is_ghost()
- &&
- cell->neighbor_child_on_subface (face, sf)->is_ghost()))
- continue;
-
- Assert (cell->face(face)->child(sf)
- ->n_active_fe_indices() == 1,
- ExcInternalError());
-
- const unsigned int subface_fe_index
- = cell->face(face)->child(sf)->nth_active_fe_index(0);
- const FiniteElement<dim,spacedim> &subface_fe
- = dof_handler.get_fe()[subface_fe_index];
+ if (cell->neighbor_child_on_subface (face, sf)->is_artificial()
+ ||
+ (dim==2 && cell->is_ghost()
+ &&
+ cell->neighbor_child_on_subface (face, sf)->is_ghost()))
+ continue;
+
+ Assert (cell->face(face)->child(sf)
+ ->n_active_fe_indices() == 1,
+ ExcInternalError());
+
+ const unsigned int subface_fe_index
+ = cell->face(face)->child(sf)->nth_active_fe_index(0);
+ const FiniteElement<dim,spacedim> &subface_fe
+ = dof_handler.get_fe()[subface_fe_index];
// first get the
// interpolation
// matrix from the
// subface to the
// virtual dofs
- Assert (dominating_fe.dofs_per_face <=
- subface_fe.dofs_per_face,
- ExcInternalError());
- ensure_existence_of_subface_matrix
- (dominating_fe,
- subface_fe,
- sf,
- subface_interpolation_matrices
- [dominating_fe_index][subface_fe_index][sf]);
-
- const FullMatrix<double> &restrict_subface_to_virtual
- = *(subface_interpolation_matrices
- [dominating_fe_index][subface_fe_index][sf]);
-
- constraint_matrix.reinit (subface_fe.dofs_per_face,
- dominating_fe.dofs_per_face);
-
- restrict_subface_to_virtual
- .mmult (constraint_matrix,
- restrict_mother_to_virtual_master_inv);
-
- slave_dofs.resize (subface_fe.dofs_per_face);
- cell->face(face)->child(sf)->get_dof_indices (slave_dofs,
- subface_fe_index);
-
- filter_constraints (master_dofs,
- slave_dofs,
- constraint_matrix,
- constraints);
- }
-
- break;
- }
-
- case FiniteElementDomination::no_requirements:
+ Assert (dominating_fe.dofs_per_face <=
+ subface_fe.dofs_per_face,
+ ExcInternalError());
+ ensure_existence_of_subface_matrix
+ (dominating_fe,
+ subface_fe,
+ sf,
+ subface_interpolation_matrices
+ [dominating_fe_index][subface_fe_index][sf]);
+
+ const FullMatrix<double> &restrict_subface_to_virtual
+ = *(subface_interpolation_matrices
+ [dominating_fe_index][subface_fe_index][sf]);
+
+ constraint_matrix.reinit (subface_fe.dofs_per_face,
+ dominating_fe.dofs_per_face);
+
+ restrict_subface_to_virtual
+ .mmult (constraint_matrix,
+ restrict_mother_to_virtual_master_inv);
+
+ slave_dofs.resize (subface_fe.dofs_per_face);
+ cell->face(face)->child(sf)->get_dof_indices (slave_dofs,
+ subface_fe_index);
+
+ filter_constraints (master_dofs,
+ slave_dofs,
+ constraint_matrix,
+ constraints);
+ }
+
+ break;
+ }
+
+ case FiniteElementDomination::no_requirements:
// there
// are no
// continuity
// elements. record
// no
// constraints
- break;
+ break;
- default:
+ default:
// we shouldn't get here
- Assert (false, ExcInternalError());
- }
- }
- else
- {
+ Assert (false, ExcInternalError());
+ }
+ }
+ else
+ {
// this face has no
// children, but it
// could still be that
// it is shared by two
// cells that use a
// different fe index
- Assert (cell->face(face)
- ->fe_index_is_active(cell->active_fe_index()) == true,
- ExcInternalError());
+ Assert (cell->face(face)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcInternalError());
// see if there is a
// neighbor that is
// cells may not have
// an active_fe_index
// set, etc
- if (!cell->at_boundary(face)
- &&
- cell->neighbor(face)->is_artificial())
- continue;
+ if (!cell->at_boundary(face)
+ &&
+ cell->neighbor(face)->is_artificial())
+ continue;
// Only if there is
// a neighbor with
// and the same h-level,
// some action has
// to be taken.
- if ((DoFHandlerSupportsDifferentFEs<DH>::value == true)
- &&
- !cell->face(face)->at_boundary ()
- &&
- (cell->neighbor(face)->active_fe_index () !=
- cell->active_fe_index ())
- &&
- (!cell->face(face)->has_children() &&
- !cell->neighbor_is_coarser(face) ))
- {
- const typename DH::cell_iterator neighbor = cell->neighbor (face);
+ if ((DoFHandlerSupportsDifferentFEs<DH>::value == true)
+ &&
+ !cell->face(face)->at_boundary ()
+ &&
+ (cell->neighbor(face)->active_fe_index () !=
+ cell->active_fe_index ())
+ &&
+ (!cell->face(face)->has_children() &&
+ !cell->neighbor_is_coarser(face) ))
+ {
+ const typename DH::cell_iterator neighbor = cell->neighbor (face);
// see which side of the
// face we have to
// constrain
- switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ()))
- {
- case FiniteElementDomination::this_element_dominates:
- {
+ switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ()))
+ {
+ case FiniteElementDomination::this_element_dominates:
+ {
// Get DoFs on
// dominating and
// dominated side of
// the face
- master_dofs.resize (cell->get_fe().dofs_per_face);
- cell->face(face)->get_dof_indices (master_dofs,
- cell->active_fe_index ());
+ master_dofs.resize (cell->get_fe().dofs_per_face);
+ cell->face(face)->get_dof_indices (master_dofs,
+ cell->active_fe_index ());
- slave_dofs.resize (neighbor->get_fe().dofs_per_face);
- cell->face(face)->get_dof_indices (slave_dofs,
- neighbor->active_fe_index ());
+ slave_dofs.resize (neighbor->get_fe().dofs_per_face);
+ cell->face(face)->get_dof_indices (slave_dofs,
+ neighbor->active_fe_index ());
// break if the n_master_dofs == 0,
// because we are attempting to
// constrain to an element that has
// has no face dofs
- if(master_dofs.size() == 0) break;
+ if(master_dofs.size() == 0) break;
// make sure
// the element
// for this
// face are
// available
- ensure_existence_of_face_matrix
- (cell->get_fe(),
- neighbor->get_fe(),
- face_interpolation_matrices
- [cell->active_fe_index()][neighbor->active_fe_index()]);
+ ensure_existence_of_face_matrix
+ (cell->get_fe(),
+ neighbor->get_fe(),
+ face_interpolation_matrices
+ [cell->active_fe_index()][neighbor->active_fe_index()]);
// Add constraints to global constraint
// matrix.
- filter_constraints (master_dofs,
- slave_dofs,
- *(face_interpolation_matrices
- [cell->active_fe_index()]
- [neighbor->active_fe_index()]),
- constraints);
-
- break;
- }
-
- case FiniteElementDomination::other_element_dominates:
- {
+ filter_constraints (master_dofs,
+ slave_dofs,
+ *(face_interpolation_matrices
+ [cell->active_fe_index()]
+ [neighbor->active_fe_index()]),
+ constraints);
+
+ break;
+ }
+
+ case FiniteElementDomination::other_element_dominates:
+ {
// we don't do anything
// here since we will
// come back to this
// we will fall into
// the first case
// clause above
- break;
- }
+ break;
+ }
- case FiniteElementDomination::either_element_can_dominate:
- {
+ case FiniteElementDomination::either_element_can_dominate:
+ {
// it appears as if
// neither element has
// any constraints on
// finite elements
// actually do have
// dofs on the face
- if ((cell->get_fe().dofs_per_face != 0)
- ||
- (cell->neighbor(face)->get_fe().dofs_per_face != 0))
- {
- Assert (cell->get_fe().dofs_per_face
- ==
- cell->neighbor(face)->get_fe().dofs_per_face,
- ExcNotImplemented());
+ if ((cell->get_fe().dofs_per_face != 0)
+ ||
+ (cell->neighbor(face)->get_fe().dofs_per_face != 0))
+ {
+ Assert (cell->get_fe().dofs_per_face
+ ==
+ cell->neighbor(face)->get_fe().dofs_per_face,
+ ExcNotImplemented());
// (ab)use the master
// and slave dofs
// arrays for a
// moment here
- master_dofs.resize (cell->get_fe().dofs_per_face);
- cell->face(face)
- ->get_dof_indices (master_dofs,
- cell->active_fe_index ());
+ master_dofs.resize (cell->get_fe().dofs_per_face);
+ cell->face(face)
+ ->get_dof_indices (master_dofs,
+ cell->active_fe_index ());
- slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face);
- cell->face(face)
- ->get_dof_indices (slave_dofs,
- cell->neighbor(face)->active_fe_index ());
+ slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face);
+ cell->face(face)
+ ->get_dof_indices (slave_dofs,
+ cell->neighbor(face)->active_fe_index ());
- for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
- AssertDimension (master_dofs[i], slave_dofs[i]);
- }
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
+ AssertDimension (master_dofs[i], slave_dofs[i]);
+ }
- break;
- }
+ break;
+ }
- case FiniteElementDomination::neither_element_dominates:
- {
+ case FiniteElementDomination::neither_element_dominates:
+ {
// we don't presently
// know what exactly to
// do here. it isn't quite
// following statement
// and see what exactly
// is going on
- Assert (false, ExcNotImplemented());
- break;
- }
+ Assert (false, ExcNotImplemented());
+ break;
+ }
- case FiniteElementDomination::no_requirements:
- {
+ case FiniteElementDomination::no_requirements:
+ {
// nothing to do here
- break;
- }
+ break;
+ }
- default:
+ default:
// we shouldn't get
// here
- Assert (false, ExcInternalError());
- }
- }
- }
+ Assert (false, ExcInternalError());
+ }
+ }
+ }
}
}
template <class DH>
void
make_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints)
+ ConstraintMatrix &constraints)
{
// Decide whether to use the
// new or old make_hanging_node_constraints
// for the moment.
if (dof_handler.get_fe().hp_constraints_are_implemented ())
internal::
- make_hp_hanging_node_constraints (dof_handler,
- constraints);
+ make_hp_hanging_node_constraints (dof_handler,
+ constraints);
else
internal::
- make_oldstyle_hanging_node_constraints (dof_handler,
- constraints,
- dealii::internal::int2type<DH::dimension>());
+ make_oldstyle_hanging_node_constraints (dof_handler,
+ constraints,
+ dealii::internal::int2type<DH::dimension>());
}
ExcMessage ("Faces for periodicity constraints must be on the boundary"));
- // In the case that both faces have
- // children, we loop over all children
- // and applu make_periodicty_constrains
- // recursively:
+ // In the case that both faces have
+ // children, we loop over all children
+ // and applu make_periodicty_constrains
+ // recursively:
if (face_1->has_children() && face_2->has_children()) {
Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
}
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]) {
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)
{
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);
}
}
// component can be assigned. Then, we sort
// them to the first selected component of the
// vector system.
- //
- // the output array dofs_by_component
- // lists for each dof the corresponding
- // vector component. if the DoFHandler is
- // based on a parallel distributed
- // triangulation then the output array is
- // index by
- // dof.locally_owned_dofs().index_within_set(indices[i])
+ //
+ // the output array dofs_by_component
+ // lists for each dof the corresponding
+ // vector component. if the DoFHandler is
+ // based on a parallel distributed
+ // triangulation then the output array is
+ // index by
+ // dof.locally_owned_dofs().index_within_set(indices[i])
template <class DH>
inline
void
extract_dofs_by_component (const DH &dof,
- const std::vector<bool> &component_select,
- const bool sort_by_blocks,
- std::vector<unsigned char> &dofs_by_component)
+ const std::vector<bool> &component_select,
+ const bool sort_by_blocks,
+ std::vector<unsigned char> &dofs_by_component)
{
const dealii::hp::FECollection<DH::dimension,DH::space_dimension>
- fe_collection (dof.get_fe());
+ fe_collection (dof.get_fe());
Assert (fe_collection.n_components() < 256, ExcNotImplemented());
Assert (dofs_by_component.size() == dof.n_locally_owned_dofs(),
- ExcDimensionMismatch(dofs_by_component.size(),
- dof.n_locally_owned_dofs()));
+ ExcDimensionMismatch(dofs_by_component.size(),
+ dof.n_locally_owned_dofs()));
// next set up a table for the degrees
// of freedom on each of the cells
// (regardless of the fact whether it
// is listed in the component_select
// argument or not)
- //
- // for each element 'f' of the
- // FECollection,
- // local_component_association[f][d]
- // then returns the vector component
- // that degree of freedom 'd' belongs
- // to (or, in the case of
- // sort_by_blocks, the block that it
- // corresponds to)
+ //
+ // for each element 'f' of the
+ // FECollection,
+ // local_component_association[f][d]
+ // then returns the vector component
+ // that degree of freedom 'd' belongs
+ // to (or, in the case of
+ // sort_by_blocks, the block that it
+ // corresponds to)
std::vector<std::vector<unsigned char> > local_component_association
- (fe_collection.size());
+ (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_component_association[f].resize(fe.dofs_per_cell);
if (sort_by_blocks == true)
- // compute the block each
- // local dof belongs to
+ // 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
+ // 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] =
// 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());
+ {
+ 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;
- }
- }
- }
+ 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
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];
- }
+ 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];
+ }
}
}
AssertDimension (dof_data.size(), dof_handler.n_dofs());
AssertIndexRange (component, n_components(dof_handler));
Assert (fe_is_primitive(dof_handler) == true,
- typename FiniteElement<dim>::ExcFENotPrimitive());
+ typename FiniteElement<dim>::ExcFENotPrimitive());
// store a flag whether we should care
// about different components. this is
dof_data = 0;
else
{
- std::vector<unsigned char> component_dofs (dof_handler.n_locally_owned_dofs());
- std::vector<bool> component_mask (dof_handler.get_fe().n_components(),
- false);
- component_mask[component] = true;
- internal::extract_dofs_by_component (dof_handler, component_mask,
- false, component_dofs);
-
- for (unsigned int i=0; i<dof_data.size(); ++i)
- if (component_dofs[i] == static_cast<unsigned char>(component))
- dof_data(i) = 0;
+ std::vector<unsigned char> component_dofs (dof_handler.n_locally_owned_dofs());
+ std::vector<bool> component_mask (dof_handler.get_fe().n_components(),
+ false);
+ component_mask[component] = true;
+ internal::extract_dofs_by_component (dof_handler, component_mask,
+ false, component_dofs);
+
+ for (unsigned int i=0; i<dof_data.size(); ++i)
+ if (component_dofs[i] == static_cast<unsigned char>(component))
+ dof_data(i) = 0;
}
// count how often we have added a value
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;
dof_indices.reserve (max_dofs_per_cell(dof_handler));
for (unsigned int present_cell = 0; cell!=endc; ++cell, ++present_cell)
{
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- dof_indices.resize (dofs_per_cell);
- cell->get_dof_indices (dof_indices);
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
// consider this dof only if it
// is the right component. if there
// is only one component, short cut
// the test
- if (!consider_components ||
- (cell->get_fe().system_to_component_index(i).first == component))
- {
+ if (!consider_components ||
+ (cell->get_fe().system_to_component_index(i).first == component))
+ {
// sum up contribution of the
// present_cell to this dof
- dof_data(dof_indices[i]) += cell_data(present_cell);
+ dof_data(dof_indices[i]) += cell_data(present_cell);
// note that we added another
// summand
- ++touch_count[dof_indices[i]];
- }
+ ++touch_count[dof_indices[i]];
+ }
}
// compute the mean value on all the
// at least once. this needs not be
// the case if the vector has more than
// one component
- Assert (consider_components || (touch_count[i]!=0),
- ExcInternalError());
- if (touch_count[i] != 0)
- dof_data(i) /= touch_count[i];
+ Assert (consider_components || (touch_count[i]!=0),
+ ExcInternalError());
+ if (touch_count[i] != 0)
+ dof_data(i) /= touch_count[i];
}
}
if (count_by_blocks == true)
{
- Assert(component_select.size() == fe.n_blocks(),
- ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
+ Assert(component_select.size() == fe.n_blocks(),
+ ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
}
else
{
- Assert(component_select.size() == n_components(dof),
- ExcDimensionMismatch(component_select.size(), n_components(dof)));
+ Assert(component_select.size() == n_components(dof),
+ ExcDimensionMismatch(component_select.size(), n_components(dof)));
}
Assert(selected_dofs.size() == dof.n_locally_owned_dofs(),
- ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs()));
+ ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs()));
// two special cases: no component
// is selected, and all components
// are selected; both rather
// stupid, but easy to catch
if (std::count (component_select.begin(), component_select.end(), true)
- == 0)
+ == 0)
{
- std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false);
- return;
+ std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false);
+ return;
}
else if (std::count (component_select.begin(), component_select.end(), true)
- == static_cast<signed int>(component_select.size()))
+ == static_cast<signed int>(component_select.size()))
{
- std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), true);
- return;
+ std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), true);
+ return;
}
// job.
std::vector<unsigned char> dofs_by_component (dof.n_locally_owned_dofs());
internal::extract_dofs_by_component (dof, component_select, count_by_blocks,
- dofs_by_component);
+ dofs_by_component);
for (unsigned int i=0; i<dof.n_locally_owned_dofs(); ++i)
if (component_select[dofs_by_component[i]] == true)
- selected_dofs[i] = true;
+ selected_dofs[i] = true;
}
if (count_by_blocks == true)
{
- Assert(component_select.size() == fe.n_blocks(),
- ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
+ Assert(component_select.size() == fe.n_blocks(),
+ ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
}
else
{
- Assert(component_select.size() == n_components(dof),
- ExcDimensionMismatch(component_select.size(), n_components(dof)));
+ Assert(component_select.size() == n_components(dof),
+ ExcDimensionMismatch(component_select.size(), n_components(dof)));
}
Assert(selected_dofs.size() == dof.n_dofs(),
- ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs()));
+ ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs()));
// two special cases: no component
// is selected, and all components
// are selected; both rather
// stupid, but easy to catch
if (std::count (component_select.begin(), component_select.end(), true)
- == 0)
+ == 0)
{
- std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
- return;
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
+ return;
};
if (std::count (component_select.begin(), component_select.end(), true)
- == static_cast<signed int>(component_select.size()))
+ == static_cast<signed int>(component_select.size()))
{
- std::fill_n (selected_dofs.begin(), dof.n_dofs(), true);
- return;
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(), true);
+ return;
};
// job.
std::vector<unsigned char> dofs_by_component (dof.n_dofs());
internal::extract_dofs_by_component (dof, component_select, count_by_blocks,
- dofs_by_component);
+ dofs_by_component);
for (unsigned int i=0; i<dof.n_dofs(); ++i)
if (component_select[dofs_by_component[i]] == true)
- selected_dofs[i] = true;
+ selected_dofs[i] = true;
}
if (count_by_blocks == true)
{
- Assert(component_select.size() == fe.n_blocks(),
- ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
+ 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_select.size() == fe.n_components(),
+ ExcDimensionMismatch(component_select.size(), fe.n_components()));
}
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
// are selected, both rather
// stupid, but easy to catch
if (std::count (component_select.begin(), component_select.end(), true)
- == 0)
+ == 0)
{
- std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
- return;
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
+ return;
};
if (std::count (component_select.begin(), component_select.end(), true)
- == static_cast<signed int>(component_select.size()))
+ == static_cast<signed int>(component_select.size()))
{
- std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true);
- return;
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true);
+ return;
};
// preset all values by false
std::vector<bool> local_selected_dofs (fe.dofs_per_cell, false);
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
if (count_by_blocks == true)
- local_selected_dofs[i]
- = component_select[fe.system_to_block_index(i).first];
+ local_selected_dofs[i]
+ = component_select[fe.system_to_block_index(i).first];
else
- if (fe.is_primitive(i))
- local_selected_dofs[i]
- = component_select[fe.system_to_component_index(i).first];
- else
+ if (fe.is_primitive(i))
+ local_selected_dofs[i]
+ = component_select[fe.system_to_component_index(i).first];
+ else
// if this shape function is
// not primitive, then we have
// to work harder. we have to
// the base element to which
// the local dof with index i
// belongs
- {
- unsigned int first_comp = 0;
- const unsigned int this_base = fe.system_to_base_index(i).first.first;
- const unsigned int this_multiplicity
- = fe.system_to_base_index(i).first.second;
-
- for (unsigned int b=0; b<this_base; ++b)
- first_comp += fe.base_element(b).n_components() *
- fe.element_multiplicity(b);
- for (unsigned int m=0; m<this_multiplicity; ++m)
- first_comp += fe.base_element(this_base).n_components();
- const unsigned int end_comp = first_comp +
- fe.base_element(this_base).n_components();
-
- Assert (first_comp < fe.n_components(), ExcInternalError());
- Assert (end_comp <= fe.n_components(), ExcInternalError());
+ {
+ unsigned int first_comp = 0;
+ const unsigned int this_base = fe.system_to_base_index(i).first.first;
+ const unsigned int this_multiplicity
+ = fe.system_to_base_index(i).first.second;
+
+ for (unsigned int b=0; b<this_base; ++b)
+ first_comp += fe.base_element(b).n_components() *
+ fe.element_multiplicity(b);
+ for (unsigned int m=0; m<this_multiplicity; ++m)
+ first_comp += fe.base_element(this_base).n_components();
+ const unsigned int end_comp = first_comp +
+ fe.base_element(this_base).n_components();
+
+ Assert (first_comp < fe.n_components(), ExcInternalError());
+ Assert (end_comp <= fe.n_components(), ExcInternalError());
// now check whether any of
// the components in between
// is set
- for (unsigned int c=first_comp; c<end_comp; ++c)
- if (component_select[c] == true)
- {
- local_selected_dofs[i] = true;
- break;
- }
- }
+ for (unsigned int c=first_comp; c<end_comp; ++c)
+ if (component_select[c] == true)
+ {
+ local_selected_dofs[i] = true;
+ break;
+ }
+ }
// then loop over all cells and do
// work
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];
}
}
extract_boundary_dofs (dof_handler, component_select,
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());
// 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)
+ == true)
selected_dofs.add_index (face_dof_indices[i]);
}
else // not primitive
{
const unsigned int first_nonzero_comp
- = (std::find (fe.get_nonzero_components(cell_index).begin(),
- fe.get_nonzero_components(cell_index).end(),
- true)
- -
- fe.get_nonzero_components(cell_index).begin());
+ = (std::find (fe.get_nonzero_components(cell_index).begin(),
+ fe.get_nonzero_components(cell_index).end(),
+ true)
+ -
+ fe.get_nonzero_components(cell_index).begin());
Assert (first_nonzero_comp < fe.n_components(),
ExcInternalError());
template <class DH>
void
extract_dofs_with_support_on_boundary (const DH &dof_handler,
- const std::vector<bool> &component_select,
- std::vector<bool> &selected_dofs,
- const std::set<types::boundary_id> &boundary_indicators)
+ const std::vector<bool> &component_select,
+ std::vector<bool> &selected_dofs,
+ const std::set<types::boundary_id> &boundary_indicators)
{
AssertDimension (component_select.size(), n_components(dof_handler));
Assert (boundary_indicators.find (numbers::internal_face_boundary_id) == boundary_indicators.end(),
- ExcInvalidBoundaryIndicator());
+ ExcInvalidBoundaryIndicator());
// let's see whether we have to
// check for certain boundary
// component is selected, or all
const bool check_vector_component
= (component_select != std::vector<bool>(component_select.size(),
- true));
+ true));
// clear and reset array by default
// values
// face which we will be visiting
// sooner or later
for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
- cell!=dof_handler.end(); ++cell)
+ cell!=dof_handler.end(); ++cell)
for (unsigned int face=0;
- face<GeometryInfo<DH::dimension>::faces_per_cell; ++face)
- if (cell->at_boundary(face))
- if (! check_boundary_indicator ||
- (boundary_indicators.find (cell->face(face)->boundary_indicator())
- != boundary_indicators.end()))
- {
- const FiniteElement<DH::dimension, DH::space_dimension> &fe = cell->get_fe();
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- cell_dof_indices.resize (dofs_per_cell);
- cell->get_dof_indices (cell_dof_indices);
-
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- if (fe.has_support_on_face(i,face))
- {
- if (!check_vector_component)
- selected_dofs[cell_dof_indices[i]] = true;
- else
+ face<GeometryInfo<DH::dimension>::faces_per_cell; ++face)
+ if (cell->at_boundary(face))
+ if (! check_boundary_indicator ||
+ (boundary_indicators.find (cell->face(face)->boundary_indicator())
+ != boundary_indicators.end()))
+ {
+ const FiniteElement<DH::dimension, DH::space_dimension> &fe = cell->get_fe();
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ cell_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (cell_dof_indices);
+
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ if (fe.has_support_on_face(i,face))
+ {
+ if (!check_vector_component)
+ selected_dofs[cell_dof_indices[i]] = true;
+ else
// check for
// component is
// required. somewhat
// but use usual
// convention (see
// docs)
- {
- if (fe.is_primitive (i))
- selected_dofs[cell_dof_indices[i]]
- = (component_select[fe.system_to_component_index(i).first]
- == true);
- else // not primitive
- {
- const unsigned int first_nonzero_comp
- = (std::find (fe.get_nonzero_components(i).begin(),
- fe.get_nonzero_components(i).end(),
- true)
- -
- fe.get_nonzero_components(i).begin());
- Assert (first_nonzero_comp < fe.n_components(),
- ExcInternalError());
-
- selected_dofs[cell_dof_indices[i]]
- = (component_select[first_nonzero_comp]
- == true);
- }
- }
- }
- }
+ {
+ if (fe.is_primitive (i))
+ selected_dofs[cell_dof_indices[i]]
+ = (component_select[fe.system_to_component_index(i).first]
+ == true);
+ else // not primitive
+ {
+ const unsigned int first_nonzero_comp
+ = (std::find (fe.get_nonzero_components(i).begin(),
+ fe.get_nonzero_components(i).end(),
+ true)
+ -
+ fe.get_nonzero_components(i).begin());
+ Assert (first_nonzero_comp < fe.n_components(),
+ ExcInternalError());
+
+ selected_dofs[cell_dof_indices[i]]
+ = (component_select[first_nonzero_comp]
+ == true);
+ }
+ }
+ }
+ }
}
{
template <int spacedim>
void extract_hanging_node_dofs (const dealii::DoFHandler<1,spacedim> &dof_handler,
- std::vector<bool> &selected_dofs)
+ std::vector<bool> &selected_dofs)
{
- Assert(selected_dofs.size() == dof_handler.n_dofs(),
- ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
+ Assert(selected_dofs.size() == dof_handler.n_dofs(),
+ ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
// preset all values by false
- std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+ std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
// there are no hanging nodes in 1d
}
template <int spacedim>
void extract_hanging_node_dofs (const dealii::DoFHandler<2,spacedim> &dof_handler,
- std::vector<bool> &selected_dofs)
+ std::vector<bool> &selected_dofs)
{
- const unsigned int dim = 2;
+ const unsigned int dim = 2;
- Assert(selected_dofs.size() == dof_handler.n_dofs(),
- ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
+ Assert(selected_dofs.size() == dof_handler.n_dofs(),
+ ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
// preset all values by false
- std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+ std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
- const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
+ const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
// this function is similar to the
// make_sparsity_pattern function,
// see there for more information
- typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->has_children())
- {
- const typename dealii::DoFHandler<dim,spacedim>::line_iterator
- line = cell->face(face);
-
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true;
-
- for (unsigned int child=0; child<2; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- selected_dofs[line->child(child)->dof_index(dof)] = true;
- }
+ typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->has_children())
+ {
+ const typename dealii::DoFHandler<dim,spacedim>::line_iterator
+ line = cell->face(face);
+
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true;
+
+ for (unsigned int child=0; child<2; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ selected_dofs[line->child(child)->dof_index(dof)] = true;
+ }
}
template <int spacedim>
void extract_hanging_node_dofs (const dealii::DoFHandler<3,spacedim> &dof_handler,
- std::vector<bool> &selected_dofs)
+ std::vector<bool> &selected_dofs)
{
- const unsigned int dim = 3;
+ const unsigned int dim = 3;
- Assert(selected_dofs.size() == dof_handler.n_dofs(),
- ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
+ Assert(selected_dofs.size() == dof_handler.n_dofs(),
+ ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
// preset all values by false
- std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+ std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
- const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
+ const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
// this function is similar to the
// make_sparsity_pattern function,
// see there for more information
- typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell->face(f)->has_children())
- {
- const typename dealii::DoFHandler<dim,spacedim>::face_iterator
- face = cell->face(f);
+ typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->face(f)->has_children())
+ {
+ const typename dealii::DoFHandler<dim,spacedim>::face_iterator
+ face = cell->face(f);
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true;
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true;
// dof numbers on the centers of
// the lines bounding this face
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true;
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true;
// next the dofs on the lines interior
// to the face; the order of these
// lines is laid down in the
// FiniteElement class documentation
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- selected_dofs[face->child(0)->line(1)->dof_index(dof)] = true;
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- selected_dofs[face->child(1)->line(2)->dof_index(dof)] = true;
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- selected_dofs[face->child(2)->line(3)->dof_index(dof)] = true;
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- selected_dofs[face->child(3)->line(0)->dof_index(dof)] = true;
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ selected_dofs[face->child(0)->line(1)->dof_index(dof)] = true;
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ selected_dofs[face->child(1)->line(2)->dof_index(dof)] = true;
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ selected_dofs[face->child(2)->line(3)->dof_index(dof)] = true;
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ selected_dofs[face->child(3)->line(0)->dof_index(dof)] = true;
// dofs on the bordering lines
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int child=0; child<2; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true;
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int child=0; child<2; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true;
// finally, for the dofs interior
// to the four child faces
- for (unsigned int child=0; child<4; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
- selected_dofs[face->child(child)->dof_index(dof)] = true;
- }
+ for (unsigned int child=0; child<4; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+ selected_dofs[face->child(child)->dof_index(dof)] = true;
+ }
}
}
}
void
extract_hanging_node_dofs (const DoFHandler<dim,spacedim> &dof_handler,
- std::vector<bool> &selected_dofs)
+ std::vector<bool> &selected_dofs)
{
internal::extract_hanging_node_dofs (dof_handler,
- selected_dofs);
+ selected_dofs);
}
template <class DH>
void
extract_subdomain_dofs (const DH &dof_handler,
- const types::subdomain_id subdomain_id,
- std::vector<bool> &selected_dofs)
+ const types::subdomain_id subdomain_id,
+ std::vector<bool> &selected_dofs)
{
Assert(selected_dofs.size() == dof_handler.n_dofs(),
- ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
+ ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
// preset all values by false
std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
endc = dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->subdomain_id() == subdomain_id)
- {
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- local_dof_indices.resize (dofs_per_cell);
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- selected_dofs[local_dof_indices[i]] = true;
- };
+ {
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ selected_dofs[local_dof_indices[i]] = true;
+ };
}
template <class DH>
void
extract_locally_owned_dofs (const DH & dof_handler,
- IndexSet & dof_set)
+ IndexSet & dof_set)
{
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
template <class DH>
void
extract_locally_active_dofs (const DH & dof_handler,
- IndexSet & dof_set)
+ IndexSet & dof_set)
{
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
std::set<unsigned int> 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);
-
- 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);
- }
+ {
+ dof_indices.resize(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(dof_indices);
+
+ 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);
+ }
dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
template <class DH>
void
extract_locally_relevant_dofs (const DH & dof_handler,
- IndexSet & dof_set)
+ IndexSet & dof_set)
{
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
std::set<unsigned int> 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);
-
- 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);
- }
+ {
+ dof_indices.resize(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(dof_indices);
+
+ 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);
+ }
dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
template <class DH>
void
extract_constant_modes (const DH &dof_handler,
- const std::vector<bool> &component_select,
- std::vector<std::vector<bool> > &constant_modes)
+ const std::vector<bool> &component_select,
+ std::vector<std::vector<bool> > &constant_modes)
{
const unsigned int n_components = dof_handler.get_fe().n_components();
Assert (n_components == component_select.size(),
- ExcDimensionMismatch(n_components,
- component_select.size()));
+ ExcDimensionMismatch(n_components,
+ component_select.size()));
std::vector<unsigned int> localized_component (n_components,
- numbers::invalid_unsigned_int);
+ numbers::invalid_unsigned_int);
unsigned int n_components_selected = 0;
for (unsigned int i=0; i<n_components; ++i)
if (component_select[i] == true)
- localized_component[i] = n_components_selected++;
+ localized_component[i] = n_components_selected++;
std::vector<unsigned char> dofs_by_component (dof_handler.n_locally_owned_dofs());
internal::extract_dofs_by_component (dof_handler, component_select, false,
- dofs_by_component);
+ dofs_by_component);
unsigned int n_selected_dofs = 0;
for (unsigned int i=0; i<n_components; ++i)
if (component_select[i] == true)
- n_selected_dofs += std::count (dofs_by_component.begin(),
- dofs_by_component.end(), i);
+ n_selected_dofs += std::count (dofs_by_component.begin(),
+ dofs_by_component.end(), i);
// First count the number of dofs
// in the current component.
constant_modes.resize (n_components_selected, std::vector<bool>(n_selected_dofs,
- false));
+ false));
std::vector<unsigned int> component_list (n_components, 0);
for (unsigned int d=0; d<n_components; ++d)
component_list[d] = component_select[d];
unsigned int counter = 0;
for (unsigned int i=0; i<dof_handler.n_locally_owned_dofs(); ++i)
if (component_select[dofs_by_component[i]])
- {
- constant_modes[localized_component[dofs_by_component[i]]][counter] = true;
- ++counter;
- }
+ {
+ constant_modes[localized_component[dofs_by_component[i]]][counter] = true;
+ ++counter;
+ }
}
template <class DH>
void
get_active_fe_indices (const DH &dof_handler,
- std::vector<unsigned int> &active_fe_indices)
+ std::vector<unsigned int> &active_fe_indices)
{
AssertDimension (active_fe_indices.size(), dof_handler.get_tria().n_active_cells());
template <class DH>
void
get_subdomain_association (const DH &dof_handler,
- std::vector<types::subdomain_id> &subdomain_association)
+ std::vector<types::subdomain_id> &subdomain_association)
{
// if the Triangulation is distributed, the
// only thing we can usefully ask is for
// its locally owned subdomain
Assert ((dynamic_cast<const parallel::distributed::
- Triangulation<DH::dimension,DH::space_dimension>*>
- (&dof_handler.get_tria()) == 0),
- 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."));
+ Triangulation<DH::dimension,DH::space_dimension>*>
+ (&dof_handler.get_tria()) == 0),
+ 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."));
Assert(subdomain_association.size() == dof_handler.n_dofs(),
- ExcDimensionMismatch(subdomain_association.size(),
- dof_handler.n_dofs()));
+ ExcDimensionMismatch(subdomain_association.size(),
+ dof_handler.n_dofs()));
// 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;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
- Assert (cell->is_artificial() == false,
- ExcMessage ("You can't call this function for meshes that "
- "have artificial cells."));
+ Assert (cell->is_artificial() == false,
+ ExcMessage ("You can't call this function for meshes that "
+ "have artificial cells."));
- const types::subdomain_id subdomain_id = cell->subdomain_id();
- const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- local_dof_indices.resize (dofs_per_cell);
- cell->get_dof_indices (local_dof_indices);
+ const types::subdomain_id subdomain_id = cell->subdomain_id();
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
// set subdomain ids. if dofs
// already have their values
// one, where we take "random"
// to be "once this way once
// that way"
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- if (subdomain_association[local_dof_indices[i]] ==
- numbers::invalid_unsigned_int)
- subdomain_association[local_dof_indices[i]] = subdomain_id;
- else
- {
- if (coin_flip == true)
- subdomain_association[local_dof_indices[i]] = subdomain_id;
- coin_flip = !coin_flip;
- }
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (subdomain_association[local_dof_indices[i]] ==
+ numbers::invalid_unsigned_int)
+ subdomain_association[local_dof_indices[i]] = subdomain_id;
+ else
+ {
+ if (coin_flip == true)
+ subdomain_association[local_dof_indices[i]] = subdomain_id;
+ coin_flip = !coin_flip;
+ }
}
Assert (std::find (subdomain_association.begin(),
- subdomain_association.end(),
- types::invalid_subdomain_id)
- == subdomain_association.end(),
- ExcInternalError());
+ subdomain_association.end(),
+ types::invalid_subdomain_id)
+ == subdomain_association.end(),
+ ExcInternalError());
}
template <class DH>
unsigned int
count_dofs_with_subdomain_association (const DH &dof_handler,
- const types::subdomain_id subdomain)
+ const types::subdomain_id subdomain)
{
std::vector<types::subdomain_id> subdomain_association (dof_handler.n_dofs());
get_subdomain_association (dof_handler, subdomain_association);
return std::count (subdomain_association.begin(),
- subdomain_association.end(),
- subdomain);
+ subdomain_association.end(),
+ subdomain);
}
template <class DH>
IndexSet
dof_indices_with_subdomain_association (const DH &dof_handler,
- const types::subdomain_id subdomain)
+ const types::subdomain_id subdomain)
{
// If we have a distributed::Triangulation only
||
(subdomain == dof_handler.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."));
IndexSet index_set (dof_handler.n_dofs());
std::vector<unsigned int> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // first generate an unsorted list of all
- // indices which we fill from the back. could
- // also insert them directly into the
- // IndexSet, but that inserts indices in the
- // middle, which is an O(n^2) algorithm and
- // hence too expensive. Could also use
- // std::set, but that is in general more
- // expensive than a vector
+ // first generate an unsorted list of all
+ // indices which we fill from the back. could
+ // also insert them directly into the
+ // IndexSet, but that inserts indices in the
+ // middle, which is an O(n^2) algorithm and
+ // hence too expensive. Could also use
+ // std::set, but that is in general more
+ // expensive than a vector
std::vector<unsigned int> subdomain_indices;
typename DH::active_cell_iterator
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 <class DH>
void
count_dofs_with_subdomain_association (const DH &dof_handler,
- const types::subdomain_id subdomain,
- std::vector<unsigned int> &n_dofs_on_subdomain)
+ const types::subdomain_id subdomain,
+ std::vector<unsigned int> &n_dofs_on_subdomain)
{
Assert (n_dofs_on_subdomain.size() == dof_handler.get_fe().n_components(),
- ExcDimensionMismatch (n_dofs_on_subdomain.size(),
- dof_handler.get_fe().n_components()));
+ ExcDimensionMismatch (n_dofs_on_subdomain.size(),
+ dof_handler.get_fe().n_components()));
std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
// in debug mode, make sure that there are
{
bool found = false;
for (typename Triangulation<DH::dimension,DH::space_dimension>::active_cell_iterator
- cell=dof_handler.get_tria().begin_active();
- cell!=dof_handler.get_tria().end(); ++cell)
- if (cell->subdomain_id() == subdomain)
- {
- found = true;
- break;
- }
+ cell=dof_handler.get_tria().begin_active();
+ cell!=dof_handler.get_tria().end(); ++cell)
+ if (cell->subdomain_id() == subdomain)
+ {
+ found = true;
+ break;
+ }
Assert (found == true,
- ExcMessage ("There are no cells for the given subdomain!"));
+ ExcMessage ("There are no cells for the given subdomain!"));
}
#endif
std::vector<unsigned char> component_association (dof_handler.n_dofs());
internal::extract_dofs_by_component (dof_handler, std::vector<bool>(), false,
- component_association);
+ component_association);
for (unsigned int c=0; c<dof_handler.get_fe().n_components(); ++c)
{
- for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
- if ((subdomain_association[i] == subdomain) &&
- (component_association[i] == static_cast<unsigned char>(c)))
- ++n_dofs_on_subdomain[c];
+ for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+ if ((subdomain_association[i] == subdomain) &&
+ (component_association[i] == static_cast<unsigned char>(c)))
+ ++n_dofs_on_subdomain[c];
}
}
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<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<unsigned int> &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
// to copy the result to those components
- if (!base.is_primitive() && !only_once)
- for (unsigned int dd=1;dd<d;++dd)
- dofs_per_component[target_component[component-d+dd]] =
- dofs_per_component[target_component[component-d]];
- }
- }
- }
+ if (!base.is_primitive() && !only_once)
+ for (unsigned int dd=1;dd<d;++dd)
+ dofs_per_component[target_component[component-d+dd]] =
+ dofs_per_component[target_component[component-d]];
+ }
+ }
+ }
}
template <int dim, int spacedim>
std::vector<unsigned int> &dofs_per_component,
unsigned int &component)
{
- // assert that all elements in the collection have the same
- // structure (base elements and multiplicity, components per base
- // element) and then simply call the function above
+ // assert that all elements in the collection have the same
+ // structure (base elements and multiplicity, components per base
+ // element) and then simply call the function above
for (unsigned int fe=1; fe<fe_collection.size(); ++fe)
- {
- Assert (fe_collection[fe].n_components() == fe_collection[0].n_components(),
- 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)
{
// vector as identity.
if (target_component.size()==0)
{
- target_component.resize(n_components);
- for (unsigned int i=0; i<n_components; ++i)
- target_component[i] = i;
+ target_component.resize(n_components);
+ for (unsigned int i=0; i<n_components; ++i)
+ target_component[i] = i;
}
else
Assert (target_component.size()==n_components,
- ExcDimensionMismatch(target_component.size(),
- n_components));
+ ExcDimensionMismatch(target_component.size(),
+ n_components));
const unsigned int max_component
= *std::max_element (target_component.begin(),
- target_component.end());
+ target_component.end());
const unsigned int n_target_components = max_component + 1;
AssertDimension (dofs_per_component.size(), n_target_components);
// computations
if (n_components == 1)
{
- dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
- return;
+ dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
+ return;
}
// separately. do so in parallel
std::vector<unsigned char> dofs_by_component (dof_handler.n_locally_owned_dofs());
internal::extract_dofs_by_component (dof_handler, std::vector<bool>(), false,
- dofs_by_component);
+ dofs_by_component);
// next count what we got
unsigned int component = 0;
internal::resolve_components(dof_handler.get_fe(),
dofs_by_component, target_component,
- only_once, dofs_per_component, component);
+ only_once, dofs_per_component, component);
Assert (n_components == component, ExcInternalError());
// finally sanity check. this is
// is actually primitive, so
// exclude other elements from this
Assert ((internal::all_elements_are_primitive(dof_handler.get_fe()) == false)
- ||
- (std::accumulate (dofs_per_component.begin(),
- dofs_per_component.end(), 0U)
- == dof_handler.n_locally_owned_dofs()),
- ExcInternalError());
+ ||
+ (std::accumulate (dofs_per_component.begin(),
+ dofs_per_component.end(), 0U)
+ == dof_handler.n_locally_owned_dofs()),
+ ExcInternalError());
// reduce information from all CPUs
#if defined(DEAL_II_USE_P4EST) && defined(DEAL_II_COMPILER_SUPPORTS_MPI)
const unsigned int spacedim = DH::space_dimension;
if (const parallel::distributed::Triangulation<dim,spacedim> * tria
- = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
- (&dof_handler.get_tria())))
+ = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+ (&dof_handler.get_tria())))
{
- std::vector<unsigned int> local_dof_count = dofs_per_component;
+ std::vector<unsigned int> local_dof_count = dofs_per_component;
- MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components,
- MPI_UNSIGNED, MPI_SUM, tria->get_communicator());
+ MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components,
+ MPI_UNSIGNED, MPI_SUM, tria->get_communicator());
}
#endif
}
template <int dim, int spacedim>
void
count_dofs_per_component (const DoFHandler<dim,spacedim> &dof_handler,
- std::vector<unsigned int> &dofs_per_component,
- std::vector<unsigned int> target_component)
+ std::vector<unsigned int> &dofs_per_component,
+ std::vector<unsigned int> target_component)
{
count_dofs_per_component (dof_handler, dofs_per_component,
- false, target_component);
+ false, target_component);
}
template <int dim, int spacedim>
void
compute_intergrid_weights_3 (
- const dealii::DoFHandler<dim,spacedim> &coarse_grid,
- const unsigned int coarse_component,
- const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
- const std::vector<dealii::Vector<double> > ¶meter_dofs,
- const std::vector<int> &weight_mapping,
- std::vector<std::map<unsigned int, float> > &weights,
- const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &begin,
- const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &end)
+ const dealii::DoFHandler<dim,spacedim> &coarse_grid,
+ const unsigned int coarse_component,
+ const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
+ const std::vector<dealii::Vector<double> > ¶meter_dofs,
+ const std::vector<int> &weight_mapping,
+ std::vector<std::map<unsigned int, float> > &weights,
+ const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &begin,
+ const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &end)
{
// aliases to the finite elements
// used by the dof handlers:
- const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe();
+ const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe();
// for each cell on the parameter grid:
// find out which degrees of freedom on the
// 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<unsigned int> 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);
- 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
- cell->get_dof_indices (parameter_dof_indices);
+ cell->get_dof_indices (parameter_dof_indices);
// loop over all dofs on this
// cell and check whether they
// are interesting for us
- for (unsigned int local_dof=0;
- local_dof<coarse_fe.dofs_per_cell;
- ++local_dof)
- if (coarse_fe.system_to_component_index(local_dof).first
- ==
- coarse_component)
- {
+ for (unsigned int local_dof=0;
+ local_dof<coarse_fe.dofs_per_cell;
+ ++local_dof)
+ if (coarse_fe.system_to_component_index(local_dof).first
+ ==
+ coarse_component)
+ {
// the how-many-th
// parameter is this on
// this cell?
- const unsigned int local_parameter_dof
- = coarse_fe.system_to_component_index(local_dof).second;
+ const unsigned int local_parameter_dof
+ = coarse_fe.system_to_component_index(local_dof).second;
- global_parameter_representation = 0;
+ global_parameter_representation = 0;
// distribute the representation of
// @p{local_parameter_dof} on the
// parameter grid cell @p{cell} to
// the global data space
- coarse_to_fine_grid_map[cell]->
- set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
- global_parameter_representation);
+ coarse_to_fine_grid_map[cell]->
+ set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
+ global_parameter_representation);
// now that we've got the global
// representation of each parameter
// dof, we've only got to clobber the
// different intergrid
// weights, have only one
// mutex for all of them
- static Threads::ThreadMutex mutex;
- Threads::ThreadMutex::ScopedLock lock (mutex);
- for (unsigned int i=0; i<global_parameter_representation.size(); ++i)
+ static Threads::ThreadMutex mutex;
+ Threads::ThreadMutex::ScopedLock lock (mutex);
+ for (unsigned int i=0; i<global_parameter_representation.size(); ++i)
// set this weight if it belongs
// to a parameter dof.
- if (weight_mapping[i] != -1)
- {
+ if (weight_mapping[i] != -1)
+ {
// only overwrite old
// value if not by
// zero
- if (global_parameter_representation(i) != 0)
- {
- const unsigned int wi = parameter_dof_indices[local_dof],
- wj = weight_mapping[i];
- weights[wi][wj] = global_parameter_representation(i);
- };
- }
- else
- Assert (global_parameter_representation(i) == 0,
- ExcInternalError());
- }
- }
+ if (global_parameter_representation(i) != 0)
+ {
+ const unsigned int wi = parameter_dof_indices[local_dof],
+ wj = weight_mapping[i];
+ weights[wi][wj] = global_parameter_representation(i);
+ };
+ }
+ else
+ Assert (global_parameter_representation(i) == 0,
+ ExcInternalError());
+ }
+ }
}
template <int dim, int spacedim>
void
compute_intergrid_weights_2 (
- const dealii::DoFHandler<dim,spacedim> &coarse_grid,
- const unsigned int coarse_component,
- const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
- const std::vector<dealii::Vector<double> > ¶meter_dofs,
- const std::vector<int> &weight_mapping,
- std::vector<std::map<unsigned int,float> > &weights)
+ const dealii::DoFHandler<dim,spacedim> &coarse_grid,
+ const unsigned int coarse_component,
+ const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
+ const std::vector<dealii::Vector<double> > ¶meter_dofs,
+ const std::vector<int> &weight_mapping,
+ std::vector<std::map<unsigned int,float> > &weights)
{
// simply distribute the range of
// cells to different threads
- typedef typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator active_cell_iterator;
- std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
- cell_intervals = Threads::split_range<active_cell_iterator> (coarse_grid.begin_active(),
- coarse_grid.end(),
- multithread_info.n_default_threads);
+ typedef typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
+ cell_intervals = Threads::split_range<active_cell_iterator> (coarse_grid.begin_active(),
+ coarse_grid.end(),
+ multithread_info.n_default_threads);
//TODO: use WorkStream here
- Threads::TaskGroup<> tasks;
- void (*fun_ptr) (const dealii::DoFHandler<dim,spacedim> &,
- const unsigned int ,
- const InterGridMap<dealii::DoFHandler<dim,spacedim> > &,
- const std::vector<dealii::Vector<double> > &,
- const std::vector<int> &,
- std::vector<std::map<unsigned int, float> > &,
- const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &,
- const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &)
- = &compute_intergrid_weights_3<dim>;
- for (unsigned int i=0; i<multithread_info.n_default_threads; ++i)
- tasks += Threads::new_task (fun_ptr,
- coarse_grid, coarse_component,
- coarse_to_fine_grid_map, parameter_dofs,
- weight_mapping, weights,
- cell_intervals[i].first,
- cell_intervals[i].second);
+ Threads::TaskGroup<> tasks;
+ void (*fun_ptr) (const dealii::DoFHandler<dim,spacedim> &,
+ const unsigned int ,
+ const InterGridMap<dealii::DoFHandler<dim,spacedim> > &,
+ const std::vector<dealii::Vector<double> > &,
+ const std::vector<int> &,
+ std::vector<std::map<unsigned int, float> > &,
+ const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &,
+ const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &)
+ = &compute_intergrid_weights_3<dim>;
+ for (unsigned int i=0; i<multithread_info.n_default_threads; ++i)
+ tasks += Threads::new_task (fun_ptr,
+ coarse_grid, coarse_component,
+ coarse_to_fine_grid_map, parameter_dofs,
+ weight_mapping, weights,
+ cell_intervals[i].first,
+ cell_intervals[i].second);
// wait for the tasks to finish
- tasks.join_all ();
+ tasks.join_all ();
}
template <int dim, int spacedim>
unsigned int
compute_intergrid_weights_1 (
- const dealii::DoFHandler<dim,spacedim> &coarse_grid,
- const unsigned int coarse_component,
- const dealii::DoFHandler<dim,spacedim> &fine_grid,
- const unsigned int fine_component,
- const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
- std::vector<std::map<unsigned int, float> > &weights,
- std::vector<int> &weight_mapping)
+ const dealii::DoFHandler<dim,spacedim> &coarse_grid,
+ const unsigned int coarse_component,
+ const dealii::DoFHandler<dim,spacedim> &fine_grid,
+ const unsigned int fine_component,
+ const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
+ std::vector<std::map<unsigned int, float> > &weights,
+ std::vector<int> &weight_mapping)
{
// aliases to the finite elements
// used by the dof handlers:
- const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe(),
- &fine_fe = fine_grid.get_fe();
+ const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe(),
+ &fine_fe = fine_grid.get_fe();
// global numbers of dofs
- const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
- n_fine_dofs = fine_grid.n_dofs();
+ const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
+ n_fine_dofs = fine_grid.n_dofs();
// local numbers of dofs
- const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell;
+ const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell;
// alias the number of dofs per
// cell belonging to the
// coarse_component which is to be
// the restriction of the fine
// grid:
- const unsigned int coarse_dofs_per_cell_component
- = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell;
+ const unsigned int coarse_dofs_per_cell_component
+ = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell;
// Try to find out whether the
// grids stem from the same coarse
// grid. This is a rather crude
// test, but better than nothing
- Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0),
- ExcGridsDontMatch());
+ Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0),
+ ExcGridsDontMatch());
// check whether the map correlates
// the right objects
- Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
- ExcGridsDontMatch ());
- Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
- ExcGridsDontMatch ());
+ Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
+ ExcGridsDontMatch ());
+ Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
+ ExcGridsDontMatch ());
// check whether component numbers
AssertIndexRange (fine_component, fine_fe.n_components());
// check whether respective finite
// elements are equal
- Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first)
- ==
- fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first),
- ExcFiniteElementsDontMatch());
+ Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first)
+ ==
+ fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first),
+ ExcFiniteElementsDontMatch());
#ifdef DEBUG
// if in debug mode, check whether
// the coarse grid is indeed
// coarser everywhere than the fine
// grid
- for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
- cell=coarse_grid.begin_active();
- cell != coarse_grid.end(); ++cell)
- Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(),
- ExcGridNotCoarser());
+ for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+ cell=coarse_grid.begin_active();
+ cell != coarse_grid.end(); ++cell)
+ Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(),
+ ExcGridNotCoarser());
#endif
// degree of freedom of the
// coarse-grid variable in the
// fine-grid element
- std::vector<dealii::Vector<double> >
- parameter_dofs (coarse_dofs_per_cell_component,
- dealii::Vector<double>(fine_dofs_per_cell));
+ std::vector<dealii::Vector<double> >
+ parameter_dofs (coarse_dofs_per_cell_component,
+ dealii::Vector<double>(fine_dofs_per_cell));
// for each coarse dof: find its
// position within the fine element
// and set this value to one in the
// respective vector (all other values
// are zero by construction)
- for (unsigned int local_coarse_dof=0;
- local_coarse_dof<coarse_dofs_per_cell_component;
- ++local_coarse_dof)
- for (unsigned int fine_dof=0; fine_dof<fine_fe.dofs_per_cell; ++fine_dof)
- if (fine_fe.system_to_component_index(fine_dof)
- ==
- std::make_pair (fine_component, local_coarse_dof))
- {
- parameter_dofs[local_coarse_dof](fine_dof) = 1.;
- break;
- };
+ for (unsigned int local_coarse_dof=0;
+ local_coarse_dof<coarse_dofs_per_cell_component;
+ ++local_coarse_dof)
+ for (unsigned int fine_dof=0; fine_dof<fine_fe.dofs_per_cell; ++fine_dof)
+ if (fine_fe.system_to_component_index(fine_dof)
+ ==
+ std::make_pair (fine_component, local_coarse_dof))
+ {
+ parameter_dofs[local_coarse_dof](fine_dof) = 1.;
+ break;
+ };
// find out how many DoFs there are
// on the grids belonging to the
// components we want to match
- unsigned int n_parameters_on_fine_grid=0;
- if (true)
- {
+ unsigned int n_parameters_on_fine_grid=0;
+ if (true)
+ {
// have a flag for each dof on
// the fine grid and set it
// to true if this is an
// interesting dof. finally count
// how many true's there
- 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<unsigned int> 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);
-
- weight_mapping.clear ();
- weight_mapping.resize (n_fine_dofs, -1);
-
- 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)
+ weights.clear ();
+ weights.resize (n_coarse_dofs);
+
+ weight_mapping.clear ();
+ weight_mapping.resize (n_fine_dofs, -1);
+
+ 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 this DoF is a
// parameter dof and has
// not yet been numbered,
// then do so
- if ((fine_fe.system_to_component_index(i).first == fine_component) &&
- (weight_mapping[local_dof_indices[i]] == -1))
- {
- weight_mapping[local_dof_indices[i]] = next_free_index;
- ++next_free_index;
- };
- };
+ if ((fine_fe.system_to_component_index(i).first == fine_component) &&
+ (weight_mapping[local_dof_indices[i]] == -1))
+ {
+ weight_mapping[local_dof_indices[i]] = next_free_index;
+ ++next_free_index;
+ };
+ };
- Assert (next_free_index == n_parameters_on_fine_grid,
- ExcInternalError());
- };
+ Assert (next_free_index == n_parameters_on_fine_grid,
+ ExcInternalError());
+ };
// for each cell on the parameter grid:
// you want to read more
// information on the algorithm
// used.
- compute_intergrid_weights_2 (coarse_grid, coarse_component,
- coarse_to_fine_grid_map, parameter_dofs,
- weight_mapping, weights);
+ compute_intergrid_weights_2 (coarse_grid, coarse_component,
+ coarse_to_fine_grid_map, parameter_dofs,
+ weight_mapping, weights);
// ok, now we have all weights for each
// a dof belongs to, but this at
// least tests some errors
#ifdef DEBUG
- for (unsigned int col=0; col<n_parameters_on_fine_grid; ++col)
- {
- double sum=0;
- for (unsigned int row=0; row<n_coarse_dofs; ++row)
- if (weights[row].find(col) != weights[row].end())
- sum += weights[row][col];
- Assert ((std::fabs(sum-1) < 1.e-12) ||
- ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
- };
+ for (unsigned int col=0; col<n_parameters_on_fine_grid; ++col)
+ {
+ double sum=0;
+ for (unsigned int row=0; row<n_coarse_dofs; ++row)
+ if (weights[row].find(col) != weights[row].end())
+ sum += weights[row][col];
+ Assert ((std::fabs(sum-1) < 1.e-12) ||
+ ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
+ };
#endif
- return n_parameters_on_fine_grid;
+ return n_parameters_on_fine_grid;
}
const unsigned int n_parameters_on_fine_grid
= internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
- fine_grid, fine_component,
- coarse_to_fine_grid_map,
- weights, weight_mapping);
+ fine_grid, fine_component,
+ coarse_to_fine_grid_map,
+ weights, weight_mapping);
// global numbers of dofs
const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
- n_fine_dofs = fine_grid.n_dofs();
+ n_fine_dofs = fine_grid.n_dofs();
// get an array in which we store
std::vector<bool> coarse_dof_is_parameter (coarse_grid.n_dofs());
if (true)
{
- std::vector<bool> mask (coarse_grid.get_fe().n_components(),
- false);
- mask[coarse_component] = true;
- extract_dofs (coarse_grid, mask, coarse_dof_is_parameter);
+ std::vector<bool> mask (coarse_grid.get_fe().n_components(),
+ false);
+ mask[coarse_component] = true;
+ extract_dofs (coarse_grid, mask, coarse_dof_is_parameter);
};
// now we know that the weights in
// possibly others)
std::vector<int> representants(n_coarse_dofs, -1);
for (unsigned int parameter_dof=0; parameter_dof<n_coarse_dofs;
- ++parameter_dof)
+ ++parameter_dof)
if (coarse_dof_is_parameter[parameter_dof] == true)
- {
+ {
// if this is the line of a
// parameter dof on the
// coarse grid, then it
// should have at least one
// dependent node on the fine
// grid
- Assert (weights[parameter_dof].size() > 0, ExcInternalError());
+ Assert (weights[parameter_dof].size() > 0, ExcInternalError());
// find the column where the
// representant is mentioned
- std::map<unsigned int,float>::const_iterator i = weights[parameter_dof].begin();
- for (; i!=weights[parameter_dof].end(); ++i)
- if (i->second == 1)
- break;
- Assert (i!=weights[parameter_dof].end(), ExcInternalError());
- const unsigned int column = i->first;
+ std::map<unsigned int,float>::const_iterator i = weights[parameter_dof].begin();
+ for (; i!=weights[parameter_dof].end(); ++i)
+ if (i->second == 1)
+ break;
+ Assert (i!=weights[parameter_dof].end(), ExcInternalError());
+ const unsigned int column = i->first;
// now we know in which column of
// weights the representant is, but
// we don't know its global index. get
// it using the inverse operation of
// the weight_mapping
- unsigned int global_dof=0;
- for (; global_dof<weight_mapping.size(); ++global_dof)
- if (weight_mapping[global_dof] == static_cast<int>(column))
- break;
- Assert (global_dof < weight_mapping.size(), ExcInternalError());
+ unsigned int global_dof=0;
+ for (; global_dof<weight_mapping.size(); ++global_dof)
+ if (weight_mapping[global_dof] == static_cast<int>(column))
+ break;
+ Assert (global_dof < weight_mapping.size(), ExcInternalError());
// now enter the representants global
// index into our list
- representants[parameter_dof] = global_dof;
- }
+ representants[parameter_dof] = global_dof;
+ }
else
- {
+ {
// consistency check: if this
// is no parameter dof on the
// coarse grid, then the
// respective row must be
// empty!
- Assert (weights[parameter_dof].size() == 0, ExcInternalError());
- };
+ Assert (weights[parameter_dof].size() == 0, ExcInternalError());
+ };
// dof, then we consider this dof
// to be unconstrained. otherwise,
// all other dofs are constrained
- {
- const unsigned int col = weight_mapping[global_dof];
- Assert (col < n_parameters_on_fine_grid, ExcInternalError());
-
- unsigned int first_used_row=0;
-
- {
- Assert (weights.size() > 0, ExcInternalError());
- std::map<unsigned int,float>::const_iterator
- col_entry = weights[0].end();
- for (; first_used_row<n_coarse_dofs; ++first_used_row)
- {
- col_entry = weights[first_used_row].find(col);
- if (col_entry != weights[first_used_row].end())
- break;
- }
-
- Assert (col_entry != weights[first_used_row].end(), ExcInternalError());
-
- if ((col_entry->second == 1) &&
- (representants[first_used_row] == static_cast<int>(global_dof)))
+ {
+ const unsigned int col = weight_mapping[global_dof];
+ Assert (col < n_parameters_on_fine_grid, ExcInternalError());
+
+ unsigned int first_used_row=0;
+
+ {
+ Assert (weights.size() > 0, ExcInternalError());
+ std::map<unsigned int,float>::const_iterator
+ col_entry = weights[0].end();
+ for (; first_used_row<n_coarse_dofs; ++first_used_row)
+ {
+ col_entry = weights[first_used_row].find(col);
+ if (col_entry != weights[first_used_row].end())
+ break;
+ }
+
+ Assert (col_entry != weights[first_used_row].end(), ExcInternalError());
+
+ if ((col_entry->second == 1) &&
+ (representants[first_used_row] == static_cast<int>(global_dof)))
// dof unconstrained or
// constrained to itself
// (in case this cell is
// mapped to itself, rather
// than to children of
// itself)
- continue;
- }
+ continue;
+ }
// 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));
+ };
+
+ constraints.add_entries (global_dof, constraint_line);
+ };
}
std::vector<int> weight_mapping;
internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
- fine_grid, fine_component,
- coarse_to_fine_grid_map,
- weights, weight_mapping);
+ fine_grid, fine_component,
+ coarse_to_fine_grid_map,
+ weights, weight_mapping);
// now compute the requested
// representation
const unsigned int n_global_parm_dofs
= std::count_if (weight_mapping.begin(), weight_mapping.end(),
- std::bind2nd (std::not_equal_to<int> (), -1));
+ std::bind2nd (std::not_equal_to<int> (), -1));
// first construct the inverse
// mapping of weight_mapping
std::vector<unsigned int> inverse_weight_mapping (n_global_parm_dofs,
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ DoFHandler<dim,spacedim>::invalid_dof_index);
for (unsigned int i=0; i<weight_mapping.size(); ++i)
{
- const unsigned int parameter_dof = weight_mapping[i];
+ const unsigned int parameter_dof = weight_mapping[i];
// if this global dof is a
// parameter
- if (parameter_dof != numbers::invalid_unsigned_int)
- {
- Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
- Assert ((inverse_weight_mapping[parameter_dof] == DoFHandler<dim,spacedim>::invalid_dof_index),
- ExcInternalError());
-
- inverse_weight_mapping[parameter_dof] = i;
- };
+ if (parameter_dof != numbers::invalid_unsigned_int)
+ {
+ Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
+ Assert ((inverse_weight_mapping[parameter_dof] == DoFHandler<dim,spacedim>::invalid_dof_index),
+ ExcInternalError());
+
+ inverse_weight_mapping[parameter_dof] = i;
+ };
};
// next copy over weights array
const unsigned int n_coarse_dofs = coarse_grid.n_dofs();
for (unsigned int i=0; i<n_coarse_dofs; ++i)
{
- std::map<unsigned int, float>::const_iterator j = weights[i].begin();
- for (; j!=weights[i].end(); ++j)
- {
- const unsigned int p = inverse_weight_mapping[j->first];
- Assert (p<n_rows, ExcInternalError());
-
- transfer_representation[p][i] = j->second;
- };
+ std::map<unsigned int, float>::const_iterator j = weights[i].begin();
+ for (; j!=weights[i].end(); ++j)
+ {
+ const unsigned int p = inverse_weight_mapping[j->first];
+ Assert (p<n_rows, ExcInternalError());
+
+ transfer_representation[p][i] = j->second;
+ };
};
}
template <class DH>
void
map_dof_to_boundary_indices (const DH &dof_handler,
- std::vector<unsigned int> &mapping)
+ std::vector<unsigned int> &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;
dofs_on_face.reserve (max_dofs_per_face(dof_handler));
// face which we will be visiting
// sooner or later
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
for (unsigned int f=0; f<GeometryInfo<DH::dimension>::faces_per_cell; ++f)
- if (cell->at_boundary(f))
- {
- const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
- dofs_on_face.resize (dofs_per_face);
- cell->face(f)->get_dof_indices (dofs_on_face,
- cell->active_fe_index());
- for (unsigned int i=0; i<dofs_per_face; ++i)
- if (mapping[dofs_on_face[i]] == DH::invalid_dof_index)
- mapping[dofs_on_face[i]] = next_boundary_index++;
- }
+ if (cell->at_boundary(f))
+ {
+ const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
+ dofs_on_face.resize (dofs_per_face);
+ cell->face(f)->get_dof_indices (dofs_on_face,
+ cell->active_fe_index());
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ if (mapping[dofs_on_face[i]] == DH::invalid_dof_index)
+ mapping[dofs_on_face[i]] = next_boundary_index++;
+ }
AssertDimension (next_boundary_index, dof_handler.n_boundary_dofs());
}
{
Assert (&dof_handler.get_fe() != 0, ExcNoFESelected());
Assert (boundary_indicators.find (numbers::internal_face_boundary_id) == boundary_indicators.end(),
- ExcInvalidBoundaryIndicator());
+ ExcInvalidBoundaryIndicator());
mapping.clear ();
mapping.insert (mapping.end(), dof_handler.n_dofs(),
- DH::invalid_dof_index);
+ DH::invalid_dof_index);
// return if there is nothing to do
if (boundary_indicators.size() == 0)
unsigned int next_boundary_index = 0;
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ endc = dof_handler.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_face.resize (dofs_per_face);
- cell->face(f)->get_dof_indices (dofs_on_face, cell->active_fe_index());
- for (unsigned int i=0; i<dofs_per_face; ++i)
- if (mapping[dofs_on_face[i]] == DH::invalid_dof_index)
- mapping[dofs_on_face[i]] = next_boundary_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_face.resize (dofs_per_face);
+ cell->face(f)->get_dof_indices (dofs_on_face, cell->active_fe_index());
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ if (mapping[dofs_on_face[i]] == DH::invalid_dof_index)
+ mapping[dofs_on_face[i]] = next_boundary_index++;
+ }
AssertDimension (next_boundary_index,
- dof_handler.n_boundary_dofs (boundary_indicators));
+ dof_handler.n_boundary_dofs (boundary_indicators));
}
namespace internal
const DH &dof_handler,
std::map<unsigned int,Point<DH::space_dimension> > &support_points)
{
- const unsigned int dim = DH::dimension;
- const unsigned int spacedim = DH::space_dimension;
+ 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;
+ 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()));
- }
+ 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<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();
+ // 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();
- 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<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<unsigned int,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);
for (unsigned int i=0;i<fe.n_components();++i)
{
- const unsigned int ib = fe.component_to_block_index(i);
- for (unsigned int j=0;j<fe.n_components();++j)
- {
- const unsigned int jb = fe.component_to_block_index(j);
- tables_by_block[0](ib,jb) |= table(i,j);
- }
+ const unsigned int ib = fe.component_to_block_index(i);
+ for (unsigned int j=0;j<fe.n_components();++j)
+ {
+ const unsigned int jb = fe.component_to_block_index(j);
+ tables_by_block[0](ib,jb) |= table(i,j);
+ }
}
}
for (unsigned int f=0;f<fe_collection.size();++f)
{
- const FiniteElement<dim,spacedim>& fe = fe_collection[f];
-
- const unsigned int nb = fe.n_blocks();
- tables_by_block[f].reinit(nb, nb);
- tables_by_block[f].fill(none);
- for (unsigned int i=0;i<fe.n_components();++i)
- {
- const unsigned int ib = fe.component_to_block_index(i);
- for (unsigned int j=0;j<fe.n_components();++j)
- {
- const unsigned int jb = fe.component_to_block_index(j);
- tables_by_block[f](ib,jb) |= table(i,j);
- }
- }
+ const FiniteElement<dim,spacedim>& fe = fe_collection[f];
+
+ const unsigned int nb = fe.n_blocks();
+ tables_by_block[f].reinit(nb, nb);
+ tables_by_block[f].fill(none);
+ for (unsigned int i=0;i<fe.n_components();++i)
+ {
+ const unsigned int ib = fe.component_to_block_index(i);
+ for (unsigned int j=0;j<fe.n_components();++j)
+ {
+ const unsigned int jb = fe.component_to_block_index(j);
+ tables_by_block[f](ib,jb) |= table(i,j);
+ }
+ }
}
}
template <int dim, int spacedim, template <int,int> class DH>
void
make_zero_boundary_constraints (const DH<dim, spacedim> &dof,
- ConstraintMatrix &zero_boundary_constraints,
- const std::vector<bool> &component_mask_)
+ ConstraintMatrix &zero_boundary_constraints,
+ const std::vector<bool> &component_mask_)
{
Assert ((component_mask_.size() == 0) ||
- (component_mask_.size() == dof.get_fe().n_components()),
- ExcMessage ("The number of components in the mask has to be either "
- "zero or equal to the number of components in the finite "
- "element."));
+ (component_mask_.size() == dof.get_fe().n_components()),
+ ExcMessage ("The number of components in the mask has to be either "
+ "zero or equal to the number of components in the finite "
+ "element."));
const unsigned int n_components = DoFTools::n_components (dof);
// 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_);
+ std::vector<bool> (n_components, true) :
+ component_mask_);
Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
- VectorTools::ExcNoComponentSelected());
+ VectorTools::ExcNoComponentSelected());
// a field to store the indices
std::vector<unsigned int> face_dofs;
endc = dof.end();
for (; cell!=endc; ++cell)
for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- {
- const FiniteElement<dim,spacedim> &fe = cell->get_fe();
+ ++face_no)
+ {
+ const FiniteElement<dim,spacedim> &fe = cell->get_fe();
- typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
+ typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
- // if face is on the boundary
- if (face->at_boundary ())
- {
+ // if face is on the boundary
+ if (face->at_boundary ())
+ {
// get indices and physical
// location on this face
- face_dofs.resize (fe.dofs_per_face);
- face->get_dof_indices (face_dofs, cell->active_fe_index());
+ face_dofs.resize (fe.dofs_per_face);
+ face->get_dof_indices (face_dofs, cell->active_fe_index());
// enter those dofs into the list
// that match the component
// signature.
- for (unsigned int i=0; i<face_dofs.size(); ++i)
- {
+ for (unsigned int i=0; i<face_dofs.size(); ++i)
+ {
// Find out if a dof
// has a contribution
// in this component,
// and if so, add it
// to the list
- const std::vector<bool> &nonzero_component_array
- = cell->get_fe().get_nonzero_components (i);
- bool nonzero = false;
- for (unsigned int c=0; c<n_components; ++c)
- if (nonzero_component_array[c] && component_mask[c])
- {
- nonzero = true;
- break;
- }
-
- if (nonzero)
- zero_boundary_constraints.add_line (face_dofs[i]);
- }
- }
- }
+ const std::vector<bool> &nonzero_component_array
+ = cell->get_fe().get_nonzero_components (i);
+ bool nonzero = false;
+ for (unsigned int c=0; c<n_components; ++c)
+ if (nonzero_component_array[c] && component_mask[c])
+ {
+ nonzero = true;
+ break;
+ }
+
+ if (nonzero)
+ zero_boundary_constraints.add_line (face_dofs[i]);
+ }
+ }
+ }
}
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,
- unsigned int offset)
+ Sparsity& block_list,
+ const DH& dof_handler,
+ const unsigned int level,
+ const std::vector<bool>& selected_dofs,
+ unsigned int 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>