55 template <
int dim,
int spacedim,
typename number>
60 const bool keep_constrained_dofs,
81 "For distributed Triangulation objects and associated "
82 "DoFHandler objects, asking for any subdomain other than the "
83 "locally owned one does not make sense."));
86 std::vector<types::global_dof_index> dofs_on_this_cell;
95 cell->is_locally_owned())
97 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
98 dofs_on_this_cell.resize(dofs_per_cell);
99 cell->get_dof_indices(dofs_on_this_cell);
106 keep_constrained_dofs);
112 template <
int dim,
int spacedim,
typename number>
118 const bool keep_constrained_dofs,
145 "For distributed Triangulation objects and associated "
146 "DoFHandler objects, asking for any subdomain other than the "
147 "locally owned one does not make sense."));
153 const std::vector<Table<2, Coupling>> dof_mask
158 std::vector<Table<2, bool>> bool_dof_mask(fe_collection.
size());
159 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
161 bool_dof_mask[f].reinit(
163 fe_collection[f].n_dofs_per_cell()));
164 bool_dof_mask[f].fill(
false);
165 for (
unsigned int i = 0; i < fe_collection[f].n_dofs_per_cell(); ++i)
166 for (
unsigned int j = 0; j < fe_collection[f].n_dofs_per_cell(); ++j)
167 if (dof_mask[f](i, j) !=
none)
168 bool_dof_mask[f](i, j) =
true;
171 std::vector<types::global_dof_index> dofs_on_this_cell(
180 cell->is_locally_owned())
183 const unsigned int dofs_per_cell =
184 fe_collection[
fe_index].n_dofs_per_cell();
186 dofs_on_this_cell.resize(dofs_per_cell);
187 cell->get_dof_indices(dofs_on_this_cell);
195 keep_constrained_dofs,
202 template <
int dim,
int spacedim>
229 ExcMessage(
"This function can only be used with with parallel "
230 "Triangulations when the Triangulations are equal."));
236 std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
239 #ifdef DEAL_II_WITH_MPI
250 Assert(this_subdomain_id ==
257 [=](
const std::pair<cell_iterator, cell_iterator> &pair) {
258 return pair.first->subdomain_id() != this_subdomain_id ||
259 pair.second->subdomain_id() != this_subdomain_id;
265 for (
const auto &cell_pair : cell_list)
267 const cell_iterator cell_row = cell_pair.first;
268 const cell_iterator cell_col = cell_pair.second;
270 if (cell_row->is_active() && cell_col->is_active())
272 const unsigned int dofs_per_cell_row =
273 cell_row->get_fe().n_dofs_per_cell();
274 const unsigned int dofs_per_cell_col =
275 cell_col->get_fe().n_dofs_per_cell();
276 std::vector<types::global_dof_index> local_dof_indices_row(
278 std::vector<types::global_dof_index> local_dof_indices_col(
280 cell_row->get_dof_indices(local_dof_indices_row);
281 cell_col->get_dof_indices(local_dof_indices_col);
282 for (
const auto &dof : local_dof_indices_row)
286 else if (cell_row->has_children())
291 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
293 for (
unsigned int i = 0; i < child_cells.size(); ++i)
296 cell_row_child = child_cells[i];
297 const unsigned int dofs_per_cell_row =
299 const unsigned int dofs_per_cell_col =
300 cell_col->get_fe().n_dofs_per_cell();
301 std::vector<types::global_dof_index> local_dof_indices_row(
303 std::vector<types::global_dof_index> local_dof_indices_col(
305 cell_row_child->get_dof_indices(local_dof_indices_row);
306 cell_col->get_dof_indices(local_dof_indices_col);
307 for (
const auto &dof : local_dof_indices_row)
317 GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
319 for (
unsigned int i = 0; i < child_cells.size(); ++i)
322 cell_col_child = child_cells[i];
323 const unsigned int dofs_per_cell_row =
325 const unsigned int dofs_per_cell_col =
327 std::vector<types::global_dof_index> local_dof_indices_row(
329 std::vector<types::global_dof_index> local_dof_indices_col(
331 cell_row->get_dof_indices(local_dof_indices_row);
332 cell_col_child->get_dof_indices(local_dof_indices_col);
333 for (
const auto &dof : local_dof_indices_row)
343 template <
int dim,
int spacedim>
347 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
354 std::map<types::boundary_id, const Function<spacedim, double> *>
356 boundary_ids[0] =
nullptr;
357 boundary_ids[1] =
nullptr;
358 make_boundary_sparsity_pattern<dim, spacedim>(dof,
360 dof_to_boundary_mapping,
372 if (sparsity.
n_rows() != 0)
382 std::vector<types::global_dof_index> dofs_on_this_face;
384 std::vector<types::global_dof_index> cols;
392 for (
const unsigned int f : cell->face_indices())
393 if (cell->at_boundary(f))
395 const unsigned int dofs_per_face =
396 cell->get_fe().n_dofs_per_face(f);
397 dofs_on_this_face.resize(dofs_per_face);
398 cell->face(f)->get_dof_indices(dofs_on_this_face,
399 cell->active_fe_index());
403 for (
const auto &dof : dofs_on_this_face)
404 cols.push_back(dof_to_boundary_mapping[dof]);
408 std::sort(cols.begin(), cols.end());
409 for (
const auto &dof : dofs_on_this_face)
418 template <
int dim,
int spacedim,
typename number>
424 const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
430 for (
unsigned int direction = 0; direction < 2; ++direction)
433 if (boundary_ids.find(direction) == boundary_ids.end())
440 while (!cell->at_boundary(direction))
441 cell = cell->neighbor(direction);
442 while (!cell->is_active())
443 cell = cell->child(direction);
445 const unsigned int dofs_per_vertex =
447 std::vector<types::global_dof_index> boundary_dof_boundary_indices(
451 for (
unsigned int i = 0; i < dofs_per_vertex; ++i)
452 boundary_dof_boundary_indices[i] =
453 dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
455 std::sort(boundary_dof_boundary_indices.begin(),
456 boundary_dof_boundary_indices.end());
457 for (
const auto &dof : boundary_dof_boundary_indices)
473 &(dof.
get_fe())) !=
nullptr);
486 if (sparsity.
n_rows() != 0)
496 std::vector<types::global_dof_index> dofs_on_this_face;
498 std::vector<types::global_dof_index> cols;
501 for (
const unsigned int f : cell->face_indices())
502 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
505 const unsigned int dofs_per_face =
506 cell->get_fe().n_dofs_per_face(f);
507 dofs_on_this_face.resize(dofs_per_face);
508 cell->face(f)->get_dof_indices(dofs_on_this_face,
509 cell->active_fe_index());
513 for (
const auto &dof : dofs_on_this_face)
514 cols.push_back(dof_to_boundary_mapping[dof]);
517 std::sort(cols.begin(), cols.end());
518 for (
const auto &dof : dofs_on_this_face)
527 template <
int dim,
int spacedim,
typename number>
532 const bool keep_constrained_dofs,
554 "For distributed Triangulation objects and associated "
555 "DoFHandler objects, asking for any subdomain other than the "
556 "locally owned one does not make sense."));
559 std::vector<types::global_dof_index> dofs_on_this_cell;
560 std::vector<types::global_dof_index> dofs_on_other_cell;
575 cell->is_locally_owned())
577 const unsigned int n_dofs_on_this_cell =
578 cell->get_fe().n_dofs_per_cell();
579 dofs_on_this_cell.resize(n_dofs_on_this_cell);
580 cell->get_dof_indices(dofs_on_this_cell);
587 keep_constrained_dofs);
589 for (
const unsigned int face : cell->face_indices())
593 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
594 if (!cell->at_boundary(face) || periodic_neighbor)
597 neighbor = cell->neighbor_or_periodic_neighbor(face);
604 while (neighbor->has_children())
605 neighbor = neighbor->child(face == 0 ? 1 : 0);
607 if (neighbor->has_children())
609 for (
unsigned int sub_nr = 0;
610 sub_nr != cell_face->n_active_descendants();
614 level_cell_iterator sub_neighbor =
616 cell->periodic_neighbor_child_on_subface(
618 cell->neighbor_child_on_subface(face, sub_nr);
620 const unsigned int n_dofs_on_neighbor =
622 dofs_on_other_cell.resize(n_dofs_on_neighbor);
623 sub_neighbor->get_dof_indices(dofs_on_other_cell);
629 keep_constrained_dofs);
634 keep_constrained_dofs);
638 if (sub_neighbor->subdomain_id() !=
639 cell->subdomain_id())
643 keep_constrained_dofs);
650 if ((!periodic_neighbor &&
651 cell->neighbor_is_coarser(face)) ||
652 (periodic_neighbor &&
653 cell->periodic_neighbor_is_coarser(face)))
654 if (neighbor->subdomain_id() == cell->subdomain_id())
657 const unsigned int n_dofs_on_neighbor =
659 dofs_on_other_cell.resize(n_dofs_on_neighbor);
661 neighbor->get_dof_indices(dofs_on_other_cell);
667 keep_constrained_dofs);
673 if (!cell->neighbor_or_periodic_neighbor(face)
675 (neighbor->subdomain_id() != cell->subdomain_id()))
681 keep_constrained_dofs);
682 if (neighbor->subdomain_id() != cell->subdomain_id())
686 keep_constrained_dofs);
696 template <
int dim,
int spacedim>
705 template <
int dim,
int spacedim>
722 for (
unsigned int i = 0; i < n_dofs; ++i)
724 const unsigned int ii =
730 for (
unsigned int j = 0; j < n_dofs; ++j)
732 const unsigned int jj =
738 dof_couplings(i, j) = component_couplings(ii, jj);
741 return dof_couplings;
746 template <
int dim,
int spacedim>
747 std::vector<Table<2, Coupling>>
752 std::vector<Table<2, Coupling>> return_value(fe.
size());
753 for (
unsigned int i = 0; i < fe.
size(); ++i)
767 template <
typename Iterator,
typename Iterator2>
770 const Iterator &cell,
771 const unsigned int face_no,
772 const Iterator2 &neighbor,
773 const unsigned int neighbor_face_no,
775 const std::vector<types::global_dof_index> &dofs_on_this_cell,
776 std::vector<types::global_dof_index> &dofs_on_other_cell,
780 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
781 neighbor->get_dof_indices(dofs_on_other_cell);
785 boost::container::small_vector<unsigned int, 64>
786 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
787 boost::container::small_vector<bool, 64> support_on_face_i(
788 neighbor->get_fe().n_dofs_per_cell());
789 boost::container::small_vector<bool, 64> support_on_face_e(
790 neighbor->get_fe().n_dofs_per_cell());
791 for (
unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
793 component_indices_neighbor[j] =
794 (neighbor->get_fe().is_primitive(j) ?
795 neighbor->get_fe().system_to_component_index(j).first :
797 .get_nonzero_components(j)
798 .first_selected_component());
799 support_on_face_i[j] =
800 neighbor->get_fe().has_support_on_face(j, face_no);
801 support_on_face_e[j] =
802 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
808 for (
int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
810 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
812 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
813 for (
unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
815 const unsigned int ii =
816 (fe.is_primitive(i) ?
817 fe.system_to_component_index(i).first :
818 fe.get_nonzero_components(i).first_selected_component());
821 const bool i_non_zero_i =
822 fe.has_support_on_face(i,
823 (f == 0 ? face_no : neighbor_face_no));
825 for (
unsigned int j = 0;
826 j < neighbor->get_fe().n_dofs_per_cell();
829 const bool j_non_zero_e = support_on_face_e[j];
830 const unsigned int jj = component_indices_neighbor[j];
832 Assert(jj < neighbor->get_fe().n_components(),
835 if ((flux_mask(ii, jj) ==
always) ||
836 (flux_mask(ii, jj) ==
nonzero && i_non_zero_i &&
838 cell_entries.emplace_back(dofs_i[i],
839 dofs_on_other_cell[j]);
840 if ((flux_mask(jj, ii) ==
always) ||
841 (flux_mask(jj, ii) ==
nonzero && j_non_zero_e &&
843 cell_entries.emplace_back(dofs_on_other_cell[j],
854 template <
int dim,
int spacedim,
typename number>
860 const bool keep_constrained_dofs,
866 const unsigned int)> &face_has_flux_coupling)
868 std::vector<types::global_dof_index> rows;
873 const ::hp::FECollection<dim, spacedim> &fe =
876 std::vector<types::global_dof_index> dofs_on_this_cell(
878 std::vector<types::global_dof_index> dofs_on_other_cell(
881 const unsigned int n_components = fe.n_components();
891 for (
unsigned int c1 = 0; c1 < n_components; ++c1)
892 for (
unsigned int c2 = 0; c2 < n_components; ++c2)
893 if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
894 int_and_flux_mask(c1, c2) =
always;
898 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
903 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
904 for (
unsigned int f = 0; f < fe.size(); ++f)
906 bool_int_and_flux_dof_mask[f].reinit(
908 fe[f].n_dofs_per_cell()));
909 bool_int_and_flux_dof_mask[f].fill(
false);
910 for (
unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
911 for (
unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
912 if (int_and_flux_dof_mask[f](i, j) != none)
913 bool_int_and_flux_dof_mask[f](i, j) =
true;
920 cell->is_locally_owned())
922 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
923 cell->get_dof_indices(dofs_on_this_cell);
931 keep_constrained_dofs,
932 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
935 for (
const unsigned int face : cell->face_indices())
937 const bool periodic_neighbor =
938 cell->has_periodic_neighbor(face);
940 if ((!cell->at_boundary(face)) || periodic_neighbor)
943 neighbor = cell->neighbor_or_periodic_neighbor(face);
949 if (neighbor->level() == cell->level() &&
950 neighbor->index() > cell->index() &&
951 neighbor->is_active() && neighbor->is_locally_owned())
963 if (neighbor->level() != cell->level() &&
964 ((!periodic_neighbor &&
965 !cell->neighbor_is_coarser(face)) ||
966 (periodic_neighbor &&
967 !cell->periodic_neighbor_is_coarser(face))) &&
968 neighbor->is_locally_owned())
971 if (!face_has_flux_coupling(cell, face))
974 const unsigned int neighbor_face_no =
976 cell->periodic_neighbor_face_no(face) :
977 cell->neighbor_face_no(face);
986 while (neighbor->has_children())
987 neighbor = neighbor->child(face == 0 ? 1 : 0);
989 if (neighbor->has_children())
991 for (
unsigned int sub_nr = 0;
992 sub_nr != cell->face(face)->n_children();
996 level_cell_iterator sub_neighbor =
998 cell->periodic_neighbor_child_on_subface(
1000 cell->neighbor_child_on_subface(face,
1002 add_cell_entries(cell,
1013 add_cell_entries(cell,
1024 cell_entries.clear();
1033 template <
int dim,
int spacedim>
1043 const bool keep_constrained_dofs =
true;
1048 keep_constrained_dofs,
1052 internal::always_couple_on_faces<dim, spacedim>);
1057 template <
int dim,
int spacedim,
typename number>
1063 const bool keep_constrained_dofs,
1067 const std::function<
1069 const unsigned int)> &face_has_flux_coupling)
1082 Assert(int_mask.n_rows() == n_comp,
1084 Assert(int_mask.n_cols() == n_comp,
1086 Assert(flux_mask.n_rows() == n_comp,
1088 Assert(flux_mask.n_cols() == n_comp,
1101 "For distributed Triangulation objects and associated "
1102 "DoFHandler objects, asking for any subdomain other than the "
1103 "locally owned one does not make sense."));
1107 face_has_flux_coupling,
1109 "The function which specifies if a flux coupling occurs over a given "
1115 keep_constrained_dofs,
1119 face_has_flux_coupling);
1127 #include "dof_tools_sparsity.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, 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
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
cell_iterator begin(const unsigned int level=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() 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
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
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 global_dof_index
unsigned int subdomain_id
unsigned short int fe_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation