Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +0200
\(\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().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].dofs_per_cell,
172  fe_collection[f].dofs_per_cell));
173  bool_dof_mask[f].fill(false);
174  for (unsigned int i = 0; i < fe_collection[f].dofs_per_cell; ++i)
175  for (unsigned int j = 0; j < fe_collection[f].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].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().dofs_per_cell;
286  const unsigned int dofs_per_cell_col =
287  cell_col->get_fe().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().dofs_per_cell;
312  const unsigned int dofs_per_cell_col =
313  cell_col->get_fe().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().dofs_per_cell;
339  const unsigned int dofs_per_cell_col =
340  cell_col_child->get_fe().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 : GeometryInfo<dim>::face_indices())
408  if (cell->at_boundary(f))
409  {
410  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
411  dofs_on_this_face.resize(dofs_per_face);
412  cell->face(f)->get_dof_indices(dofs_on_this_face,
413  cell->active_fe_index());
414 
415  // make sparsity pattern for this cell
416  for (unsigned int i = 0; i < dofs_per_face; ++i)
417  for (unsigned int j = 0; j < dofs_per_face; ++j)
418  sparsity.add(dof_to_boundary_mapping[dofs_on_this_face[i]],
419  dof_to_boundary_mapping[dofs_on_this_face[j]]);
420  }
421  }
422 
423 
424 
425  template <int dim,
426  int spacedim,
427  typename SparsityPatternType,
428  typename number>
429  void
431  const DoFHandler<dim, spacedim> &dof,
432  const std::map<types::boundary_id, const Function<spacedim, number> *>
433  & boundary_ids,
434  const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
435  SparsityPatternType & sparsity)
436  {
437  if (dim == 1)
438  {
439  // first check left, then right boundary point
440  for (unsigned int direction = 0; direction < 2; ++direction)
441  {
442  // if this boundary is not requested, then go on with next one
443  if (boundary_ids.find(direction) == boundary_ids.end())
444  continue;
445 
446  // find active cell at that boundary: first go to left/right,
447  // then to children
449  dof.begin(0);
450  while (!cell->at_boundary(direction))
451  cell = cell->neighbor(direction);
452  while (!cell->is_active())
453  cell = cell->child(direction);
454 
455  const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
456  std::vector<types::global_dof_index> boundary_dof_boundary_indices(
457  dofs_per_vertex);
458 
459  // next get boundary mapped dof indices of boundary dofs
460  for (unsigned int i = 0; i < dofs_per_vertex; ++i)
461  boundary_dof_boundary_indices[i] =
462  dof_to_boundary_mapping[cell->vertex_dof_index(direction, i)];
463 
464  for (unsigned int i = 0; i < dofs_per_vertex; ++i)
465  sparsity.add_entries(boundary_dof_boundary_indices[i],
466  boundary_dof_boundary_indices.begin(),
467  boundary_dof_boundary_indices.end());
468  }
469  return;
470  }
471 
472  const types::global_dof_index n_dofs = dof.n_dofs();
473  (void)n_dofs;
474 
475  AssertDimension(dof_to_boundary_mapping.size(), n_dofs);
476  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
477  boundary_ids.end(),
479  Assert(sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
480  ExcDimensionMismatch(sparsity.n_rows(),
481  dof.n_boundary_dofs(boundary_ids)));
482  Assert(sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
483  ExcDimensionMismatch(sparsity.n_cols(),
484  dof.n_boundary_dofs(boundary_ids)));
485 #ifdef DEBUG
486  if (sparsity.n_rows() != 0)
487  {
488  types::global_dof_index max_element = 0;
489  for (const types::global_dof_index index : dof_to_boundary_mapping)
490  if ((index != numbers::invalid_dof_index) && (index > max_element))
491  max_element = index;
492  AssertDimension(max_element, sparsity.n_rows() - 1);
493  }
494 #endif
495 
496  std::vector<types::global_dof_index> dofs_on_this_face;
497  dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
499  cell = dof.begin_active(),
500  endc = dof.end();
501  for (; cell != endc; ++cell)
502  for (const unsigned int f : GeometryInfo<dim>::face_indices())
503  if (boundary_ids.find(cell->face(f)->boundary_id()) !=
504  boundary_ids.end())
505  {
506  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
507  dofs_on_this_face.resize(dofs_per_face);
508  cell->face(f)->get_dof_indices(dofs_on_this_face,
509  cell->active_fe_index());
510 
511  // make sparsity pattern for this cell
512  for (unsigned int i = 0; i < dofs_per_face; ++i)
513  for (unsigned int j = 0; j < dofs_per_face; ++j)
514  sparsity.add(dof_to_boundary_mapping[dofs_on_this_face[i]],
515  dof_to_boundary_mapping[dofs_on_this_face[j]]);
516  }
517  }
518 
519 
520 
521  template <int dim,
522  int spacedim,
523  typename SparsityPatternType,
524  typename number>
525  void
527  SparsityPatternType & sparsity,
528  const AffineConstraints<number> &constraints,
529  const bool keep_constrained_dofs,
531 
532  // TODO: QA: reduce the indentation level of this method..., Maier 2012
533 
534  {
535  const types::global_dof_index n_dofs = dof.n_dofs();
536  (void)n_dofs;
537 
538  AssertDimension(sparsity.n_rows(), n_dofs);
539  AssertDimension(sparsity.n_cols(), n_dofs);
540 
541  // If we have a distributed::Triangulation only allow locally_owned
542  // subdomain. Not setting a subdomain is also okay, because we skip
543  // ghost cells in the loop below.
544  Assert((dof.get_triangulation().locally_owned_subdomain() ==
546  (subdomain_id == numbers::invalid_subdomain_id) ||
547  (subdomain_id ==
548  dof.get_triangulation().locally_owned_subdomain()),
549  ExcMessage(
550  "For parallel::distributed::Triangulation objects and "
551  "associated DoF handler objects, asking for any subdomain other "
552  "than the locally owned one does not make sense."));
553 
554  std::vector<types::global_dof_index> dofs_on_this_cell;
555  std::vector<types::global_dof_index> dofs_on_other_cell;
556  dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
557  dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
559  cell = dof.begin_active(),
560  endc = dof.end();
561 
562  // TODO: in an old implementation, we used user flags before to tag
563  // faces that were already touched. this way, we could reduce the work
564  // a little bit. now, we instead add only data from one side. this
565  // should be OK, but we need to actually verify it.
566 
567  // In case we work with a distributed sparsity pattern of Trilinos
568  // type, we only have to do the work if the current cell is owned by
569  // the calling processor. Otherwise, just continue.
570  for (; cell != endc; ++cell)
571  if (((subdomain_id == numbers::invalid_subdomain_id) ||
572  (subdomain_id == cell->subdomain_id())) &&
573  cell->is_locally_owned())
574  {
575  const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell;
576  dofs_on_this_cell.resize(n_dofs_on_this_cell);
577  cell->get_dof_indices(dofs_on_this_cell);
578 
579  // make sparsity pattern for this cell. if no constraints pattern
580  // was given, then the following call acts as if simply no
581  // constraints existed
582  constraints.add_entries_local_to_global(dofs_on_this_cell,
583  sparsity,
584  keep_constrained_dofs);
585 
586  for (const unsigned int face : GeometryInfo<dim>::face_indices())
587  {
588  typename DoFHandler<dim, spacedim>::face_iterator cell_face =
589  cell->face(face);
590  const bool periodic_neighbor = cell->has_periodic_neighbor(face);
591  if (!cell->at_boundary(face) || periodic_neighbor)
592  {
594  neighbor = cell->neighbor_or_periodic_neighbor(face);
595 
596  // in 1d, we do not need to worry whether the neighbor
597  // might have children and then loop over those children.
598  // rather, we may as well go straight to the cell behind
599  // this particular cell's most terminal child
600  if (dim == 1)
601  while (neighbor->has_children())
602  neighbor = neighbor->child(face == 0 ? 1 : 0);
603 
604  if (neighbor->has_children())
605  {
606  for (unsigned int sub_nr = 0;
607  sub_nr != cell_face->number_of_children();
608  ++sub_nr)
609  {
610  const typename DoFHandler<dim, spacedim>::
611  level_cell_iterator sub_neighbor =
612  periodic_neighbor ?
613  cell->periodic_neighbor_child_on_subface(
614  face, sub_nr) :
615  cell->neighbor_child_on_subface(face, sub_nr);
616 
617  const unsigned int n_dofs_on_neighbor =
618  sub_neighbor->get_fe().dofs_per_cell;
619  dofs_on_other_cell.resize(n_dofs_on_neighbor);
620  sub_neighbor->get_dof_indices(dofs_on_other_cell);
621 
622  constraints.add_entries_local_to_global(
623  dofs_on_this_cell,
624  dofs_on_other_cell,
625  sparsity,
626  keep_constrained_dofs);
627  constraints.add_entries_local_to_global(
628  dofs_on_other_cell,
629  dofs_on_this_cell,
630  sparsity,
631  keep_constrained_dofs);
632  // only need to add this when the neighbor is not
633  // owned by the current processor, otherwise we add
634  // the entries for the neighbor there
635  if (sub_neighbor->subdomain_id() !=
636  cell->subdomain_id())
637  constraints.add_entries_local_to_global(
638  dofs_on_other_cell,
639  sparsity,
640  keep_constrained_dofs);
641  }
642  }
643  else
644  {
645  // Refinement edges are taken care of by coarser
646  // cells
647  if ((!periodic_neighbor &&
648  cell->neighbor_is_coarser(face)) ||
649  (periodic_neighbor &&
650  cell->periodic_neighbor_is_coarser(face)))
651  if (neighbor->subdomain_id() == cell->subdomain_id())
652  continue;
653 
654  const unsigned int n_dofs_on_neighbor =
655  neighbor->get_fe().dofs_per_cell;
656  dofs_on_other_cell.resize(n_dofs_on_neighbor);
657 
658  neighbor->get_dof_indices(dofs_on_other_cell);
659 
660  constraints.add_entries_local_to_global(
661  dofs_on_this_cell,
662  dofs_on_other_cell,
663  sparsity,
664  keep_constrained_dofs);
665 
666  // only need to add these in case the neighbor cell
667  // is not locally owned - otherwise, we touch each
668  // face twice and hence put the indices the other way
669  // around
670  if (!cell->neighbor_or_periodic_neighbor(face)
671  ->is_active() ||
672  (neighbor->subdomain_id() != cell->subdomain_id()))
673  {
674  constraints.add_entries_local_to_global(
675  dofs_on_other_cell,
676  dofs_on_this_cell,
677  sparsity,
678  keep_constrained_dofs);
679  if (neighbor->subdomain_id() != cell->subdomain_id())
680  constraints.add_entries_local_to_global(
681  dofs_on_other_cell,
682  sparsity,
683  keep_constrained_dofs);
684  }
685  }
686  }
687  }
688  }
689  }
690 
691 
692 
693  template <int dim, int spacedim, typename SparsityPatternType>
694  void
696  SparsityPatternType & sparsity)
697  {
699  make_flux_sparsity_pattern(dof, sparsity, dummy);
700  }
701 
702  template <int dim, int spacedim>
706  const Table<2, Coupling> & component_couplings)
707  {
708  Assert(component_couplings.n_rows() == fe.n_components(),
709  ExcDimensionMismatch(component_couplings.n_rows(),
710  fe.n_components()));
711  Assert(component_couplings.n_cols() == fe.n_components(),
712  ExcDimensionMismatch(component_couplings.n_cols(),
713  fe.n_components()));
714 
715  const unsigned int n_dofs = fe.dofs_per_cell;
716 
717  Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
718 
719  for (unsigned int i = 0; i < n_dofs; ++i)
720  {
721  const unsigned int ii =
722  (fe.is_primitive(i) ?
723  fe.system_to_component_index(i).first :
725  Assert(ii < fe.n_components(), ExcInternalError());
726 
727  for (unsigned int j = 0; j < n_dofs; ++j)
728  {
729  const unsigned int jj =
730  (fe.is_primitive(j) ?
731  fe.system_to_component_index(j).first :
733  Assert(jj < fe.n_components(), ExcInternalError());
734 
735  dof_couplings(i, j) = component_couplings(ii, jj);
736  }
737  }
738  return dof_couplings;
739  }
740 
741 
742 
743  template <int dim, int spacedim>
744  std::vector<Table<2, Coupling>>
747  const Table<2, Coupling> & component_couplings)
748  {
749  std::vector<Table<2, Coupling>> return_value(fe.size());
750  for (unsigned int i = 0; i < fe.size(); ++i)
751  return_value[i] =
752  dof_couplings_from_component_couplings(fe[i], component_couplings);
753 
754  return return_value;
755  }
756 
757 
758 
759  namespace internal
760  {
761  namespace
762  {
763  // implementation of the same function in namespace DoFTools for
764  // non-hp DoFHandlers
765  template <int dim,
766  int spacedim,
767  typename SparsityPatternType,
768  typename number>
769  void
771  const DoFHandler<dim, spacedim> &dof,
772  SparsityPatternType & sparsity,
773  const AffineConstraints<number> &constraints,
774  const bool keep_constrained_dofs,
775  const Table<2, Coupling> & int_mask,
776  const Table<2, Coupling> & flux_mask,
778  const std::function<
780  const unsigned int)> &face_has_flux_coupling)
781  {
782  if (dof.hp_capability_enabled == false)
783  {
784  const FiniteElement<dim, spacedim> &fe = dof.get_fe();
785 
786  std::vector<types::global_dof_index> dofs_on_this_cell(
787  fe.dofs_per_cell);
788  std::vector<types::global_dof_index> dofs_on_other_cell(
789  fe.dofs_per_cell);
790 
791  const Table<2, Coupling>
792  int_dof_mask =
794  flux_dof_mask =
796 
797  Table<2, bool> support_on_face(fe.dofs_per_cell,
799  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
800  for (const unsigned int f : GeometryInfo<dim>::face_indices())
801  support_on_face(i, f) = fe.has_support_on_face(i, f);
802 
803  // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
804  // to constraints.add_entries_local_to_global()
805  Table<2, bool> bool_int_dof_mask(fe.dofs_per_cell,
806  fe.dofs_per_cell);
807  bool_int_dof_mask.fill(false);
808  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
809  for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
810  if (int_dof_mask(i, j) != none)
811  bool_int_dof_mask(i, j) = true;
812 
814  cell = dof.begin_active(),
815  endc = dof.end();
816  for (; cell != endc; ++cell)
817  if (((subdomain_id == numbers::invalid_subdomain_id) ||
818  (subdomain_id == cell->subdomain_id())) &&
819  cell->is_locally_owned())
820  {
821  cell->get_dof_indices(dofs_on_this_cell);
822  // make sparsity pattern for this cell
823  constraints.add_entries_local_to_global(dofs_on_this_cell,
824  sparsity,
825  keep_constrained_dofs,
826  bool_int_dof_mask);
827  // Loop over all interior neighbors
828  for (const unsigned int face_n :
830  {
832  cell_face = cell->face(face_n);
833 
834  const bool periodic_neighbor =
835  cell->has_periodic_neighbor(face_n);
836 
837  if (cell->at_boundary(face_n) && (!periodic_neighbor))
838  {
839  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
840  {
841  const bool i_non_zero_i =
842  support_on_face(i, face_n);
843  for (unsigned int j = 0; j < fe.dofs_per_cell;
844  ++j)
845  {
846  const bool j_non_zero_i =
847  support_on_face(j, face_n);
848 
849  if (flux_dof_mask(i, j) == always ||
850  (flux_dof_mask(i, j) == nonzero &&
851  i_non_zero_i && j_non_zero_i))
852  sparsity.add(dofs_on_this_cell[i],
853  dofs_on_this_cell[j]);
854  }
855  }
856  }
857  else
858  {
859  if (!face_has_flux_coupling(cell, face_n))
860  continue;
861 
862  typename DoFHandler<dim,
863  spacedim>::level_cell_iterator
864  neighbor =
865  cell->neighbor_or_periodic_neighbor(face_n);
866  // If the cells are on the same level (and both are
867  // active, locally-owned cells) then only add to the
868  // sparsity pattern if the current cell is 'greater'
869  // in the total ordering.
870  if (neighbor->level() == cell->level() &&
871  neighbor->index() > cell->index() &&
872  neighbor->is_active() &&
873  neighbor->is_locally_owned())
874  continue;
875  // If we are more refined then the neighbor, then we
876  // will automatically find the active neighbor cell
877  // when we call 'neighbor (face_n)' above. The
878  // opposite is not true; if the neighbor is more
879  // refined then the call 'neighbor (face_n)' will
880  // *not* return an active cell. Hence, only add things
881  // to the sparsity pattern if (when the levels are
882  // different) the neighbor is coarser than the current
883  // cell.
884  //
885  // Like above, do not use this optimization if the
886  // neighbor is not locally owned.
887  if (neighbor->level() != cell->level() &&
888  ((!periodic_neighbor &&
889  !cell->neighbor_is_coarser(face_n)) ||
890  (periodic_neighbor &&
891  !cell->periodic_neighbor_is_coarser(face_n))) &&
892  neighbor->is_locally_owned())
893  continue; // (the neighbor is finer)
894 
895  const unsigned int neighbor_face_n =
896  periodic_neighbor ?
897  cell->periodic_neighbor_face_no(face_n) :
898  cell->neighbor_face_no(face_n);
899 
900 
901  // In 1D, go straight to the cell behind this
902  // particular cell's most terminal cell. This makes us
903  // skip the if (neighbor->has_children()) section
904  // below. We need to do this since we otherwise
905  // iterate over the children of the face, which are
906  // always 0 in 1D.
907  if (dim == 1)
908  while (neighbor->has_children())
909  neighbor = neighbor->child(face_n == 0 ? 1 : 0);
910 
911  if (neighbor->has_children())
912  {
913  for (unsigned int sub_nr = 0;
914  sub_nr != cell_face->n_children();
915  ++sub_nr)
916  {
917  const typename DoFHandler<dim, spacedim>::
918  level_cell_iterator sub_neighbor =
919  periodic_neighbor ?
920  cell
921  ->periodic_neighbor_child_on_subface(
922  face_n, sub_nr) :
923  cell->neighbor_child_on_subface(face_n,
924  sub_nr);
925 
926  sub_neighbor->get_dof_indices(
927  dofs_on_other_cell);
928  for (unsigned int i = 0; i < fe.dofs_per_cell;
929  ++i)
930  {
931  const bool i_non_zero_i =
932  support_on_face(i, face_n);
933  const bool i_non_zero_e =
934  support_on_face(i, neighbor_face_n);
935  for (unsigned int j = 0;
936  j < fe.dofs_per_cell;
937  ++j)
938  {
939  const bool j_non_zero_i =
940  support_on_face(j, face_n);
941  const bool j_non_zero_e =
942  support_on_face(j, neighbor_face_n);
943 
944  if (flux_dof_mask(i, j) == always)
945  {
946  sparsity.add(
947  dofs_on_this_cell[i],
948  dofs_on_other_cell[j]);
949  sparsity.add(
950  dofs_on_other_cell[i],
951  dofs_on_this_cell[j]);
952  sparsity.add(
953  dofs_on_this_cell[i],
954  dofs_on_this_cell[j]);
955  sparsity.add(
956  dofs_on_other_cell[i],
957  dofs_on_other_cell[j]);
958  }
959  else if (flux_dof_mask(i, j) ==
960  nonzero)
961  {
962  if (i_non_zero_i && j_non_zero_e)
963  sparsity.add(
964  dofs_on_this_cell[i],
965  dofs_on_other_cell[j]);
966  if (i_non_zero_e && j_non_zero_i)
967  sparsity.add(
968  dofs_on_other_cell[i],
969  dofs_on_this_cell[j]);
970  if (i_non_zero_i && j_non_zero_i)
971  sparsity.add(
972  dofs_on_this_cell[i],
973  dofs_on_this_cell[j]);
974  if (i_non_zero_e && j_non_zero_e)
975  sparsity.add(
976  dofs_on_other_cell[i],
977  dofs_on_other_cell[j]);
978  }
979 
980  if (flux_dof_mask(j, i) == always)
981  {
982  sparsity.add(
983  dofs_on_this_cell[j],
984  dofs_on_other_cell[i]);
985  sparsity.add(
986  dofs_on_other_cell[j],
987  dofs_on_this_cell[i]);
988  sparsity.add(
989  dofs_on_this_cell[j],
990  dofs_on_this_cell[i]);
991  sparsity.add(
992  dofs_on_other_cell[j],
993  dofs_on_other_cell[i]);
994  }
995  else if (flux_dof_mask(j, i) ==
996  nonzero)
997  {
998  if (j_non_zero_i && i_non_zero_e)
999  sparsity.add(
1000  dofs_on_this_cell[j],
1001  dofs_on_other_cell[i]);
1002  if (j_non_zero_e && i_non_zero_i)
1003  sparsity.add(
1004  dofs_on_other_cell[j],
1005  dofs_on_this_cell[i]);
1006  if (j_non_zero_i && i_non_zero_i)
1007  sparsity.add(
1008  dofs_on_this_cell[j],
1009  dofs_on_this_cell[i]);
1010  if (j_non_zero_e && i_non_zero_e)
1011  sparsity.add(
1012  dofs_on_other_cell[j],
1013  dofs_on_other_cell[i]);
1014  }
1015  }
1016  }
1017  }
1018  }
1019  else
1020  {
1021  neighbor->get_dof_indices(dofs_on_other_cell);
1022  for (unsigned int i = 0; i < fe.dofs_per_cell;
1023  ++i)
1024  {
1025  const bool i_non_zero_i =
1026  support_on_face(i, face_n);
1027  const bool i_non_zero_e =
1028  support_on_face(i, neighbor_face_n);
1029  for (unsigned int j = 0; j < fe.dofs_per_cell;
1030  ++j)
1031  {
1032  const bool j_non_zero_i =
1033  support_on_face(j, face_n);
1034  const bool j_non_zero_e =
1035  support_on_face(j, neighbor_face_n);
1036  if (flux_dof_mask(i, j) == always)
1037  {
1038  sparsity.add(dofs_on_this_cell[i],
1039  dofs_on_other_cell[j]);
1040  sparsity.add(dofs_on_other_cell[i],
1041  dofs_on_this_cell[j]);
1042  sparsity.add(dofs_on_this_cell[i],
1043  dofs_on_this_cell[j]);
1044  sparsity.add(dofs_on_other_cell[i],
1045  dofs_on_other_cell[j]);
1046  }
1047  if (flux_dof_mask(i, j) == nonzero)
1048  {
1049  if (i_non_zero_i && j_non_zero_e)
1050  sparsity.add(dofs_on_this_cell[i],
1051  dofs_on_other_cell[j]);
1052  if (i_non_zero_e && j_non_zero_i)
1053  sparsity.add(dofs_on_other_cell[i],
1054  dofs_on_this_cell[j]);
1055  if (i_non_zero_i && j_non_zero_i)
1056  sparsity.add(dofs_on_this_cell[i],
1057  dofs_on_this_cell[j]);
1058  if (i_non_zero_e && j_non_zero_e)
1059  sparsity.add(dofs_on_other_cell[i],
1060  dofs_on_other_cell[j]);
1061  }
1062 
1063  if (flux_dof_mask(j, i) == always)
1064  {
1065  sparsity.add(dofs_on_this_cell[j],
1066  dofs_on_other_cell[i]);
1067  sparsity.add(dofs_on_other_cell[j],
1068  dofs_on_this_cell[i]);
1069  sparsity.add(dofs_on_this_cell[j],
1070  dofs_on_this_cell[i]);
1071  sparsity.add(dofs_on_other_cell[j],
1072  dofs_on_other_cell[i]);
1073  }
1074  if (flux_dof_mask(j, i) == nonzero)
1075  {
1076  if (j_non_zero_i && i_non_zero_e)
1077  sparsity.add(dofs_on_this_cell[j],
1078  dofs_on_other_cell[i]);
1079  if (j_non_zero_e && i_non_zero_i)
1080  sparsity.add(dofs_on_other_cell[j],
1081  dofs_on_this_cell[i]);
1082  if (j_non_zero_i && i_non_zero_i)
1083  sparsity.add(dofs_on_this_cell[j],
1084  dofs_on_this_cell[i]);
1085  if (j_non_zero_e && i_non_zero_e)
1086  sparsity.add(dofs_on_other_cell[j],
1087  dofs_on_other_cell[i]);
1088  }
1089  }
1090  }
1091  }
1092  }
1093  }
1094  }
1095  }
1096  else
1097  {
1098  // while the implementation above is quite optimized and caches a
1099  // lot of data (see e.g. the int/flux_dof_mask tables), this is no
1100  // longer practical for the hp version since we would have to have
1101  // it for all combinations of elements in the hp::FECollection.
1102  // consequently, the implementation here is simpler and probably
1103  // less efficient but at least readable...
1104 
1105  const ::hp::FECollection<dim, spacedim> &fe =
1106  dof.get_fe_collection();
1107 
1108  std::vector<types::global_dof_index> dofs_on_this_cell(
1109  dof.get_fe_collection().max_dofs_per_cell());
1110  std::vector<types::global_dof_index> dofs_on_other_cell(
1111  dof.get_fe_collection().max_dofs_per_cell());
1112 
1113  const unsigned int n_components = fe.n_components();
1114  AssertDimension(int_mask.size(0), n_components);
1115  AssertDimension(int_mask.size(1), n_components);
1116  AssertDimension(flux_mask.size(0), n_components);
1117  AssertDimension(flux_mask.size(1), n_components);
1118 
1119  // note that we also need to set the respective entries if flux_mask
1120  // says so. this is necessary since we need to consider all degrees
1121  // of freedom on a cell for interior faces.
1122  Table<2, Coupling> int_and_flux_mask(n_components, n_components);
1123  for (unsigned int c1 = 0; c1 < n_components; ++c1)
1124  for (unsigned int c2 = 0; c2 < n_components; ++c2)
1125  if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
1126  int_and_flux_mask(c1, c2) = always;
1127 
1128  std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
1129  dof_couplings_from_component_couplings(fe, int_and_flux_mask);
1130 
1131  // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
1132  // can pass it to constraints.add_entries_local_to_global()
1133  std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
1134  for (unsigned int f = 0; f < fe.size(); ++f)
1135  {
1136  bool_int_and_flux_dof_mask[f].reinit(
1137  TableIndices<2>(fe[f].dofs_per_cell, fe[f].dofs_per_cell));
1138  bool_int_and_flux_dof_mask[f].fill(false);
1139  for (unsigned int i = 0; i < fe[f].dofs_per_cell; ++i)
1140  for (unsigned int j = 0; j < fe[f].dofs_per_cell; ++j)
1141  if (int_and_flux_dof_mask[f](i, j) != none)
1142  bool_int_and_flux_dof_mask[f](i, j) = true;
1143  }
1144 
1145 
1146  typename ::DoFHandler<dim, spacedim>::active_cell_iterator
1147  cell = dof.begin_active(),
1148  endc = dof.end();
1149  for (; cell != endc; ++cell)
1150  if (((subdomain_id == numbers::invalid_subdomain_id) ||
1151  (subdomain_id == cell->subdomain_id())) &&
1152  cell->is_locally_owned())
1153  {
1154  dofs_on_this_cell.resize(cell->get_fe().dofs_per_cell);
1155  cell->get_dof_indices(dofs_on_this_cell);
1156 
1157  // make sparsity pattern for this cell also taking into
1158  // account the couplings due to face contributions on the same
1159  // cell
1160  constraints.add_entries_local_to_global(
1161  dofs_on_this_cell,
1162  sparsity,
1163  keep_constrained_dofs,
1164  bool_int_and_flux_dof_mask[cell->active_fe_index()]);
1165 
1166  // Loop over interior faces
1167  for (const unsigned int face :
1169  {
1170  const typename ::DoFHandler<dim,
1171  spacedim>::face_iterator
1172  cell_face = cell->face(face);
1173 
1174  const bool periodic_neighbor =
1175  cell->has_periodic_neighbor(face);
1176 
1177  if ((!cell->at_boundary(face)) || periodic_neighbor)
1178  {
1179  typename ::DoFHandler<dim, spacedim>::
1180  level_cell_iterator neighbor =
1181  cell->neighbor_or_periodic_neighbor(face);
1182 
1183  if (!face_has_flux_coupling(cell, face))
1184  continue;
1185 
1186  // Like the non-hp case: If the cells are on the same
1187  // level (and both are active, locally-owned cells)
1188  // then only add to the sparsity pattern if the
1189  // current cell is 'greater' in the total ordering.
1190  if (neighbor->level() == cell->level() &&
1191  neighbor->index() > cell->index() &&
1192  neighbor->is_active() &&
1193  neighbor->is_locally_owned())
1194  continue;
1195  // Again, like the non-hp case: If we are more refined
1196  // then the neighbor, then we will automatically find
1197  // the active neighbor cell when we call 'neighbor
1198  // (face)' above. The opposite is not true; if the
1199  // neighbor is more refined then the call 'neighbor
1200  // (face)' will *not* return an active cell. Hence,
1201  // only add things to the sparsity pattern if (when
1202  // the levels are different) the neighbor is coarser
1203  // than the current cell.
1204  //
1205  // Like above, do not use this optimization if the
1206  // neighbor is not locally owned.
1207  if (neighbor->level() != cell->level() &&
1208  ((!periodic_neighbor &&
1209  !cell->neighbor_is_coarser(face)) ||
1210  (periodic_neighbor &&
1211  !cell->periodic_neighbor_is_coarser(face))) &&
1212  neighbor->is_locally_owned())
1213  continue; // (the neighbor is finer)
1214 
1215  // In 1D, go straight to the cell behind this
1216  // particular cell's most terminal cell. This makes us
1217  // skip the if (neighbor->has_children()) section
1218  // below. We need to do this since we otherwise
1219  // iterate over the children of the face, which are
1220  // always 0 in 1D.
1221  if (dim == 1)
1222  while (neighbor->has_children())
1223  neighbor = neighbor->child(face == 0 ? 1 : 0);
1224 
1225  if (neighbor->has_children())
1226  {
1227  for (unsigned int sub_nr = 0;
1228  sub_nr != cell_face->n_children();
1229  ++sub_nr)
1230  {
1231  const typename ::DoFHandler<dim,
1232  spacedim>::
1233  level_cell_iterator sub_neighbor =
1234  periodic_neighbor ?
1235  cell
1236  ->periodic_neighbor_child_on_subface(
1237  face, sub_nr) :
1238  cell->neighbor_child_on_subface(face,
1239  sub_nr);
1240 
1241  dofs_on_other_cell.resize(
1242  sub_neighbor->get_fe().dofs_per_cell);
1243  sub_neighbor->get_dof_indices(
1244  dofs_on_other_cell);
1245  for (unsigned int i = 0;
1246  i < cell->get_fe().dofs_per_cell;
1247  ++i)
1248  {
1249  const unsigned int ii =
1250  (cell->get_fe().is_primitive(i) ?
1251  cell->get_fe()
1252  .system_to_component_index(i)
1253  .first :
1254  cell->get_fe()
1255  .get_nonzero_components(i)
1256  .first_selected_component());
1257 
1258  Assert(ii < cell->get_fe().n_components(),
1259  ExcInternalError());
1260 
1261  for (unsigned int j = 0;
1262  j <
1263  sub_neighbor->get_fe().dofs_per_cell;
1264  ++j)
1265  {
1266  const unsigned int jj =
1267  (sub_neighbor->get_fe()
1268  .is_primitive(j) ?
1269  sub_neighbor->get_fe()
1270  .system_to_component_index(j)
1271  .first :
1272  sub_neighbor->get_fe()
1273  .get_nonzero_components(j)
1274  .first_selected_component());
1275 
1276  Assert(jj < sub_neighbor->get_fe()
1277  .n_components(),
1278  ExcInternalError());
1279 
1280  if ((flux_mask(ii, jj) == always) ||
1281  (flux_mask(ii, jj) == nonzero))
1282  {
1283  sparsity.add(
1284  dofs_on_this_cell[i],
1285  dofs_on_other_cell[j]);
1286  }
1287 
1288  if ((flux_mask(jj, ii) == always) ||
1289  (flux_mask(jj, ii) == nonzero))
1290  {
1291  sparsity.add(
1292  dofs_on_other_cell[j],
1293  dofs_on_this_cell[i]);
1294  }
1295  }
1296  }
1297  }
1298  }
1299  else
1300  {
1301  dofs_on_other_cell.resize(
1302  neighbor->get_fe().dofs_per_cell);
1303  neighbor->get_dof_indices(dofs_on_other_cell);
1304  for (unsigned int i = 0;
1305  i < cell->get_fe().dofs_per_cell;
1306  ++i)
1307  {
1308  const unsigned int ii =
1309  (cell->get_fe().is_primitive(i) ?
1310  cell->get_fe()
1311  .system_to_component_index(i)
1312  .first :
1313  cell->get_fe()
1314  .get_nonzero_components(i)
1315  .first_selected_component());
1316 
1317  Assert(ii < cell->get_fe().n_components(),
1318  ExcInternalError());
1319 
1320  for (unsigned int j = 0;
1321  j < neighbor->get_fe().dofs_per_cell;
1322  ++j)
1323  {
1324  const unsigned int jj =
1325  (neighbor->get_fe().is_primitive(j) ?
1326  neighbor->get_fe()
1327  .system_to_component_index(j)
1328  .first :
1329  neighbor->get_fe()
1330  .get_nonzero_components(j)
1331  .first_selected_component());
1332 
1333  Assert(
1334  jj < neighbor->get_fe().n_components(),
1335  ExcInternalError());
1336 
1337  if ((flux_mask(ii, jj) == always) ||
1338  (flux_mask(ii, jj) == nonzero))
1339  {
1340  sparsity.add(dofs_on_this_cell[i],
1341  dofs_on_other_cell[j]);
1342  }
1343 
1344  if ((flux_mask(jj, ii) == always) ||
1345  (flux_mask(jj, ii) == nonzero))
1346  {
1347  sparsity.add(dofs_on_other_cell[j],
1348  dofs_on_this_cell[i]);
1349  }
1350  }
1351  }
1352  }
1353  }
1354  }
1355  }
1356  }
1357  }
1358  } // namespace
1359 
1360  } // namespace internal
1361 
1362 
1363 
1364  template <int dim, int spacedim, typename SparsityPatternType>
1365  void
1367  SparsityPatternType & sparsity,
1368  const Table<2, Coupling> & int_mask,
1369  const Table<2, Coupling> & flux_mask,
1371  {
1373 
1374  const bool keep_constrained_dofs = true;
1375 
1377  sparsity,
1378  dummy,
1379  keep_constrained_dofs,
1380  int_mask,
1381  flux_mask,
1382  subdomain_id,
1383  internal::always_couple_on_faces<dim, spacedim>);
1384  }
1385 
1386 
1387 
1388  template <int dim,
1389  int spacedim,
1390  typename SparsityPatternType,
1391  typename number>
1392  void
1394  const DoFHandler<dim, spacedim> &dof,
1395  SparsityPatternType & sparsity,
1396  const AffineConstraints<number> &constraints,
1397  const bool keep_constrained_dofs,
1398  const Table<2, Coupling> & int_mask,
1399  const Table<2, Coupling> & flux_mask,
1401  const std::function<
1403  const unsigned int)> &face_has_flux_coupling)
1404  {
1405  // do the error checking and frame code here, and then pass on to more
1406  // specialized functions in the internal namespace
1407  const types::global_dof_index n_dofs = dof.n_dofs();
1408  (void)n_dofs;
1409  const unsigned int n_comp = dof.get_fe(0).n_components();
1410  (void)n_comp;
1411 
1412  Assert(sparsity.n_rows() == n_dofs,
1413  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1414  Assert(sparsity.n_cols() == n_dofs,
1415  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1416  Assert(int_mask.n_rows() == n_comp,
1417  ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1418  Assert(int_mask.n_cols() == n_comp,
1419  ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1420  Assert(flux_mask.n_rows() == n_comp,
1421  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1422  Assert(flux_mask.n_cols() == n_comp,
1423  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1424 
1425  // If we have a distributed::Triangulation only allow locally_owned
1426  // subdomain. Not setting a subdomain is also okay, because we skip
1427  // ghost cells in the loop below.
1428  Assert((dof.get_triangulation().locally_owned_subdomain() ==
1430  (subdomain_id == numbers::invalid_subdomain_id) ||
1431  (subdomain_id ==
1432  dof.get_triangulation().locally_owned_subdomain()),
1433  ExcMessage(
1434  "For parallel::distributed::Triangulation objects and "
1435  "associated DoF handler objects, asking for any subdomain other "
1436  "than the locally owned one does not make sense."));
1437 
1438  Assert(
1439  face_has_flux_coupling,
1440  ExcMessage(
1441  "The function which specifies if a flux coupling occurs over a given "
1442  "face is empty."));
1443 
1445  sparsity,
1446  constraints,
1447  keep_constrained_dofs,
1448  int_mask,
1449  flux_mask,
1450  subdomain_id,
1451  face_has_flux_coupling);
1452  }
1453 
1454 } // end of namespace DoFTools
1455 
1456 
1457 // --------------------------------------------------- explicit instantiations
1458 
1459 #include "dof_tools_sparsity.inst"
1460 
1461 
1462 
const bool hp_capability_enabled
Definition: dof_handler.h:425
const Triangulation< dim, spacedim > & get_triangulation() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
const types::subdomain_id invalid_subdomain_id
Definition: types.h:285
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:3253
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:389
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
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:1403
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:3233
types::global_dof_index n_boundary_dofs() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
types::global_dof_index n_dofs() const
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:407
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
const unsigned int dofs_per_cell
Definition: fe_base.h:280
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:3062
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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)
Definition: table.h:687
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:371
const types::boundary_id internal_face_boundary_id
Definition: types.h:250
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:343
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:206
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
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:1078
static ::ExceptionBase & ExcInternalError()