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