54 template <
int dim,
int spacedim,
typename number>
59 const bool keep_constrained_dofs,
80 "For distributed Triangulation objects and associated "
81 "DoFHandler objects, asking for any subdomain other than the "
82 "locally owned one does not make sense."));
86 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.size());
87 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
89 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
92 std::vector<types::global_dof_index> dofs_on_this_cell;
100 (subdomain_id == cell->subdomain_id())) &&
101 cell->is_locally_owned())
103 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
104 dofs_on_this_cell.resize(dofs_per_cell);
105 cell->get_dof_indices(dofs_on_this_cell);
111 if (fe_dof_mask[fe_index].empty())
114 keep_constrained_dofs);
118 keep_constrained_dofs,
119 fe_dof_mask[fe_index]);
125 template <
int dim,
int spacedim,
typename number>
131 const bool keep_constrained_dofs,
158 "For distributed Triangulation objects and associated "
159 "DoFHandler objects, asking for any subdomain other than the "
160 "locally owned one does not make sense."));
166 const std::vector<Table<2, Coupling>> dof_mask
169 std::vector<Table<2, bool>> fe_dof_mask(fe_collection.
size());
170 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
172 fe_dof_mask[f] = fe_collection[f].get_local_dof_sparsity_pattern();
177 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.
size());
178 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
180 bool_dof_mask[f].reinit(
182 fe_collection[f].n_dofs_per_cell()));
183 bool_dof_mask[f].fill(
false);
184 for (
unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
185 for (
unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
186 if (dof_mask[f](i, j) !=
none &&
187 (fe_dof_mask[f].empty() || fe_dof_mask[f](i, j)))
188 bool_dof_mask[f](i, j) =
true;
191 std::vector<types::global_dof_index> dofs_on_this_cell(
199 (subdomain_id == cell->subdomain_id())) &&
200 cell->is_locally_owned())
203 const unsigned int dofs_per_cell =
204 fe_collection[fe_index].n_dofs_per_cell();
206 dofs_on_this_cell.resize(dofs_per_cell);
207 cell->get_dof_indices(dofs_on_this_cell);
215 keep_constrained_dofs,
216 bool_dof_mask[fe_index]);
222 template <
int dim,
int spacedim>
249 ExcMessage(
"This function can only be used with with parallel "
250 "Triangulations when the Triangulations are equal."));
256 std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
259#ifdef DEAL_II_WITH_MPI
270 Assert(this_subdomain_id ==
277 [=](
const std::pair<cell_iterator, cell_iterator> &pair) {
278 return pair.first->subdomain_id() != this_subdomain_id ||
279 pair.second->subdomain_id() != this_subdomain_id;
285 for (
const auto &cell_pair : cell_list)
287 const cell_iterator cell_row = cell_pair.first;
288 const cell_iterator cell_col = cell_pair.second;
290 if (cell_row->is_active() && cell_col->is_active())
292 const unsigned int dofs_per_cell_row =
293 cell_row->get_fe().n_dofs_per_cell();
294 const unsigned int dofs_per_cell_col =
295 cell_col->get_fe().n_dofs_per_cell();
296 std::vector<types::global_dof_index> local_dof_indices_row(
298 std::vector<types::global_dof_index> local_dof_indices_col(
300 cell_row->get_dof_indices(local_dof_indices_row);
301 cell_col->get_dof_indices(local_dof_indices_col);
302 for (
const auto &dof : local_dof_indices_row)
306 else if (cell_row->has_children())
311 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
313 for (
unsigned int i = 0; i < child_cells.size(); ++i)
316 cell_row_child = child_cells[i];
317 const unsigned int dofs_per_cell_row =
319 const unsigned int dofs_per_cell_col =
320 cell_col->get_fe().n_dofs_per_cell();
321 std::vector<types::global_dof_index> local_dof_indices_row(
323 std::vector<types::global_dof_index> local_dof_indices_col(
325 cell_row_child->get_dof_indices(local_dof_indices_row);
326 cell_col->get_dof_indices(local_dof_indices_col);
327 for (
const auto &dof : local_dof_indices_row)
337 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
339 for (
unsigned int i = 0; i < child_cells.size(); ++i)
342 cell_col_child = child_cells[i];
343 const unsigned int dofs_per_cell_row =
344 cell_row->get_fe().n_dofs_per_cell();
345 const unsigned int dofs_per_cell_col =
347 std::vector<types::global_dof_index> local_dof_indices_row(
349 std::vector<types::global_dof_index> local_dof_indices_col(
351 cell_row->get_dof_indices(local_dof_indices_row);
352 cell_col_child->get_dof_indices(local_dof_indices_col);
353 for (
const auto &dof : local_dof_indices_row)
363 template <
int dim,
int spacedim>
367 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
374 std::map<types::boundary_id, const Function<spacedim, double> *>
376 boundary_ids[0] =
nullptr;
377 boundary_ids[1] =
nullptr;
378 make_boundary_sparsity_pattern<dim, spacedim>(dof,
380 dof_to_boundary_mapping,
392 if (sparsity.
n_rows() != 0)
402 std::vector<types::global_dof_index> dofs_on_this_face;
404 std::vector<types::global_dof_index> cols;
412 for (
const unsigned int f : cell->face_indices())
413 if (cell->at_boundary(f))
415 const unsigned int dofs_per_face =
416 cell->get_fe().n_dofs_per_face(f);
417 dofs_on_this_face.resize(dofs_per_face);
418 cell->face(f)->get_dof_indices(dofs_on_this_face,
419 cell->active_fe_index());
423 for (
const auto &dof : dofs_on_this_face)
424 cols.push_back(dof_to_boundary_mapping[dof]);
428 std::sort(cols.begin(), cols.end());
429 for (
const auto &dof : dofs_on_this_face)
438 template <
int dim,
int spacedim,
typename number>
444 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
450 for (
unsigned int direction = 0; direction < 2; ++direction)
453 if (boundary_ids.find(direction) == boundary_ids.end())
460 while (!cell->at_boundary(direction))
461 cell = cell->neighbor(direction);
462 while (!cell->is_active())
463 cell = cell->child(direction);
465 const unsigned int dofs_per_vertex =
467 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
471 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
472 boundary_dof_boundary_indices[i] =
473 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
475 std::sort(boundary_dof_boundary_indices.begin(),
476 boundary_dof_boundary_indices.end());
477 for (
const auto &dof : boundary_dof_boundary_indices)
493 &(dof.
get_fe())) !=
nullptr);
506 if (sparsity.
n_rows() != 0)
516 std::vector<types::global_dof_index> dofs_on_this_face;
518 std::vector<types::global_dof_index> cols;
521 for (
const unsigned int f : cell->face_indices())
522 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
525 const unsigned int dofs_per_face =
526 cell->get_fe().n_dofs_per_face(f);
527 dofs_on_this_face.resize(dofs_per_face);
528 cell->face(f)->get_dof_indices(dofs_on_this_face,
529 cell->active_fe_index());
533 for (
const auto &dof : dofs_on_this_face)
534 cols.push_back(dof_to_boundary_mapping[dof]);
537 std::sort(cols.begin(), cols.end());
538 for (
const auto &dof : dofs_on_this_face)
547 template <
int dim,
int spacedim,
typename number>
552 const bool keep_constrained_dofs,
574 "For distributed Triangulation objects and associated "
575 "DoFHandler objects, asking for any subdomain other than the "
576 "locally owned one does not make sense."));
579 std::vector<types::global_dof_index> dofs_on_this_cell;
580 std::vector<types::global_dof_index> dofs_on_other_cell;
594 (subdomain_id == cell->subdomain_id())) &&
595 cell->is_locally_owned())
597 const unsigned int n_dofs_on_this_cell =
598 cell->get_fe().n_dofs_per_cell();
599 dofs_on_this_cell.resize(n_dofs_on_this_cell);
600 cell->get_dof_indices(dofs_on_this_cell);
607 keep_constrained_dofs);
609 for (
const unsigned int face : cell->face_indices())
613 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
614 if (!cell->at_boundary(face) || periodic_neighbor)
617 neighbor = cell->neighbor_or_periodic_neighbor(face);
624 while (neighbor->has_children())
625 neighbor = neighbor->child(face == 0 ? 1 : 0);
627 if (neighbor->has_children())
629 for (
unsigned int sub_nr = 0;
630 sub_nr != cell_face->n_active_descendants();
634 level_cell_iterator sub_neighbor =
636 cell->periodic_neighbor_child_on_subface(
638 cell->neighbor_child_on_subface(face, sub_nr);
640 const unsigned int n_dofs_on_neighbor =
642 dofs_on_other_cell.resize(n_dofs_on_neighbor);
643 sub_neighbor->get_dof_indices(dofs_on_other_cell);
649 keep_constrained_dofs);
654 keep_constrained_dofs);
658 if (sub_neighbor->subdomain_id() !=
659 cell->subdomain_id())
663 keep_constrained_dofs);
670 if ((!periodic_neighbor &&
671 cell->neighbor_is_coarser(face)) ||
672 (periodic_neighbor &&
673 cell->periodic_neighbor_is_coarser(face)))
674 if (neighbor->subdomain_id() == cell->subdomain_id())
677 const unsigned int n_dofs_on_neighbor =
679 dofs_on_other_cell.resize(n_dofs_on_neighbor);
681 neighbor->get_dof_indices(dofs_on_other_cell);
687 keep_constrained_dofs);
693 if (!cell->neighbor_or_periodic_neighbor(face)
695 (neighbor->subdomain_id() != cell->subdomain_id()))
701 keep_constrained_dofs);
702 if (neighbor->subdomain_id() != cell->subdomain_id())
706 keep_constrained_dofs);
716 template <
int dim,
int spacedim>
725 template <
int dim,
int spacedim>
742 for (
unsigned int i = 0; i < n_dofs; ++i)
744 const unsigned int ii =
750 for (
unsigned int j = 0; j < n_dofs; ++j)
752 const unsigned int jj =
758 dof_couplings(i, j) = component_couplings(ii, jj);
761 return dof_couplings;
766 template <
int dim,
int spacedim>
767 std::vector<Table<2, Coupling>>
772 std::vector<Table<2, Coupling>> return_value(fe.
size());
773 for (
unsigned int i = 0; i < fe.
size(); ++i)
787 template <
typename Iterator,
typename Iterator2>
790 const Iterator &cell,
791 const unsigned int face_no,
792 const Iterator2 &neighbor,
793 const unsigned int neighbor_face_no,
795 const std::vector<types::global_dof_index> &dofs_on_this_cell,
796 std::vector<types::global_dof_index> &dofs_on_other_cell,
800 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
801 neighbor->get_dof_indices(dofs_on_other_cell);
805 boost::container::small_vector<unsigned int, 64>
806 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
807 boost::container::small_vector<bool, 64> support_on_face_i(
808 neighbor->get_fe().n_dofs_per_cell());
809 boost::container::small_vector<bool, 64> support_on_face_e(
810 neighbor->get_fe().n_dofs_per_cell());
811 for (
unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
813 component_indices_neighbor[j] =
814 (neighbor->get_fe().is_primitive(j) ?
815 neighbor->get_fe().system_to_component_index(j).first :
817 .get_nonzero_components(j)
818 .first_selected_component());
819 support_on_face_i[j] =
820 neighbor->get_fe().has_support_on_face(j, face_no);
821 support_on_face_e[j] =
822 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
828 for (
int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
830 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
832 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
833 for (
unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
835 const unsigned int ii =
836 (fe.is_primitive(i) ?
837 fe.system_to_component_index(i).first :
838 fe.get_nonzero_components(i).first_selected_component());
841 const bool i_non_zero_i =
842 fe.has_support_on_face(i,
843 (f == 0 ? face_no : neighbor_face_no));
845 for (
unsigned int j = 0;
846 j < neighbor->get_fe().n_dofs_per_cell();
849 const bool j_non_zero_e = support_on_face_e[j];
850 const unsigned int jj = component_indices_neighbor[j];
852 Assert(jj < neighbor->get_fe().n_components(),
855 if ((flux_mask(ii, jj) ==
always) ||
856 (flux_mask(ii, jj) ==
nonzero && i_non_zero_i &&
858 cell_entries.emplace_back(dofs_i[i],
859 dofs_on_other_cell[j]);
860 if ((flux_mask(jj, ii) ==
always) ||
861 (flux_mask(jj, ii) ==
nonzero && j_non_zero_e &&
863 cell_entries.emplace_back(dofs_on_other_cell[j],
874 template <
int dim,
int spacedim,
typename number>
880 const bool keep_constrained_dofs,
886 const unsigned int)> &face_has_flux_coupling)
888 std::vector<types::global_dof_index> rows;
893 const ::hp::FECollection<dim, spacedim> &fe =
896 std::vector<types::global_dof_index> dofs_on_this_cell(
898 std::vector<types::global_dof_index> dofs_on_other_cell(
901 const unsigned int n_components = fe.n_components();
911 for (
unsigned int c1 = 0; c1 < n_components; ++c1)
912 for (
unsigned int c2 = 0; c2 < n_components; ++c2)
913 if (int_mask(c1, c2) !=
none || flux_mask(c1, c2) !=
none)
914 int_and_flux_mask(c1, c2) =
always;
918 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
923 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
924 for (
unsigned int f = 0; f < fe.size(); ++f)
926 bool_int_and_flux_dof_mask[f].reinit(
928 fe[f].n_dofs_per_cell()));
929 bool_int_and_flux_dof_mask[f].fill(
false);
930 for (
unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
931 for (
unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
932 if (int_and_flux_dof_mask[f](i, j) !=
none)
933 bool_int_and_flux_dof_mask[f](i, j) =
true;
937 for (
const auto &cell : dof.active_cell_iterators())
940 cell->is_locally_owned())
942 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
943 cell->get_dof_indices(dofs_on_this_cell);
951 keep_constrained_dofs,
952 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
955 for (
const unsigned int face : cell->face_indices())
957 const bool periodic_neighbor =
958 cell->has_periodic_neighbor(face);
960 if ((!cell->at_boundary(face)) || periodic_neighbor)
963 neighbor = cell->neighbor_or_periodic_neighbor(face);
969 if (neighbor->level() == cell->level() &&
970 neighbor->index() > cell->index() &&
971 neighbor->is_active() && neighbor->is_locally_owned())
983 if (neighbor->level() != cell->level() &&
984 ((!periodic_neighbor &&
985 !cell->neighbor_is_coarser(face)) ||
986 (periodic_neighbor &&
987 !cell->periodic_neighbor_is_coarser(face))) &&
988 neighbor->is_locally_owned())
991 if (!face_has_flux_coupling(cell, face))
994 const unsigned int neighbor_face_no =
996 cell->periodic_neighbor_face_no(face) :
997 cell->neighbor_face_no(face);
1006 while (neighbor->has_children())
1007 neighbor = neighbor->child(face == 0 ? 1 : 0);
1009 if (neighbor->has_children())
1011 for (
unsigned int sub_nr = 0;
1012 sub_nr != cell->face(face)->n_children();
1016 level_cell_iterator sub_neighbor =
1018 cell->periodic_neighbor_child_on_subface(
1020 cell->neighbor_child_on_subface(face,
1022 add_cell_entries(cell,
1033 add_cell_entries(cell,
1044 cell_entries.clear();
1053 template <
int dim,
int spacedim>
1063 const bool keep_constrained_dofs =
true;
1068 keep_constrained_dofs,
1072 internal::always_couple_on_faces<dim, spacedim>);
1077 template <
int dim,
int spacedim,
typename number>
1083 const bool keep_constrained_dofs,
1087 const std::function<
1089 const unsigned int)> &face_has_flux_coupling)
1102 Assert(int_mask.n_rows() == n_comp,
1104 Assert(int_mask.n_cols() == n_comp,
1106 Assert(flux_mask.n_rows() == n_comp,
1108 Assert(flux_mask.n_cols() == n_comp,
1121 "For distributed Triangulation objects and associated "
1122 "DoFHandler objects, asking for any subdomain other than the "
1123 "locally owned one does not make sense."));
1127 face_has_flux_coupling,
1129 "The function which specifies if a flux coupling occurs over a given "
1132 internal::make_flux_sparsity_pattern(dof,
1135 keep_constrained_dofs,
1139 face_has_flux_coupling);
1147#include "dof_tools_sparsity.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
typename LevelSelector::cell_iterator level_cell_iterator
cell_iterator begin(const unsigned int level=0) const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
types::global_dof_index size_type
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false)=0
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int size() const
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
void make_boundary_sparsity_pattern(const DoFHandler< dim, spacedim > &dof, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPatternBase &sparsity_pattern)
const types::boundary_id internal_face_boundary_id
const types::subdomain_id invalid_subdomain_id
const types::global_dof_index invalid_dof_index
unsigned int subdomain_id
unsigned short int fe_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation