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