deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40:00+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2013 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16#include <deal.II/base/table.h>
19
22
26
27#include <deal.II/fe/fe.h>
29#include <deal.II/fe/fe_tools.h>
31
34#include <deal.II/grid/tria.h>
36
40
43#include <deal.II/lac/vector.h>
44
45#include <algorithm>
46#include <numeric>
47
49
50
51
52namespace DoFTools
53{
54 template <int dim, int spacedim, typename number>
55 void
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<
75 &dof.get_triangulation()))
76 {
77 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
78 (subdomain_id == triangulation->locally_owned_subdomain()),
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
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<
139 &dof.get_triangulation()))
140 {
141 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
142 (subdomain_id == triangulation->locally_owned_subdomain()),
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
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 ||
245 &dof_col.get_triangulation()) != nullptr)
246 {
247 const types::subdomain_id this_subdomain_id =
249 Assert(this_subdomain_id ==
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<
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 {
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<
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 {
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
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
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
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(),
470
471 const bool fe_is_hermite = (dynamic_cast<const FE_Hermite<dim, spacedim> *>(
472 &(dof.get_fe())) != nullptr);
473
474 Assert(fe_is_hermite ||
475 sparsity.n_rows() == dof.n_boundary_dofs(boundary_ids),
476 ExcDimensionMismatch(sparsity.n_rows(),
477 dof.n_boundary_dofs(boundary_ids)));
478 Assert(fe_is_hermite ||
479 sparsity.n_cols() == dof.n_boundary_dofs(boundary_ids),
480 ExcDimensionMismatch(sparsity.n_cols(),
481 dof.n_boundary_dofs(boundary_ids)));
482 (void)fe_is_hermite;
483
484#ifdef DEBUG
485 if (sparsity.n_rows() != 0)
486 {
487 types::global_dof_index max_element = 0;
488 for (const types::global_dof_index index : dof_to_boundary_mapping)
489 if ((index != numbers::invalid_dof_index) && (index > max_element))
490 max_element = index;
491 AssertDimension(max_element, sparsity.n_rows() - 1);
492 }
493#endif
494
495 std::vector<types::global_dof_index> dofs_on_this_face;
496 dofs_on_this_face.reserve(dof.get_fe_collection().max_dofs_per_face());
497 std::vector<types::global_dof_index> cols;
498
499 for (const auto &cell : dof.active_cell_iterators())
500 for (const unsigned int f : cell->face_indices())
501 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
502 boundary_ids.end())
503 {
504 const unsigned int dofs_per_face =
505 cell->get_fe().n_dofs_per_face(f);
506 dofs_on_this_face.resize(dofs_per_face);
507 cell->face(f)->get_dof_indices(dofs_on_this_face,
508 cell->active_fe_index());
509
510 // make sparsity pattern for this cell
511 cols.clear();
512 for (const auto &dof : dofs_on_this_face)
513 cols.push_back(dof_to_boundary_mapping[dof]);
514
515 // Like the other one: sort once.
516 std::sort(cols.begin(), cols.end());
517 for (const auto &dof : dofs_on_this_face)
518 sparsity.add_row_entries(dof_to_boundary_mapping[dof],
519 make_array_view(cols),
520 true);
521 }
522 }
523
524
525
526 template <int dim, int spacedim, typename number>
527 void
529 SparsityPatternBase &sparsity,
530 const AffineConstraints<number> &constraints,
531 const bool keep_constrained_dofs,
532 const types::subdomain_id subdomain_id)
533
534 // TODO: QA: reduce the indentation level of this method..., Maier 2012
535
536 {
537 const types::global_dof_index n_dofs = dof.n_dofs();
538 (void)n_dofs;
539
540 AssertDimension(sparsity.n_rows(), n_dofs);
541 AssertDimension(sparsity.n_cols(), n_dofs);
542
543 // If we have a distributed Triangulation only allow locally_owned
544 // subdomain. Not setting a subdomain is also okay, because we skip
545 // ghost cells in the loop below.
546 if (const auto *triangulation = dynamic_cast<
548 &dof.get_triangulation()))
549 {
550 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
551 (subdomain_id == triangulation->locally_owned_subdomain()),
553 "For distributed Triangulation objects and associated "
554 "DoFHandler objects, asking for any subdomain other than the "
555 "locally owned one does not make sense."));
556 }
557
558 std::vector<types::global_dof_index> dofs_on_this_cell;
559 std::vector<types::global_dof_index> dofs_on_other_cell;
560 dofs_on_this_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
561 dofs_on_other_cell.reserve(dof.get_fe_collection().max_dofs_per_cell());
562
563 // TODO: in an old implementation, we used user flags before to tag
564 // faces that were already touched. this way, we could reduce the work
565 // a little bit. now, we instead add only data from one side. this
566 // should be OK, but we need to actually verify it.
567
568 // In case we work with a distributed sparsity pattern of Trilinos
569 // type, we only have to do the work if the current cell is owned by
570 // the calling processor. Otherwise, just continue.
571 for (const auto &cell : dof.active_cell_iterators())
572 if (((subdomain_id == numbers::invalid_subdomain_id) ||
573 (subdomain_id == cell->subdomain_id())) &&
574 cell->is_locally_owned())
575 {
576 const unsigned int n_dofs_on_this_cell =
577 cell->get_fe().n_dofs_per_cell();
578 dofs_on_this_cell.resize(n_dofs_on_this_cell);
579 cell->get_dof_indices(dofs_on_this_cell);
580
581 // make sparsity pattern for this cell. if no constraints pattern
582 // was given, then the following call acts as if simply no
583 // constraints existed
584 constraints.add_entries_local_to_global(dofs_on_this_cell,
585 sparsity,
586 keep_constrained_dofs);
587
588 for (const unsigned int face : cell->face_indices())
589 {
591 cell->face(face);
592 const bool periodic_neighbor = cell->has_periodic_neighbor(face);
593 if (!cell->at_boundary(face) || periodic_neighbor)
594 {
596 neighbor = cell->neighbor_or_periodic_neighbor(face);
597
598 // in 1d, we do not need to worry whether the neighbor
599 // might have children and then loop over those children.
600 // rather, we may as well go straight to the cell behind
601 // this particular cell's most terminal child
602 if (dim == 1)
603 while (neighbor->has_children())
604 neighbor = neighbor->child(face == 0 ? 1 : 0);
605
606 if (neighbor->has_children())
607 {
608 for (unsigned int sub_nr = 0;
609 sub_nr != cell_face->n_active_descendants();
610 ++sub_nr)
611 {
612 const typename DoFHandler<dim, spacedim>::
613 level_cell_iterator sub_neighbor =
614 periodic_neighbor ?
615 cell->periodic_neighbor_child_on_subface(
616 face, sub_nr) :
617 cell->neighbor_child_on_subface(face, sub_nr);
618
619 const unsigned int n_dofs_on_neighbor =
620 sub_neighbor->get_fe().n_dofs_per_cell();
621 dofs_on_other_cell.resize(n_dofs_on_neighbor);
622 sub_neighbor->get_dof_indices(dofs_on_other_cell);
623
624 constraints.add_entries_local_to_global(
625 dofs_on_this_cell,
626 dofs_on_other_cell,
627 sparsity,
628 keep_constrained_dofs);
629 constraints.add_entries_local_to_global(
630 dofs_on_other_cell,
631 dofs_on_this_cell,
632 sparsity,
633 keep_constrained_dofs);
634 // only need to add this when the neighbor is not
635 // owned by the current processor, otherwise we add
636 // the entries for the neighbor there
637 if (sub_neighbor->subdomain_id() !=
638 cell->subdomain_id())
639 constraints.add_entries_local_to_global(
640 dofs_on_other_cell,
641 sparsity,
642 keep_constrained_dofs);
643 }
644 }
645 else
646 {
647 // Refinement edges are taken care of by coarser
648 // cells
649 if ((!periodic_neighbor &&
650 cell->neighbor_is_coarser(face)) ||
651 (periodic_neighbor &&
652 cell->periodic_neighbor_is_coarser(face)))
653 if (neighbor->subdomain_id() == cell->subdomain_id())
654 continue;
655
656 const unsigned int n_dofs_on_neighbor =
657 neighbor->get_fe().n_dofs_per_cell();
658 dofs_on_other_cell.resize(n_dofs_on_neighbor);
659
660 neighbor->get_dof_indices(dofs_on_other_cell);
661
662 constraints.add_entries_local_to_global(
663 dofs_on_this_cell,
664 dofs_on_other_cell,
665 sparsity,
666 keep_constrained_dofs);
667
668 // only need to add these in case the neighbor cell
669 // is not locally owned - otherwise, we touch each
670 // face twice and hence put the indices the other way
671 // around
672 if (!cell->neighbor_or_periodic_neighbor(face)
673 ->is_active() ||
674 (neighbor->subdomain_id() != cell->subdomain_id()))
675 {
676 constraints.add_entries_local_to_global(
677 dofs_on_other_cell,
678 dofs_on_this_cell,
679 sparsity,
680 keep_constrained_dofs);
681 if (neighbor->subdomain_id() != cell->subdomain_id())
682 constraints.add_entries_local_to_global(
683 dofs_on_other_cell,
684 sparsity,
685 keep_constrained_dofs);
686 }
687 }
688 }
689 }
690 }
691 }
692
693
694
695 template <int dim, int spacedim>
696 void
703
704 template <int dim, int spacedim>
708 const Table<2, Coupling> &component_couplings)
709 {
710 Assert(component_couplings.n_rows() == fe.n_components(),
711 ExcDimensionMismatch(component_couplings.n_rows(),
712 fe.n_components()));
713 Assert(component_couplings.n_cols() == fe.n_components(),
714 ExcDimensionMismatch(component_couplings.n_cols(),
715 fe.n_components()));
716
717 const unsigned int n_dofs = fe.n_dofs_per_cell();
718
719 Table<2, Coupling> dof_couplings(n_dofs, n_dofs);
720
721 for (unsigned int i = 0; i < n_dofs; ++i)
722 {
723 const unsigned int ii =
724 (fe.is_primitive(i) ?
725 fe.system_to_component_index(i).first :
728
729 for (unsigned int j = 0; j < n_dofs; ++j)
730 {
731 const unsigned int jj =
732 (fe.is_primitive(j) ?
733 fe.system_to_component_index(j).first :
736
737 dof_couplings(i, j) = component_couplings(ii, jj);
738 }
739 }
740 return dof_couplings;
741 }
742
743
744
745 template <int dim, int spacedim>
746 std::vector<Table<2, Coupling>>
749 const Table<2, Coupling> &component_couplings)
750 {
751 std::vector<Table<2, Coupling>> return_value(fe.size());
752 for (unsigned int i = 0; i < fe.size(); ++i)
753 return_value[i] =
754 dof_couplings_from_component_couplings(fe[i], component_couplings);
755
756 return return_value;
757 }
758
759
760
761 namespace internal
762 {
763 namespace
764 {
765 // helper function
766 template <typename Iterator, typename Iterator2>
767 void
768 add_cell_entries(
769 const Iterator &cell,
770 const unsigned int face_no,
771 const Iterator2 &neighbor,
772 const unsigned int neighbor_face_no,
773 const Table<2, Coupling> &flux_mask,
774 const std::vector<types::global_dof_index> &dofs_on_this_cell,
775 std::vector<types::global_dof_index> &dofs_on_other_cell,
776 std::vector<std::pair<SparsityPatternBase::size_type,
777 SparsityPatternBase::size_type>> &cell_entries)
778 {
779 dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell());
780 neighbor->get_dof_indices(dofs_on_other_cell);
781
782 // Keep expensive data structures in separate vectors for inner j loop
783 // in separate vectors
784 boost::container::small_vector<unsigned int, 64>
785 component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell());
786 boost::container::small_vector<bool, 64> support_on_face_i(
787 neighbor->get_fe().n_dofs_per_cell());
788 boost::container::small_vector<bool, 64> support_on_face_e(
789 neighbor->get_fe().n_dofs_per_cell());
790 for (unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j)
791 {
792 component_indices_neighbor[j] =
793 (neighbor->get_fe().is_primitive(j) ?
794 neighbor->get_fe().system_to_component_index(j).first :
795 neighbor->get_fe()
796 .get_nonzero_components(j)
797 .first_selected_component());
798 support_on_face_i[j] =
799 neighbor->get_fe().has_support_on_face(j, face_no);
800 support_on_face_e[j] =
801 neighbor->get_fe().has_support_on_face(j, neighbor_face_no);
802 }
803
804 // For the parallel setting, must include also the diagonal
805 // neighbor-neighbor coupling, otherwise those get included on the
806 // other cell
807 for (int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); ++f)
808 {
809 const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe();
810 const auto &dofs_i =
811 (f == 0) ? dofs_on_this_cell : dofs_on_other_cell;
812 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
813 {
814 const unsigned int ii =
815 (fe.is_primitive(i) ?
816 fe.system_to_component_index(i).first :
817 fe.get_nonzero_components(i).first_selected_component());
818
819 Assert(ii < fe.n_components(), ExcInternalError());
820 const bool i_non_zero_i =
821 fe.has_support_on_face(i,
822 (f == 0 ? face_no : neighbor_face_no));
823
824 for (unsigned int j = 0;
825 j < neighbor->get_fe().n_dofs_per_cell();
826 ++j)
827 {
828 const bool j_non_zero_e = support_on_face_e[j];
829 const unsigned int jj = component_indices_neighbor[j];
830
831 Assert(jj < neighbor->get_fe().n_components(),
833
834 if ((flux_mask(ii, jj) == always) ||
835 (flux_mask(ii, jj) == nonzero && i_non_zero_i &&
836 j_non_zero_e))
837 cell_entries.emplace_back(dofs_i[i],
838 dofs_on_other_cell[j]);
839 if ((flux_mask(jj, ii) == always) ||
840 (flux_mask(jj, ii) == nonzero && j_non_zero_e &&
841 i_non_zero_i))
842 cell_entries.emplace_back(dofs_on_other_cell[j],
843 dofs_i[i]);
844 }
845 }
846 }
847 }
848
849
850
851 // implementation of the same function in namespace DoFTools for
852 // non-hp-DoFHandlers
853 template <int dim, int spacedim, typename number>
854 void
856 const DoFHandler<dim, spacedim> &dof,
857 SparsityPatternBase &sparsity,
858 const AffineConstraints<number> &constraints,
859 const bool keep_constrained_dofs,
860 const Table<2, Coupling> &int_mask,
861 const Table<2, Coupling> &flux_mask,
862 const types::subdomain_id subdomain_id,
863 const std::function<
865 const unsigned int)> &face_has_flux_coupling)
866 {
867 std::vector<types::global_dof_index> rows;
868 std::vector<std::pair<SparsityPatternBase::size_type,
870 cell_entries;
871
872 const ::hp::FECollection<dim, spacedim> &fe =
873 dof.get_fe_collection();
874
875 std::vector<types::global_dof_index> dofs_on_this_cell(
877 std::vector<types::global_dof_index> dofs_on_other_cell(
879
880 const unsigned int n_components = fe.n_components();
881 AssertDimension(int_mask.size(0), n_components);
882 AssertDimension(int_mask.size(1), n_components);
883 AssertDimension(flux_mask.size(0), n_components);
884 AssertDimension(flux_mask.size(1), n_components);
885
886 // note that we also need to set the respective entries if flux_mask
887 // says so. this is necessary since we need to consider all degrees
888 // of freedom on a cell for interior faces.
889 Table<2, Coupling> int_and_flux_mask(n_components, n_components);
890 for (unsigned int c1 = 0; c1 < n_components; ++c1)
891 for (unsigned int c2 = 0; c2 < n_components; ++c2)
892 if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none)
893 int_and_flux_mask(c1, c2) = always;
894
895 // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
896 // to constraints.add_entries_local_to_global()
897 std::vector<Table<2, Coupling>> int_and_flux_dof_mask =
898 dof_couplings_from_component_couplings(fe, int_and_flux_mask);
899
900 // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
901 // can pass it to constraints.add_entries_local_to_global()
902 std::vector<Table<2, bool>> bool_int_and_flux_dof_mask(fe.size());
903 for (unsigned int f = 0; f < fe.size(); ++f)
904 {
905 bool_int_and_flux_dof_mask[f].reinit(
906 TableIndices<2>(fe[f].n_dofs_per_cell(),
907 fe[f].n_dofs_per_cell()));
908 bool_int_and_flux_dof_mask[f].fill(false);
909 for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
910 for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j)
911 if (int_and_flux_dof_mask[f](i, j) != none)
912 bool_int_and_flux_dof_mask[f](i, j) = true;
913 }
914
915
916 for (const auto &cell : dof.active_cell_iterators())
918 (subdomain_id == cell->subdomain_id())) &&
919 cell->is_locally_owned())
920 {
921 dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell());
922 cell->get_dof_indices(dofs_on_this_cell);
923
924 // make sparsity pattern for this cell also taking into
925 // account the couplings due to face contributions on the same
926 // cell
927 constraints.add_entries_local_to_global(
928 dofs_on_this_cell,
929 sparsity,
930 keep_constrained_dofs,
931 bool_int_and_flux_dof_mask[cell->active_fe_index()]);
932
933 // Loop over interior faces
934 for (const unsigned int face : cell->face_indices())
935 {
936 const bool periodic_neighbor =
937 cell->has_periodic_neighbor(face);
938
939 if ((!cell->at_boundary(face)) || periodic_neighbor)
940 {
942 neighbor = cell->neighbor_or_periodic_neighbor(face);
943
944 // If the cells are on the same level (and both are
945 // active, locally-owned cells) then only add to the
946 // sparsity pattern if the current cell is 'greater' in
947 // the total ordering.
948 if (neighbor->level() == cell->level() &&
949 neighbor->index() > cell->index() &&
950 neighbor->is_active() && neighbor->is_locally_owned())
951 continue;
952
953 // If we are more refined then the neighbor, then we
954 // will automatically find the active neighbor cell when
955 // we call 'neighbor (face)' above. The opposite is not
956 // true; if the neighbor is more refined then the call
957 // 'neighbor (face)' will *not* return an active
958 // cell. Hence, only add things to the sparsity pattern
959 // if (when the levels are different) the neighbor is
960 // coarser than the current cell, except in the case
961 // when the neighbor is not locally owned.
962 if (neighbor->level() != cell->level() &&
963 ((!periodic_neighbor &&
964 !cell->neighbor_is_coarser(face)) ||
965 (periodic_neighbor &&
966 !cell->periodic_neighbor_is_coarser(face))) &&
967 neighbor->is_locally_owned())
968 continue; // (the neighbor is finer)
969
970 if (!face_has_flux_coupling(cell, face))
971 continue;
972
973 const unsigned int neighbor_face_no =
974 periodic_neighbor ?
975 cell->periodic_neighbor_face_no(face) :
976 cell->neighbor_face_no(face);
977
978 // In 1d, go straight to the cell behind this
979 // particular cell's most terminal cell. This makes us
980 // skip the if (neighbor->has_children()) section
981 // below. We need to do this since we otherwise
982 // iterate over the children of the face, which are
983 // always 0 in 1d.
984 if (dim == 1)
985 while (neighbor->has_children())
986 neighbor = neighbor->child(face == 0 ? 1 : 0);
987
988 if (neighbor->has_children())
989 {
990 for (unsigned int sub_nr = 0;
991 sub_nr != cell->face(face)->n_children();
992 ++sub_nr)
993 {
994 const typename DoFHandler<dim, spacedim>::
995 level_cell_iterator sub_neighbor =
996 periodic_neighbor ?
997 cell->periodic_neighbor_child_on_subface(
998 face, sub_nr) :
999 cell->neighbor_child_on_subface(face,
1000 sub_nr);
1001 add_cell_entries(cell,
1002 face,
1003 sub_neighbor,
1004 neighbor_face_no,
1005 flux_mask,
1006 dofs_on_this_cell,
1007 dofs_on_other_cell,
1008 cell_entries);
1009 }
1010 }
1011 else
1012 add_cell_entries(cell,
1013 face,
1014 neighbor,
1015 neighbor_face_no,
1016 flux_mask,
1017 dofs_on_this_cell,
1018 dofs_on_other_cell,
1019 cell_entries);
1020 }
1021 }
1022 sparsity.add_entries(make_array_view(cell_entries));
1023 cell_entries.clear();
1024 }
1025 }
1026 } // namespace
1027
1028 } // namespace internal
1029
1030
1031
1032 template <int dim, int spacedim>
1033 void
1035 SparsityPatternBase &sparsity,
1036 const Table<2, Coupling> &int_mask,
1037 const Table<2, Coupling> &flux_mask,
1038 const types::subdomain_id subdomain_id)
1039 {
1041
1042 const bool keep_constrained_dofs = true;
1043
1045 sparsity,
1046 dummy,
1047 keep_constrained_dofs,
1048 int_mask,
1049 flux_mask,
1050 subdomain_id,
1051 internal::always_couple_on_faces<dim, spacedim>);
1052 }
1053
1054
1055
1056 template <int dim, int spacedim, typename number>
1057 void
1059 const DoFHandler<dim, spacedim> &dof,
1060 SparsityPatternBase &sparsity,
1061 const AffineConstraints<number> &constraints,
1062 const bool keep_constrained_dofs,
1063 const Table<2, Coupling> &int_mask,
1064 const Table<2, Coupling> &flux_mask,
1065 const types::subdomain_id subdomain_id,
1066 const std::function<
1068 const unsigned int)> &face_has_flux_coupling)
1069 {
1070 // do the error checking and frame code here, and then pass on to more
1071 // specialized functions in the internal namespace
1072 const types::global_dof_index n_dofs = dof.n_dofs();
1073 (void)n_dofs;
1074 const unsigned int n_comp = dof.get_fe(0).n_components();
1075 (void)n_comp;
1076
1077 Assert(sparsity.n_rows() == n_dofs,
1078 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1079 Assert(sparsity.n_cols() == n_dofs,
1080 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1081 Assert(int_mask.n_rows() == n_comp,
1082 ExcDimensionMismatch(int_mask.n_rows(), n_comp));
1083 Assert(int_mask.n_cols() == n_comp,
1084 ExcDimensionMismatch(int_mask.n_cols(), n_comp));
1085 Assert(flux_mask.n_rows() == n_comp,
1086 ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
1087 Assert(flux_mask.n_cols() == n_comp,
1088 ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
1089
1090 // If we have a distributed Triangulation only allow locally_owned
1091 // subdomain. Not setting a subdomain is also okay, because we skip
1092 // ghost cells in the loop below.
1093 if (const auto *triangulation = dynamic_cast<
1095 &dof.get_triangulation()))
1096 {
1097 Assert((subdomain_id == numbers::invalid_subdomain_id) ||
1098 (subdomain_id == triangulation->locally_owned_subdomain()),
1099 ExcMessage(
1100 "For distributed Triangulation objects and associated "
1101 "DoFHandler objects, asking for any subdomain other than the "
1102 "locally owned one does not make sense."));
1103 }
1104
1105 Assert(
1106 face_has_flux_coupling,
1107 ExcMessage(
1108 "The function which specifies if a flux coupling occurs over a given "
1109 "face is empty."));
1110
1111 internal::make_flux_sparsity_pattern(dof,
1112 sparsity,
1113 constraints,
1114 keep_constrained_dofs,
1115 int_mask,
1116 flux_mask,
1117 subdomain_id,
1118 face_has_flux_coupling);
1119 }
1120
1121} // end of namespace DoFTools
1122
1123
1124// --------------------------------------------------- explicit instantiations
1125
1126#include "dof_tools_sparsity.inst"
1127
1128
1129
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
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
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
typename LevelSelector::cell_iterator level_cell_iterator
cell_iterator begin(const unsigned int level=0) const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
types::global_dof_index size_type
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false)=0
size_type n_cols() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int size() const
Definition collection.h:315
unsigned int max_dofs_per_face() const
unsigned int max_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &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)
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id invalid_subdomain_id
Definition types.h:341
const types::global_dof_index invalid_dof_index
Definition types.h:252
unsigned int subdomain_id
Definition types.h:43
unsigned short int fe_index
Definition types.h:59
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation