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