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