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