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