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."));
85 std::vector<types::global_dof_index> dofs_on_this_cell;
93 (subdomain_id == cell->subdomain_id())) &&
94 cell->is_locally_owned())
96 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
97 dofs_on_this_cell.resize(dofs_per_cell);
98 cell->get_dof_indices(dofs_on_this_cell);
105 keep_constrained_dofs);
111 template <
int dim,
int spacedim,
typename number>
117 const bool keep_constrained_dofs,
144 "For distributed Triangulation objects and associated "
145 "DoFHandler objects, asking for any subdomain other than the "
146 "locally owned one does not make sense."));
152 const std::vector<Table<2, Coupling>> dof_mask
157 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.
size());
158 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
160 bool_dof_mask[f].reinit(
162 fe_collection[f].n_dofs_per_cell()));
163 bool_dof_mask[f].fill(
false);
164 for (
unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
165 for (
unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
166 if (dof_mask[f](i, j) !=
none)
167 bool_dof_mask[f](i, j) =
true;
170 std::vector<types::global_dof_index> dofs_on_this_cell(
178 (subdomain_id == cell->subdomain_id())) &&
179 cell->is_locally_owned())
182 const unsigned int dofs_per_cell =
183 fe_collection[fe_index].n_dofs_per_cell();
185 dofs_on_this_cell.resize(dofs_per_cell);
186 cell->get_dof_indices(dofs_on_this_cell);
194 keep_constrained_dofs,
195 bool_dof_mask[fe_index]);
201 template <
int dim,
int spacedim>
228 ExcMessage(
"This function can only be used with with parallel "
229 "Triangulations when the Triangulations are equal."));
235 std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
238#ifdef DEAL_II_WITH_MPI
249 Assert(this_subdomain_id ==
256 [=](
const std::pair<cell_iterator, cell_iterator> &pair) {
257 return pair.first->subdomain_id() != this_subdomain_id ||
258 pair.second->subdomain_id() != this_subdomain_id;
264 for (
const auto &cell_pair : cell_list)
266 const cell_iterator cell_row = cell_pair.first;
267 const cell_iterator cell_col = cell_pair.second;
269 if (cell_row->is_active() && cell_col->is_active())
271 const unsigned int dofs_per_cell_row =
272 cell_row->get_fe().n_dofs_per_cell();
273 const unsigned int dofs_per_cell_col =
274 cell_col->get_fe().n_dofs_per_cell();
275 std::vector<types::global_dof_index> local_dof_indices_row(
277 std::vector<types::global_dof_index> local_dof_indices_col(
279 cell_row->get_dof_indices(local_dof_indices_row);
280 cell_col->get_dof_indices(local_dof_indices_col);
281 for (
const auto &dof : local_dof_indices_row)
285 else if (cell_row->has_children())
290 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
292 for (
unsigned int i = 0; i < child_cells.size(); ++i)
295 cell_row_child = child_cells[i];
296 const unsigned int dofs_per_cell_row =
298 const unsigned int dofs_per_cell_col =
299 cell_col->get_fe().n_dofs_per_cell();
300 std::vector<types::global_dof_index> local_dof_indices_row(
302 std::vector<types::global_dof_index> local_dof_indices_col(
304 cell_row_child->get_dof_indices(local_dof_indices_row);
305 cell_col->get_dof_indices(local_dof_indices_col);
306 for (
const auto &dof : local_dof_indices_row)
316 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
318 for (
unsigned int i = 0; i < child_cells.size(); ++i)
321 cell_col_child = child_cells[i];
322 const unsigned int dofs_per_cell_row =
323 cell_row->get_fe().n_dofs_per_cell();
324 const unsigned int dofs_per_cell_col =
326 std::vector<types::global_dof_index> local_dof_indices_row(
328 std::vector<types::global_dof_index> local_dof_indices_col(
330 cell_row->get_dof_indices(local_dof_indices_row);
331 cell_col_child->get_dof_indices(local_dof_indices_col);
332 for (
const auto &dof : local_dof_indices_row)
342 template <
int dim,
int spacedim>
346 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
353 std::map<types::boundary_id, const Function<spacedim, double> *>
355 boundary_ids[0] =
nullptr;
356 boundary_ids[1] =
nullptr;
357 make_boundary_sparsity_pattern<dim, spacedim>(dof,
359 dof_to_boundary_mapping,
371 if (sparsity.
n_rows() != 0)
381 std::vector<types::global_dof_index> dofs_on_this_face;
383 std::vector<types::global_dof_index> cols;
391 for (
const unsigned int f : cell->face_indices())
392 if (cell->at_boundary(f))
394 const unsigned int dofs_per_face =
395 cell->get_fe().n_dofs_per_face(f);
396 dofs_on_this_face.resize(dofs_per_face);
397 cell->face(f)->get_dof_indices(dofs_on_this_face,
398 cell->active_fe_index());
402 for (
const auto &dof : dofs_on_this_face)
403 cols.push_back(dof_to_boundary_mapping[dof]);
407 std::sort(cols.begin(), cols.end());
408 for (
const auto &dof : dofs_on_this_face)
417 template <
int dim,
int spacedim,
typename number>
423 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
429 for (
unsigned int direction = 0; direction < 2; ++direction)
432 if (boundary_ids.find(direction) == boundary_ids.end())
439 while (!cell->at_boundary(direction))
440 cell = cell->neighbor(direction);
441 while (!cell->is_active())
442 cell = cell->child(direction);
444 const unsigned int dofs_per_vertex =
446 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
450 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
451 boundary_dof_boundary_indices[i] =
452 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
454 std::sort(boundary_dof_boundary_indices.begin(),
455 boundary_dof_boundary_indices.end());
456 for (
const auto &dof : boundary_dof_boundary_indices)
472 &(dof.
get_fe())) !=
nullptr);
485 if (sparsity.
n_rows() != 0)
495 std::vector<types::global_dof_index> dofs_on_this_face;
497 std::vector<types::global_dof_index> cols;
500 for (
const unsigned int f : cell->face_indices())
501 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
504 const unsigned int dofs_per_face =
505 cell->get_fe().n_dofs_per_face(f);
506 dofs_on_this_face.resize(dofs_per_face);
507 cell->face(f)->get_dof_indices(dofs_on_this_face,
508 cell->active_fe_index());
512 for (
const auto &dof : dofs_on_this_face)
513 cols.push_back(dof_to_boundary_mapping[dof]);
516 std::sort(cols.begin(), cols.end());
517 for (
const auto &dof : dofs_on_this_face)
526 template <
int dim,
int spacedim,
typename number>
531 const bool keep_constrained_dofs,
553 "For distributed Triangulation objects and associated "
554 "DoFHandler objects, asking for any subdomain other than the "
555 "locally owned one does not make sense."));
558 std::vector<types::global_dof_index> dofs_on_this_cell;
559 std::vector<types::global_dof_index> dofs_on_other_cell;
573 (subdomain_id == cell->subdomain_id())) &&
574 cell->is_locally_owned())
576 const unsigned int n_dofs_on_this_cell =
577 cell->get_fe().n_dofs_per_cell();
578 dofs_on_this_cell.resize(n_dofs_on_this_cell);
579 cell->get_dof_indices(dofs_on_this_cell);
586 keep_constrained_dofs);
588 for (
const unsigned int face : cell->face_indices())
592 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
593 if (!cell->at_boundary(face) || periodic_neighbor)
596 neighbor = cell->neighbor_or_periodic_neighbor(face);
603 while (neighbor->has_children())
604 neighbor = neighbor->child(face == 0 ? 1 : 0);
606 if (neighbor->has_children())
608 for (
unsigned int sub_nr = 0;
609 sub_nr != cell_face->n_active_descendants();
613 level_cell_iterator sub_neighbor =
615 cell->periodic_neighbor_child_on_subface(
617 cell->neighbor_child_on_subface(face, sub_nr);
619 const unsigned int n_dofs_on_neighbor =
621 dofs_on_other_cell.resize(n_dofs_on_neighbor);
622 sub_neighbor->get_dof_indices(dofs_on_other_cell);
628 keep_constrained_dofs);
633 keep_constrained_dofs);
637 if (sub_neighbor->subdomain_id() !=
638 cell->subdomain_id())
642 keep_constrained_dofs);
649 if ((!periodic_neighbor &&
650 cell->neighbor_is_coarser(face)) ||
651 (periodic_neighbor &&
652 cell->periodic_neighbor_is_coarser(face)))
653 if (neighbor->subdomain_id() == cell->subdomain_id())
656 const unsigned int n_dofs_on_neighbor =
658 dofs_on_other_cell.resize(n_dofs_on_neighbor);
660 neighbor->get_dof_indices(dofs_on_other_cell);
666 keep_constrained_dofs);
672 if (!cell->neighbor_or_periodic_neighbor(face)
674 (neighbor->subdomain_id() != cell->subdomain_id()))
680 keep_constrained_dofs);
681 if (neighbor->subdomain_id() != cell->subdomain_id())
685 keep_constrained_dofs);
695 template <
int dim,
int spacedim>
704 template <
int dim,
int spacedim>
721 for (
unsigned int i = 0; i < n_dofs; ++i)
723 const unsigned int ii =
729 for (
unsigned int j = 0; j < n_dofs; ++j)
731 const unsigned int jj =
737 dof_couplings(i, j) = component_couplings(ii, jj);
740 return dof_couplings;
745 template <
int dim,
int spacedim>
746 std::vector<Table<2, Coupling>>
751 std::vector<Table<2, Coupling>> return_value(fe.
size());
752 for (
unsigned int i = 0; i < fe.
size(); ++i)
766 template <
typename Iterator,
typename Iterator2>
769 const Iterator &cell,
770 const unsigned int face_no,
771 const Iterator2 &neighbor,
772 const unsigned int neighbor_face_no,
774 const std::vector<types::global_dof_index> &dofs_on_this_cell,
775 std::vector<types::global_dof_index> &dofs_on_other_cell,
779 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
780 neighbor->get_dof_indices(dofs_on_other_cell);
784 boost::container::small_vector<unsigned int, 64>
785 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
786 boost::container::small_vector<bool, 64> support_on_face_i(
787 neighbor->get_fe().n_dofs_per_cell());
788 boost::container::small_vector<bool, 64> support_on_face_e(
789 neighbor->get_fe().n_dofs_per_cell());
790 for (
unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
792 component_indices_neighbor[j] =
793 (neighbor->get_fe().is_primitive(j) ?
794 neighbor->get_fe().system_to_component_index(j).first :
796 .get_nonzero_components(j)
797 .first_selected_component());
798 support_on_face_i[j] =
799 neighbor->get_fe().has_support_on_face(j, face_no);
800 support_on_face_e[j] =
801 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
807 for (
int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
809 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
811 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
812 for (
unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
814 const unsigned int ii =
815 (fe.is_primitive(i) ?
816 fe.system_to_component_index(i).first :
817 fe.get_nonzero_components(i).first_selected_component());
820 const bool i_non_zero_i =
821 fe.has_support_on_face(i,
822 (f == 0 ? face_no : neighbor_face_no));
824 for (
unsigned int j = 0;
825 j < neighbor->get_fe().n_dofs_per_cell();
828 const bool j_non_zero_e = support_on_face_e[j];
829 const unsigned int jj = component_indices_neighbor[j];
831 Assert(jj < neighbor->get_fe().n_components(),
834 if ((flux_mask(ii, jj) ==
always) ||
835 (flux_mask(ii, jj) ==
nonzero && i_non_zero_i &&
837 cell_entries.emplace_back(dofs_i[i],
838 dofs_on_other_cell[j]);
839 if ((flux_mask(jj, ii) ==
always) ||
840 (flux_mask(jj, ii) ==
nonzero && j_non_zero_e &&
842 cell_entries.emplace_back(dofs_on_other_cell[j],
853 template <
int dim,
int spacedim,
typename number>
859 const bool keep_constrained_dofs,
865 const unsigned int)> &face_has_flux_coupling)
867 std::vector<types::global_dof_index> rows;
872 const ::hp::FECollection<dim, spacedim> &fe =
875 std::vector<types::global_dof_index> dofs_on_this_cell(
877 std::vector<types::global_dof_index> dofs_on_other_cell(
880 const unsigned int n_components = fe.n_components();
890 for (
unsigned int c1 = 0; c1 < n_components; ++c1)
891 for (
unsigned int c2 = 0; c2 < n_components; ++c2)
892 if (int_mask(c1, c2) !=
none || flux_mask(c1, c2) !=
none)
893 int_and_flux_mask(c1, c2) =
always;
897 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
902 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
903 for (
unsigned int f = 0; f < fe.size(); ++f)
905 bool_int_and_flux_dof_mask[f].reinit(
907 fe[f].n_dofs_per_cell()));
908 bool_int_and_flux_dof_mask[f].fill(
false);
909 for (
unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
910 for (
unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
911 if (int_and_flux_dof_mask[f](i, j) !=
none)
912 bool_int_and_flux_dof_mask[f](i, j) =
true;
916 for (
const auto &cell : dof.active_cell_iterators())
919 cell->is_locally_owned())
921 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
922 cell->get_dof_indices(dofs_on_this_cell);
930 keep_constrained_dofs,
931 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
934 for (
const unsigned int face : cell->face_indices())
936 const bool periodic_neighbor =
937 cell->has_periodic_neighbor(face);
939 if ((!cell->at_boundary(face)) || periodic_neighbor)
942 neighbor = cell->neighbor_or_periodic_neighbor(face);
948 if (neighbor->level() == cell->level() &&
949 neighbor->index() > cell->index() &&
950 neighbor->is_active() && neighbor->is_locally_owned())
962 if (neighbor->level() != cell->level() &&
963 ((!periodic_neighbor &&
964 !cell->neighbor_is_coarser(face)) ||
965 (periodic_neighbor &&
966 !cell->periodic_neighbor_is_coarser(face))) &&
967 neighbor->is_locally_owned())
970 if (!face_has_flux_coupling(cell, face))
973 const unsigned int neighbor_face_no =
975 cell->periodic_neighbor_face_no(face) :
976 cell->neighbor_face_no(face);
985 while (neighbor->has_children())
986 neighbor = neighbor->child(face == 0 ? 1 : 0);
988 if (neighbor->has_children())
990 for (
unsigned int sub_nr = 0;
991 sub_nr != cell->face(face)->n_children();
995 level_cell_iterator sub_neighbor =
997 cell->periodic_neighbor_child_on_subface(
999 cell->neighbor_child_on_subface(face,
1001 add_cell_entries(cell,
1012 add_cell_entries(cell,
1023 cell_entries.clear();
1032 template <
int dim,
int spacedim>
1042 const bool keep_constrained_dofs =
true;
1047 keep_constrained_dofs,
1051 internal::always_couple_on_faces<dim, spacedim>);
1056 template <
int dim,
int spacedim,
typename number>
1062 const bool keep_constrained_dofs,
1066 const std::function<
1068 const unsigned int)> &face_has_flux_coupling)
1081 Assert(int_mask.n_rows() == n_comp,
1083 Assert(int_mask.n_cols() == n_comp,
1085 Assert(flux_mask.n_rows() == n_comp,
1087 Assert(flux_mask.n_cols() == n_comp,
1100 "For distributed Triangulation objects and associated "
1101 "DoFHandler objects, asking for any subdomain other than the "
1102 "locally owned one does not make sense."));
1106 face_has_flux_coupling,
1108 "The function which specifies if a flux coupling occurs over a given "
1111 internal::make_flux_sparsity_pattern(dof,
1114 keep_constrained_dofs,
1118 face_has_flux_coupling);
1126#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