Reference documentation for deal.II version GIT 2896a7e638 2023-05-31 13:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_tools_sparsity.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/table.h>
19 #include <deal.II/base/utilities.h>
20 
23 
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 
34 #include <deal.II/grid/tria.h>
36 
38 #include <deal.II/hp/fe_values.h>
40 
43 #include <deal.II/lac/vector.h>
44 
45 #include <algorithm>
46 #include <numeric>
47 
49 
50 
51 
52 namespace DoFTools
53 {
54  template <int dim, int spacedim, typename number>
55  void
57  SparsityPatternBase & sparsity,
58  const AffineConstraints<number> &constraints,
59  const bool keep_constrained_dofs,
61  {
62  const types::global_dof_index n_dofs = dof.n_dofs();
63  (void)n_dofs;
64 
65  Assert(sparsity.n_rows() == n_dofs,
66  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
67  Assert(sparsity.n_cols() == n_dofs,
68  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
69 
70  // If we have a distributed Triangulation only allow locally_owned
71  // subdomain. Not setting a subdomain is also okay, because we skip
72  // ghost cells in the loop below.
73  if (const auto *triangulation = dynamic_cast<
75  &dof.get_triangulation()))
76  {
78  (subdomain_id == triangulation->locally_owned_subdomain()),
79  ExcMessage(
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."));
83  }
84 
85  std::vector<types::global_dof_index> dofs_on_this_cell;
86  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
87 
88  // In case we work with a distributed sparsity pattern of Trilinos
89  // type, we only have to do the work if the current cell is owned by
90  // the calling processor. Otherwise, just continue.
91  for (const auto &cell : dof.active_cell_iterators())
93  (subdomain_id == cell->subdomain_id())) &&
94  cell->is_locally_owned())
95  {
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);
99 
100  // make sparsity pattern for this cell. if no constraints pattern
101  // was given, then the following call acts as if simply no
102  // constraints existed
103  constraints.add_entries_local_to_global(dofs_on_this_cell,
104  sparsity,
105  keep_constrained_dofs);
106  }
107  }
108 
109 
110 
111  template <int dim, int spacedim, typename number>
112  void
114  const Table<2, Coupling> & couplings,
115  SparsityPatternBase & sparsity,
116  const AffineConstraints<number> &constraints,
117  const bool keep_constrained_dofs,
119  {
120  const types::global_dof_index n_dofs = dof.n_dofs();
121  (void)n_dofs;
122 
123  Assert(sparsity.n_rows() == n_dofs,
124  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
125  Assert(sparsity.n_cols() == n_dofs,
126  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
127  Assert(couplings.n_rows() == dof.get_fe(0).n_components(),
128  ExcDimensionMismatch(couplings.n_rows(),
129  dof.get_fe(0).n_components()));
130  Assert(couplings.n_cols() == dof.get_fe(0).n_components(),
131  ExcDimensionMismatch(couplings.n_cols(),
132  dof.get_fe(0).n_components()));
133 
134  // If we have a distributed Triangulation only allow locally_owned
135  // subdomain. Not setting a subdomain is also okay, because we skip
136  // ghost cells in the loop below.
137  if (const auto *triangulation = dynamic_cast<
139  &dof.get_triangulation()))
140  {
142  (subdomain_id == triangulation->locally_owned_subdomain()),
143  ExcMessage(
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."));
147  }
148 
149  const hp::FECollection<dim, spacedim> &fe_collection =
150  dof.get_fe_collection();
151 
152  const std::vector<Table<2, Coupling>> dof_mask //(fe_collection.size())
153  = dof_couplings_from_component_couplings(fe_collection, couplings);
154 
155  // Convert the dof_mask to bool_dof_mask so we can pass it
156  // to constraints.add_entries_local_to_global()
157  std::vector<Table<2, bool>> bool_dof_mask(fe_collection.size());
158  for (unsigned int f = 0; f < fe_collection.size(); ++f)
159  {
160  bool_dof_mask[f].reinit(
161  TableIndices<2>(fe_collection[f].n_dofs_per_cell(),
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;
168  }
169 
170  std::vector<types::global_dof_index> dofs_on_this_cell(
171  fe_collection.max_dofs_per_cell());
172 
173  // In case we work with a distributed sparsity pattern of Trilinos
174  // type, we only have to do the work if the current cell is owned by
175  // the calling processor. Otherwise, just continue.
176  for (const auto &cell : dof.active_cell_iterators())
178  (subdomain_id == cell->subdomain_id())) &&
179  cell->is_locally_owned())
180  {
181  const types::fe_index fe_index = cell->active_fe_index();
182  const unsigned int dofs_per_cell =
183  fe_collection[fe_index].n_dofs_per_cell();
184 
185  dofs_on_this_cell.resize(dofs_per_cell);
186  cell->get_dof_indices(dofs_on_this_cell);
187 
188 
189  // make sparsity pattern for this cell. if no constraints pattern
190  // was given, then the following call acts as if simply no
191  // constraints existed
192  constraints.add_entries_local_to_global(dofs_on_this_cell,
193  sparsity,
194  keep_constrained_dofs,
195  bool_dof_mask[fe_index]);
196  }
197  }
198 
199 
200 
201  template <int dim, int spacedim>
202  void
204  const DoFHandler<dim, spacedim> &dof_col,
205  SparsityPatternBase & sparsity)
206  {
207  const types::global_dof_index n_dofs_row = dof_row.n_dofs();
208  const types::global_dof_index n_dofs_col = dof_col.n_dofs();
209  (void)n_dofs_row;
210  (void)n_dofs_col;
211 
212  Assert(sparsity.n_rows() == n_dofs_row,
213  ExcDimensionMismatch(sparsity.n_rows(), n_dofs_row));
214  Assert(sparsity.n_cols() == n_dofs_col,
215  ExcDimensionMismatch(sparsity.n_cols(), n_dofs_col));
216 
217  // It doesn't make sense to assemble sparsity patterns when the
218  // Triangulations are both parallel (i.e., different cells are assigned to
219  // different processors) and unequal: no processor will be responsible for
220  // assembling coupling terms between dofs on a cell owned by one processor
221  // and dofs on a cell owned by a different processor.
222  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
223  &dof_row.get_triangulation()) != nullptr ||
224  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
225  &dof_col.get_triangulation()) != nullptr)
226  {
227  Assert(&dof_row.get_triangulation() == &dof_col.get_triangulation(),
228  ExcMessage("This function can only be used with with parallel "
229  "Triangulations when the Triangulations are equal."));
230  }
231 
232  // TODO: Looks like wasteful memory management here
233 
234  using cell_iterator = typename DoFHandler<dim, spacedim>::cell_iterator;
235  std::list<std::pair<cell_iterator, cell_iterator>> cell_list =
236  GridTools::get_finest_common_cells(dof_row, dof_col);
237 
238 #ifdef DEAL_II_WITH_MPI
239  // get_finest_common_cells returns all cells (locally owned and otherwise)
240  // for shared::Tria, but we don't want to assemble on cells that are not
241  // locally owned so remove them
242  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
243  &dof_row.get_triangulation()) != nullptr ||
245  &dof_col.get_triangulation()) != nullptr)
246  {
247  const types::subdomain_id this_subdomain_id =
249  Assert(this_subdomain_id ==
251  ExcInternalError());
252  cell_list.erase(
253  std::remove_if(
254  cell_list.begin(),
255  cell_list.end(),
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;
259  }),
260  cell_list.end());
261  }
262 #endif
263 
264  for (const auto &cell_pair : cell_list)
265  {
266  const cell_iterator cell_row = cell_pair.first;
267  const cell_iterator cell_col = cell_pair.second;
268 
269  if (cell_row->is_active() && cell_col->is_active())
270  {
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(
276  dofs_per_cell_row);
277  std::vector<types::global_dof_index> local_dof_indices_col(
278  dofs_per_cell_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)
282  sparsity.add_row_entries(dof,
283  make_array_view(local_dof_indices_col));
284  }
285  else if (cell_row->has_children())
286  {
287  const std::vector<
289  child_cells =
290  GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
291  cell_row);
292  for (unsigned int i = 0; i < child_cells.size(); ++i)
293  {
295  cell_row_child = child_cells[i];
296  const unsigned int dofs_per_cell_row =
297  cell_row_child->get_fe().n_dofs_per_cell();
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(
301  dofs_per_cell_row);
302  std::vector<types::global_dof_index> local_dof_indices_col(
303  dofs_per_cell_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)
307  sparsity.add_row_entries(
308  dof, make_array_view(local_dof_indices_col));
309  }
310  }
311  else
312  {
313  std::vector<
315  child_cells =
316  GridTools::get_active_child_cells<DoFHandler<dim, spacedim>>(
317  cell_col);
318  for (unsigned int i = 0; i < child_cells.size(); ++i)
319  {
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 =
325  cell_col_child->get_fe().n_dofs_per_cell();
326  std::vector<types::global_dof_index> local_dof_indices_row(
327  dofs_per_cell_row);
328  std::vector<types::global_dof_index> local_dof_indices_col(
329  dofs_per_cell_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)
333  sparsity.add_row_entries(
334  dof, make_array_view(local_dof_indices_col));
335  }
336  }
337  }
338  }
339 
340 
341 
342  template <int dim, int spacedim>
343  void
345  const DoFHandler<dim, spacedim> & dof,
346  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
347  SparsityPatternBase & sparsity)
348  {
349  if (dim == 1)
350  {
351  // there are only 2 boundary indicators in 1d, so it is no
352  // performance problem to call the other function
353  std::map<types::boundary_id, const Function<spacedim, double> *>
354  boundary_ids;
355  boundary_ids[0] = nullptr;
356  boundary_ids[1] = nullptr;
357  make_boundary_sparsity_pattern<dim, spacedim>(dof,
358  boundary_ids,
359  dof_to_boundary_mapping,
360  sparsity);
361  return;
362  }
363 
364  const types::global_dof_index n_dofs = dof.n_dofs();
365  (void)n_dofs;
366 
367  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
368  AssertDimension(sparsity.n_rows(), dof.n_boundary_dofs());
369  AssertDimension(sparsity.n_cols(), dof.n_boundary_dofs());
370 #ifdef DEBUG
371  if (sparsity.n_rows() != 0)
372  {
373  types::global_dof_index max_element = 0;
374  for (const types::global_dof_index index : dof_to_boundary_mapping)
375  if ((index != numbers::invalid_dof_index) && (index > max_element))
376  max_element = index;
377  AssertDimension(max_element, sparsity.n_rows() - 1);
378  }
379 #endif
380 
381  std::vector<types::global_dof_index> dofs_on_this_face;
382  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
383  std::vector<types::global_dof_index> cols;
384 
385  // loop over all faces to check whether they are at a boundary. note
386  // that we need not take special care of single lines (using
387  // @p{cell->has_boundary_lines}), since we do not support boundaries of
388  // dimension dim-2, and so every boundary line is also part of a
389  // boundary face.
390  for (const auto &cell : dof.active_cell_iterators())
391  for (const unsigned int f : cell->face_indices())
392  if (cell->at_boundary(f))
393  {
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());
399 
400  // make sparsity pattern for this cell
401  cols.clear();
402  for (const auto &dof : dofs_on_this_face)
403  cols.push_back(dof_to_boundary_mapping[dof]);
404  // We are not guaranteed that the mapping to a second index space
405  // is increasing so sort here to use the faster add_row_entries()
406  // path
407  std::sort(cols.begin(), cols.end());
408  for (const auto &dof : dofs_on_this_face)
409  sparsity.add_row_entries(dof_to_boundary_mapping[dof],
410  make_array_view(cols),
411  true);
412  }
413  }
414 
415 
416 
417  template <int dim, int spacedim, typename number>
418  void
420  const DoFHandler<dim, spacedim> &dof,
421  const std::map<types::boundary_id, const Function<spacedim, number> *>
422  & boundary_ids,
423  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
424  SparsityPatternBase & sparsity)
425  {
426  if (dim == 1)
427  {
428  // first check left, then right boundary point
429  for (unsigned int direction = 0; direction < 2; ++direction)
430  {
431  // if this boundary is not requested, then go on with next one
432  if (boundary_ids.find(direction) == boundary_ids.end())
433  continue;
434 
435  // find active cell at that boundary: first go to left/right,
436  // then to children
438  dof.begin(0);
439  while (!cell->at_boundary(direction))
440  cell = cell->neighbor(direction);
441  while (!cell->is_active())
442  cell = cell->child(direction);
443 
444  const unsigned int dofs_per_vertex =
445  cell->get_fe().n_dofs_per_vertex();
446  std::vector<types::global_dof_index> boundary_dof_boundary_indices(
447  dofs_per_vertex);
448 
449  // next get boundary mapped dof indices of boundary dofs
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)];
453 
454  std::sort(boundary_dof_boundary_indices.begin(),
455  boundary_dof_boundary_indices.end());
456  for (const auto &dof : boundary_dof_boundary_indices)
457  sparsity.add_row_entries(
458  dof, make_array_view(boundary_dof_boundary_indices), true);
459  }
460  return;
461  }
462 
463  const types::global_dof_index n_dofs = dof.n_dofs();
464  (void)n_dofs;
465 
466  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
467  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
468  boundary_ids.end(),
470  Assert(sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
471  ExcDimensionMismatch(sparsity.n_rows(),
472  dof.n_boundary_dofs(boundary_ids)));
473  Assert(sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
474  ExcDimensionMismatch(sparsity.n_cols(),
475  dof.n_boundary_dofs(boundary_ids)));
476 #ifdef DEBUG
477  if (sparsity.n_rows() != 0)
478  {
479  types::global_dof_index max_element = 0;
480  for (const types::global_dof_index index : dof_to_boundary_mapping)
481  if ((index != numbers::invalid_dof_index) && (index > max_element))
482  max_element = index;
483  AssertDimension(max_element, sparsity.n_rows() - 1);
484  }
485 #endif
486 
487  std::vector<types::global_dof_index> dofs_on_this_face;
488  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
489  std::vector<types::global_dof_index> cols;
490 
491  for (const auto &cell : dof.active_cell_iterators())
492  for (const unsigned int f : cell->face_indices())
493  if (boundary_ids.find(cell->face(f)->boundary_id()) !=
494  boundary_ids.end())
495  {
496  const unsigned int dofs_per_face =
497  cell->get_fe().n_dofs_per_face(f);
498  dofs_on_this_face.resize(dofs_per_face);
499  cell->face(f)->get_dof_indices(dofs_on_this_face,
500  cell->active_fe_index());
501 
502  // make sparsity pattern for this cell
503  cols.clear();
504  for (const auto &dof : dofs_on_this_face)
505  cols.push_back(dof_to_boundary_mapping[dof]);
506 
507  // Like the other one: sort once.
508  std::sort(cols.begin(), cols.end());
509  for (const auto &dof : dofs_on_this_face)
510  sparsity.add_row_entries(dof_to_boundary_mapping[dof],
511  make_array_view(cols),
512  true);
513  }
514  }
515 
516 
517 
518  template <int dim, int spacedim, typename number>
519  void
521  SparsityPatternBase & sparsity,
522  const AffineConstraints<number> &constraints,
523  const bool keep_constrained_dofs,
525 
526  // TODO: QA: reduce the indentation level of this method..., Maier 2012
527 
528  {
529  const types::global_dof_index n_dofs = dof.n_dofs();
530  (void)n_dofs;
531 
532  AssertDimension(sparsity.n_rows(), n_dofs);
533  AssertDimension(sparsity.n_cols(), n_dofs);
534 
535  // If we have a distributed Triangulation only allow locally_owned
536  // subdomain. Not setting a subdomain is also okay, because we skip
537  // ghost cells in the loop below.
538  if (const auto *triangulation = dynamic_cast<
540  &dof.get_triangulation()))
541  {
543  (subdomain_id == triangulation->locally_owned_subdomain()),
544  ExcMessage(
545  "For distributed Triangulation objects and associated "
546  "DoFHandler objects, asking for any subdomain other than the "
547  "locally owned one does not make sense."));
548  }
549 
550  std::vector<types::global_dof_index> dofs_on_this_cell;
551  std::vector<types::global_dof_index> dofs_on_other_cell;
552  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
553  dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
554 
555  // TODO: in an old implementation, we used user flags before to tag
556  // faces that were already touched. this way, we could reduce the work
557  // a little bit. now, we instead add only data from one side. this
558  // should be OK, but we need to actually verify it.
559 
560  // In case we work with a distributed sparsity pattern of Trilinos
561  // type, we only have to do the work if the current cell is owned by
562  // the calling processor. Otherwise, just continue.
563  for (const auto &cell : dof.active_cell_iterators())
565  (subdomain_id == cell->subdomain_id())) &&
566  cell->is_locally_owned())
567  {
568  const unsigned int n_dofs_on_this_cell =
569  cell->get_fe().n_dofs_per_cell();
570  dofs_on_this_cell.resize(n_dofs_on_this_cell);
571  cell->get_dof_indices(dofs_on_this_cell);
572 
573  // make sparsity pattern for this cell. if no constraints pattern
574  // was given, then the following call acts as if simply no
575  // constraints existed
576  constraints.add_entries_local_to_global(dofs_on_this_cell,
577  sparsity,
578  keep_constrained_dofs);
579 
580  for (const unsigned int face : cell->face_indices())
581  {
582  typename DoFHandler<dim, spacedim>::face_iterator cell_face =
583  cell->face(face);
584  const bool periodic_neighbor = cell->has_periodic_neighbor(face);
585  if (!cell->at_boundary(face) || periodic_neighbor)
586  {
588  neighbor = cell->neighbor_or_periodic_neighbor(face);
589 
590  // in 1d, we do not need to worry whether the neighbor
591  // might have children and then loop over those children.
592  // rather, we may as well go straight to the cell behind
593  // this particular cell's most terminal child
594  if (dim == 1)
595  while (neighbor->has_children())
596  neighbor = neighbor->child(face == 0 ? 1 : 0);
597 
598  if (neighbor->has_children())
599  {
600  for (unsigned int sub_nr = 0;
601  sub_nr != cell_face->n_active_descendants();
602  ++sub_nr)
603  {
604  const typename DoFHandler<dim, spacedim>::
605  level_cell_iterator sub_neighbor =
606  periodic_neighbor ?
607  cell->periodic_neighbor_child_on_subface(
608  face, sub_nr) :
609  cell->neighbor_child_on_subface(face, sub_nr);
610 
611  const unsigned int n_dofs_on_neighbor =
612  sub_neighbor->get_fe().n_dofs_per_cell();
613  dofs_on_other_cell.resize(n_dofs_on_neighbor);
614  sub_neighbor->get_dof_indices(dofs_on_other_cell);
615 
616  constraints.add_entries_local_to_global(
617  dofs_on_this_cell,
618  dofs_on_other_cell,
619  sparsity,
620  keep_constrained_dofs);
621  constraints.add_entries_local_to_global(
622  dofs_on_other_cell,
623  dofs_on_this_cell,
624  sparsity,
625  keep_constrained_dofs);
626  // only need to add this when the neighbor is not
627  // owned by the current processor, otherwise we add
628  // the entries for the neighbor there
629  if (sub_neighbor->subdomain_id() !=
630  cell->subdomain_id())
631  constraints.add_entries_local_to_global(
632  dofs_on_other_cell,
633  sparsity,
634  keep_constrained_dofs);
635  }
636  }
637  else
638  {
639  // Refinement edges are taken care of by coarser
640  // cells
641  if ((!periodic_neighbor &&
642  cell->neighbor_is_coarser(face)) ||
643  (periodic_neighbor &&
644  cell->periodic_neighbor_is_coarser(face)))
645  if (neighbor->subdomain_id() == cell->subdomain_id())
646  continue;
647 
648  const unsigned int n_dofs_on_neighbor =
649  neighbor->get_fe().n_dofs_per_cell();
650  dofs_on_other_cell.resize(n_dofs_on_neighbor);
651 
652  neighbor->get_dof_indices(dofs_on_other_cell);
653 
654  constraints.add_entries_local_to_global(
655  dofs_on_this_cell,
656  dofs_on_other_cell,
657  sparsity,
658  keep_constrained_dofs);
659 
660  // only need to add these in case the neighbor cell
661  // is not locally owned - otherwise, we touch each
662  // face twice and hence put the indices the other way
663  // around
664  if (!cell->neighbor_or_periodic_neighbor(face)
665  ->is_active() ||
666  (neighbor->subdomain_id() != cell->subdomain_id()))
667  {
668  constraints.add_entries_local_to_global(
669  dofs_on_other_cell,
670  dofs_on_this_cell,
671  sparsity,
672  keep_constrained_dofs);
673  if (neighbor->subdomain_id() != cell->subdomain_id())
674  constraints.add_entries_local_to_global(
675  dofs_on_other_cell,
676  sparsity,
677  keep_constrained_dofs);
678  }
679  }
680  }
681  }
682  }
683  }
684 
685 
686 
687  template <int dim, int spacedim>
688  void
690  SparsityPatternBase & sparsity)
691  {
693  make_flux_sparsity_pattern(dof, sparsity, dummy);
694  }
695 
696  template <int dim, int spacedim>
700  const Table<2, Coupling> & component_couplings)
701  {
702  Assert(component_couplings.n_rows() == fe.n_components(),
703  ExcDimensionMismatch(component_couplings.n_rows(),
704  fe.n_components()));
705  Assert(component_couplings.n_cols() == fe.n_components(),
706  ExcDimensionMismatch(component_couplings.n_cols(),
707  fe.n_components()));
708 
709  const unsigned int n_dofs = fe.n_dofs_per_cell();
710 
711  Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
712 
713  for (unsigned int i = 0; i < n_dofs; ++i)
714  {
715  const unsigned int ii =
716  (fe.is_primitive(i) ?
717  fe.system_to_component_index(i).first :
719  Assert(ii < fe.n_components(), ExcInternalError());
720 
721  for (unsigned int j = 0; j < n_dofs; ++j)
722  {
723  const unsigned int jj =
724  (fe.is_primitive(j) ?
725  fe.system_to_component_index(j).first :
727  Assert(jj < fe.n_components(), ExcInternalError());
728 
729  dof_couplings(i, j) = component_couplings(ii, jj);
730  }
731  }
732  return dof_couplings;
733  }
734 
735 
736 
737  template <int dim, int spacedim>
738  std::vector<Table<2, Coupling>>
741  const Table<2, Coupling> & component_couplings)
742  {
743  std::vector<Table<2, Coupling>> return_value(fe.size());
744  for (unsigned int i = 0; i < fe.size(); ++i)
745  return_value[i] =
746  dof_couplings_from_component_couplings(fe[i], component_couplings);
747 
748  return return_value;
749  }
750 
751 
752 
753  namespace internal
754  {
755  namespace
756  {
757  // implementation of the same function in namespace DoFTools for
758  // non-hp-DoFHandlers
759  template <int dim, int spacedim, typename number>
760  void
762  const DoFHandler<dim, spacedim> &dof,
763  SparsityPatternBase & sparsity,
764  const AffineConstraints<number> &constraints,
765  const bool keep_constrained_dofs,
766  const Table<2, Coupling> & int_mask,
767  const Table<2, Coupling> & flux_mask,
769  const std::function<
771  const unsigned int)> &face_has_flux_coupling)
772  {
773  std::vector<types::global_dof_index> rows;
774  std::vector<std::pair<SparsityPatternBase::size_type,
776  cell_entries;
777 
778  if (dof.has_hp_capabilities() == false)
779  {
780  const FiniteElement<dim, spacedim> &fe = dof.get_fe();
781 
782  std::vector<types::global_dof_index> dofs_on_this_cell(
783  fe.n_dofs_per_cell());
784  std::vector<types::global_dof_index> dofs_on_other_cell(
785  fe.n_dofs_per_cell());
786 
787  const Table<2, Coupling>
788  int_dof_mask =
790  flux_dof_mask =
792 
793  Table<2, bool> support_on_face(fe.n_dofs_per_cell(),
795  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
796  for (const unsigned int f : GeometryInfo<dim>::face_indices())
797  support_on_face(i, f) = fe.has_support_on_face(i, f);
798 
799  // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
800  // to constraints.add_entries_local_to_global()
801  Table<2, bool> bool_int_dof_mask(fe.n_dofs_per_cell(),
802  fe.n_dofs_per_cell());
803  bool_int_dof_mask.fill(false);
804  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
805  for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
806  if (int_dof_mask(i, j) != none)
807  bool_int_dof_mask(i, j) = true;
808 
809  for (const auto &cell : dof.active_cell_iterators())
811  (subdomain_id == cell->subdomain_id())) &&
812  cell->is_locally_owned())
813  {
814  cell->get_dof_indices(dofs_on_this_cell);
815  // make sparsity pattern for this cell
816  constraints.add_entries_local_to_global(dofs_on_this_cell,
817  sparsity,
818  keep_constrained_dofs,
819  bool_int_dof_mask);
820  // Loop over all interior neighbors
821  for (const unsigned int face_n : cell->face_indices())
822  {
824  cell_face = cell->face(face_n);
825 
826  const bool periodic_neighbor =
827  cell->has_periodic_neighbor(face_n);
828 
829  if (cell->at_boundary(face_n) && (!periodic_neighbor))
830  {
831  for (unsigned int i = 0; i < fe.n_dofs_per_cell();
832  ++i)
833  {
834  const bool i_non_zero_i =
835  support_on_face(i, face_n);
836  for (unsigned int j = 0; j < fe.n_dofs_per_cell();
837  ++j)
838  {
839  const bool j_non_zero_i =
840  support_on_face(j, face_n);
841 
842  if (flux_dof_mask(i, j) == always ||
843  (flux_dof_mask(i, j) == nonzero &&
844  i_non_zero_i && j_non_zero_i))
845  cell_entries.emplace_back(
846  dofs_on_this_cell[i],
847  dofs_on_this_cell[j]);
848  }
849  }
850  sparsity.add_entries(make_array_view(cell_entries));
851  cell_entries.clear();
852  }
853  else
854  {
855  typename DoFHandler<dim,
856  spacedim>::level_cell_iterator
857  neighbor =
858  cell->neighbor_or_periodic_neighbor(face_n);
859  // If the cells are on the same level (and both are
860  // active, locally-owned cells) then only add to the
861  // sparsity pattern if the current cell is 'greater'
862  // in the total ordering.
863  if (neighbor->level() == cell->level() &&
864  neighbor->index() > cell->index() &&
865  neighbor->is_active() &&
866  neighbor->is_locally_owned())
867  continue;
868  // If we are more refined then the neighbor, then we
869  // will automatically find the active neighbor cell
870  // when we call 'neighbor (face_n)' above. The
871  // opposite is not true; if the neighbor is more
872  // refined then the call 'neighbor (face_n)' will
873  // *not* return an active cell. Hence, only add things
874  // to the sparsity pattern if (when the levels are
875  // different) the neighbor is coarser than the current
876  // cell.
877  //
878  // Like above, do not use this optimization if the
879  // neighbor is not locally owned.
880  if (neighbor->level() != cell->level() &&
881  ((!periodic_neighbor &&
882  !cell->neighbor_is_coarser(face_n)) ||
883  (periodic_neighbor &&
884  !cell->periodic_neighbor_is_coarser(face_n))) &&
885  neighbor->is_locally_owned())
886  continue; // (the neighbor is finer)
887 
888  if (!face_has_flux_coupling(cell, face_n))
889  continue;
890 
891 
892  const unsigned int neighbor_face_n =
893  periodic_neighbor ?
894  cell->periodic_neighbor_face_no(face_n) :
895  cell->neighbor_face_no(face_n);
896 
897 
898  // In 1d, go straight to the cell behind this
899  // particular cell's most terminal cell. This makes us
900  // skip the if (neighbor->has_children()) section
901  // below. We need to do this since we otherwise
902  // iterate over the children of the face, which are
903  // always 0 in 1d.
904  if (dim == 1)
905  while (neighbor->has_children())
906  neighbor = neighbor->child(face_n == 0 ? 1 : 0);
907 
908  if (neighbor->has_children())
909  {
910  for (unsigned int sub_nr = 0;
911  sub_nr != cell_face->n_children();
912  ++sub_nr)
913  {
914  const typename DoFHandler<dim, spacedim>::
915  level_cell_iterator sub_neighbor =
916  periodic_neighbor ?
917  cell
918  ->periodic_neighbor_child_on_subface(
919  face_n, sub_nr) :
920  cell->neighbor_child_on_subface(face_n,
921  sub_nr);
922 
923  sub_neighbor->get_dof_indices(
924  dofs_on_other_cell);
925  for (unsigned int i = 0;
926  i < fe.n_dofs_per_cell();
927  ++i)
928  {
929  const bool i_non_zero_i =
930  support_on_face(i, face_n);
931  const bool i_non_zero_e =
932  support_on_face(i, neighbor_face_n);
933  for (unsigned int j = 0;
934  j < fe.n_dofs_per_cell();
935  ++j)
936  {
937  const bool j_non_zero_i =
938  support_on_face(j, face_n);
939  const bool j_non_zero_e =
940  support_on_face(j, neighbor_face_n);
941 
942  if (flux_dof_mask(i, j) == always)
943  {
944  cell_entries.emplace_back(
945  dofs_on_this_cell[i],
946  dofs_on_other_cell[j]);
947  cell_entries.emplace_back(
948  dofs_on_other_cell[i],
949  dofs_on_this_cell[j]);
950  cell_entries.emplace_back(
951  dofs_on_this_cell[i],
952  dofs_on_this_cell[j]);
953  cell_entries.emplace_back(
954  dofs_on_other_cell[i],
955  dofs_on_other_cell[j]);
956  }
957  else if (flux_dof_mask(i, j) ==
958  nonzero)
959  {
960  if (i_non_zero_i && j_non_zero_e)
961  cell_entries.emplace_back(
962  dofs_on_this_cell[i],
963  dofs_on_other_cell[j]);
964  if (i_non_zero_e && j_non_zero_i)
965  cell_entries.emplace_back(
966  dofs_on_other_cell[i],
967  dofs_on_this_cell[j]);
968  if (i_non_zero_i && j_non_zero_i)
969  cell_entries.emplace_back(
970  dofs_on_this_cell[i],
971  dofs_on_this_cell[j]);
972  if (i_non_zero_e && j_non_zero_e)
973  cell_entries.emplace_back(
974  dofs_on_other_cell[i],
975  dofs_on_other_cell[j]);
976  }
977 
978  if (flux_dof_mask(j, i) == always)
979  {
980  cell_entries.emplace_back(
981  dofs_on_this_cell[j],
982  dofs_on_other_cell[i]);
983  cell_entries.emplace_back(
984  dofs_on_other_cell[j],
985  dofs_on_this_cell[i]);
986  cell_entries.emplace_back(
987  dofs_on_this_cell[j],
988  dofs_on_this_cell[i]);
989  cell_entries.emplace_back(
990  dofs_on_other_cell[j],
991  dofs_on_other_cell[i]);
992  }
993  else if (flux_dof_mask(j, i) ==
994  nonzero)
995  {
996  if (j_non_zero_i && i_non_zero_e)
997  cell_entries.emplace_back(
998  dofs_on_this_cell[j],
999  dofs_on_other_cell[i]);
1000  if (j_non_zero_e && i_non_zero_i)
1001  cell_entries.emplace_back(
1002  dofs_on_other_cell[j],
1003  dofs_on_this_cell[i]);
1004  if (j_non_zero_i && i_non_zero_i)
1005  cell_entries.emplace_back(
1006  dofs_on_this_cell[j],
1007  dofs_on_this_cell[i]);
1008  if (j_non_zero_e && i_non_zero_e)
1009  cell_entries.emplace_back(
1010  dofs_on_other_cell[j],
1011  dofs_on_other_cell[i]);
1012  }
1013  }
1014  }
1015  }
1016  sparsity.add_entries(
1017  make_array_view(cell_entries));
1018  cell_entries.clear();
1019  }
1020  else
1021  {
1022  neighbor->get_dof_indices(dofs_on_other_cell);
1023  for (unsigned int i = 0; i < fe.n_dofs_per_cell();
1024  ++i)
1025  {
1026  const bool i_non_zero_i =
1027  support_on_face(i, face_n);
1028  const bool i_non_zero_e =
1029  support_on_face(i, neighbor_face_n);
1030  for (unsigned int j = 0;
1031  j < fe.n_dofs_per_cell();
1032  ++j)
1033  {
1034  const bool j_non_zero_i =
1035  support_on_face(j, face_n);
1036  const bool j_non_zero_e =
1037  support_on_face(j, neighbor_face_n);
1038  if (flux_dof_mask(i, j) == always)
1039  {
1040  cell_entries.emplace_back(
1041  dofs_on_this_cell[i],
1042  dofs_on_other_cell[j]);
1043  cell_entries.emplace_back(
1044  dofs_on_other_cell[i],
1045  dofs_on_this_cell[j]);
1046  cell_entries.emplace_back(
1047  dofs_on_this_cell[i],
1048  dofs_on_this_cell[j]);
1049  cell_entries.emplace_back(
1050  dofs_on_other_cell[i],
1051  dofs_on_other_cell[j]);
1052  }
1053  if (flux_dof_mask(i, j) == nonzero)
1054  {
1055  if (i_non_zero_i && j_non_zero_e)
1056  cell_entries.emplace_back(
1057  dofs_on_this_cell[i],
1058  dofs_on_other_cell[j]);
1059  if (i_non_zero_e && j_non_zero_i)
1060  cell_entries.emplace_back(
1061  dofs_on_other_cell[i],
1062  dofs_on_this_cell[j]);
1063  if (i_non_zero_i && j_non_zero_i)
1064  cell_entries.emplace_back(
1065  dofs_on_this_cell[i],
1066  dofs_on_this_cell[j]);
1067  if (i_non_zero_e && j_non_zero_e)
1068  cell_entries.emplace_back(
1069  dofs_on_other_cell[i],
1070  dofs_on_other_cell[j]);
1071  }
1072 
1073  if (flux_dof_mask(j, i) == always)
1074  {
1075  cell_entries.emplace_back(
1076  dofs_on_this_cell[j],
1077  dofs_on_other_cell[i]);
1078  cell_entries.emplace_back(
1079  dofs_on_other_cell[j],
1080  dofs_on_this_cell[i]);
1081  cell_entries.emplace_back(
1082  dofs_on_this_cell[j],
1083  dofs_on_this_cell[i]);
1084  cell_entries.emplace_back(
1085  dofs_on_other_cell[j],
1086  dofs_on_other_cell[i]);
1087  }
1088  if (flux_dof_mask(j, i) == nonzero)
1089  {
1090  if (j_non_zero_i && i_non_zero_e)
1091  cell_entries.emplace_back(
1092  dofs_on_this_cell[j],
1093  dofs_on_other_cell[i]);
1094  if (j_non_zero_e && i_non_zero_i)
1095  cell_entries.emplace_back(
1096  dofs_on_other_cell[j],
1097  dofs_on_this_cell[i]);
1098  if (j_non_zero_i && i_non_zero_i)
1099  cell_entries.emplace_back(
1100  dofs_on_this_cell[j],
1101  dofs_on_this_cell[i]);
1102  if (j_non_zero_e && i_non_zero_e)
1103  cell_entries.emplace_back(
1104  dofs_on_other_cell[j],
1105  dofs_on_other_cell[i]);
1106  }
1107  }
1108  }
1109  sparsity.add_entries(
1110  make_array_view(cell_entries));
1111  cell_entries.clear();
1112  }
1113  }
1114  }
1115  }
1116  }
1117  else
1118  {
1119  // while the implementation above is quite optimized and caches a
1120  // lot of data (see e.g. the int/flux_dof_mask tables), this is no
1121  // longer practical for the hp-version since we would have to have
1122  // it for all combinations of elements in the hp::FECollection.
1123  // consequently, the implementation here is simpler and probably
1124  // less efficient but at least readable...
1125 
1126  const ::hp::FECollection<dim, spacedim> &fe =
1127  dof.get_fe_collection();
1128 
1129  std::vector<types::global_dof_index> dofs_on_this_cell(
1131  std::vector<types::global_dof_index> dofs_on_other_cell(
1133 
1134  const unsigned int n_components = fe.n_components();
1135  AssertDimension(int_mask.size(0), n_components);
1136  AssertDimension(int_mask.size(1), n_components);
1137  AssertDimension(flux_mask.size(0), n_components);
1138  AssertDimension(flux_mask.size(1), n_components);
1139 
1140  // note that we also need to set the respective entries if flux_mask
1141  // says so. this is necessary since we need to consider all degrees
1142  // of freedom on a cell for interior faces.
1143  Table<2, Coupling> int_and_flux_mask(n_components, n_components);
1144  for (unsigned int c1 = 0; c1 < n_components; ++c1)
1145  for (unsigned int c2 = 0; c2 < n_components; ++c2)
1146  if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
1147  int_and_flux_mask(c1, c2) = always;
1148 
1149  std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
1150  dof_couplings_from_component_couplings(fe, int_and_flux_mask);
1151 
1152  // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
1153  // can pass it to constraints.add_entries_local_to_global()
1154  std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
1155  for (unsigned int f = 0; f < fe.size(); ++f)
1156  {
1157  bool_int_and_flux_dof_mask[f].reinit(
1158  TableIndices<2>(fe[f].n_dofs_per_cell(),
1159  fe[f].n_dofs_per_cell()));
1160  bool_int_and_flux_dof_mask[f].fill(false);
1161  for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
1162  for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
1163  if (int_and_flux_dof_mask[f](i, j) != none)
1164  bool_int_and_flux_dof_mask[f](i, j) = true;
1165  }
1166 
1167 
1168  for (const auto &cell : dof.active_cell_iterators())
1170  (subdomain_id == cell->subdomain_id())) &&
1171  cell->is_locally_owned())
1172  {
1173  dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
1174  cell->get_dof_indices(dofs_on_this_cell);
1175 
1176  // make sparsity pattern for this cell also taking into
1177  // account the couplings due to face contributions on the same
1178  // cell
1179  constraints.add_entries_local_to_global(
1180  dofs_on_this_cell,
1181  sparsity,
1182  keep_constrained_dofs,
1183  bool_int_and_flux_dof_mask[cell->active_fe_index()]);
1184 
1185  // Loop over interior faces
1186  for (const unsigned int face : cell->face_indices())
1187  {
1189  cell_face = cell->face(face);
1190 
1191  const bool periodic_neighbor =
1192  cell->has_periodic_neighbor(face);
1193 
1194  if ((!cell->at_boundary(face)) || periodic_neighbor)
1195  {
1196  typename DoFHandler<dim,
1197  spacedim>::level_cell_iterator
1198  neighbor =
1199  cell->neighbor_or_periodic_neighbor(face);
1200 
1201  // Like the non-hp-case: If the cells are on the same
1202  // level (and both are active, locally-owned cells)
1203  // then only add to the sparsity pattern if the
1204  // current cell is 'greater' in the total ordering.
1205  if (neighbor->level() == cell->level() &&
1206  neighbor->index() > cell->index() &&
1207  neighbor->is_active() &&
1208  neighbor->is_locally_owned())
1209  continue;
1210  // Again, like the non-hp-case: If we are more refined
1211  // then the neighbor, then we will automatically find
1212  // the active neighbor cell when we call 'neighbor
1213  // (face)' above. The opposite is not true; if the
1214  // neighbor is more refined then the call 'neighbor
1215  // (face)' will *not* return an active cell. Hence,
1216  // only add things to the sparsity pattern if (when
1217  // the levels are different) the neighbor is coarser
1218  // than the current cell.
1219  //
1220  // Like above, do not use this optimization if the
1221  // neighbor is not locally owned.
1222  if (neighbor->level() != cell->level() &&
1223  ((!periodic_neighbor &&
1224  !cell->neighbor_is_coarser(face)) ||
1225  (periodic_neighbor &&
1226  !cell->periodic_neighbor_is_coarser(face))) &&
1227  neighbor->is_locally_owned())
1228  continue; // (the neighbor is finer)
1229 
1230 
1231  if (!face_has_flux_coupling(cell, face))
1232  continue;
1233 
1234  // In 1d, go straight to the cell behind this
1235  // particular cell's most terminal cell. This makes us
1236  // skip the if (neighbor->has_children()) section
1237  // below. We need to do this since we otherwise
1238  // iterate over the children of the face, which are
1239  // always 0 in 1d.
1240  if (dim == 1)
1241  while (neighbor->has_children())
1242  neighbor = neighbor->child(face == 0 ? 1 : 0);
1243 
1244  if (neighbor->has_children())
1245  {
1246  for (unsigned int sub_nr = 0;
1247  sub_nr != cell_face->n_children();
1248  ++sub_nr)
1249  {
1250  const typename DoFHandler<dim, spacedim>::
1251  level_cell_iterator sub_neighbor =
1252  periodic_neighbor ?
1253  cell
1254  ->periodic_neighbor_child_on_subface(
1255  face, sub_nr) :
1256  cell->neighbor_child_on_subface(face,
1257  sub_nr);
1258 
1259  dofs_on_other_cell.resize(
1260  sub_neighbor->get_fe().n_dofs_per_cell());
1261  sub_neighbor->get_dof_indices(
1262  dofs_on_other_cell);
1263  for (unsigned int i = 0;
1264  i < cell->get_fe().n_dofs_per_cell();
1265  ++i)
1266  {
1267  const unsigned int ii =
1268  (cell->get_fe().is_primitive(i) ?
1269  cell->get_fe()
1270  .system_to_component_index(i)
1271  .first :
1272  cell->get_fe()
1273  .get_nonzero_components(i)
1274  .first_selected_component());
1275 
1276  Assert(ii < cell->get_fe().n_components(),
1277  ExcInternalError());
1278 
1279  for (unsigned int j = 0;
1280  j < sub_neighbor->get_fe()
1281  .n_dofs_per_cell();
1282  ++j)
1283  {
1284  const unsigned int jj =
1285  (sub_neighbor->get_fe()
1286  .is_primitive(j) ?
1287  sub_neighbor->get_fe()
1289  .first :
1290  sub_neighbor->get_fe()
1293 
1294  Assert(jj < sub_neighbor->get_fe()
1295  .n_components(),
1296  ExcInternalError());
1297 
1298  if ((flux_mask(ii, jj) == always) ||
1299  (flux_mask(ii, jj) == nonzero))
1300  {
1301  cell_entries.emplace_back(
1302  dofs_on_this_cell[i],
1303  dofs_on_other_cell[j]);
1304  }
1305 
1306  if ((flux_mask(jj, ii) == always) ||
1307  (flux_mask(jj, ii) == nonzero))
1308  {
1309  cell_entries.emplace_back(
1310  dofs_on_other_cell[j],
1311  dofs_on_this_cell[i]);
1312  }
1313  }
1314  }
1315  }
1316  }
1317  else
1318  {
1319  dofs_on_other_cell.resize(
1320  neighbor->get_fe().n_dofs_per_cell());
1321  neighbor->get_dof_indices(dofs_on_other_cell);
1322  for (unsigned int i = 0;
1323  i < cell->get_fe().n_dofs_per_cell();
1324  ++i)
1325  {
1326  const unsigned int ii =
1327  (cell->get_fe().is_primitive(i) ?
1328  cell->get_fe()
1329  .system_to_component_index(i)
1330  .first :
1331  cell->get_fe()
1332  .get_nonzero_components(i)
1333  .first_selected_component());
1334 
1335  Assert(ii < cell->get_fe().n_components(),
1336  ExcInternalError());
1337 
1338  for (unsigned int j = 0;
1339  j < neighbor->get_fe().n_dofs_per_cell();
1340  ++j)
1341  {
1342  const unsigned int jj =
1343  (neighbor->get_fe().is_primitive(j) ?
1344  neighbor->get_fe()
1345  .system_to_component_index(j)
1346  .first :
1347  neighbor->get_fe()
1348  .get_nonzero_components(j)
1349  .first_selected_component());
1350 
1351  Assert(
1352  jj < neighbor->get_fe().n_components(),
1353  ExcInternalError());
1354 
1355  if ((flux_mask(ii, jj) == always) ||
1356  (flux_mask(ii, jj) == nonzero))
1357  {
1358  cell_entries.emplace_back(
1359  dofs_on_this_cell[i],
1360  dofs_on_other_cell[j]);
1361  }
1362 
1363  if ((flux_mask(jj, ii) == always) ||
1364  (flux_mask(jj, ii) == nonzero))
1365  {
1366  cell_entries.emplace_back(
1367  dofs_on_other_cell[j],
1368  dofs_on_this_cell[i]);
1369  }
1370  }
1371  }
1372  }
1373  }
1374  }
1375  sparsity.add_entries(make_array_view(cell_entries));
1376  cell_entries.clear();
1377  }
1378  }
1379  }
1380  } // namespace
1381 
1382  } // namespace internal
1383 
1384 
1385 
1386  template <int dim, int spacedim>
1387  void
1389  SparsityPatternBase & sparsity,
1390  const Table<2, Coupling> & int_mask,
1391  const Table<2, Coupling> & flux_mask,
1393  {
1395 
1396  const bool keep_constrained_dofs = true;
1397 
1399  sparsity,
1400  dummy,
1401  keep_constrained_dofs,
1402  int_mask,
1403  flux_mask,
1404  subdomain_id,
1405  internal::always_couple_on_faces<dim, spacedim>);
1406  }
1407 
1408 
1409 
1410  template <int dim, int spacedim, typename number>
1411  void
1413  const DoFHandler<dim, spacedim> &dof,
1414  SparsityPatternBase & sparsity,
1415  const AffineConstraints<number> &constraints,
1416  const bool keep_constrained_dofs,
1417  const Table<2, Coupling> & int_mask,
1418  const Table<2, Coupling> & flux_mask,
1420  const std::function<
1422  const unsigned int)> &face_has_flux_coupling)
1423  {
1424  // do the error checking and frame code here, and then pass on to more
1425  // specialized functions in the internal namespace
1426  const types::global_dof_index n_dofs = dof.n_dofs();
1427  (void)n_dofs;
1428  const unsigned int n_comp = dof.get_fe(0).n_components();
1429  (void)n_comp;
1430 
1431  Assert(sparsity.n_rows() == n_dofs,
1432  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1433  Assert(sparsity.n_cols() == n_dofs,
1434  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1435  Assert(int_mask.n_rows() == n_comp,
1436  ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1437  Assert(int_mask.n_cols() == n_comp,
1438  ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1439  Assert(flux_mask.n_rows() == n_comp,
1440  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1441  Assert(flux_mask.n_cols() == n_comp,
1442  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1443 
1444  // If we have a distributed Triangulation only allow locally_owned
1445  // subdomain. Not setting a subdomain is also okay, because we skip
1446  // ghost cells in the loop below.
1447  if (const auto *triangulation = dynamic_cast<
1449  &dof.get_triangulation()))
1450  {
1452  (subdomain_id == triangulation->locally_owned_subdomain()),
1453  ExcMessage(
1454  "For distributed Triangulation objects and associated "
1455  "DoFHandler objects, asking for any subdomain other than the "
1456  "locally owned one does not make sense."));
1457  }
1458 
1459  Assert(
1460  face_has_flux_coupling,
1461  ExcMessage(
1462  "The function which specifies if a flux coupling occurs over a given "
1463  "face is empty."));
1464 
1466  sparsity,
1467  constraints,
1468  keep_constrained_dofs,
1469  int_mask,
1470  flux_mask,
1471  subdomain_id,
1472  face_has_flux_coupling);
1473  }
1474 
1475 } // end of namespace DoFTools
1476 
1477 
1478 // --------------------------------------------------- explicit instantiations
1479 
1480 #include "dof_tools_sparsity.inst"
1481 
1482 
1483 
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:704
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
bool has_hp_capabilities() 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
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
types::global_dof_index size_type
size_type n_rows() const
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
size_type n_cols() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int size() const
Definition: collection.h:264
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
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)
Definition: exceptions.h:1614
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
const types::boundary_id internal_face_boundary_id
Definition: types.h:277
const types::subdomain_id invalid_subdomain_id
Definition: types.h:298
const types::global_dof_index invalid_dof_index
Definition: types.h:233
unsigned int global_dof_index
Definition: types.h:82
unsigned int subdomain_id
Definition: types.h:44
unsigned short int fe_index
Definition: types.h:60
unsigned int boundary_id
Definition: types.h:141
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation