Reference documentation for deal.II version 9.4.1
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
dof_handler.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2021 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/config.h>
17
18#include <deal.II/base/geometry_info.h>
19#include <deal.II/base/memory_consumption.h>
20#include <deal.II/base/mpi.templates.h>
21
22#include <deal.II/distributed/cell_data_transfer.templates.h>
23#include <deal.II/distributed/fully_distributed_tria.h>
24#include <deal.II/distributed/shared_tria.h>
25#include <deal.II/distributed/tria.h>
26
27#include <deal.II/dofs/dof_handler.h>
28#include <deal.II/dofs/dof_handler_policy.h>
29
30#include <deal.II/grid/filtered_iterator.h>
31#include <deal.II/grid/grid_tools.h>
32#include <deal.II/grid/tria.h>
33#include <deal.II/grid/tria_accessor.h>
34#include <deal.II/grid/tria_iterator.h>
35#include <deal.II/grid/tria_levels.h>
36
37#include <algorithm>
38#include <memory>
39#include <set>
40#include <unordered_set>
41
42DEAL_II_NAMESPACE_OPEN
43
44template <int dim, int spacedim>
45const unsigned int DoFHandler<dim, spacedim>::default_fe_index;
46
47
48template <int dim, int spacedim>
49const typename DoFHandler<dim, spacedim>::active_fe_index_type
50 DoFHandler<dim, spacedim>::invalid_active_fe_index;
51
52
53namespace internal
54{
55 template <int dim, int spacedim>
56 std::string
57 policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
58 PolicyBase<dim, spacedim> &policy)
59 {
60 std::string policy_name;
61 if (dynamic_cast<const typename ::internal::DoFHandlerImplementation::
62 Policy::Sequential<dim, spacedim> *>(&policy) ||
63 dynamic_cast<const typename ::internal::DoFHandlerImplementation::
64 Policy::Sequential<dim, spacedim> *>(&policy))
65 policy_name = "Policy::Sequential<";
66 else if (dynamic_cast<
67 const typename ::internal::DoFHandlerImplementation::
68 Policy::ParallelDistributed<dim, spacedim> *>(&policy) ||
69 dynamic_cast<
70 const typename ::internal::DoFHandlerImplementation::
71 Policy::ParallelDistributed<dim, spacedim> *>(&policy))
72 policy_name = "Policy::ParallelDistributed<";
73 else if (dynamic_cast<
74 const typename ::internal::DoFHandlerImplementation::
75 Policy::ParallelShared<dim, spacedim> *>(&policy) ||
76 dynamic_cast<
77 const typename ::internal::DoFHandlerImplementation::
78 Policy::ParallelShared<dim, spacedim> *>(&policy))
79 policy_name = "Policy::ParallelShared<";
80 else
81 AssertThrow(false, ExcNotImplemented());
82 policy_name += Utilities::int_to_string(dim) + "," +
83 Utilities::int_to_string(spacedim) + ">";
84 return policy_name;
85 }
86
87
88 namespace DoFHandlerImplementation
89 {
94 struct Implementation
95 {
100 template <int spacedim>
101 static unsigned int
102 max_couplings_between_dofs(const DoFHandler<1, spacedim> &dof_handler)
103 {
104 return std::min(static_cast<types::global_dof_index>(
105 3 * dof_handler.fe_collection.max_dofs_per_vertex() +
106 2 * dof_handler.fe_collection.max_dofs_per_line()),
107 dof_handler.n_dofs());
108 }
109
110 template <int spacedim>
111 static unsigned int
112 max_couplings_between_dofs(const DoFHandler<2, spacedim> &dof_handler)
113 {
114 // get these numbers by drawing pictures
115 // and counting...
116 // example:
117 // | | |
118 // --x-----x--x--X--
119 // | | | |
120 // | x--x--x
121 // | | | |
122 // --x--x--*--x--x--
123 // | | | |
124 // x--x--x |
125 // | | | |
126 // --X--x--x-----x--
127 // | | |
128 // x = vertices connected with center vertex *;
129 // = total of 19
130 // (the X vertices are connected with * if
131 // the vertices adjacent to X are hanging
132 // nodes)
133 // count lines -> 28 (don't forget to count
134 // mother and children separately!)
135 types::global_dof_index max_couplings;
136 switch (dof_handler.tria->max_adjacent_cells())
137 {
138 case 4:
139 max_couplings =
140 19 * dof_handler.fe_collection.max_dofs_per_vertex() +
141 28 * dof_handler.fe_collection.max_dofs_per_line() +
142 8 * dof_handler.fe_collection.max_dofs_per_quad();
143 break;
144 case 5:
145 max_couplings =
146 21 * dof_handler.fe_collection.max_dofs_per_vertex() +
147 31 * dof_handler.fe_collection.max_dofs_per_line() +
148 9 * dof_handler.fe_collection.max_dofs_per_quad();
149 break;
150 case 6:
151 max_couplings =
152 28 * dof_handler.fe_collection.max_dofs_per_vertex() +
153 42 * dof_handler.fe_collection.max_dofs_per_line() +
154 12 * dof_handler.fe_collection.max_dofs_per_quad();
155 break;
156 case 7:
157 max_couplings =
158 30 * dof_handler.fe_collection.max_dofs_per_vertex() +
159 45 * dof_handler.fe_collection.max_dofs_per_line() +
160 13 * dof_handler.fe_collection.max_dofs_per_quad();
161 break;
162 case 8:
163 max_couplings =
164 37 * dof_handler.fe_collection.max_dofs_per_vertex() +
165 56 * dof_handler.fe_collection.max_dofs_per_line() +
166 16 * dof_handler.fe_collection.max_dofs_per_quad();
167 break;
168
169 // the following numbers are not based on actual counting but by
170 // extrapolating the number sequences from the previous ones (for
171 // example, for n_dofs_per_vertex(), the sequence above is 19, 21,
172 // 28, 30, 37, and is continued as follows):
173 case 9:
174 max_couplings =
175 39 * dof_handler.fe_collection.max_dofs_per_vertex() +
176 59 * dof_handler.fe_collection.max_dofs_per_line() +
177 17 * dof_handler.fe_collection.max_dofs_per_quad();
178 break;
179 case 10:
180 max_couplings =
181 46 * dof_handler.fe_collection.max_dofs_per_vertex() +
182 70 * dof_handler.fe_collection.max_dofs_per_line() +
183 20 * dof_handler.fe_collection.max_dofs_per_quad();
184 break;
185 case 11:
186 max_couplings =
187 48 * dof_handler.fe_collection.max_dofs_per_vertex() +
188 73 * dof_handler.fe_collection.max_dofs_per_line() +
189 21 * dof_handler.fe_collection.max_dofs_per_quad();
190 break;
191 case 12:
192 max_couplings =
193 55 * dof_handler.fe_collection.max_dofs_per_vertex() +
194 84 * dof_handler.fe_collection.max_dofs_per_line() +
195 24 * dof_handler.fe_collection.max_dofs_per_quad();
196 break;
197 case 13:
198 max_couplings =
199 57 * dof_handler.fe_collection.max_dofs_per_vertex() +
200 87 * dof_handler.fe_collection.max_dofs_per_line() +
201 25 * dof_handler.fe_collection.max_dofs_per_quad();
202 break;
203 case 14:
204 max_couplings =
205 63 * dof_handler.fe_collection.max_dofs_per_vertex() +
206 98 * dof_handler.fe_collection.max_dofs_per_line() +
207 28 * dof_handler.fe_collection.max_dofs_per_quad();
208 break;
209 case 15:
210 max_couplings =
211 65 * dof_handler.fe_collection.max_dofs_per_vertex() +
212 103 * dof_handler.fe_collection.max_dofs_per_line() +
213 29 * dof_handler.fe_collection.max_dofs_per_quad();
214 break;
215 case 16:
216 max_couplings =
217 72 * dof_handler.fe_collection.max_dofs_per_vertex() +
218 114 * dof_handler.fe_collection.max_dofs_per_line() +
219 32 * dof_handler.fe_collection.max_dofs_per_quad();
220 break;
221
222 default:
223 Assert(false, ExcNotImplemented());
224 max_couplings = 0;
225 }
226 return std::min(max_couplings, dof_handler.n_dofs());
227 }
228
229 template <int spacedim>
230 static unsigned int
231 max_couplings_between_dofs(const DoFHandler<3, spacedim> &dof_handler)
232 {
233 // TODO:[?] Invent significantly better estimates than the ones in this
234 // function
235
236 // doing the same thing here is a rather complicated thing, compared
237 // to the 2d case, since it is hard to draw pictures with several
238 // refined hexahedra :-) so I presently only give a coarse
239 // estimate for the case that at most 8 hexes meet at each vertex
240 //
241 // can anyone give better estimate here?
242 const unsigned int max_adjacent_cells =
243 dof_handler.tria->max_adjacent_cells();
244
245 types::global_dof_index max_couplings;
246 if (max_adjacent_cells <= 8)
247 max_couplings =
248 7 * 7 * 7 * dof_handler.fe_collection.max_dofs_per_vertex() +
249 7 * 6 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_line() +
250 9 * 4 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_quad() +
251 27 * dof_handler.fe_collection.max_dofs_per_hex();
252 else
253 {
254 Assert(false, ExcNotImplemented());
255 max_couplings = 0;
256 }
257
258 return std::min(max_couplings, dof_handler.n_dofs());
259 }
260
265 template <int dim, int spacedim>
266 static void
267 reset_to_empty_objects(DoFHandler<dim, spacedim> &dof_handler)
268 {
269 dof_handler.object_dof_indices.clear();
270 dof_handler.object_dof_indices.resize(dof_handler.tria->n_levels());
271 dof_handler.object_dof_indices.shrink_to_fit();
272
273 dof_handler.object_dof_ptr.clear();
274 dof_handler.object_dof_ptr.resize(dof_handler.tria->n_levels());
275 dof_handler.object_dof_ptr.shrink_to_fit();
276
277 dof_handler.cell_dof_cache_indices.clear();
278 dof_handler.cell_dof_cache_indices.resize(dof_handler.tria->n_levels());
279 dof_handler.cell_dof_cache_indices.shrink_to_fit();
280
281 dof_handler.cell_dof_cache_ptr.clear();
282 dof_handler.cell_dof_cache_ptr.resize(dof_handler.tria->n_levels());
283 dof_handler.cell_dof_cache_ptr.shrink_to_fit();
284 }
285
289 template <int dim, int spacedim>
290 static void
291 reserve_cells(DoFHandler<dim, spacedim> &dof_handler,
292 const unsigned int n_inner_dofs_per_cell)
293 {
294 for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
295 {
296 // 1) object_dof_indices
297 dof_handler.object_dof_ptr[i][dim].assign(
298 dof_handler.tria->n_raw_cells(i) + 1, 0);
299
300 for (const auto &cell :
301 dof_handler.tria->cell_iterators_on_level(i))
302 if (cell->is_active() && !cell->is_artificial())
303 dof_handler.object_dof_ptr[i][dim][cell->index() + 1] =
304 n_inner_dofs_per_cell;
305
306 for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
307 dof_handler.object_dof_ptr[i][dim][j + 1] +=
308 dof_handler.object_dof_ptr[i][dim][j];
309
310 dof_handler.object_dof_indices[i][dim].resize(
311 dof_handler.object_dof_ptr[i][dim].back(),
312 numbers::invalid_dof_index);
313
314 // 2) cell_dof_cache_indices
315 dof_handler.cell_dof_cache_ptr[i].assign(
316 dof_handler.tria->n_raw_cells(i) + 1, 0);
317
318 for (const auto &cell :
319 dof_handler.tria->cell_iterators_on_level(i))
320 if (cell->is_active() && !cell->is_artificial())
321 dof_handler.cell_dof_cache_ptr[i][cell->index() + 1] =
322 dof_handler.get_fe().n_dofs_per_cell();
323
324 for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
325 dof_handler.cell_dof_cache_ptr[i][j + 1] +=
326 dof_handler.cell_dof_cache_ptr[i][j];
327
328 dof_handler.cell_dof_cache_indices[i].resize(
329 dof_handler.cell_dof_cache_ptr[i].back(),
330 numbers::invalid_dof_index);
331 }
332 }
333
338 template <int dim, int spacedim, typename T>
339 static void
340 reserve_subentities(DoFHandler<dim, spacedim> &dof_handler,
341 const unsigned int structdim,
342 const unsigned int n_raw_entities,
343 const T & cell_process)
344 {
345 if (dof_handler.tria->n_cells() == 0)
346 return;
347
348 dof_handler.object_dof_ptr[0][structdim].assign(n_raw_entities + 1, -1);
349 // determine for each entity the number of dofs
350 for (const auto &cell : dof_handler.tria->cell_iterators())
351 if (cell->is_active() && !cell->is_artificial())
352 cell_process(
353 cell,
354 [&](const unsigned int n_dofs_per_entity,
355 const unsigned int index) {
356 auto &n_dofs_per_entity_target =
357 dof_handler.object_dof_ptr[0][structdim][index + 1];
358
359 // make sure that either the entity has not been visited or
360 // the entity has the same number of dofs assigned
361 Assert((n_dofs_per_entity_target ==
362 static_cast<
363 typename DoFHandler<dim, spacedim>::offset_type>(
364 -1) ||
365 n_dofs_per_entity_target == n_dofs_per_entity),
366 ExcNotImplemented());
367
368 n_dofs_per_entity_target = n_dofs_per_entity;
369 });
370
371 // convert the absolute numbers to CRS
372 dof_handler.object_dof_ptr[0][structdim][0] = 0;
373 for (unsigned int i = 1; i < n_raw_entities + 1; ++i)
374 {
375 if (dof_handler.object_dof_ptr[0][structdim][i] ==
376 static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
377 -1))
378 dof_handler.object_dof_ptr[0][structdim][i] =
379 dof_handler.object_dof_ptr[0][structdim][i - 1];
380 else
381 dof_handler.object_dof_ptr[0][structdim][i] +=
382 dof_handler.object_dof_ptr[0][structdim][i - 1];
383 }
384
385 // allocate memory for indices
386 dof_handler.object_dof_indices[0][structdim].resize(
387 dof_handler.object_dof_ptr[0][structdim].back(),
388 numbers::invalid_dof_index);
389 }
390
397 template <int dim, int spacedim>
398 static void
399 reserve_space(DoFHandler<dim, spacedim> &dof_handler)
400 {
401 reset_to_empty_objects(dof_handler);
402
403 const auto &fe = dof_handler.get_fe();
404
405 // cell
406 reserve_cells(dof_handler,
407 dim == 1 ? fe.n_dofs_per_line() :
408 (dim == 2 ? fe.n_dofs_per_quad(0) :
409 fe.n_dofs_per_hex()));
410
411 // vertices
412 reserve_subentities(dof_handler,
413 0,
414 dof_handler.tria->n_vertices(),
415 [&](const auto &cell, const auto &process) {
416 for (const auto vertex_index :
417 cell->vertex_indices())
418 process(fe.n_dofs_per_vertex(),
419 cell->vertex_index(vertex_index));
420 });
421
422 // lines
423 if (dim == 2 || dim == 3)
424 reserve_subentities(dof_handler,
425 1,
426 dof_handler.tria->n_raw_lines(),
427 [&](const auto &cell, const auto &process) {
428 for (const auto line_index :
429 cell->line_indices())
430 process(fe.n_dofs_per_line(),
431 cell->line(line_index)->index());
432 });
433
434 // quads
435 if (dim == 3)
436 reserve_subentities(dof_handler,
437 2,
438 dof_handler.tria->n_raw_quads(),
439 [&](const auto &cell, const auto &process) {
440 for (const auto face_index :
441 cell->face_indices())
442 process(fe.n_dofs_per_quad(face_index),
443 cell->face(face_index)->index());
444 });
445 }
446
447 template <int spacedim>
448 static void
449 reserve_space_mg(DoFHandler<1, spacedim> &dof_handler)
450 {
451 Assert(dof_handler.get_triangulation().n_levels() > 0,
452 ExcMessage("Invalid triangulation"));
453 dof_handler.clear_mg_space();
454
455 const ::Triangulation<1, spacedim> &tria =
456 dof_handler.get_triangulation();
457 const unsigned int dofs_per_line =
458 dof_handler.get_fe().n_dofs_per_line();
459 const unsigned int n_levels = tria.n_levels();
460
461 for (unsigned int i = 0; i < n_levels; ++i)
462 {
463 dof_handler.mg_levels.emplace_back(
464 new internal::DoFHandlerImplementation::DoFLevel<1>);
465 dof_handler.mg_levels.back()->dof_object.dofs =
466 std::vector<types::global_dof_index>(tria.n_raw_lines(i) *
467 dofs_per_line,
468 numbers::invalid_dof_index);
469 }
470
471 const unsigned int n_vertices = tria.n_vertices();
472
473 dof_handler.mg_vertex_dofs.resize(n_vertices);
474
475 std::vector<unsigned int> max_level(n_vertices, 0);
476 std::vector<unsigned int> min_level(n_vertices, n_levels);
477
478 for (typename ::Triangulation<1, spacedim>::cell_iterator cell =
479 tria.begin();
480 cell != tria.end();
481 ++cell)
482 {
483 const unsigned int level = cell->level();
484
485 for (const auto vertex : cell->vertex_indices())
486 {
487 const unsigned int vertex_index = cell->vertex_index(vertex);
488
489 if (min_level[vertex_index] > level)
490 min_level[vertex_index] = level;
491
492 if (max_level[vertex_index] < level)
493 max_level[vertex_index] = level;
494 }
495 }
496
497 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
498 if (tria.vertex_used(vertex))
499 {
500 Assert(min_level[vertex] < n_levels, ExcInternalError());
501 Assert(max_level[vertex] >= min_level[vertex],
502 ExcInternalError());
503 dof_handler.mg_vertex_dofs[vertex].init(
504 min_level[vertex],
505 max_level[vertex],
506 dof_handler.get_fe().n_dofs_per_vertex());
507 }
508
509 else
510 {
511 Assert(min_level[vertex] == n_levels, ExcInternalError());
512 Assert(max_level[vertex] == 0, ExcInternalError());
513 dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
514 }
515 }
516
517 template <int spacedim>
518 static void
519 reserve_space_mg(DoFHandler<2, spacedim> &dof_handler)
520 {
521 Assert(dof_handler.get_triangulation().n_levels() > 0,
522 ExcMessage("Invalid triangulation"));
523 dof_handler.clear_mg_space();
524
525 const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe();
526 const ::Triangulation<2, spacedim> &tria =
527 dof_handler.get_triangulation();
528 const unsigned int n_levels = tria.n_levels();
529
530 for (unsigned int i = 0; i < n_levels; ++i)
531 {
532 dof_handler.mg_levels.emplace_back(
533 std::make_unique<
534 internal::DoFHandlerImplementation::DoFLevel<2>>());
535 dof_handler.mg_levels.back()->dof_object.dofs =
536 std::vector<types::global_dof_index>(
537 tria.n_raw_quads(i) *
538 fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/),
539 numbers::invalid_dof_index);
540 }
541
542 dof_handler.mg_faces =
543 std::make_unique<internal::DoFHandlerImplementation::DoFFaces<2>>();
544 dof_handler.mg_faces->lines.dofs =
545 std::vector<types::global_dof_index>(tria.n_raw_lines() *
546 fe.n_dofs_per_line(),
547 numbers::invalid_dof_index);
548
549 const unsigned int n_vertices = tria.n_vertices();
550
551 dof_handler.mg_vertex_dofs.resize(n_vertices);
552
553 std::vector<unsigned int> max_level(n_vertices, 0);
554 std::vector<unsigned int> min_level(n_vertices, n_levels);
555
556 for (typename ::Triangulation<2, spacedim>::cell_iterator cell =
557 tria.begin();
558 cell != tria.end();
559 ++cell)
560 {
561 const unsigned int level = cell->level();
562
563 for (const auto vertex : cell->vertex_indices())
564 {
565 const unsigned int vertex_index = cell->vertex_index(vertex);
566
567 if (min_level[vertex_index] > level)
568 min_level[vertex_index] = level;
569
570 if (max_level[vertex_index] < level)
571 max_level[vertex_index] = level;
572 }
573 }
574
575 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
576 if (tria.vertex_used(vertex))
577 {
578 Assert(min_level[vertex] < n_levels, ExcInternalError());
579 Assert(max_level[vertex] >= min_level[vertex],
580 ExcInternalError());
581 dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
582 max_level[vertex],
583 fe.n_dofs_per_vertex());
584 }
585
586 else
587 {
588 Assert(min_level[vertex] == n_levels, ExcInternalError());
589 Assert(max_level[vertex] == 0, ExcInternalError());
590 dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
591 }
592 }
593
594 template <int spacedim>
595 static void
596 reserve_space_mg(DoFHandler<3, spacedim> &dof_handler)
597 {
598 Assert(dof_handler.get_triangulation().n_levels() > 0,
599 ExcMessage("Invalid triangulation"));
600 dof_handler.clear_mg_space();
601
602 const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe();
603 const ::Triangulation<3, spacedim> &tria =
604 dof_handler.get_triangulation();
605 const unsigned int n_levels = tria.n_levels();
606
607 for (unsigned int i = 0; i < n_levels; ++i)
608 {
609 dof_handler.mg_levels.emplace_back(
610 std::make_unique<
611 internal::DoFHandlerImplementation::DoFLevel<3>>());
612 dof_handler.mg_levels.back()->dof_object.dofs =
613 std::vector<types::global_dof_index>(tria.n_raw_hexs(i) *
614 fe.n_dofs_per_hex(),
615 numbers::invalid_dof_index);
616 }
617
618 dof_handler.mg_faces =
619 std::make_unique<internal::DoFHandlerImplementation::DoFFaces<3>>();
620 dof_handler.mg_faces->lines.dofs =
621 std::vector<types::global_dof_index>(tria.n_raw_lines() *
622 fe.n_dofs_per_line(),
623 numbers::invalid_dof_index);
624
625 // TODO: the implementation makes the assumption that all faces have the
626 // same number of dofs
627 AssertDimension(fe.n_unique_faces(), 1);
628 dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
629 tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
630 numbers::invalid_dof_index);
631
632 const unsigned int n_vertices = tria.n_vertices();
633
634 dof_handler.mg_vertex_dofs.resize(n_vertices);
635
636 std::vector<unsigned int> max_level(n_vertices, 0);
637 std::vector<unsigned int> min_level(n_vertices, n_levels);
638
639 for (typename ::Triangulation<3, spacedim>::cell_iterator cell =
640 tria.begin();
641 cell != tria.end();
642 ++cell)
643 {
644 const unsigned int level = cell->level();
645
646 for (const auto vertex : cell->vertex_indices())
647 {
648 const unsigned int vertex_index = cell->vertex_index(vertex);
649
650 if (min_level[vertex_index] > level)
651 min_level[vertex_index] = level;
652
653 if (max_level[vertex_index] < level)
654 max_level[vertex_index] = level;
655 }
656 }
657
658 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
659 if (tria.vertex_used(vertex))
660 {
661 Assert(min_level[vertex] < n_levels, ExcInternalError());
662 Assert(max_level[vertex] >= min_level[vertex],
663 ExcInternalError());
664 dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
665 max_level[vertex],
666 fe.n_dofs_per_vertex());
667 }
668
669 else
670 {
671 Assert(min_level[vertex] == n_levels, ExcInternalError());
672 Assert(max_level[vertex] == 0, ExcInternalError());
673 dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
674 }
675 }
676 };
677 } // namespace DoFHandlerImplementation
678
679
680
681 namespace hp
682 {
683 namespace DoFHandlerImplementation
684 {
689 struct Implementation
690 {
696 template <int dim, int spacedim>
697 static void
698 ensure_absence_of_future_fe_indices(
699 DoFHandler<dim, spacedim> &dof_handler)
700 {
701 (void)dof_handler;
702 for (const auto &cell : dof_handler.active_cell_iterators())
703 if (cell->is_locally_owned())
704 Assert(
705 !cell->future_fe_index_set(),
706 ExcMessage(
707 "There shouldn't be any cells flagged for p-adaptation when partitioning."));
708 }
709
710
711
716 template <int dim, int spacedim>
717 static void
718 reserve_space_vertices(DoFHandler<dim, spacedim> &dof_handler)
719 {
720 // The final step in all of the reserve_space() functions is to set
721 // up vertex dof information. since vertices are sequentially
722 // numbered, what we do first is to set up an array in which
723 // we record whether a vertex is associated with any of the
724 // given fe's, by setting a bit. in a later step, we then
725 // actually allocate memory for the required dofs
726 //
727 // in the following, we only need to consider vertices that are
728 // adjacent to either a locally owned or a ghost cell; we never
729 // store anything on vertices that are only surrounded by
730 // artificial cells. so figure out that subset of vertices
731 // first
732 std::vector<bool> locally_used_vertices(
733 dof_handler.tria->n_vertices(), false);
734 for (const auto &cell : dof_handler.active_cell_iterators())
735 if (!cell->is_artificial())
736 for (const auto v : cell->vertex_indices())
737 locally_used_vertices[cell->vertex_index(v)] = true;
738
739 std::vector<std::vector<bool>> vertex_fe_association(
740 dof_handler.fe_collection.size(),
741 std::vector<bool>(dof_handler.tria->n_vertices(), false));
742
743 for (const auto &cell : dof_handler.active_cell_iterators())
744 if (!cell->is_artificial())
745 for (const auto v : cell->vertex_indices())
746 vertex_fe_association[cell->active_fe_index()]
747 [cell->vertex_index(v)] = true;
748
749 // in debug mode, make sure that each vertex is associated
750 // with at least one FE (note that except for unused
751 // vertices, all vertices are actually active). this is of
752 // course only true for vertices that are part of either
753 // ghost or locally owned cells
754#ifdef DEBUG
755 for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
756 if (locally_used_vertices[v] == true)
757 if (dof_handler.tria->vertex_used(v) == true)
758 {
759 unsigned int fe = 0;
760 for (; fe < dof_handler.fe_collection.size(); ++fe)
761 if (vertex_fe_association[fe][v] == true)
762 break;
763 Assert(fe != dof_handler.fe_collection.size(),
764 ExcInternalError());
765 }
766#endif
767
768 const unsigned int d = 0;
769 const unsigned int l = 0;
770
771 dof_handler.hp_object_fe_ptr[d].clear();
772 dof_handler.hp_object_fe_indices[d].clear();
773 dof_handler.object_dof_ptr[l][d].clear();
774 dof_handler.object_dof_indices[l][d].clear();
775
776 dof_handler.hp_object_fe_ptr[d].reserve(
777 dof_handler.tria->n_vertices() + 1);
778
779 unsigned int vertex_slots_needed = 0;
780 unsigned int fe_slots_needed = 0;
781
782 for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
783 {
784 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
785
786 if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
787 {
788 for (unsigned int fe = 0;
789 fe < dof_handler.fe_collection.size();
790 ++fe)
791 if (vertex_fe_association[fe][v] == true)
792 {
793 fe_slots_needed++;
794 vertex_slots_needed +=
795 dof_handler.get_fe(fe).n_dofs_per_vertex();
796 }
797 }
798 }
799
800 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
801
802 dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
803 dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
804
805 dof_handler.object_dof_indices[l][d].reserve(vertex_slots_needed);
806
807 for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
808 if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
809 {
810 for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
811 ++fe)
812 if (vertex_fe_association[fe][v] == true)
813 {
814 dof_handler.hp_object_fe_indices[d].push_back(fe);
815 dof_handler.object_dof_ptr[l][d].push_back(
816 dof_handler.object_dof_indices[l][d].size());
817
818 for (unsigned int i = 0;
819 i < dof_handler.get_fe(fe).n_dofs_per_vertex();
820 i++)
821 dof_handler.object_dof_indices[l][d].push_back(
822 numbers::invalid_dof_index);
823 }
824 }
825
826
827 dof_handler.object_dof_ptr[l][d].push_back(
828 dof_handler.object_dof_indices[l][d].size());
829
830 AssertDimension(vertex_slots_needed,
831 dof_handler.object_dof_indices[l][d].size());
832 AssertDimension(fe_slots_needed,
833 dof_handler.hp_object_fe_indices[d].size());
834 AssertDimension(fe_slots_needed + 1,
835 dof_handler.object_dof_ptr[l][d].size());
836 AssertDimension(dof_handler.tria->n_vertices() + 1,
837 dof_handler.hp_object_fe_ptr[d].size());
838
839 dof_handler.object_dof_indices[l][d].assign(
840 vertex_slots_needed, numbers::invalid_dof_index);
841 }
842
843
844
849 template <int dim, int spacedim>
850 static void
851 reserve_space_cells(DoFHandler<dim, spacedim> &dof_handler)
852 {
853 (void)dof_handler;
854 // count how much space we need on each level for the cell
855 // dofs and set the dof_*_offsets data. initially set the
856 // latter to an invalid index, and only later set it to
857 // something reasonable for active dof_handler.cells
858 //
859 // note that for dof_handler.cells, the situation is simpler
860 // than for other (lower dimensional) objects since exactly
861 // one finite element is used for it
862 for (unsigned int level = 0; level < dof_handler.tria->n_levels();
863 ++level)
864 {
865 dof_handler.object_dof_ptr[level][dim] =
866 std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
867 dof_handler.tria->n_raw_cells(level),
868 static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
869 -1));
870 dof_handler.cell_dof_cache_ptr[level] =
871 std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
872 dof_handler.tria->n_raw_cells(level),
873 static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
874 -1));
875
876 types::global_dof_index next_free_dof = 0;
877 types::global_dof_index cache_size = 0;
878
879 for (auto cell :
880 dof_handler.active_cell_iterators_on_level(level))
881 if (cell->is_active() && !cell->is_artificial())
882 {
883 dof_handler.object_dof_ptr[level][dim][cell->index()] =
884 next_free_dof;
885 next_free_dof +=
886 cell->get_fe().template n_dofs_per_object<dim>();
887
888 dof_handler.cell_dof_cache_ptr[level][cell->index()] =
889 cache_size;
890 cache_size += cell->get_fe().n_dofs_per_cell();
891 }
892
893 dof_handler.object_dof_indices[level][dim] =
894 std::vector<types::global_dof_index>(
895 next_free_dof, numbers::invalid_dof_index);
896 dof_handler.cell_dof_cache_indices[level] =
897 std::vector<types::global_dof_index>(
898 cache_size, numbers::invalid_dof_index);
899 }
900 }
901
902
903
908 template <int dim, int spacedim>
909 static void
910 reserve_space_faces(DoFHandler<dim, spacedim> &dof_handler)
911 {
912 // FACE DOFS
913 //
914 // Count face dofs, then allocate as much space
915 // as we need and prime the linked list for faces (see the
916 // description in hp::DoFLevel) with the indices we will
917 // need. Note that our task is more complicated than for the
918 // cell case above since two adjacent cells may have different
919 // active FE indices, in which case we need to allocate
920 // *two* sets of face dofs for the same face. But they don't
921 // *have* to be different, and so we need to prepare for this
922 // as well.
923 //
924 // The way we do things is that we loop over all active
925 // cells (these are the only ones that have DoFs
926 // anyway) and all their faces. We note in the
927 // user flags whether we have previously visited a face and
928 // if so skip it (consequently, we have to save and later
929 // restore the face flags)
930 {
931 std::vector<bool> saved_face_user_flags;
932 switch (dim)
933 {
934 case 2:
935 {
936 const_cast<::Triangulation<dim, spacedim> &>(
937 *dof_handler.tria)
938 .save_user_flags_line(saved_face_user_flags);
939 const_cast<::Triangulation<dim, spacedim> &>(
940 *dof_handler.tria)
941 .clear_user_flags_line();
942
943 break;
944 }
945
946 case 3:
947 {
948 const_cast<::Triangulation<dim, spacedim> &>(
949 *dof_handler.tria)
950 .save_user_flags_quad(saved_face_user_flags);
951 const_cast<::Triangulation<dim, spacedim> &>(
952 *dof_handler.tria)
953 .clear_user_flags_quad();
954
955 break;
956 }
957
958 default:
959 Assert(false, ExcNotImplemented());
960 }
961
962 const unsigned int d = dim - 1;
963 const unsigned int l = 0;
964
965 dof_handler.hp_object_fe_ptr[d].clear();
966 dof_handler.hp_object_fe_indices[d].clear();
967 dof_handler.object_dof_ptr[l][d].clear();
968 dof_handler.object_dof_indices[l][d].clear();
969
970 dof_handler.hp_object_fe_ptr[d].resize(
971 dof_handler.tria->n_raw_faces() + 1);
972
973 // An array to hold how many slots (see the hp::DoFLevel
974 // class) we will have to store on each level
975 unsigned int n_face_slots = 0;
976
977 for (const auto &cell : dof_handler.active_cell_iterators())
978 if (!cell->is_artificial())
979 for (const auto face : cell->face_indices())
980 if (cell->face(face)->user_flag_set() == false)
981 {
982 unsigned int fe_slots_needed = 0;
983
984 if (cell->at_boundary(face) ||
985 cell->face(face)->has_children() ||
986 cell->neighbor_is_coarser(face) ||
987 (!cell->at_boundary(face) &&
988 cell->neighbor(face)->is_artificial()) ||
989 (!cell->at_boundary(face) &&
990 !cell->neighbor(face)->is_artificial() &&
991 (cell->active_fe_index() ==
992 cell->neighbor(face)->active_fe_index())))
993 {
994 fe_slots_needed = 1;
995 n_face_slots +=
996 dof_handler.get_fe(cell->active_fe_index())
997 .template n_dofs_per_object<dim - 1>(face);
998 }
999 else
1000 {
1001 fe_slots_needed = 2;
1002 n_face_slots +=
1003 dof_handler.get_fe(cell->active_fe_index())
1004 .template n_dofs_per_object<dim - 1>(face) +
1005 dof_handler
1006 .get_fe(cell->neighbor(face)->active_fe_index())
1007 .template n_dofs_per_object<dim - 1>(
1008 cell->neighbor_face_no(face));
1009 }
1010
1011 // mark this face as visited
1012 cell->face(face)->set_user_flag();
1013
1014 dof_handler
1015 .hp_object_fe_ptr[d][cell->face(face)->index() + 1] =
1016 fe_slots_needed;
1017 }
1018
1019 for (unsigned int i = 1; i < dof_handler.hp_object_fe_ptr[d].size();
1020 i++)
1021 dof_handler.hp_object_fe_ptr[d][i] +=
1022 dof_handler.hp_object_fe_ptr[d][i - 1];
1023
1024
1025 dof_handler.hp_object_fe_indices[d].resize(
1026 dof_handler.hp_object_fe_ptr[d].back());
1027 dof_handler.object_dof_ptr[l][d].resize(
1028 dof_handler.hp_object_fe_ptr[d].back() + 1);
1029
1030 dof_handler.object_dof_indices[l][d].reserve(n_face_slots);
1031
1032
1033 // With the memory now allocated, loop over the
1034 // dof_handler cells again and prime the _offset values as
1035 // well as the fe_index fields
1036 switch (dim)
1037 {
1038 case 2:
1039 {
1040 const_cast<::Triangulation<dim, spacedim> &>(
1041 *dof_handler.tria)
1042 .clear_user_flags_line();
1043
1044 break;
1045 }
1046
1047 case 3:
1048 {
1049 const_cast<::Triangulation<dim, spacedim> &>(
1050 *dof_handler.tria)
1051 .clear_user_flags_quad();
1052
1053 break;
1054 }
1055
1056 default:
1057 Assert(false, ExcNotImplemented());
1058 }
1059
1060 for (const auto &cell : dof_handler.active_cell_iterators())
1061 if (!cell->is_artificial())
1062 for (const auto face : cell->face_indices())
1063 if (!cell->face(face)->user_flag_set())
1064 {
1065 // Same decision tree as before
1066 if (cell->at_boundary(face) ||
1067 cell->face(face)->has_children() ||
1068 cell->neighbor_is_coarser(face) ||
1069 (!cell->at_boundary(face) &&
1070 cell->neighbor(face)->is_artificial()) ||
1071 (!cell->at_boundary(face) &&
1072 !cell->neighbor(face)->is_artificial() &&
1073 (cell->active_fe_index() ==
1074 cell->neighbor(face)->active_fe_index())))
1075 {
1076 const unsigned int fe = cell->active_fe_index();
1077 const unsigned int n_dofs =
1078 dof_handler.get_fe(fe)
1079 .template n_dofs_per_object<dim - 1>(face);
1080 const unsigned int offset =
1081 dof_handler
1082 .hp_object_fe_ptr[d][cell->face(face)->index()];
1083
1084 dof_handler.hp_object_fe_indices[d][offset] = fe;
1085 dof_handler.object_dof_ptr[l][d][offset + 1] = n_dofs;
1086
1087 for (unsigned int i = 0; i < n_dofs; ++i)
1088 dof_handler.object_dof_indices[l][d].push_back(
1089 numbers::invalid_dof_index);
1090 }
1091 else
1092 {
1093 unsigned int fe_1 = cell->active_fe_index();
1094 unsigned int face_no_1 = face;
1095 unsigned int fe_2 =
1096 cell->neighbor(face)->active_fe_index();
1097 unsigned int face_no_2 = cell->neighbor_face_no(face);
1098
1099 if (fe_2 < fe_1)
1100 {
1101 std::swap(fe_1, fe_2);
1102 std::swap(face_no_1, face_no_2);
1103 }
1104
1105 const unsigned int n_dofs_1 =
1106 dof_handler.get_fe(fe_1)
1107 .template n_dofs_per_object<dim - 1>(face_no_1);
1108
1109 const unsigned int n_dofs_2 =
1110 dof_handler.get_fe(fe_2)
1111 .template n_dofs_per_object<dim - 1>(face_no_2);
1112
1113 const unsigned int offset =
1114 dof_handler
1115 .hp_object_fe_ptr[d][cell->face(face)->index()];
1116
1117 dof_handler.hp_object_fe_indices[d].push_back(
1118 cell->active_fe_index());
1119 dof_handler.object_dof_ptr[l][d].push_back(
1120 dof_handler.object_dof_indices[l][d].size());
1121
1122 dof_handler.hp_object_fe_indices[d][offset + 0] =
1123 fe_1;
1124 dof_handler.hp_object_fe_indices[d][offset + 1] =
1125 fe_2;
1126 dof_handler.object_dof_ptr[l][d][offset + 1] =
1127 n_dofs_1;
1128 dof_handler.object_dof_ptr[l][d][offset + 2] =
1129 n_dofs_2;
1130
1131
1132 for (unsigned int i = 0; i < n_dofs_1 + n_dofs_2; ++i)
1133 dof_handler.object_dof_indices[l][d].push_back(
1134 numbers::invalid_dof_index);
1135 }
1136
1137 // mark this face as visited
1138 cell->face(face)->set_user_flag();
1139 }
1140
1141 for (unsigned int i = 1;
1142 i < dof_handler.object_dof_ptr[l][d].size();
1143 i++)
1144 dof_handler.object_dof_ptr[l][d][i] +=
1145 dof_handler.object_dof_ptr[l][d][i - 1];
1146
1147 // at the end, restore the user flags for the faces
1148 switch (dim)
1149 {
1150 case 2:
1151 {
1152 const_cast<::Triangulation<dim, spacedim> &>(
1153 *dof_handler.tria)
1154 .load_user_flags_line(saved_face_user_flags);
1155
1156 break;
1157 }
1158
1159 case 3:
1160 {
1161 const_cast<::Triangulation<dim, spacedim> &>(
1162 *dof_handler.tria)
1163 .load_user_flags_quad(saved_face_user_flags);
1164
1165 break;
1166 }
1167
1168 default:
1169 Assert(false, ExcNotImplemented());
1170 }
1171 }
1172 }
1173
1174
1175
1182 template <int spacedim>
1183 static void
1184 reserve_space(::DoFHandler<1, spacedim> &dof_handler)
1185 {
1186 Assert(dof_handler.fe_collection.size() > 0,
1187 (typename ::DoFHandler<1, spacedim>::ExcNoFESelected()));
1188 Assert(dof_handler.tria->n_levels() > 0,
1189 ExcMessage("The current Triangulation must not be empty."));
1190 Assert(dof_handler.tria->n_levels() ==
1191 dof_handler.hp_cell_future_fe_indices.size(),
1192 ExcInternalError());
1193
1194 internal::DoFHandlerImplementation::Implementation::
1195 reset_to_empty_objects(dof_handler);
1196
1197 Threads::TaskGroup<> tasks;
1198 tasks +=
1199 Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
1200 tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
1201 dof_handler);
1202 tasks.join_all();
1203 }
1204
1205
1206
1207 template <int spacedim>
1208 static void
1209 reserve_space(::DoFHandler<2, spacedim> &dof_handler)
1210 {
1211 Assert(dof_handler.fe_collection.size() > 0,
1212 (typename ::DoFHandler<1, spacedim>::ExcNoFESelected()));
1213 Assert(dof_handler.tria->n_levels() > 0,
1214 ExcMessage("The current Triangulation must not be empty."));
1215 Assert(dof_handler.tria->n_levels() ==
1216 dof_handler.hp_cell_future_fe_indices.size(),
1217 ExcInternalError());
1218
1219 internal::DoFHandlerImplementation::Implementation::
1220 reset_to_empty_objects(dof_handler);
1221
1222 Threads::TaskGroup<> tasks;
1223 tasks +=
1224 Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
1225 tasks +=
1226 Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
1227 tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
1228 dof_handler);
1229 tasks.join_all();
1230 }
1231
1232
1233
1234 template <int spacedim>
1235 static void
1236 reserve_space(::DoFHandler<3, spacedim> &dof_handler)
1237 {
1238 Assert(dof_handler.fe_collection.size() > 0,
1239 (typename ::DoFHandler<1, spacedim>::ExcNoFESelected()));
1240 Assert(dof_handler.tria->n_levels() > 0,
1241 ExcMessage("The current Triangulation must not be empty."));
1242 Assert(dof_handler.tria->n_levels() ==
1243 dof_handler.hp_cell_future_fe_indices.size(),
1244 ExcInternalError());
1245
1246 internal::DoFHandlerImplementation::Implementation::
1247 reset_to_empty_objects(dof_handler);
1248
1249 Threads::TaskGroup<> tasks;
1250 tasks +=
1251 Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
1252 tasks +=
1253 Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
1254 tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
1255 dof_handler);
1256
1257 // While the tasks above are running, we can turn to line dofs
1258
1259 // the situation here is pretty much like with vertices:
1260 // there can be an arbitrary number of finite elements
1261 // associated with each line.
1262 //
1263 // the algorithm we use is somewhat similar to what we do in
1264 // reserve_space_vertices()
1265 {
1266 // what we do first is to set up an array in which we
1267 // record whether a line is associated with any of the
1268 // given fe's, by setting a bit. in a later step, we
1269 // then actually allocate memory for the required dofs
1270 std::vector<std::vector<bool>> line_fe_association(
1271 dof_handler.fe_collection.size(),
1272 std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
1273
1274 for (const auto &cell : dof_handler.active_cell_iterators())
1275 if (!cell->is_artificial())
1276 for (const auto l : cell->line_indices())
1277 line_fe_association[cell->active_fe_index()]
1278 [cell->line_index(l)] = true;
1279
1280 // first check which of the lines is used at all,
1281 // i.e. is associated with a finite element. we do this
1282 // since not all lines may actually be used, in which
1283 // case we do not have to allocate any memory at all
1284 std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
1285 false);
1286 for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1287 ++line)
1288 for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1289 ++fe)
1290 if (line_fe_association[fe][line] == true)
1291 {
1292 line_is_used[line] = true;
1293 break;
1294 }
1295
1296
1297
1298 const unsigned int d = 1;
1299 const unsigned int l = 0;
1300
1301 dof_handler.hp_object_fe_ptr[d].clear();
1302 dof_handler.hp_object_fe_indices[d].clear();
1303 dof_handler.object_dof_ptr[l][d].clear();
1304 dof_handler.object_dof_indices[l][d].clear();
1305
1306 dof_handler.hp_object_fe_ptr[d].reserve(
1307 dof_handler.tria->n_raw_lines() + 1);
1308
1309 unsigned int line_slots_needed = 0;
1310 unsigned int fe_slots_needed = 0;
1311
1312 for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1313 ++line)
1314 {
1315 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1316
1317 if (line_is_used[line] == true)
1318 {
1319 for (unsigned int fe = 0;
1320 fe < dof_handler.fe_collection.size();
1321 ++fe)
1322 if (line_fe_association[fe][line] == true)
1323 {
1324 fe_slots_needed++;
1325 line_slots_needed +=
1326 dof_handler.get_fe(fe).n_dofs_per_line();
1327 }
1328 }
1329 }
1330
1331 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1332
1333 // make sure that all entries have been set
1334 AssertDimension(dof_handler.hp_object_fe_ptr[d].size(),
1335 dof_handler.tria->n_raw_lines() + 1);
1336
1337 dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1338 dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1339
1340 dof_handler.object_dof_indices[l][d].reserve(line_slots_needed);
1341
1342 for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1343 ++line)
1344 if (line_is_used[line] == true)
1345 {
1346 for (unsigned int fe = 0;
1347 fe < dof_handler.fe_collection.size();
1348 ++fe)
1349 if (line_fe_association[fe][line] == true)
1350 {
1351 dof_handler.hp_object_fe_indices[d].push_back(fe);
1352 dof_handler.object_dof_ptr[l][d].push_back(
1353 dof_handler.object_dof_indices[l][d].size());
1354
1355 for (unsigned int i = 0;
1356 i < dof_handler.get_fe(fe).n_dofs_per_line();
1357 i++)
1358 dof_handler.object_dof_indices[l][d].push_back(
1359 numbers::invalid_dof_index);
1360 }
1361 }
1362
1363 dof_handler.object_dof_ptr[l][d].push_back(
1364 dof_handler.object_dof_indices[l][d].size());
1365
1366 // make sure that all entries have been set
1367 AssertDimension(dof_handler.hp_object_fe_indices[d].size(),
1368 fe_slots_needed);
1369 AssertDimension(dof_handler.object_dof_ptr[l][d].size(),
1370 fe_slots_needed + 1);
1371 AssertDimension(dof_handler.object_dof_indices[l][d].size(),
1372 line_slots_needed);
1373 }
1374
1375 // Ensure that everything is done at this point.
1376 tasks.join_all();
1377 }
1378
1379
1380
1392 template <int dim, int spacedim>
1393 static void
1394 communicate_active_fe_indices(DoFHandler<dim, spacedim> &dof_handler)
1395 {
1396 Assert(
1397 dof_handler.hp_capability_enabled == true,
1398 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
1399
1400 using active_fe_index_type =
1401 typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
1402
1403 if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1404 dynamic_cast<
1405 const ::parallel::shared::Triangulation<dim, spacedim>
1406 *>(&dof_handler.get_triangulation()))
1407 {
1408 // we have a shared triangulation. in this case, every processor
1409 // knows about all cells, but every processor only has knowledge
1410 // about the active FE index on the cells it owns.
1411 //
1412 // we can create a complete set of active FE indices by letting
1413 // every processor create a vector of indices for all cells,
1414 // filling only those on the cells it owns and setting the indices
1415 // on the other cells to zero. then we add all of these vectors
1416 // up, and because every vector entry has exactly one processor
1417 // that owns it, the sum is correct
1418 std::vector<active_fe_index_type> active_fe_indices(
1419 tr->n_active_cells(), 0u);
1420 for (const auto &cell : dof_handler.active_cell_iterators())
1421 if (cell->is_locally_owned())
1422 active_fe_indices[cell->active_cell_index()] =
1423 cell->active_fe_index();
1424
1425 Utilities::MPI::sum(active_fe_indices,
1426 tr->get_communicator(),
1427 active_fe_indices);
1428
1429 // now go back and fill the active FE index on all other
1430 // cells. we would like to call cell->set_active_fe_index(),
1431 // but that function does not allow setting these indices on
1432 // non-locally_owned cells. so we have to work around the
1433 // issue a little bit by accessing the underlying data
1434 // structures directly
1435 for (const auto &cell : dof_handler.active_cell_iterators())
1436 if (!cell->is_locally_owned())
1437 dof_handler
1438 .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1439 active_fe_indices[cell->active_cell_index()];
1440 }
1441 else if (const ::parallel::
1442 DistributedTriangulationBase<dim, spacedim> *tr =
1443 dynamic_cast<
1444 const ::parallel::
1445 DistributedTriangulationBase<dim, spacedim> *>(
1446 &dof_handler.get_triangulation()))
1447 {
1448 // For completely distributed meshes, use the function that is
1449 // able to move data from locally owned cells on one processor to
1450 // the corresponding ghost cells on others. To this end, we need
1451 // to have functions that can pack and unpack the data we want to
1452 // transport -- namely, the single unsigned int active_fe_index
1453 // objects
1454 auto pack =
1455 [](const typename ::DoFHandler<dim, spacedim>::
1456 active_cell_iterator &cell) -> active_fe_index_type {
1457 return cell->active_fe_index();
1458 };
1459
1460 auto unpack =
1461 [&dof_handler](
1462 const typename ::DoFHandler<dim, spacedim>::
1463 active_cell_iterator & cell,
1464 const active_fe_index_type active_fe_index) -> void {
1465 // we would like to say
1466 // cell->set_active_fe_index(active_fe_index);
1467 // but this is not allowed on cells that are not
1468 // locally owned, and we are on a ghost cell
1469 dof_handler
1470 .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1471 active_fe_index;
1472 };
1473
1474 GridTools::exchange_cell_data_to_ghosts<
1475 active_fe_index_type,
1476 ::DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1477 }
1478 else
1479 {
1480 // a sequential triangulation. there is nothing we need to do here
1481 Assert(
1482 (dynamic_cast<
1483 const ::parallel::TriangulationBase<dim, spacedim> *>(
1484 &dof_handler.get_triangulation()) == nullptr),
1485 ExcInternalError());
1486 }
1487 }
1488
1489
1490
1504 template <int dim, int spacedim>
1505 static void
1506 communicate_future_fe_indices(DoFHandler<dim, spacedim> &dof_handler)
1507 {
1508 Assert(
1509 dof_handler.hp_capability_enabled == true,
1510 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
1511
1512 using active_fe_index_type =
1513 typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
1514
1515 if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1516 dynamic_cast<
1517 const ::parallel::shared::Triangulation<dim, spacedim>
1518 *>(&dof_handler.get_triangulation()))
1519 {
1520 std::vector<active_fe_index_type> future_fe_indices(
1521 tr->n_active_cells(), 0u);
1522 for (const auto &cell : dof_handler.active_cell_iterators() |
1523 IteratorFilters::LocallyOwnedCell())
1524 future_fe_indices[cell->active_cell_index()] =
1525 dof_handler
1526 .hp_cell_future_fe_indices[cell->level()][cell->index()];
1527
1528 Utilities::MPI::sum(future_fe_indices,
1529 tr->get_communicator(),
1530 future_fe_indices);
1531
1532 for (const auto &cell : dof_handler.active_cell_iterators())
1533 if (!cell->is_locally_owned())
1534 dof_handler
1535 .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1536 future_fe_indices[cell->active_cell_index()];
1537 }
1538 else if (const ::parallel::
1539 DistributedTriangulationBase<dim, spacedim> *tr =
1540 dynamic_cast<
1541 const ::parallel::
1542 DistributedTriangulationBase<dim, spacedim> *>(
1543 &dof_handler.get_triangulation()))
1544 {
1545 auto pack =
1546 [&dof_handler](
1547 const typename ::DoFHandler<dim, spacedim>::
1548 active_cell_iterator &cell) -> active_fe_index_type {
1549 return dof_handler
1550 .hp_cell_future_fe_indices[cell->level()][cell->index()];
1551 };
1552
1553 auto unpack =
1554 [&dof_handler](
1555 const typename ::DoFHandler<dim, spacedim>::
1556 active_cell_iterator & cell,
1557 const active_fe_index_type future_fe_index) -> void {
1558 dof_handler
1559 .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1560 future_fe_index;
1561 };
1562
1563 GridTools::exchange_cell_data_to_ghosts<
1564 active_fe_index_type,
1565 ::DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1566 }
1567 else
1568 {
1569 Assert(
1570 (dynamic_cast<
1571 const ::parallel::TriangulationBase<dim, spacedim> *>(
1572 &dof_handler.get_triangulation()) == nullptr),
1573 ExcInternalError());
1574 }
1575 }
1576
1577
1578
1599 template <int dim, int spacedim>
1600 static void
1601 collect_fe_indices_on_cells_to_be_refined(
1602 DoFHandler<dim, spacedim> &dof_handler)
1603 {
1604 const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1605
1606 for (const auto &cell : dof_handler.active_cell_iterators())
1607 if (cell->is_locally_owned())
1608 {
1609 if (cell->refine_flag_set())
1610 {
1611 // Store the active FE index of each cell that will be
1612 // refined to and distribute it later on its children.
1613 // Pick their future index if flagged for p-refinement.
1614 fe_transfer->refined_cells_fe_index.insert(
1615 {cell, cell->future_fe_index()});
1616 }
1617 else if (cell->coarsen_flag_set())
1618 {
1619 // From all cells that will be coarsened, determine their
1620 // parent and calculate its proper active FE index, so that
1621 // it can be set after refinement. But first, check if that
1622 // particular cell has a parent at all.
1623 Assert(cell->level() > 0, ExcInternalError());
1624 const auto &parent = cell->parent();
1625
1626 // Check if the active FE index for the current cell has
1627 // been determined already.
1628 if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1629 fe_transfer->coarsened_cells_fe_index.end())
1630 {
1631 // Find a suitable active FE index for the parent cell
1632 // based on the 'least dominant finite element' of its
1633 // children. Consider the childrens' hypothetical future
1634 // index when they have been flagged for p-refinement.
1635#ifdef DEBUG
1636 for (const auto &child : parent->child_iterators())
1637 Assert(child->is_active() &&
1638 child->coarsen_flag_set(),
1639 typename ::Triangulation<
1640 dim>::ExcInconsistentCoarseningFlags());
1641#endif
1642
1643 const unsigned int fe_index = ::internal::hp::
1644 DoFHandlerImplementation::Implementation::
1645 dominated_future_fe_on_children<dim, spacedim>(
1646 parent);
1647
1648 fe_transfer->coarsened_cells_fe_index.insert(
1649 {parent, fe_index});
1650 }
1651 }
1652 else
1653 {
1654 // No h-refinement is scheduled for this cell.
1655 // However, it may have p-refinement indicators, so we
1656 // choose a new active FE index based on its flags.
1657 if (cell->future_fe_index_set() == true)
1658 fe_transfer->persisting_cells_fe_index.insert(
1659 {cell, cell->future_fe_index()});
1660 }
1661 }
1662 }
1663
1664
1665
1670 template <int dim, int spacedim>
1671 static void
1672 distribute_fe_indices_on_refined_cells(
1673 DoFHandler<dim, spacedim> &dof_handler)
1674 {
1675 const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1676
1677 // Set active FE indices on persisting cells.
1678 for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1679 {
1680 const auto &cell = persist.first;
1681
1682 if (cell->is_locally_owned())
1683 {
1684 Assert(cell->is_active(), ExcInternalError());
1685 cell->set_active_fe_index(persist.second);
1686 }
1687 }
1688
1689 // Distribute active FE indices from all refined cells on their
1690 // respective children.
1691 for (const auto &refine : fe_transfer->refined_cells_fe_index)
1692 {
1693 const auto &parent = refine.first;
1694
1695 for (const auto &child : parent->child_iterators())
1696 if (child->is_locally_owned())
1697 {
1698 Assert(child->is_active(), ExcInternalError());
1699 child->set_active_fe_index(refine.second);
1700 }
1701 }
1702
1703 // Set active FE indices on coarsened cells that have been determined
1704 // before the actual coarsening happened.
1705 for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1706 {
1707 const auto &cell = coarsen.first;
1708
1709 if (cell->is_locally_owned())
1710 {
1711 Assert(cell->is_active(), ExcInternalError());
1712 cell->set_active_fe_index(coarsen.second);
1713 }
1714 }
1715 }
1716
1717
1728 template <int dim, int spacedim>
1729 static unsigned int
1730 determine_fe_from_children(
1731 const typename Triangulation<dim, spacedim>::cell_iterator &,
1732 const std::vector<unsigned int> & children_fe_indices,
1733 const ::hp::FECollection<dim, spacedim> &fe_collection)
1734 {
1735 Assert(!children_fe_indices.empty(), ExcInternalError());
1736
1737 // convert vector to set
1738 const std::set<unsigned int> children_fe_indices_set(
1739 children_fe_indices.begin(), children_fe_indices.end());
1740
1741 const unsigned int dominated_fe_index =
1742 fe_collection.find_dominated_fe_extended(children_fe_indices_set,
1743 /*codim=*/0);
1744
1745 Assert(dominated_fe_index != numbers::invalid_unsigned_int,
1746 ExcNoDominatedFiniteElementOnChildren());
1747
1748 return dominated_fe_index;
1749 }
1750
1751
1759 template <int dim, int spacedim>
1760 static unsigned int
1761 dominated_future_fe_on_children(
1762 const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1763 {
1764 Assert(
1765 !parent->is_active(),
1766 ExcMessage(
1767 "You ask for information on children of this cell which is only "
1768 "available for active cells. This cell has no children."));
1769
1770 const auto &dof_handler = parent->get_dof_handler();
1771 Assert(
1772 dof_handler.has_hp_capabilities(),
1773 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
1774
1775 std::set<unsigned int> future_fe_indices_children;
1776 for (const auto &child : parent->child_iterators())
1777 {
1778 Assert(
1779 child->is_active(),
1780 ExcMessage(
1781 "You ask for information on children of this cell which is only "
1782 "available for active cells. One of its children is not active."));
1783
1784 // Ghost siblings might occur on parallel::shared::Triangulation
1785 // objects. The public interface does not allow to access future
1786 // FE indices on ghost cells. However, we need this information
1787 // here and thus call the internal function that does not check
1788 // for cell ownership. This requires that future FE indices have
1789 // been communicated prior to calling this function.
1790 const unsigned int future_fe_index_child =
1791 ::internal::DoFCellAccessorImplementation::
1792 Implementation::future_fe_index<dim, spacedim, false>(*child);
1793
1794 future_fe_indices_children.insert(future_fe_index_child);
1795 }
1796 Assert(!future_fe_indices_children.empty(), ExcInternalError());
1797
1798 const unsigned int future_fe_index =
1799 dof_handler.fe_collection.find_dominated_fe_extended(
1800 future_fe_indices_children,
1801 /*codim=*/0);
1802
1803 Assert(future_fe_index != numbers::invalid_unsigned_int,
1804 ExcNoDominatedFiniteElementOnChildren());
1805
1806 return future_fe_index;
1807 }
1808 };
1809
1810
1811
1815 template <int dim, int spacedim>
1816 void
1817 communicate_future_fe_indices(DoFHandler<dim, spacedim> &dof_handler)
1818 {
1819 Implementation::communicate_future_fe_indices<dim, spacedim>(
1820 dof_handler);
1821 }
1822
1823
1824
1828 template <int dim, int spacedim>
1829 unsigned int
1830 dominated_future_fe_on_children(
1831 const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1832 {
1833 return Implementation::dominated_future_fe_on_children<dim, spacedim>(
1834 parent);
1835 }
1836 } // namespace DoFHandlerImplementation
1837 } // namespace hp
1838} // namespace internal
1839
1840
1841
1842template <int dim, int spacedim>
1843DoFHandler<dim, spacedim>::DoFHandler()
1844 : hp_capability_enabled(true)
1845 , tria(nullptr, typeid(*this).name())
1846 , mg_faces(nullptr)
1847{}
1848
1849
1850
1851template <int dim, int spacedim>
1852DoFHandler<dim, spacedim>::DoFHandler(const Triangulation<dim, spacedim> &tria)
1853 : DoFHandler()
1854{
1855 reinit(tria);
1856}
1857
1858
1859
1860template <int dim, int spacedim>
1861DoFHandler<dim, spacedim>::~DoFHandler()
1862{
1863 // unsubscribe all attachments to signals of the underlying triangulation
1864 for (auto &connection : this->tria_listeners)
1865 connection.disconnect();
1866 this->tria_listeners.clear();
1867
1868 for (auto &connection : this->tria_listeners_for_transfer)
1869 connection.disconnect();
1870 this->tria_listeners_for_transfer.clear();
1871
1872 // release allocated memory
1873 // virtual functions called in constructors and destructors never use the
1874 // override in a derived class
1875 // for clarity be explicit on which function is called
1876 DoFHandler<dim, spacedim>::clear();
1877
1878 // also release the policy. this needs to happen before the
1879 // current object disappears because the policy objects
1880 // store references to the DoFhandler object they work on
1881 this->policy.reset();
1882}
1883
1884
1885
1886template <int dim, int spacedim>
1887void
1888DoFHandler<dim, spacedim>::initialize(const Triangulation<dim, spacedim> &tria,
1889 const FiniteElement<dim, spacedim> &fe)
1890{
1891 this->initialize(tria, hp::FECollection<dim, spacedim>(fe));
1892}
1893
1894
1895
1896template <int dim, int spacedim>
1897void
1898DoFHandler<dim, spacedim>::initialize(const Triangulation<dim, spacedim> &tria,
1899 const hp::FECollection<dim, spacedim> &fe)
1900{
1901 this->reinit(tria);
1902 this->distribute_dofs(fe);
1903}
1904
1905
1906
1907template <int dim, int spacedim>
1908void
1909DoFHandler<dim, spacedim>::reinit(const Triangulation<dim, spacedim> &tria)
1910{
1911 //
1912 // call destructor
1913 //
1914 // remove association with old triangulation
1915 for (auto &connection : this->tria_listeners)
1916 connection.disconnect();
1917 this->tria_listeners.clear();
1918
1919 for (auto &connection : this->tria_listeners_for_transfer)
1920 connection.disconnect();
1921 this->tria_listeners_for_transfer.clear();
1922
1923 // release allocated memory and policy
1924 DoFHandler<dim, spacedim>::clear();
1925 this->policy.reset();
1926
1927 // reset the finite element collection
1928 this->fe_collection = hp::FECollection<dim, spacedim>();
1929
1930 //
1931 // call constructor
1932 //
1933 // establish connection to new triangulation
1934 this->tria = &tria;
1935 this->setup_policy();
1936
1937 // start in hp-mode and let distribute_dofs toggle it if necessary
1938 hp_capability_enabled = true;
1939 this->connect_to_triangulation_signals();
1940 this->create_active_fe_table();
1941}
1942
1943
1944
1945/*------------------------ Cell iterator functions ------------------------*/
1946
1947template <int dim, int spacedim>
1948typename DoFHandler<dim, spacedim>::cell_iterator
1949DoFHandler<dim, spacedim>::begin(const unsigned int level) const
1950{
1951 typename Triangulation<dim, spacedim>::cell_iterator cell =
1952 this->get_triangulation().begin(level);
1953 if (cell == this->get_triangulation().end(level))
1954 return end(level);
1955 return cell_iterator(*cell, this);
1956}
1957
1958
1959
1960template <int dim, int spacedim>
1961typename DoFHandler<dim, spacedim>::active_cell_iterator
1962DoFHandler<dim, spacedim>::begin_active(const unsigned int level) const
1963{
1964 // level is checked in begin
1965 cell_iterator i = begin(level);
1966 if (i.state() != IteratorState::valid)
1967 return i;
1968 while (i->has_children())
1969 if ((++i).state() != IteratorState::valid)
1970 return i;
1971 return i;
1972}
1973
1974
1975
1976template <int dim, int spacedim>
1977typename DoFHandler<dim, spacedim>::cell_iterator
1978DoFHandler<dim, spacedim>::end() const
1979{
1980 return cell_iterator(&this->get_triangulation(), -1, -1, this);
1981}
1982
1983
1984
1985template <int dim, int spacedim>
1986typename DoFHandler<dim, spacedim>::cell_iterator
1987DoFHandler<dim, spacedim>::end(const unsigned int level) const
1988{
1989 typename Triangulation<dim, spacedim>::cell_iterator cell =
1990 this->get_triangulation().end(level);
1991 if (cell.state() != IteratorState::valid)
1992 return end();
1993 return cell_iterator(*cell, this);
1994}
1995
1996
1997
1998template <int dim, int spacedim>
1999typename DoFHandler<dim, spacedim>::active_cell_iterator
2000DoFHandler<dim, spacedim>::end_active(const unsigned int level) const
2001{
2002 typename Triangulation<dim, spacedim>::cell_iterator cell =
2003 this->get_triangulation().end_active(level);
2004 if (cell.state() != IteratorState::valid)
2005 return active_cell_iterator(end());
2006 return active_cell_iterator(*cell, this);
2007}
2008
2009
2010
2011template <int dim, int spacedim>
2012typename DoFHandler<dim, spacedim>::level_cell_iterator
2013DoFHandler<dim, spacedim>::begin_mg(const unsigned int level) const
2014{
2015 Assert(this->has_level_dofs(),
2016 ExcMessage("You can only iterate over mg "
2017 "levels if mg dofs got distributed."));
2018 typename Triangulation<dim, spacedim>::cell_iterator cell =
2019 this->get_triangulation().begin(level);
2020 if (cell == this->get_triangulation().end(level))
2021 return end_mg(level);
2022 return level_cell_iterator(*cell, this);
2023}
2024
2025
2026
2027template <int dim, int spacedim>
2028typename DoFHandler<dim, spacedim>::level_cell_iterator
2029DoFHandler<dim, spacedim>::end_mg(const unsigned int level) const
2030{
2031 Assert(this->has_level_dofs(),
2032 ExcMessage("You can only iterate over mg "
2033 "levels if mg dofs got distributed."));
2034 typename Triangulation<dim, spacedim>::cell_iterator cell =
2035 this->get_triangulation().end(level);
2036 if (cell.state() != IteratorState::valid)
2037 return end();
2038 return level_cell_iterator(*cell, this);
2039}
2040
2041
2042
2043template <int dim, int spacedim>
2044typename DoFHandler<dim, spacedim>::level_cell_iterator
2045DoFHandler<dim, spacedim>::end_mg() const
2046{
2047 return level_cell_iterator(&this->get_triangulation(), -1, -1, this);
2048}
2049
2050
2051
2052template <int dim, int spacedim>
2053IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>
2054DoFHandler<dim, spacedim>::cell_iterators() const
2055{
2056 return IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>(
2057 begin(), end());
2058}
2059
2060
2061
2062template <int dim, int spacedim>
2063IteratorRange<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2064DoFHandler<dim, spacedim>::active_cell_iterators() const
2065{
2066 return IteratorRange<
2067 typename DoFHandler<dim, spacedim>::active_cell_iterator>(begin_active(),
2068 end());
2069}
2070
2071
2072
2073template <int dim, int spacedim>
2074IteratorRange<typename DoFHandler<dim, spacedim>::level_cell_iterator>
2075DoFHandler<dim, spacedim>::mg_cell_iterators() const
2076{
2077 return IteratorRange<typename DoFHandler<dim, spacedim>::level_cell_iterator>(
2078 begin_mg(), end_mg());
2079}
2080
2081
2082
2083template <int dim, int spacedim>
2084IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>
2085DoFHandler<dim, spacedim>::cell_iterators_on_level(
2086 const unsigned int level) const
2087{
2088 return IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>(
2089 begin(level), end(level));
2090}
2091
2092
2093
2094template <int dim, int spacedim>
2095IteratorRange<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2096DoFHandler<dim, spacedim>::active_cell_iterators_on_level(
2097 const unsigned int level) const
2098{
2099 return IteratorRange<
2100 typename DoFHandler<dim, spacedim>::active_cell_iterator>(
2101 begin_active(level), end_active(level));
2102}
2103
2104
2105
2106template <int dim, int spacedim>
2107IteratorRange<typename DoFHandler<dim, spacedim>::level_cell_iterator>
2108DoFHandler<dim, spacedim>::mg_cell_iterators_on_level(
2109 const unsigned int level) const
2110{
2111 return IteratorRange<typename DoFHandler<dim, spacedim>::level_cell_iterator>(
2112 begin_mg(level), end_mg(level));
2113}
2114
2115
2116
2117//---------------------------------------------------------------------------
2118
2119
2120
2121template <int dim, int spacedim>
2122types::global_dof_index
2123DoFHandler<dim, spacedim>::n_boundary_dofs() const
2124{
2125 Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2126 ExcNotImplementedWithHP());
2127
2128 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2129
2130 std::unordered_set<types::global_dof_index> boundary_dofs;
2131 std::vector<types::global_dof_index> dofs_on_face;
2132 dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2133
2134 const IndexSet &owned_dofs = locally_owned_dofs();
2135
2136 // loop over all faces to check whether they are at a
2137 // boundary. note that we need not take special care of single
2138 // lines in 3d (using @p{cell->has_boundary_lines}), since we do
2139 // not support boundaries of dimension dim-2, and so every
2140 // boundary line is also part of a boundary face.
2141 for (const auto &cell : this->active_cell_iterators())
2142 if (cell->is_locally_owned() && cell->at_boundary())
2143 {
2144 for (const auto iface : cell->face_indices())
2145 {
2146 const auto face = cell->face(iface);
2147 if (face->at_boundary())
2148 {
2149 const unsigned int dofs_per_face =
2150 cell->get_fe().n_dofs_per_face(iface);
2151 dofs_on_face.resize(dofs_per_face);
2152
2153 face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2154 for (unsigned int i = 0; i < dofs_per_face; ++i)
2155 {
2156 const unsigned int global_idof_index = dofs_on_face[i];
2157 if (owned_dofs.is_element(global_idof_index))
2158 {
2159 boundary_dofs.insert(global_idof_index);
2160 }
2161 }
2162 }
2163 }
2164 }
2165 return boundary_dofs.size();
2166}
2167
2168
2169
2170template <int dim, int spacedim>
2171types::global_dof_index
2172DoFHandler<dim, spacedim>::n_boundary_dofs(
2173 const std::set<types::boundary_id> &boundary_ids) const
2174{
2175 Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2176 ExcNotImplementedWithHP());
2177
2178 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2179 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2180 boundary_ids.end(),
2181 ExcInvalidBoundaryIndicator());
2182
2183 // same as above, but with additional checks for set of boundary
2184 // indicators
2185 std::unordered_set<types::global_dof_index> boundary_dofs;
2186 std::vector<types::global_dof_index> dofs_on_face;
2187 dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2188
2189 const IndexSet &owned_dofs = locally_owned_dofs();
2190
2191 for (const auto &cell : this->active_cell_iterators())
2192 if (cell->is_locally_owned() && cell->at_boundary())
2193 {
2194 for (const auto iface : cell->face_indices())
2195 {
2196 const auto face = cell->face(iface);
2197 const unsigned int boundary_id = face->boundary_id();
2198 if (face->at_boundary() &&
2199 (boundary_ids.find(boundary_id) != boundary_ids.end()))
2200 {
2201 const unsigned int dofs_per_face =
2202 cell->get_fe().n_dofs_per_face(iface);
2203 dofs_on_face.resize(dofs_per_face);
2204
2205 face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2206 for (unsigned int i = 0; i < dofs_per_face; ++i)
2207 {
2208 const unsigned int global_idof_index = dofs_on_face[i];
2209 if (owned_dofs.is_element(global_idof_index))
2210 {
2211 boundary_dofs.insert(global_idof_index);
2212 }
2213 }
2214 }
2215 }
2216 }
2217 return boundary_dofs.size();
2218}
2219
2220
2221
2222template <int dim, int spacedim>
2223std::size_t
2224DoFHandler<dim, spacedim>::memory_consumption() const
2225{
2226 std::size_t mem = MemoryConsumption::memory_consumption(this->tria) +
2227 MemoryConsumption::memory_consumption(this->fe_collection) +
2228 MemoryConsumption::memory_consumption(this->number_cache);
2229
2230 mem += MemoryConsumption::memory_consumption(cell_dof_cache_indices) +
2231 MemoryConsumption::memory_consumption(cell_dof_cache_ptr) +
2232 MemoryConsumption::memory_consumption(object_dof_indices) +
2233 MemoryConsumption::memory_consumption(object_dof_ptr) +
2234 MemoryConsumption::memory_consumption(hp_object_fe_indices) +
2235 MemoryConsumption::memory_consumption(hp_object_fe_ptr) +
2236 MemoryConsumption::memory_consumption(hp_cell_active_fe_indices) +
2237 MemoryConsumption::memory_consumption(hp_cell_future_fe_indices);
2238
2239
2240 if (hp_capability_enabled)
2241 {
2242 // nothing to add
2243 }
2244 else
2245 {
2246 // collect size of multigrid data structures
2247
2248 mem += MemoryConsumption::memory_consumption(this->block_info_object);
2249
2250 for (unsigned int level = 0; level < this->mg_levels.size(); ++level)
2251 mem += this->mg_levels[level]->memory_consumption();
2252
2253 if (this->mg_faces != nullptr)
2254 mem += MemoryConsumption::memory_consumption(*this->mg_faces);
2255
2256 for (unsigned int i = 0; i < this->mg_vertex_dofs.size(); ++i)
2257 mem += sizeof(MGVertexDoFs) +
2258 (1 + this->mg_vertex_dofs[i].get_finest_level() -
2259 this->mg_vertex_dofs[i].get_coarsest_level()) *
2260 sizeof(types::global_dof_index);
2261 }
2262
2263 return mem;
2264}
2265
2266
2267
2268template <int dim, int spacedim>
2269void
2270DoFHandler<dim, spacedim>::set_fe(const FiniteElement<dim, spacedim> &fe)
2271{
2272 this->set_fe(hp::FECollection<dim, spacedim>(fe));
2273}
2274
2275
2276
2277template <int dim, int spacedim>
2278void
2279DoFHandler<dim, spacedim>::set_fe(const hp::FECollection<dim, spacedim> &ff)
2280{
2281 Assert(
2282 this->tria != nullptr,
2283 ExcMessage(
2284 "You need to set the Triangulation in the DoFHandler using reinit() or "
2285 "in the constructor before you can distribute DoFs."));
2286 Assert(this->tria->n_levels() > 0,
2287 ExcMessage("The Triangulation you are using is empty!"));
2288 Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
2289
2290 // don't create a new object if the one we have is already appropriate
2291 if (this->fe_collection != ff)
2292 {
2293 this->fe_collection = hp::FECollection<dim, spacedim>(ff);
2294
2295 const bool contains_multiple_fes = (this->fe_collection.size() > 1);
2296
2297 // disable hp-mode if only a single finite element has been registered
2298 if (hp_capability_enabled && !contains_multiple_fes)
2299 {
2300 hp_capability_enabled = false;
2301
2302 // unsubscribe connections to signals that are only relevant for
2303 // hp-mode, since we only have a single element here
2304 for (auto &connection : this->tria_listeners_for_transfer)
2305 connection.disconnect();
2306 this->tria_listeners_for_transfer.clear();
2307
2308 // release active and future finite element tables
2309 this->hp_cell_active_fe_indices.clear();
2310 this->hp_cell_active_fe_indices.shrink_to_fit();
2311 this->hp_cell_future_fe_indices.clear();
2312 this->hp_cell_future_fe_indices.shrink_to_fit();
2313 }
2314
2315 // re-enabling hp-mode is not permitted since the active and future FE
2316 // tables are no longer available
2317 AssertThrow(
2318 hp_capability_enabled || !contains_multiple_fes,
2319 ExcMessage(
2320 "You cannot re-enable hp-capabilities after you registered a single "
2321 "finite element. Please create a new DoFHandler object instead."));
2322 }
2323
2324 if (hp_capability_enabled)
2325 {
2326 // make sure every processor knows the active FE indices
2327 // on both its own cells and all ghost cells
2328 ::internal::hp::DoFHandlerImplementation::Implementation::
2329 communicate_active_fe_indices(*this);
2330
2331 // make sure that the FE collection is large enough to
2332 // cover all FE indices presently in use on the mesh
2333 for (const auto &cell : this->active_cell_iterators())
2334 if (!cell->is_artificial())
2335 Assert(cell->active_fe_index() < this->fe_collection.size(),
2336 ExcInvalidFEIndex(cell->active_fe_index(),
2337 this->fe_collection.size()));
2338 }
2339}
2340
2341
2342
2343template <int dim, int spacedim>
2344void
2345DoFHandler<dim, spacedim>::distribute_dofs(
2346 const FiniteElement<dim, spacedim> &fe)
2347{
2348 this->distribute_dofs(hp::FECollection<dim, spacedim>(fe));
2349}
2350
2351
2352
2353template <int dim, int spacedim>
2354void
2355DoFHandler<dim, spacedim>::distribute_dofs(
2356 const hp::FECollection<dim, spacedim> &ff)
2357{
2358 Assert(
2359 this->tria != nullptr,
2360 ExcMessage(
2361 "You need to set the Triangulation in the DoFHandler using reinit() or "
2362 "in the constructor before you can distribute DoFs."));
2363 Assert(this->tria->n_levels() > 0,
2364 ExcMessage("The Triangulation you are using is empty!"));
2365 Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
2366
2367 //
2368 // register the new finite element collection
2369 //
2370 // don't create a new object if the one we have is identical
2371 if (this->fe_collection != ff)
2372 {
2373 this->fe_collection = hp::FECollection<dim, spacedim>(ff);
2374
2375 const bool contains_multiple_fes = (this->fe_collection.size() > 1);
2376
2377 // disable hp-mode if only a single finite element has been registered
2378 if (hp_capability_enabled && !contains_multiple_fes)
2379 {
2380 hp_capability_enabled = false;
2381
2382 // unsubscribe connections to signals that are only relevant for
2383 // hp-mode, since we only have a single element here
2384 for (auto &connection : this->tria_listeners_for_transfer)
2385 connection.disconnect();
2386 this->tria_listeners_for_transfer.clear();
2387
2388 // release active and future finite element tables
2389 this->hp_cell_active_fe_indices.clear();
2390 this->hp_cell_active_fe_indices.shrink_to_fit();
2391 this->hp_cell_future_fe_indices.clear();
2392 this->hp_cell_future_fe_indices.shrink_to_fit();
2393 }
2394
2395 // re-enabling hp-mode is not permitted since the active and future FE
2396 // tables are no longer available
2397 AssertThrow(
2398 hp_capability_enabled || !contains_multiple_fes,
2399 ExcMessage(
2400 "You cannot re-enable hp-capabilities after you registered a single "
2401 "finite element. Please call reinit() or create a new DoFHandler "
2402 "object instead."));
2403 }
2404
2405 //
2406 // enumerate all degrees of freedom
2407 //
2408 if (hp_capability_enabled)
2409 {
2410 // make sure every processor knows the active FE indices
2411 // on both its own cells and all ghost cells
2412 internal::hp::DoFHandlerImplementation::Implementation::
2413 communicate_active_fe_indices(*this);
2414
2415#ifdef DEBUG
2416 // make sure that the FE collection is large enough to
2417 // cover all FE indices presently in use on the mesh
2418 for (const auto &cell : this->active_cell_iterators())
2419 {
2420 if (!cell->is_artificial())
2421 Assert(cell->active_fe_index() < this->fe_collection.size(),
2422 ExcInvalidFEIndex(cell->active_fe_index(),
2423 this->fe_collection.size()));
2424 if (cell->is_locally_owned())
2425 Assert(cell->future_fe_index() < this->fe_collection.size(),
2426 ExcInvalidFEIndex(cell->future_fe_index(),
2427 this->fe_collection.size()));
2428 }
2429#endif
2430 }
2431
2432 {
2433 // We would like to enumerate all dofs for shared::Triangulations. If an
2434 // underlying shared::Tria allows artificial cells, we need to restore the
2435 // true cell owners temporarily.
2436 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
2437 // the current set of subdomain ids, set subdomain ids to the "true" owner
2438 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
2439 // object, and later restore these flags when it is destroyed.
2440 const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
2441 spacedim>
2442 subdomain_modifier(this->get_triangulation());
2443
2444 // Adjust size of levels to the triangulation. Note that we still have to
2445 // allocate space for all degrees of freedom on this mesh (including ghost
2446 // and cells that are entirely stored on different processors), though we
2447 // may not assign numbers to some of them (i.e. they will remain at
2448 // invalid_dof_index). We need to allocate the space because we will want
2449 // to be able to query the dof_indices on each cell, and simply be told
2450 // that we don't know them on some cell (i.e. get back invalid_dof_index)
2451 if (hp_capability_enabled)
2452 internal::hp::DoFHandlerImplementation::Implementation::reserve_space(
2453 *this);
2454 else
2455 internal::DoFHandlerImplementation::Implementation::reserve_space(*this);
2456 }
2457
2458 // hand the actual work over to the policy
2459 this->number_cache = this->policy->distribute_dofs();
2460
2461 // do some housekeeping: compress indices
2462 // if(hp_capability_enabled)
2463 // {
2464 // Threads::TaskGroup<> tg;
2465 // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2466 // tg += Threads::new_task(
2467 // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2468 // *this->levels_hp[level],
2469 // this->fe_collection);
2470 // tg.join_all();
2471 // }
2472
2473 // Initialize the block info object only if this is a sequential
2474 // triangulation. It doesn't work correctly yet if it is parallel and has not
2475 // yet been implemented for hp-mode.
2476 if (!hp_capability_enabled &&
2477 dynamic_cast<const parallel::DistributedTriangulationBase<dim, spacedim>
2478 *>(&*this->tria) == nullptr)
2479 this->block_info_object.initialize(*this, false, true);
2480}
2481
2482
2483
2484template <int dim, int spacedim>
2485void
2486DoFHandler<dim, spacedim>::distribute_mg_dofs()
2487{
2488 AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2489
2490 Assert(
2491 this->object_dof_indices.size() > 0,
2492 ExcMessage(
2493 "Distribute active DoFs using distribute_dofs() before calling distribute_mg_dofs()."));
2494
2495 Assert(
2496 ((this->tria->get_mesh_smoothing() &
2497 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) !=
2498 Triangulation<dim, spacedim>::none),
2499 ExcMessage(
2500 "The mesh smoothing requirement 'limit_level_difference_at_vertices' has to be set for using multigrid!"));
2501
2502 this->clear_mg_space();
2503
2504 internal::DoFHandlerImplementation::Implementation::reserve_space_mg(*this);
2505 this->mg_number_cache = this->policy->distribute_mg_dofs();
2506
2507 // initialize the block info object only if this is a sequential
2508 // triangulation. it doesn't work correctly yet if it is parallel
2509 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2510 &*this->tria) == nullptr)
2511 this->block_info_object.initialize(*this, true, false);
2512}
2513
2514
2515
2516template <int dim, int spacedim>
2517void
2518DoFHandler<dim, spacedim>::initialize_local_block_info()
2519{
2520 AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2521
2522 this->block_info_object.initialize_local(*this);
2523}
2524
2525
2526
2527template <int dim, int spacedim>
2528void
2529DoFHandler<dim, spacedim>::setup_policy()
2530{
2531 // decide whether we need a sequential or a parallel distributed policy
2532 if (dynamic_cast<const ::parallel::shared::Triangulation<dim, spacedim>
2533 *>(&this->get_triangulation()) != nullptr)
2534 this->policy = std::make_unique<internal::DoFHandlerImplementation::Policy::
2535 ParallelShared<dim, spacedim>>(*this);
2536 else if (dynamic_cast<
2537 const ::parallel::DistributedTriangulationBase<dim, spacedim>
2538 *>(&this->get_triangulation()) == nullptr)
2539 this->policy = std::make_unique<
2540 internal::DoFHandlerImplementation::Policy::Sequential<dim, spacedim>>(
2541 *this);
2542 else
2543 this->policy =
2544 std::make_unique<internal::DoFHandlerImplementation::Policy::
2545 ParallelDistributed<dim, spacedim>>(*this);
2546}
2547
2548
2549
2550template <int dim, int spacedim>
2551void
2552DoFHandler<dim, spacedim>::clear()
2553{
2554 // release memory
2555 this->clear_space();
2556 this->clear_mg_space();
2557}
2558
2559
2560
2561template <int dim, int spacedim>
2562void
2563DoFHandler<dim, spacedim>::clear_space()
2564{
2565 cell_dof_cache_indices.clear();
2566
2567 cell_dof_cache_ptr.clear();
2568
2569 object_dof_indices.clear();
2570
2571 object_dof_ptr.clear();
2572
2573 this->number_cache.clear();
2574
2575 this->hp_cell_active_fe_indices.clear();
2576 this->hp_cell_future_fe_indices.clear();
2577}
2578
2579
2580
2581template <int dim, int spacedim>
2582void
2583DoFHandler<dim, spacedim>::clear_mg_space()
2584{
2585 this->mg_levels.clear();
2586 this->mg_faces.reset();
2587
2588 std::vector<MGVertexDoFs> tmp;
2589
2590 std::swap(this->mg_vertex_dofs, tmp);
2591
2592 this->mg_number_cache.clear();
2593}
2594
2595
2596
2597template <int dim, int spacedim>
2598void
2599DoFHandler<dim, spacedim>::renumber_dofs(
2600 const std::vector<types::global_dof_index> &new_numbers)
2601{
2602 if (hp_capability_enabled)
2603 {
2604 Assert(this->hp_cell_future_fe_indices.size() > 0,
2605 ExcMessage(
2606 "You need to distribute DoFs before you can renumber them."));
2607
2608 AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2609
2610#ifdef DEBUG
2611 // assert that the new indices are consecutively numbered if we are
2612 // working on a single processor. this doesn't need to
2613 // hold in the case of a parallel mesh since we map the interval
2614 // [0...n_dofs()) into itself but only globally, not on each processor
2615 if (this->n_locally_owned_dofs() == this->n_dofs())
2616 {
2617 std::vector<types::global_dof_index> tmp(new_numbers);
2618 std::sort(tmp.begin(), tmp.end());
2619 std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2620 types::global_dof_index i = 0;
2621 for (; p != tmp.end(); ++p, ++i)
2622 Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2623 }
2624 else
2625 for (const auto new_number : new_numbers)
2626 Assert(new_number < this->n_dofs(),
2627 ExcMessage(
2628 "New DoF index is not less than the total number of dofs."));
2629#endif
2630
2631 // uncompress the internal storage scheme of dofs on cells so that
2632 // we can access dofs in turns. uncompress in parallel, starting
2633 // with the most expensive levels (the highest ones)
2634 //{
2635 // Threads::TaskGroup<> tg;
2636 // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2637 // tg += Threads::new_task(
2638 // &::internal::hp::DoFLevel::uncompress_data<dim, spacedim>,
2639 // *this->levels_hp[level],
2640 // this->fe_collection);
2641 // tg.join_all();
2642 //}
2643
2644 // do the renumbering
2645 this->number_cache = this->policy->renumber_dofs(new_numbers);
2646
2647 // now re-compress the dof indices
2648 //{
2649 // Threads::TaskGroup<> tg;
2650 // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2651 // tg += Threads::new_task(
2652 // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2653 // *this->levels_hp[level],
2654 // this->fe_collection);
2655 // tg.join_all();
2656 //}
2657 }
2658 else
2659 {
2660 Assert(this->object_dof_indices.size() > 0,
2661 ExcMessage(
2662 "You need to distribute DoFs before you can renumber them."));
2663
2664#ifdef DEBUG
2665 if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
2666 &*this->tria) != nullptr)
2667 {
2668 Assert(new_numbers.size() == this->n_dofs() ||
2669 new_numbers.size() == this->n_locally_owned_dofs(),
2670 ExcMessage("Incorrect size of the input array."));
2671 }
2672 else if (dynamic_cast<
2673 const parallel::DistributedTriangulationBase<dim, spacedim> *>(
2674 &*this->tria) != nullptr)
2675 {
2676 AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2677 }
2678 else
2679 {
2680 AssertDimension(new_numbers.size(), this->n_dofs());
2681 }
2682
2683 // assert that the new indices are consecutively numbered if we are
2684 // working on a single processor. this doesn't need to
2685 // hold in the case of a parallel mesh since we map the interval
2686 // [0...n_dofs()) into itself but only globally, not on each processor
2687 if (this->n_locally_owned_dofs() == this->n_dofs())
2688 {
2689 std::vector<types::global_dof_index> tmp(new_numbers);
2690 std::sort(tmp.begin(), tmp.end());
2691 std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2692 types::global_dof_index i = 0;
2693 for (; p != tmp.end(); ++p, ++i)
2694 Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2695 }
2696 else
2697 for (const auto new_number : new_numbers)
2698 Assert(new_number < this->n_dofs(),
2699 ExcMessage(
2700 "New DoF index is not less than the total number of dofs."));
2701#endif
2702
2703 this->number_cache = this->policy->renumber_dofs(new_numbers);
2704 }
2705}
2706
2707
2708
2709template <int dim, int spacedim>
2710void
2711DoFHandler<dim, spacedim>::renumber_dofs(
2712 const unsigned int level,
2713 const std::vector<types::global_dof_index> &new_numbers)
2714{
2715 AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2716
2717 Assert(
2718 this->mg_levels.size() > 0 && this->object_dof_indices.size() > 0,
2719 ExcMessage(
2720 "You need to distribute active and level DoFs before you can renumber level DoFs."));
2721 AssertIndexRange(level, this->get_triangulation().n_global_levels());
2722 AssertDimension(new_numbers.size(),
2723 this->locally_owned_mg_dofs(level).n_elements());
2724
2725#ifdef DEBUG
2726 // assert that the new indices are consecutively numbered if we are working
2727 // on a single processor. this doesn't need to hold in the case of a
2728 // parallel mesh since we map the interval [0...n_dofs(level)) into itself
2729 // but only globally, not on each processor
2730 if (this->n_locally_owned_dofs() == this->n_dofs())
2731 {
2732 std::vector<types::global_dof_index> tmp(new_numbers);
2733 std::sort(tmp.begin(), tmp.end());
2734 std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2735 types::global_dof_index i = 0;
2736 for (; p != tmp.end(); ++p, ++i)
2737 Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2738 }
2739 else
2740 for (const auto new_number : new_numbers)
2741 Assert(new_number < this->n_dofs(level),
2742 ExcMessage(
2743 "New DoF index is not less than the total number of dofs."));
2744#endif
2745
2746 this->mg_number_cache[level] =
2747 this->policy->renumber_mg_dofs(level, new_numbers);
2748}
2749
2750
2751
2752template <int dim, int spacedim>
2753unsigned int
2754DoFHandler<dim, spacedim>::max_couplings_between_boundary_dofs() const
2755{
2756 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2757
2758 switch (dim)
2759 {
2760 case 1:
2761 return this->fe_collection.max_dofs_per_vertex();
2762 case 2:
2763 return (3 * this->fe_collection.max_dofs_per_vertex() +
2764 2 * this->fe_collection.max_dofs_per_line());
2765 case 3:
2766 // we need to take refinement of one boundary face into
2767 // consideration here; in fact, this function returns what
2768 // #max_coupling_between_dofs<2> returns
2769 //
2770 // we assume here, that only four faces meet at the boundary;
2771 // this assumption is not justified and needs to be fixed some
2772 // time. fortunately, omitting it for now does no harm since
2773 // the matrix will cry foul if its requirements are not
2774 // satisfied
2775 return (19 * this->fe_collection.max_dofs_per_vertex() +
2776 28 * this->fe_collection.max_dofs_per_line() +
2777 8 * this->fe_collection.max_dofs_per_quad());
2778 default:
2779 Assert(false, ExcNotImplemented());
2780 return 0;
2781 }
2782}
2783
2784
2785
2786template <int dim, int spacedim>
2787unsigned int
2788DoFHandler<dim, spacedim>::max_couplings_between_dofs() const
2789{
2790 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2791 return internal::DoFHandlerImplementation::Implementation::
2792 max_couplings_between_dofs(*this);
2793}
2794
2795
2796
2797template <int dim, int spacedim>
2798void
2799DoFHandler<dim, spacedim>::set_active_fe_indices(
2800 const std::vector<unsigned int> &active_fe_indices)
2801{
2802 Assert(active_fe_indices.size() == this->get_triangulation().n_active_cells(),
2803 ExcDimensionMismatch(active_fe_indices.size(),
2804 this->get_triangulation().n_active_cells()));
2805
2806 this->create_active_fe_table();
2807 // we could set the values directly, since they are stored as
2808 // protected data of this object, but for simplicity we use the
2809 // cell-wise access. this way we also have to pass some debug-mode
2810 // tests which we would have to duplicate ourselves otherwise
2811 for (const auto &cell : this->active_cell_iterators())
2812 if (cell->is_locally_owned())
2813 cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
2814}
2815
2816
2817
2818template <int dim, int spacedim>
2819void
2820DoFHandler<dim, spacedim>::get_active_fe_indices(
2821 std::vector<unsigned int> &active_fe_indices) const
2822{
2823 active_fe_indices.resize(this->get_triangulation().n_active_cells());
2824
2825 // we could try to extract the values directly, since they are
2826 // stored as protected data of this object, but for simplicity we
2827 // use the cell-wise access.
2828 for (const auto &cell : this->active_cell_iterators())
2829 if (!cell->is_artificial())
2830 active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
2831}
2832
2833
2834
2835template <int dim, int spacedim>
2836void
2837DoFHandler<dim, spacedim>::connect_to_triangulation_signals()
2838{
2839 // make sure this is called during initialization in hp-mode
2840 Assert(hp_capability_enabled, ExcOnlyAvailableWithHP());
2841
2842 // connect functions to signals of the underlying triangulation
2843 this->tria_listeners.push_back(this->tria->signals.create.connect(
2844 [this]() { this->reinit(*(this->tria)); }));
2845 this->tria_listeners.push_back(
2846 this->tria->signals.clear.connect([this]() { this->clear(); }));
2847
2848 // attach corresponding callback functions dealing with the transfer of
2849 // active FE indices depending on the type of triangulation
2850 if (dynamic_cast<
2851 const ::parallel::fullydistributed::Triangulation<dim, spacedim>
2852 *>(&this->get_triangulation()))
2853 {
2854 // no transfer of active FE indices for this class
2855 }
2856 else if (dynamic_cast<
2857 const ::parallel::distributed::Triangulation<dim, spacedim>
2858 *>(&this->get_triangulation()))
2859 {
2860 // repartitioning signals
2861 this->tria_listeners_for_transfer.push_back(
2862 this->tria->signals.pre_distributed_repartition.connect([this]() {
2863 internal::hp::DoFHandlerImplementation::Implementation::
2864 ensure_absence_of_future_fe_indices<dim, spacedim>(*this);
2865 }));
2866 this->tria_listeners_for_transfer.push_back(
2867 this->tria->signals.pre_distributed_repartition.connect(
2868 [this]() { this->pre_distributed_transfer_action(); }));
2869 this->tria_listeners_for_transfer.push_back(
2870 this->tria->signals.post_distributed_repartition.connect(
2871 [this] { this->post_distributed_transfer_action(); }));
2872
2873 // refinement signals
2874 this->tria_listeners_for_transfer.push_back(
2875 this->tria->signals.post_p4est_refinement.connect(
2876 [this]() { this->pre_distributed_transfer_action(); }));
2877 this->tria_listeners_for_transfer.push_back(
2878 this->tria->signals.post_distributed_refinement.connect(
2879 [this]() { this->post_distributed_transfer_action(); }));
2880
2881 // serialization signals
2882 this->tria_listeners_for_transfer.push_back(
2883 this->tria->signals.post_distributed_save.connect(
2884 [this]() { this->active_fe_index_transfer.reset(); }));
2885 this->tria_listeners_for_transfer.push_back(
2886 this->tria->signals.post_distributed_load.connect(
2887 [this]() { this->update_active_fe_table(); }));
2888 }
2889 else if (dynamic_cast<
2890 const ::parallel::shared::Triangulation<dim, spacedim> *>(
2891 &this->get_triangulation()) != nullptr)
2892 {
2893 // partitioning signals
2894 this->tria_listeners_for_transfer.push_back(
2895 this->tria->signals.pre_partition.connect([this]() {
2896 internal::hp::DoFHandlerImplementation::Implementation::
2897 ensure_absence_of_future_fe_indices(*this);
2898 }));
2899
2900 // refinement signals
2901 this->tria_listeners_for_transfer.push_back(
2902 this->tria->signals.pre_refinement.connect([this]() {
2903 internal::hp::DoFHandlerImplementation::Implementation::
2904 communicate_future_fe_indices(*this);
2905 }));
2906 this->tria_listeners_for_transfer.push_back(
2907 this->tria->signals.pre_refinement.connect(
2908 [this] { this->pre_transfer_action(); }));
2909 this->tria_listeners_for_transfer.push_back(
2910 this->tria->signals.post_refinement.connect(
2911 [this] { this->post_transfer_action(); }));
2912 }
2913 else
2914 {
2915 // refinement signals
2916 this->tria_listeners_for_transfer.push_back(
2917 this->tria->signals.pre_refinement.connect(
2918 [this] { this->pre_transfer_action(); }));
2919 this->tria_listeners_for_transfer.push_back(
2920 this->tria->signals.post_refinement.connect(
2921 [this] { this->post_transfer_action(); }));
2922 }
2923}
2924
2925
2926
2927template <int dim, int spacedim>
2928void
2929DoFHandler<dim, spacedim>::create_active_fe_table()
2930{
2931 AssertThrow(hp_capability_enabled == true, ExcOnlyAvailableWithHP());
2932
2933
2934 // Create sufficiently many hp::DoFLevels.
2935 // while (this->levels_hp.size() < this->tria->n_levels())
2936 // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
2937
2938 this->hp_cell_active_fe_indices.resize(this->tria->n_levels());
2939 this->hp_cell_future_fe_indices.resize(this->tria->n_levels());
2940
2941 // then make sure that on each level we have the appropriate size
2942 // of active FE indices; preset them to zero, i.e. the default FE
2943 for (unsigned int level = 0; level < this->hp_cell_future_fe_indices.size();
2944 ++level)
2945 {
2946 if (this->hp_cell_active_fe_indices[level].size() == 0 &&
2947 this->hp_cell_future_fe_indices[level].size() == 0)
2948 {
2949 this->hp_cell_active_fe_indices[level].resize(
2950 this->tria->n_raw_cells(level), 0);
2951 this->hp_cell_future_fe_indices[level].resize(
2952 this->tria->n_raw_cells(level), invalid_active_fe_index);
2953 }
2954 else
2955 {
2956 // Either the active FE indices have size zero because
2957 // they were just created, or the correct size. Other
2958 // sizes indicate that something went wrong.
2959 Assert(this->hp_cell_active_fe_indices[level].size() ==
2960 this->tria->n_raw_cells(level) &&
2961 this->hp_cell_future_fe_indices[level].size() ==
2962 this->tria->n_raw_cells(level),
2963 ExcInternalError());
2964 }
2965
2966 // it may be that the previous table was compressed; in that
2967 // case, restore the correct active FE index. the fact that
2968 // this no longer matches the indices in the table is of no
2969 // importance because the current function is called at a
2970 // point where we have to recreate the dof_indices tables in
2971 // the levels anyway
2972 // this->levels_hp[level]->normalize_active_fe_indices();
2973 }
2974}
2975
2976
2977
2978template <int dim, int spacedim>
2979void
2980DoFHandler<dim, spacedim>::update_active_fe_table()
2981{
2982 // // Normally only one level is added, but if this Triangulation
2983 // // is created by copy_triangulation, it can be more than one level.
2984 // while (this->levels_hp.size() < this->tria->n_levels())
2985 // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
2986 //
2987 // // Coarsening can lead to the loss of levels. Hence remove them.
2988 // while (this->levels_hp.size() > this->tria->n_levels())
2989 // {
2990 // // drop the last element. that also releases the memory pointed to
2991 // this->levels_hp.pop_back();
2992 // }
2993
2994 this->hp_cell_active_fe_indices.resize(this->tria->n_levels());
2995 this->hp_cell_active_fe_indices.shrink_to_fit();
2996
2997 this->hp_cell_future_fe_indices.resize(this->tria->n_levels());
2998 this->hp_cell_future_fe_indices.shrink_to_fit();
2999
3000 for (unsigned int i = 0; i < this->hp_cell_future_fe_indices.size(); ++i)
3001 {
3002 // Resize active FE indices vectors. Use zero indicator to extend.
3003 this->hp_cell_active_fe_indices[i].resize(this->tria->n_raw_cells(i), 0);
3004
3005 // Resize future FE indices vectors. Make sure that all
3006 // future FE indices have been cleared after refinement happened.
3007 //
3008 // We have used future FE indices to update all active FE indices
3009 // before refinement happened, thus we are safe to clear them now.
3010 this->hp_cell_future_fe_indices[i].assign(this->tria->n_raw_cells(i),
3011 invalid_active_fe_index);
3012 }
3013}
3014
3015
3016template <int dim, int spacedim>
3017void
3018DoFHandler<dim, spacedim>::pre_transfer_action()
3019{
3020 Assert(this->active_fe_index_transfer == nullptr, ExcInternalError());
3021
3022 this->active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3023
3024 ::internal::hp::DoFHandlerImplementation::Implementation::
3025 collect_fe_indices_on_cells_to_be_refined(*this);
3026}
3027
3028
3029
3030template <int dim, int spacedim>
3031void
3032DoFHandler<dim, spacedim>::pre_distributed_transfer_action()
3033{
3034#ifndef DEAL_II_WITH_P4EST
3035 Assert(false,
3036 ExcMessage(
3037 "You are attempting to use a functionality that is only available "
3038 "if deal.II was configured to use p4est, but cmake did not find a "
3039 "valid p4est library."));
3040#else
3041 // the implementation below requires a p:d:T currently
3042 Assert(
3043 (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
3044 &this->get_triangulation()) != nullptr),
3045 ExcNotImplemented());
3046
3047 Assert(active_fe_index_transfer == nullptr, ExcInternalError());
3048
3049 active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3050
3051 // If we work on a p::d::Triangulation, we have to transfer all
3052 // active FE indices since ownership of cells may change. We will
3053 // use our p::d::CellDataTransfer member to achieve this. Further,
3054 // we prepare the values in such a way that they will correspond to
3055 // the active FE indices on the new mesh.
3056
3057 // Gather all current future FE indices.
3058 active_fe_index_transfer->active_fe_indices.resize(
3059 get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
3060
3061 for (const auto &cell : active_cell_iterators())
3062 if (cell->is_locally_owned())
3063 active_fe_index_transfer->active_fe_indices[cell->active_cell_index()] =
3064 cell->future_fe_index();
3065
3066 // Create transfer object and attach to it.
3067 const auto *distributed_tria =
3068 dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
3069 &this->get_triangulation());
3070
3071 active_fe_index_transfer->cell_data_transfer = std::make_unique<
3072 parallel::distributed::
3073 CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3074 *distributed_tria,
3075 /*transfer_variable_size_data=*/false,
3076 /*refinement_strategy=*/
3077 &::AdaptationStrategies::Refinement::
3078 preserve<dim, spacedim, unsigned int>,
3079 /*coarsening_strategy=*/
3080 [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3081 const std::vector<unsigned int> &children_fe_indices)
3082 -> unsigned int {
3083 return ::internal::hp::DoFHandlerImplementation::Implementation::
3084 determine_fe_from_children<dim, spacedim>(parent,
3085 children_fe_indices,
3086 fe_collection);
3087 });
3088
3089 active_fe_index_transfer->cell_data_transfer
3090 ->prepare_for_coarsening_and_refinement(
3091 active_fe_index_transfer->active_fe_indices);
3092#endif
3093}
3094
3095
3096
3097template <int dim, int spacedim>
3098void
3099DoFHandler<dim, spacedim>::post_transfer_action()
3100{
3101 update_active_fe_table();
3102
3103 Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3104
3105 ::internal::hp::DoFHandlerImplementation::Implementation::
3106 distribute_fe_indices_on_refined_cells(*this);
3107
3108 // We have to distribute the information about active FE indices
3109 // of all cells (including the artificial ones) on all processors,
3110 // if a parallel::shared::Triangulation has been used.
3111 ::internal::hp::DoFHandlerImplementation::Implementation::
3112 communicate_active_fe_indices(*this);
3113
3114 // Free memory.
3115 this->active_fe_index_transfer.reset();
3116}
3117
3118
3119
3120template <int dim, int spacedim>
3121void
3122DoFHandler<dim, spacedim>::post_distributed_transfer_action()
3123{
3124#ifndef DEAL_II_WITH_P4EST
3125 Assert(false, ExcInternalError());
3126#else
3127 update_active_fe_table();
3128
3129 Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3130
3131 // Unpack active FE indices.
3132 this->active_fe_index_transfer->active_fe_indices.resize(
3133 this->get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
3134 this->active_fe_index_transfer->cell_data_transfer->unpack(
3135 this->active_fe_index_transfer->active_fe_indices);
3136
3137 // Update all locally owned active FE indices.
3138 this->set_active_fe_indices(
3139 this->active_fe_index_transfer->active_fe_indices);
3140
3141 // Update active FE indices on ghost cells.
3142 ::internal::hp::DoFHandlerImplementation::Implementation::
3143 communicate_active_fe_indices(*this);
3144
3145 // Free memory.
3146 this->active_fe_index_transfer.reset();
3147#endif
3148}
3149
3150
3151
3152template <int dim, int spacedim>
3153void
3154DoFHandler<dim, spacedim>::prepare_for_serialization_of_active_fe_indices()
3155{
3156#ifndef DEAL_II_WITH_P4EST
3157 Assert(false,
3158 ExcMessage(
3159 "You are attempting to use a functionality that is only available "
3160 "if deal.II was configured to use p4est, but cmake did not find a "
3161 "valid p4est library."));
3162#else
3163 // the implementation below requires a p:d:T currently
3164 Assert(
3165 (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
3166 &this->get_triangulation()) != nullptr),
3167 ExcNotImplemented());
3168
3169 Assert(active_fe_index_transfer == nullptr, ExcInternalError());
3170
3171 active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3172
3173 // Create transfer object and attach to it.
3174 const auto *distributed_tria =
3175 dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
3176 &this->get_triangulation());
3177
3178 active_fe_index_transfer->cell_data_transfer = std::make_unique<
3179 parallel::distributed::
3180 CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3181 *distributed_tria,
3182 /*transfer_variable_size_data=*/false,
3183 /*refinement_strategy=*/
3184 &::AdaptationStrategies::Refinement::
3185 preserve<dim, spacedim, unsigned int>,
3186 /*coarsening_strategy=*/
3187 [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3188 const std::vector<unsigned int> &children_fe_indices)
3189 -> unsigned int {
3190 return ::internal::hp::DoFHandlerImplementation::Implementation::
3191 determine_fe_from_children<dim, spacedim>(parent,
3192 children_fe_indices,
3193 fe_collection);
3194 });
3195
3196 // If we work on a p::d::Triangulation, we have to transfer all
3197 // active FE indices since ownership of cells may change.
3198
3199 // Gather all current active FE indices
3200 get_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3201
3202 // Attach to transfer object
3203 active_fe_index_transfer->cell_data_transfer->prepare_for_serialization(
3204 active_fe_index_transfer->active_fe_indices);
3205#endif
3206}
3207
3208
3209
3210template <int dim, int spacedim>
3211void
3212DoFHandler<dim, spacedim>::deserialize_active_fe_indices()
3213{
3214#ifndef DEAL_II_WITH_P4EST
3215 Assert(false,
3216 ExcMessage(
3217 "You are attempting to use a functionality that is only available "
3218 "if deal.II was configured to use p4est, but cmake did not find a "
3219 "valid p4est library."));
3220#else
3221 // the implementation below requires a p:d:T currently
3222 Assert(
3223 (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
3224 &this->get_triangulation()) != nullptr),
3225 ExcNotImplemented());
3226
3227 Assert(active_fe_index_transfer == nullptr, ExcInternalError());
3228
3229 active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3230
3231 // Create transfer object and attach to it.
3232 const auto *distributed_tria =
3233 dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
3234 &this->get_triangulation());
3235
3236 active_fe_index_transfer->cell_data_transfer = std::make_unique<
3237 parallel::distributed::
3238 CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3239 *distributed_tria,
3240 /*transfer_variable_size_data=*/false,
3241 /*refinement_strategy=*/
3242 &::AdaptationStrategies::Refinement::
3243 preserve<dim, spacedim, unsigned int>,
3244 /*coarsening_strategy=*/
3245 [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3246 const std::vector<unsigned int> &children_fe_indices)
3247 -> unsigned int {
3248 return ::internal::hp::DoFHandlerImplementation::Implementation::
3249 determine_fe_from_children<dim, spacedim>(parent,
3250 children_fe_indices,
3251 fe_collection);
3252 });
3253
3254 // Unpack active FE indices.
3255 active_fe_index_transfer->active_fe_indices.resize(
3256 get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
3257 active_fe_index_transfer->cell_data_transfer->deserialize(
3258 active_fe_index_transfer->active_fe_indices);
3259
3260 // Update all locally owned active FE indices.
3261 set_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3262
3263 // Update active FE indices on ghost cells.
3264 ::internal::hp::DoFHandlerImplementation::Implementation::
3265 communicate_active_fe_indices(*this);
3266
3267 // Free memory.
3268 active_fe_index_transfer.reset();
3269#endif
3270}
3271
3272
3273
3274template <int dim, int spacedim>
3275DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs()
3276 : coarsest_level(numbers::invalid_unsigned_int)
3277 , finest_level(0)
3278{}
3279
3280
3281
3282template <int dim, int spacedim>
3283void
3284DoFHandler<dim, spacedim>::MGVertexDoFs::init(
3285 const unsigned int cl,
3286 const unsigned int fl,
3287 const unsigned int dofs_per_vertex)
3288{
3289 coarsest_level = cl;
3290 finest_level = fl;
3291
3292 if (coarsest_level <= finest_level)
3293 {
3294 const unsigned int n_levels = finest_level - coarsest_level + 1;
3295 const unsigned int n_indices = n_levels * dofs_per_vertex;
3296
3297 indices = std::make_unique<types::global_dof_index[]>(n_indices);
3298 std::fill(indices.get(),
3299 indices.get() + n_indices,
3300 numbers::invalid_dof_index);
3301 }
3302 else
3303 indices.reset();
3304}
3305
3306
3307
3308template <int dim, int spacedim>
3309unsigned int
3310DoFHandler<dim, spacedim>::MGVertexDoFs::get_coarsest_level() const
3311{
3312 return coarsest_level;
3313}
3314
3315
3316
3317template <int dim, int spacedim>
3318unsigned int
3319DoFHandler<dim, spacedim>::MGVertexDoFs::get_finest_level() const
3320{
3321 return finest_level;
3322}
3323
3324/*-------------- Explicit Instantiations -------------------------------*/
3325#include "dof_handler.inst"
3326
3327
3328
3329DEAL_II_NAMESPACE_CLOSE
DoFHandler::MGVertexDoFs::init
void init(const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex)
DoFHandler::MGVertexDoFs::get_finest_level
unsigned int get_finest_level() const
DoFHandler::MGVertexDoFs::get_coarsest_level
unsigned int get_coarsest_level() const
DoFHandler::end
cell_iterator end() const
DoFHandler::pre_distributed_transfer_action
void pre_distributed_transfer_action()
DoFHandler::post_transfer_action
void post_transfer_action()
DoFHandler::memory_consumption
virtual std::size_t memory_consumption() const
DoFHandler::hp_cell_active_fe_indices
std::vector< std::vector< active_fe_index_type > > hp_cell_active_fe_indices
Definition: dof_handler.h:1564
DoFHandler::max_couplings_between_dofs
unsigned int max_couplings_between_dofs() const
DoFHandler::mg_levels
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1584
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
DoFHandler::active_fe_index_transfer
std::unique_ptr< ActiveFEIndexTransfer > active_fe_index_transfer
Definition: dof_handler.h:1596
DoFHandler::fe_collection
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1488
DoFHandler::create_active_fe_table
void create_active_fe_table()
DoFHandler::tria
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1481
DoFHandler::set_active_fe_indices
void set_active_fe_indices(const std::vector< unsigned int > &active_fe_indices)
DoFHandler::end_mg
level_cell_iterator end_mg() const
DoFHandler::renumber_dofs
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
DoFHandler::begin_mg
level_cell_iterator begin_mg(const unsigned int level=0) const
DoFHandler::distribute_dofs
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
DoFHandler::n_boundary_dofs
types::global_dof_index n_boundary_dofs() const
DoFHandler::clear_mg_space
void clear_mg_space()
DoFHandler::connect_to_triangulation_signals
void connect_to_triangulation_signals()
DoFHandler::cell_dof_cache_indices
std::vector< std::vector< types::global_dof_index > > cell_dof_cache_indices
Definition: dof_handler.h:1519
DoFHandler::mg_vertex_dofs
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1577
DoFHandler::hp_cell_future_fe_indices
std::vector< std::vector< active_fe_index_type > > hp_cell_future_fe_indices
Definition: dof_handler.h:1571
DoFHandler::post_distributed_transfer_action
void post_distributed_transfer_action()
DoFHandler::object_dof_indices
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1533
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::pre_transfer_action
void pre_transfer_action()
DoFHandler::clear_space
void clear_space()
DoFHandler::hp_object_fe_indices
std::array< std::vector< active_fe_index_type >, dim+1 > hp_object_fe_indices
Definition: dof_handler.h:1552
DoFHandler::reinit
void reinit(const Triangulation< dim, spacedim > &tria)
DoFHandler::object_dof_ptr
std::vector< std::array< std::vector< offset_type >, dim+1 > > object_dof_ptr
Definition: dof_handler.h:1544
DoFHandler::mg_faces
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1590
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
DoFHandler::distribute_mg_dofs
void distribute_mg_dofs()
DoFHandler::end_active
active_cell_iterator end_active(const unsigned int level) const
DoFHandler::set_fe
void set_fe(const FiniteElement< dim, spacedim > &fe)
DoFHandler::hp_capability_enabled
bool hp_capability_enabled
Definition: dof_handler.h:1475
DoFHandler::initialize_local_block_info
void initialize_local_block_info()
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
DoFHandler::update_active_fe_table
void update_active_fe_table()
DoFHandler::cell_dof_cache_ptr
std::vector< std::vector< offset_type > > cell_dof_cache_ptr
Definition: dof_handler.h:1525
DoFHandler::level_cell_iterator
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:502
DoFHandler::prepare_for_serialization_of_active_fe_indices
void prepare_for_serialization_of_active_fe_indices()
DoFHandler::begin
cell_iterator begin(const unsigned int level=0) const
DoFHandler::hp_object_fe_ptr
std::array< std::vector< offset_type >, dim+1 > hp_object_fe_ptr
Definition: dof_handler.h:1557
DoFHandler::clear
void clear()
DoFHandler::get_active_fe_indices
void get_active_fe_indices(std::vector< unsigned int > &active_fe_indices) const
DoFHandler::max_couplings_between_boundary_dofs
unsigned int max_couplings_between_boundary_dofs() const
DoFHandler::setup_policy
void setup_policy()
DoFHandler::initialize
void initialize(const Triangulation< dim, spacedim > &tria, const FiniteElement< dim, spacedim > &fe)
DoFHandler::active_fe_index_type
unsigned short int active_fe_index_type
Definition: dof_handler.h:529
DoFHandler::deserialize_active_fe_indices
void deserialize_active_fe_indices()
DoFHandler::~DoFHandler
virtual ~DoFHandler() override
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_dofs_per_line
unsigned int n_dofs_per_line() const
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1767
TriaRawIterator::state
IteratorState::IteratorStates state() const
Definition: tria_iterator.h:1058
Triangulation::get_mesh_smoothing
virtual const MeshSmoothing & get_mesh_smoothing() const
Triangulation::begin
cell_iterator begin(const unsigned int level=0) const
Triangulation::n_raw_lines
unsigned int n_raw_lines() const
Triangulation::n_raw_faces
unsigned int n_raw_faces() const
Triangulation::n_levels
unsigned int n_levels() const
Triangulation::end
cell_iterator end() const
Triangulation::n_raw_cells
unsigned int n_raw_cells(const unsigned int level) const
Triangulation::max_adjacent_cells
unsigned int max_adjacent_cells() const
Triangulation::vertex_used
bool vertex_used(const unsigned int index) const
Triangulation::n_raw_quads
unsigned int n_raw_quads() const
Triangulation::n_cells
unsigned int n_cells() const
Triangulation::signals
Signals signals
Definition: tria.h:2448
Triangulation::n_vertices
unsigned int n_vertices() const
Triangulation::n_raw_hexs
unsigned int n_raw_hexs(const unsigned int level) const
hp::Collection::size
unsigned int size() const
Definition: collection.h:264
hp::FECollection::max_dofs_per_line
unsigned int max_dofs_per_line() const
Definition: fe_collection.h:978
hp::FECollection::max_dofs_per_hex
unsigned int max_dofs_per_hex() const
Definition: fe_collection.h:1008
hp::FECollection::max_dofs_per_vertex
unsigned int max_dofs_per_vertex() const
Definition: fe_collection.h:963
hp::FECollection::max_dofs_per_quad
unsigned int max_dofs_per_quad() const
Definition: fe_collection.h:993
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
dof_handler_policy.h
dof_handler.h
filtered_iterator.h
fully_distributed_tria.h
geometry_info.h
level
unsigned int level
Definition: grid_out.cc:4606
grid_tools.h
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
DoFHandler::mg_cell_iterators
IteratorRange< level_cell_iterator > mg_cell_iterators() const
DoFHandler::cell_iterators_on_level
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
DoFHandler::cell_iterators
IteratorRange< cell_iterator > cell_iterators() const
DoFHandler::active_cell_iterators_on_level
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
Triangulation::cell_iterators_on_level
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Triangulation::cell_iterators
IteratorRange< cell_iterator > cell_iterators() const
DoFHandler::mg_cell_iterators_on_level
IteratorRange< level_cell_iterator > mg_cell_iterators_on_level(const unsigned int level) const
Threads::TaskGroup::join_all
void join_all() const
Definition: thread_management.h:1726
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
internal::hp::DoFHandlerImplementation::ExcNoDominatedFiniteElementOnChildren
static ::ExceptionBase & ExcNoDominatedFiniteElementOnChildren()
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)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
GridTools::exchange_cell_data_to_ghosts
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition: thread_management.h:1530
memory_consumption.h
DoFTools::get_active_fe_indices
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1371
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:38
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:266
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:67
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:53
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
hp
Definition: hp.h:118
internal::TriangulationImplementation::n_active_cells
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13741
internal::hp::DoFHandlerImplementation::dominated_future_fe_on_children
unsigned int dominated_future_fe_on_children(const typename DoFHandler< dim, spacedim >::cell_iterator &parent)
Definition: dof_handler.cc:1830
internal::hp::DoFHandlerImplementation::communicate_future_fe_indices
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1817
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
internal::policy_to_string
std::string policy_to_string(const ::internal::DoFHandlerImplementation::Policy::PolicyBase< dim, spacedim > &policy)
Definition: dof_handler.cc:57
numbers::internal_face_boundary_id
const types::boundary_id internal_face_boundary_id
Definition: types.h:260
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:216
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5605
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
shared_tria.h
Triangulation::Signals::post_distributed_load
boost::signals2::signal< void()> post_distributed_load
Definition: tria.h:2442
Triangulation::Signals::post_distributed_refinement
boost::signals2::signal< void()> post_distributed_refinement
Definition: tria.h:2398
Triangulation::Signals::pre_refinement
boost::signals2::signal< void()> pre_refinement
Definition: tria.h:2091
Triangulation::Signals::create
boost::signals2::signal< void()> create
Definition: tria.h:2082
Triangulation::Signals::clear
boost::signals2::signal< void()> clear
Definition: tria.h:2161
Triangulation::Signals::post_refinement
boost::signals2::signal< void()> post_refinement
Definition: tria.h:2098
Triangulation::Signals::pre_partition
boost::signals2::signal< void()> pre_partition
Definition: tria.h:2106
Triangulation::Signals::pre_distributed_repartition
boost::signals2::signal< void()> pre_distributed_repartition
Definition: tria.h:2405
Triangulation::Signals::post_p4est_refinement
boost::signals2::signal< void()> post_p4est_refinement
Definition: tria.h:2388
Triangulation::Signals::post_distributed_repartition
boost::signals2::signal< void()> post_distributed_repartition
Definition: tria.h:2412
Triangulation::Signals::post_distributed_save
boost::signals2::signal< void()> post_distributed_save
Definition: tria.h:2427
internal::DoFHandlerImplementation::Implementation::reserve_space_mg
static void reserve_space_mg(DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:596
internal::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:399
internal::DoFHandlerImplementation::Implementation::reserve_subentities
static void reserve_subentities(DoFHandler< dim, spacedim > &dof_handler, const unsigned int structdim, const unsigned int n_raw_entities, const T &cell_process)
Definition: dof_handler.cc:340
internal::DoFHandlerImplementation::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:112
internal::DoFHandlerImplementation::Implementation::reset_to_empty_objects
static void reset_to_empty_objects(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:267
internal::DoFHandlerImplementation::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:231
internal::DoFHandlerImplementation::Implementation::reserve_space_mg
static void reserve_space_mg(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:449
internal::DoFHandlerImplementation::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:102
internal::DoFHandlerImplementation::Implementation::reserve_space_mg
static void reserve_space_mg(DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:519
internal::DoFHandlerImplementation::Implementation::reserve_cells
static void reserve_cells(DoFHandler< dim, spacedim > &dof_handler, const unsigned int n_inner_dofs_per_cell)
Definition: dof_handler.cc:291
internal::hp::DoFHandlerImplementation::Implementation::collect_fe_indices_on_cells_to_be_refined
static void collect_fe_indices_on_cells_to_be_refined(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1601
internal::hp::DoFHandlerImplementation::Implementation::determine_fe_from_children
static unsigned int determine_fe_from_children(const typename Triangulation< dim, spacedim >::cell_iterator &, const std::vector< unsigned int > &children_fe_indices, const ::hp::FECollection< dim, spacedim > &fe_collection)
Definition: dof_handler.cc:1730
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_cells
static void reserve_space_cells(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:851
internal::hp::DoFHandlerImplementation::Implementation::dominated_future_fe_on_children
static unsigned int dominated_future_fe_on_children(const typename DoFHandler< dim, spacedim >::cell_iterator &parent)
Definition: dof_handler.cc:1761
internal::hp::DoFHandlerImplementation::Implementation::distribute_fe_indices_on_refined_cells
static void distribute_fe_indices_on_refined_cells(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1672
internal::hp::DoFHandlerImplementation::Implementation::ensure_absence_of_future_fe_indices
static void ensure_absence_of_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:698
internal::hp::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(::DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:1184
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_faces
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:910
internal::hp::DoFHandlerImplementation::Implementation::communicate_future_fe_indices
static void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1506
internal::hp::DoFHandlerImplementation::Implementation::communicate_active_fe_indices
static void communicate_active_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1394
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_vertices
static void reserve_space_vertices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:718
internal::hp::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(::DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:1209
internal::hp::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(::DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:1236
tria_accessor.h
tria
const ::Triangulation< dim, spacedim > & tria
Definition: tria_description.cc:717
tria_iterator.h
tria_levels.h
Generated by doxygen 1.9.5