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