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