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 #ifndef dealii_dof_handler_h
17 #define dealii_dof_handler_h
18
19
20
21 #include <deal.II/base/config.h>
22
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/function.h>
25 #include <deal.II/base/index_set.h>
26 #include <deal.II/base/iterator_range.h>
27 #include <deal.II/base/smartpointer.h>
28
29 #include <deal.II/distributed/tria_base.h>
30
31 #include <deal.II/dofs/block_info.h>
32 #include <deal.II/dofs/dof_accessor.h>
33 #include <deal.II/dofs/dof_faces.h>
34 #include <deal.II/dofs/dof_iterator_selector.h>
35 #include <deal.II/dofs/dof_levels.h>
36 #include <deal.II/dofs/number_cache.h>
37
38 #include <deal.II/hp/fe_collection.h>
39
40 #include <boost/serialization/split_member.hpp>
41
42 #include <map>
43 #include <memory>
44 #include <set>
45 #include <vector>
46
47 DEAL_II_NAMESPACE_OPEN
48
49 // Forward declarations
50 #ifndef DOXYGEN
51 template <int dim, int spacedim>
52 class FiniteElement;
53 template <int dim, int spacedim>
54 class Triangulation;
55
56 namespace internal
57 {
58 namespace DoFHandlerImplementation
59 {
60 struct Implementation;
61
62 namespace Policy
63 {
64 template <int dim, int spacedim>
65 class PolicyBase;
66 struct Implementation;
67 } // namespace Policy
68 } // namespace DoFHandlerImplementation
69
70 namespace DoFAccessorImplementation
71 {
72 struct Implementation;
73 }
74
75 namespace DoFCellAccessorImplementation
76 {
77 struct Implementation;
78 }
79
80 namespace hp
81 {
82 namespace DoFHandlerImplementation
83 {
84 struct Implementation;
85 }
86 } // namespace hp
87 } // namespace internal
88
89 namespace parallel
90 {
91 namespace distributed
92 {
93 template <int dim, int spacedim, typename VectorType>
94 class CellDataTransfer;
95 }
96 } // namespace parallel
97 #endif
98
99 /**
100 * Given a triangulation and a description of a finite element, this
101 * class enumerates degrees of freedom on all vertices, edges, faces,
102 * and cells of the triangulation. As a result, it also provides a
103 * <i>basis</i> for a discrete space $V_h$ whose elements are finite
104 * element functions defined on each cell by a FiniteElement object.
105 * This class satisfies the
106 * @ref ConceptMeshType "MeshType concept"
107 * requirements.
108 *
109 * It is first used in the step-2 tutorial program.
110 *
111 * For each vertex, line, quad, etc, this class stores a list of the indices
112 * of degrees of freedom living on this object. These indices refer to the
113 * unconstrained degrees of freedom, i.e. constrained degrees of freedom are
114 * numbered in the same way as unconstrained ones, and are only later
115 * eliminated. This leads to the fact that indices in global vectors and
116 * matrices also refer to all degrees of freedom and some kind of condensation
117 * is needed to restrict the systems of equations to the unconstrained degrees
118 * of freedom only. The actual layout of storage of the indices is described
119 * in the dealii::internal::DoFHandlerImplementation::DoFLevel class
120 * documentation.
121 *
122 * The class offers iterators to traverse all cells, in much the same way as
123 * the Triangulation class does. Using the begin() and end() functions (and
124 * companions, like begin_active()), one can obtain iterators to walk over
125 * cells, and query the degree of freedom structures as well as the
126 * triangulation data. These iterators are built on top of those of the
127 * Triangulation class, but offer the additional information on degrees of
128 * freedom functionality compared to pure triangulation iterators. The order
129 * in which dof iterators are presented by the <tt>++</tt> and <tt>\--</tt>
130 * operators is the same as that for the corresponding iterators traversing
131 * the triangulation on which this DoFHandler is constructed.
132 *
133 * The <tt>spacedim</tt> parameter has to be used if one wants to solve
134 * problems on surfaces. If not specified, this parameter takes the default
135 * value <tt>=dim</tt> implying that we want to solve problems in a domain
136 * whose dimension equals the dimension of the space in which it is embedded.
137 *
138 *
139 * <h3>Distribution of indices for degrees of freedom</h3>
140 *
141 * The degrees of freedom (`dofs') are distributed on the given triangulation
142 * by the function distribute_dofs(). It gets passed a finite element object
143 * describing how many degrees of freedom are located on vertices, lines, etc.
144 * It traverses the triangulation cell by cell and numbers the dofs of that
145 * cell if not yet numbered. For non-multigrid algorithms, only active cells
146 * are considered. Active cells are defined to be those cells which have no
147 * children, i.e. they are the most refined ones.
148 *
149 * Since the triangulation is traversed starting with the cells of the
150 * coarsest active level and going to more refined levels, the lowest numbers
151 * for dofs are given to the largest cells as well as their bounding lines and
152 * vertices, with the dofs of more refined cells getting higher numbers.
153 *
154 * This numbering implies very large bandwidths of the resulting matrices and
155 * is thus vastly suboptimal for some solution algorithms. For this reason,
156 * the DoFRenumbering class offers several algorithms to reorder the dof
157 * numbering according. See there for a discussion of the implemented
158 * algorithms.
159 *
160 *
161 * <h3>Interaction with distributed meshes</h3>
162 *
163 * Upon construction, this class takes a reference to a triangulation object.
164 * In most cases, this will be a reference to an object of type Triangulation,
165 * i.e. the class that represents triangulations that entirely reside on a
166 * single processor. However, it can also be of type
167 * parallel::distributed::Triangulation (see, for example, step-32, step-40
168 * and in particular the
169 * @ref distributed
170 * module) in which case the DoFHandler object will proceed to only manage
171 * degrees of freedom on locally owned and ghost cells. This process is
172 * entirely transparent to the used.
173 *
174 *
175 * <h3>User defined renumbering schemes</h3>
176 *
177 * The DoFRenumbering class offers a number of renumbering schemes like the
178 * Cuthill-McKee scheme. Basically, the function sets up an array in which for
179 * each degree of freedom we store the new index this DoF should have after
180 * renumbering. Using this array, the renumber_dofs() function of the present
181 * class is called, which actually performs the change from old DoF indices to
182 * the ones given in the array. In some cases, however, a user may want to
183 * compute their own renumbering order; in this case, one can allocate an array
184 * with one element per degree of freedom and fill it with the number that the
185 * respective degree of freedom shall be assigned. This number may, for
186 * example, be obtained by sorting the support points of the degrees of
187 * freedom in downwind direction. Then call the
188 * <tt>renumber_dofs(vector<types::global_dof_index>)</tt> function with the
189 * array, which converts old into new degree of freedom indices.
190 *
191 *
192 * <h3>Serializing (loading or storing) DoFHandler objects</h3>
193 *
194 * Like many other classes in deal.II, the DoFHandler class can stream its
195 * contents to an archive using BOOST's serialization facilities. The data so
196 * stored can later be retrieved again from the archive to restore the
197 * contents of this object. This facility is frequently used to save the state
198 * of a program to disk for possible later resurrection, often in the context
199 * of checkpoint/restart strategies for long running computations or on
200 * computers that aren't very reliable (e.g. on very large clusters where
201 * individual nodes occasionally fail and then bring down an entire MPI job).
202 *
203 * The model for doing so is similar for the DoFHandler class as it is for the
204 * Triangulation class (see the section in the general documentation of that
205 * class). In particular, the load() function does not exactly restore the
206 * same state as was stored previously using the save() function. Rather, the
207 * function assumes that you load data into a DoFHandler object that is
208 * already associated with a triangulation that has a content that matches the
209 * one that was used when the data was saved. Likewise, the load() function
210 * assumes that the current object is already associated with a finite element
211 * object that matches the one that was associated with it when data was
212 * saved; the latter can be achieved by calling DoFHandler::distribute_dofs()
213 * using the same kind of finite element before re-loading data from the
214 * serialization archive.
215 *
216 * @ingroup dofs
217 */
218 template <int dim, int spacedim = dim>
219 class DoFHandler : public Subscriptor
220 {
221 using ActiveSelector =
222 dealii::internal::DoFHandlerImplementation::Iterators<dim, spacedim, false>;
223 using LevelSelector =
224 dealii::internal::DoFHandlerImplementation::Iterators<dim, spacedim, true>;
225
226 public:
227 /**
228 * An alias that is used to identify cell iterators in DoFHandler objects.
229 * The concept of iterators is discussed at length in the
230 * @ref Iterators "iterators documentation module".
231 *
232 * The current alias works, in essence, like the corresponding
233 * Triangulation::cell_accessor alias. However, it also makes available
234 * the member functions of DoFCellAccessor, in addition to the ones
235 * already available through the CellAccessor class.
236 *
237 * @ingroup Iterators
238 */
239 using cell_accessor = typename ActiveSelector::CellAccessor;
240
241 /**
242 * An alias that is used to identify iterators that point to faces.
243 * The concept of iterators is discussed at length in the
244 * @ref Iterators "iterators documentation module".
245 *
246 * The current alias works, in essence, like the corresponding
247 * Triangulation::face_accessor alias. However, it also makes available
248 * the member functions of DoFAccessor, in addition to the ones
249 * already available through the TriaAccessor class.
250 *
251 * @ingroup Iterators
252 */
253 using face_accessor = typename ActiveSelector::FaceAccessor;
254
255 /**
256 * An alias that defines an iterator over the (one-dimensional) lines
257 * of a mesh. In one-dimensional meshes, these are the cells of the mesh,
258 * whereas in two-dimensional meshes the lines are the faces of cells.
259 *
260 * @ingroup Iterators
261 */
262 using line_iterator = typename ActiveSelector::line_iterator;
263
264 /**
265 * An alias that allows iterating over the <i>active</i> lines, i.e.,
266 * that subset of lines that have no children. In one-dimensional meshes,
267 * these are the cells of the mesh, whereas in two-dimensional
268 * meshes the lines are the faces of cells.
269 *
270 * In two- or three-dimensional meshes, lines without children (i.e.,
271 * the active lines) are part of at least one active cell. Each such line may
272 * additionally be a child of a line of a coarser cell adjacent to a cell
273 * that is active. (This coarser neighbor would then also be active.)
274 *
275 * @ingroup Iterators
276 */
277 using active_line_iterator = typename ActiveSelector::active_line_iterator;
278
279 /**
280 * An alias that defines an iterator over the (two-dimensional) quads
281 * of a mesh. In two-dimensional meshes, these are the cells of the mesh,
282 * whereas in three-dimensional meshes the quads are the faces of cells.
283 *
284 * @ingroup Iterators
285 */
286 using quad_iterator = typename ActiveSelector::quad_iterator;
287
288 /**
289 * An alias that allows iterating over the <i>active</i> quads, i.e.,
290 * that subset of quads that have no children. In two-dimensional meshes,
291 * these are the cells of the mesh, whereas in three-dimensional
292 * meshes the quads are the faces of cells.
293 *
294 * In three-dimensional meshes, quads without children (i.e.,
295 * the active quads) are faces of at least one active cell. Each such quad may
296 * additionally be a child of a quad face of a coarser cell adjacent to a cell
297 * that is active. (This coarser neighbor would then also be active.)
298 *
299 * @ingroup Iterators
300 */
301 using active_quad_iterator = typename ActiveSelector::active_quad_iterator;
302
303 /**
304 * An alias that defines an iterator over the (three-dimensional) hexes
305 * of a mesh. This iterator only makes sense in three-dimensional meshes,
306 * where hexes are the cells of the mesh.
307 *
308 * @ingroup Iterators
309 */
310 using hex_iterator = typename ActiveSelector::hex_iterator;
311
312 /**
313 * An alias that allows iterating over the <i>active</i> hexes of a mesh.
314 * This iterator only makes sense in three-dimensional meshes,
315 * where hexes are the cells of the mesh. Consequently, in these
316 * three-dimensional meshes, this iterator is equivalent to the
317 * @p active_cell_iterator alias.
318 *
319 * @ingroup Iterators
320 */
321 using active_hex_iterator = typename ActiveSelector::active_hex_iterator;
322
323 /**
324 * An alias that is used to identify
325 * @ref GlossActive "active cell iterators".
326 * The concept of iterators is discussed at length in the
327 * @ref Iterators "iterators documentation module".
328 *
329 * The current alias identifies active cells in a DoFHandler object. While
330 * the actual data type of the alias is hidden behind a few layers of
331 * (unfortunately necessary) indirections, it is in essence
332 * TriaActiveIterator<DoFCellAccessor>. The TriaActiveIterator class works
333 * like a pointer to active objects that when you dereference it yields an
334 * object of type DoFCellAccessor. DoFCellAccessor is a class that
335 * identifies properties that are specific to cells in a DoFHandler, but it
336 * is derived (and consequently inherits) from both DoFAccessor,
337 * TriaCellAccessor and TriaAccessor that describe what you can ask of more
338 * general objects (lines, faces, as well as cells) in a triangulation and
339 * DoFHandler objects.
340 *
341 * @ingroup Iterators
342 */
343 using active_cell_iterator = typename ActiveSelector::active_cell_iterator;
344
345 /**
346 * An alias that is used to identify cell iterators. The concept of
347 * iterators is discussed at length in the
348 * @ref Iterators "iterators documentation module".
349 *
350 * The current alias identifies cells in a DoFHandler object. Some of
351 * these cells may in fact be active (see
352 * @ref GlossActive "active cell iterators")
353 * in which case they can in fact be asked for the degrees of freedom that
354 * live on them. On the other hand, if the cell is not active, any such
355 * query will result in an error. Note that this is what distinguishes this
356 * alias from the level_cell_iterator alias.
357 *
358 * While the actual data type of the alias is hidden behind a few layers
359 * of (unfortunately necessary) indirections, it is in essence
360 * TriaIterator<DoFCellAccessor>. The TriaIterator class works like a
361 * pointer to objects that when you dereference it yields an object of type
362 * DoFCellAccessor. DoFCellAccessor is a class that identifies properties
363 * that are specific to cells in a DoFHandler, but it is derived (and
364 * consequently inherits) from both DoFAccessor, TriaCellAccessor and
365 * TriaAccessor that describe what you can ask of more general objects
366 * (lines, faces, as well as cells) in a triangulation and DoFHandler
367 * objects.
368 *
369 * @ingroup Iterators
370 */
371 using cell_iterator = typename ActiveSelector::cell_iterator;
372
373 /**
374 * An alias that is used to identify iterators that point to faces.
375 * The concept of iterators is discussed at length in the
376 * @ref Iterators "iterators documentation module".
377 *
378 * While the actual data type of the alias is hidden behind a few layers
379 * of (unfortunately necessary) indirections, it is in essence
380 * TriaIterator<DoFAccessor>. The
381 * TriaIterator class works like a pointer to objects that when
382 * you dereference it yields an object of type DoFAccessor. DoFAccessor,
383 * in turn, is a class that can be used to query DoF indices on faces,
384 * but it is also derived from TriaAccessor and consequently can be used
385 * to query geometric properties such as vertices of faces, their area, etc.
386 *
387 * @ingroup Iterators
388 */
389 using face_iterator = typename ActiveSelector::face_iterator;
390
391 /**
392 * An alias that is used to identify iterators that point to active faces,
393 * i.e., to faces that have no children. Active faces must be faces of at
394 * least one active cell.
395 *
396 * Other than the "active" qualification, this alias is identical to the
397 * @p face_iterator alias. In particular, dereferencing either yields
398 * the same kind of object.
399 *
400 * @ingroup Iterators
401 */
402 using active_face_iterator = typename ActiveSelector::active_face_iterator;
403
404 using level_cell_accessor = typename LevelSelector::CellAccessor;
405 using level_face_accessor = typename LevelSelector::FaceAccessor;
406
407 using level_cell_iterator = typename LevelSelector::cell_iterator;
408 using level_face_iterator = typename LevelSelector::face_iterator;
409
410
411 /**
412 * Make the dimension available in function templates.
413 */
414 static const unsigned int dimension = dim;
415
416 /**
417 * Make the space dimension available in function templates.
418 */
419 static const unsigned int space_dimension = spacedim;
420
421 /**
422 * Boolean indicating whether or not the current DoFHander has hp
423 * capabilities.
424 */
425 const bool hp_capability_enabled;
426
427 /**
428 * The default index of the finite element to be used on a given cell.
429 */
430 static const unsigned int default_fe_index = 0;
431
432 /**
433 * Invalid index of the finite element to be used on a given cell.
434 */
435 static const unsigned int invalid_fe_index = numbers::invalid_unsigned_int;
436
437 /**
438 * The type in which we store the active FE index.
439 */
440 using active_fe_index_type = unsigned short int;
441
442 /**
443 * The type in which we store the offsets in the CRS data structures.
444 */
445 using offset_type = unsigned int;
446
447 /**
448 * Invalid active_fe_index which will be used as a default value to determine
449 * whether a future_fe_index has been set or not.
450 */
451 static const active_fe_index_type invalid_active_fe_index =
452 static_cast<active_fe_index_type>(-1);
453
454 /**
455 * Standard constructor, not initializing any data. After constructing an
456 * object with this constructor, use initialize() to make a valid
457 * DoFHandler.
458 */
459 DoFHandler(const bool hp_capability_enabled = false);
460
461 /**
462 * Constructor. Take @p tria as the triangulation to work on.
463 */
464 DoFHandler(const Triangulation<dim, spacedim> &tria,
465 const bool hp_capability_enabled = false);
466
467 /**
468 * Copy constructor. DoFHandler objects are large and expensive.
469 * They should not be copied, in particular not by accident, but
470 * rather deliberately constructed. As a consequence, this constructor
471 * is explicitly removed from the interface of this class.
472 */
473 DoFHandler(const DoFHandler &) = delete;
474
475 /**
476 * Destructor.
477 */
478 virtual ~DoFHandler() override;
479
480 /**
481 * Copy operator. DoFHandler objects are large and expensive.
482 * They should not be copied, in particular not by accident, but
483 * rather deliberately constructed. As a consequence, this operator
484 * is explicitly removed from the interface of this class.
485 */
486 DoFHandler &
487 operator=(const DoFHandler &) = delete;
488
489 /**
490 * Assign a Triangulation and a FiniteElement to the DoFHandler and compute
491 * the distribution of degrees of freedom over the mesh.
492 */
493 void
494 initialize(const Triangulation<dim, spacedim> &tria,
495 const FiniteElement<dim, spacedim> &fe);
496
497 /**
498 * Same as above but taking an hp::FECollection object.
499 */
500 void
501 initialize(const Triangulation<dim, spacedim> & tria,
502 const hp::FECollection<dim, spacedim> &fe);
503
504 /**
505 * Assign a FiniteElement @p fe to this object.
506 *
507 * @note This function makes a copy of the finite element given as
508 * argument, and stores it as a member variable. Consequently, it is
509 * possible to write code such as
510 * @code
511 * dof_handler.set_fe(FE_Q<dim>(2));
512 * @endcode
513 * You can then access the finite element later on by calling
514 * DoFHandler::get_fe(). However, it is often more convenient to
515 * keep a named finite element object as a member variable in your
516 * main class and refer to it directly whenever you need to access
517 * properties of the finite element (such as
518 * FiniteElementData::dofs_per_cell). This is what all tutorial programs do.
519 *
520 * @warning This function only sets a FiniteElement. Degrees of freedom have
521 * either not been distributed yet, or are distributed using a previously set
522 * element. In both cases, accessing degrees of freedom will lead to invalid
523 * results. To restore consistency, call distribute_dofs().
524 */
525 void
526 set_fe(const FiniteElement<dim, spacedim> &fe);
527
528 /**
529 * Same as above but taking an hp::FECollection object.
530 */
531 void
532 set_fe(const hp::FECollection<dim, spacedim> &fe);
533
534 /**
535 * Go through the triangulation and set the active FE indices of all
536 * active cells to the values given in @p active_fe_indices.
537 */
538 void
539 set_active_fe_indices(const std::vector<unsigned int> &active_fe_indices);
540
541 /**
542 * Go through the triangulation and store the active FE indices of all
543 * active cells to the vector @p active_fe_indices. This vector is
544 * resized, if necessary.
545 */
546 void
547 get_active_fe_indices(std::vector<unsigned int> &active_fe_indices) const;
548
549 /**
550 * Go through the triangulation and "distribute" the degrees of
551 * freedom needed for the given finite element. "Distributing"
552 * degrees of freedom involves allocating memory to store the
553 * indices on all entities on which degrees of freedom can be
554 * located (e.g., vertices, edges, faces, etc.) and to then enumerate
555 * all degrees of freedom. In other words, while the mesh and the
556 * finite element object by themselves simply define a finite
557 * element space $V_h$, the process of distributing degrees of
558 * freedom makes sure that there is a basis for this space and that
559 * the shape functions of this basis are enumerated in an indexable,
560 * predictable way.
561 *
562 * The exact order in which degrees of freedom on a mesh are
563 * ordered, i.e., the order in which basis functions of the finite
564 * element space are enumerated, is something that deal.II treats as
565 * an implementation detail. By and large, degrees of freedom are
566 * enumerated in the same order in which we traverse cells, but you
567 * should not rely on any specific numbering. In contrast, if you
568 * want a particular ordering, use the functions in namespace
569 * DoFRenumbering.
570 *
571 * This function is first discussed in the introduction to the
572 * step-2 tutorial program.
573 *
574 * @note This function makes a copy of the finite element given as
575 * argument, and stores it as a member variable, similarly to the above
576 * function set_fe().
577 */
578 void
579 distribute_dofs(const FiniteElement<dim, spacedim> &fe);
580
581 /**
582 * Same as above but taking an hp::FECollection object.
583 */
584 void
585 distribute_dofs(const hp::FECollection<dim, spacedim> &fe);
586
587 /**
588 * Distribute level degrees of freedom on each level for geometric
589 * multigrid. The active DoFs need to be distributed using distribute_dofs()
590 * before calling this function.
591 */
592 void
593 distribute_mg_dofs();
594
595 /**
596 * This function returns whether this DoFHandler has DoFs distributed on
597 * each multigrid level or in other words if distribute_mg_dofs() has been
598 * called.
599 */
600 bool
601 has_level_dofs() const;
602
603 /**
604 * This function returns whether this DoFHandler has active DoFs. This is
605 * equivalent to asking whether (i) distribute_dofs() has been called and
606 * (ii) the finite element for which degrees of freedom have been
607 * distributed actually has degrees of freedom (which is not the case for
608 * FE_Nothing, for example).
609 *
610 * If this object is based on a parallel::distributed::Triangulation, then
611 * the current function returns true if <i>any</i> partition of the parallel
612 * DoFHandler object has any degrees of freedom. In other words, the
613 * function returns true even if the Triangulation does not own any active
614 * cells on the current MPI process, but at least one process owns cells and
615 * at least this one process has any degrees of freedom associated with it.
616 */
617 bool
618 has_active_dofs() const;
619
620 /**
621 * After distribute_dofs() with an FESystem element, the block structure of
622 * global and level vectors is stored in a BlockInfo object accessible with
623 * block_info(). This function initializes the local block structure on each
624 * cell in the same object.
625 */
626 void
627 initialize_local_block_info();
628
629 /**
630 * Clear all data of this object.
631 */
632 void
633 clear();
634
635 /**
636 * Renumber degrees of freedom based on a list of new DoF indices for each
637 * of the degrees of freedom.
638 *
639 * This function is called by the functions in DoFRenumbering function after
640 * computing a new ordering of the degree of freedom indices. However, it
641 * can of course also be called from user code.
642 *
643 * @arg new_number This array must have a size equal to the number of
644 * degrees of freedom owned by the current processor, i.e. the size must be
645 * equal to what n_locally_owned_dofs() returns. If only one processor
646 * participates in storing the current mesh, then this equals the total
647 * number of degrees of freedom, i.e. the result of n_dofs(). The contents
648 * of this array are the new global indices for each freedom listed in the
649 * IndexSet returned by locally_owned_dofs(). In the case of a sequential
650 * mesh this means that the array is a list of new indices for each of the
651 * degrees of freedom on the current mesh. In the case that we have a
652 * parallel::shared::Triangulation or
653 * parallel::distributed::Triangulation underlying this DoFHandler object,
654 * the array is a list of new indices for all the locally owned degrees of
655 * freedom, enumerated in the same order as the currently locally owned
656 * DoFs. In other words, assume that degree of freedom <code>i</code> is
657 * currently locally owned, then
658 * <code>new_numbers[locally_owned_dofs().index_within_set(i)]</code>
659 * returns the new global DoF index of <code>i</code>. Since the IndexSet of
660 * locally_owned_dofs() is complete in the sequential case, the latter
661 * convention for the content of the array reduces to the former in the case
662 * that only one processor participates in the mesh.
663 *
664 * @note While it follows from the above, it may be surprising to know that
665 * the <i>number</i> of locally owned degrees of freedom in a parallel
666 * computation is an invariant
667 * under renumbering, even if the <i>indices</i> associated with these
668 * locally owned degrees of freedom are not. At a fundamental level,
669 * this invariant exists because the <i>decision</i> whether a degree of
670 * freedom is locally owned or not has nothing to do with that
671 * degree of freedom's (old or new) index. Indeed, degrees of freedom
672 * are locally owned if they are on a locally owned cell and not on
673 * an interface between cells where the neighboring cell has a lower
674 * @ref GlossSubdomainId "subdomain id". Since both of these conditions
675 * are independent of the index associated with the DoF, a locally
676 * owned degree of freedom will also be locally owned after renumbering.
677 * On the other hand, properties such as whether the set of indices
678 * of locally owned DoFs forms a contiguous range or not
679 * (i.e., whether the locally_owned_dofs() returns an IndexSet object
680 * for which IndexSet::is_contiguous() returns @p true) are of
681 * course affected by the exact renumbering performed here. For example,
682 * while the initial numbering of DoF indices done in distribute_dofs()
683 * yields a contiguous numbering, the renumbering performed by
684 * DoFRenumbering::component_wise() will, in general, not yield
685 * contiguous locally owned DoF indices.
686 */
687 void
688 renumber_dofs(const std::vector<types::global_dof_index> &new_numbers);
689
690 /**
691 * The same function as above, but renumber the degrees of freedom of a
692 * single level of a multigrid hierarchy.
693 */
694 void
695 renumber_dofs(const unsigned int level,
696 const std::vector<types::global_dof_index> &new_numbers);
697
698 /**
699 * Return the maximum number of degrees of freedom a degree of freedom in
700 * the given triangulation with the given finite element may couple with.
701 * This is the maximum number of entries per line in the system matrix; this
702 * information can therefore be used upon construction of the
703 * SparsityPattern object.
704 *
705 * The returned number is not really the maximum number but an estimate
706 * based on the finite element and the maximum number of cells meeting at a
707 * vertex. The number holds for the constrained matrix as well.
708 *
709 * The determination of the number of couplings can be done by simple
710 * picture drawing. An example can be found in the implementation of this
711 * function.
712 *
713 * @note This function is most often used to determine the maximal row
714 * length for sparsity patterns. Unfortunately, while the estimates returned
715 * by this function are rather accurate in 1d and 2d, they are often
716 * significantly too high in 3d, leading the SparsityPattern class to
717 * allocate much too much memory in some cases. Unless someone comes around
718 * to improving the present function for 3d, there is not very much one can
719 * do about these cases. The typical way to work around this problem is to
720 * use an intermediate compressed sparsity pattern that only allocates
721 * memory on demand. Refer to the step-2 and step-11 example programs on how
722 * to do this. The problem is also discussed in the documentation of the
723 * module on
724 * @ref Sparsity.
725 */
726 unsigned int
727 max_couplings_between_dofs() const;
728
729 /**
730 * Return the number of degrees of freedom located on the boundary another
731 * dof on the boundary can couple with.
732 *
733 * The number is the same as for max_couplings_between_dofs() in one
734 * dimension less.
735 *
736 * @note The same applies to this function as to max_couplings_per_dofs() as
737 * regards the performance of this function. Think about one of the dynamic
738 * sparsity pattern classes instead (see
739 * @ref Sparsity).
740 */
741 unsigned int
742 max_couplings_between_boundary_dofs() const;
743
744 /*--------------------------------------*/
745
746 /**
747 * @name Cell iterator functions
748 */
749
750 /*
751 * @{
752 */
753
754 /**
755 * Iterator to the first used cell on level @p level.
756 */
757 cell_iterator
758 begin(const unsigned int level = 0) const;
759
760 /**
761 * Iterator to the first active cell on level @p level. If the given level
762 * does not contain any active cells (i.e., all cells on this level are
763 * further refined), then this function returns
764 * <code>end_active(level)</code> so that loops of the kind
765 * @code
766 * for (cell=dof_handler.begin_active(level);
767 * cell!=dof_handler.end_active(level);
768 * ++cell)
769 * {
770 * ...
771 * }
772 * @endcode
773 * have zero iterations, as may be expected if there are no active cells on
774 * this level.
775 */
776 active_cell_iterator
777 begin_active(const unsigned int level = 0) const;
778
779 /**
780 * Iterator past the end; this iterator serves for comparisons of iterators
781 * with past-the-end or before-the-beginning states.
782 */
783 cell_iterator
784 end() const;
785
786 /**
787 * Return an iterator which is the first iterator not on the given level. If
788 * @p level is the last level, then this returns <tt>end()</tt>.
789 */
790 cell_iterator
791 end(const unsigned int level) const;
792
793 /**
794 * Return an active iterator which is the first active iterator not on the
795 * given level. If @p level is the last level, then this returns
796 * <tt>end()</tt>.
797 */
798 active_cell_iterator
799 end_active(const unsigned int level) const;
800
801 /**
802 * Iterator to the first used cell on level @p level. This returns a
803 * level_cell_iterator that returns level dofs when dof_indices() is called.
804 */
805 level_cell_iterator
806 begin_mg(const unsigned int level = 0) const;
807
808 /**
809 * Iterator past the last cell on level @p level. This returns a
810 * level_cell_iterator that returns level dofs when dof_indices() is called.
811 */
812 level_cell_iterator
813 end_mg(const unsigned int level) const;
814
815 /**
816 * Iterator past the end; this iterator serves for comparisons of iterators
817 * with past-the-end or before-the-beginning states.
818 */
819 level_cell_iterator
820 end_mg() const;
821
822 /**
823 * @name Cell iterator functions returning ranges of iterators
824 */
825
826 /**
827 * Return an iterator range that contains all cells (active or not) that
828 * make up this DoFHandler. Such a range is useful to initialize range-based
829 * for loops as supported by C++11. See the example in the documentation of
830 * active_cell_iterators().
831 *
832 * @return The half open range <code>[this->begin(), this->end())</code>
833 *
834 * @ingroup CPP11
835 */
836 IteratorRange<cell_iterator>
837 cell_iterators() const;
838
839 /**
840 * Return an iterator range that contains all active cells that make up this
841 * DoFHandler. Such a range is useful to initialize range-based for loops as
842 * supported by C++11, see also
843 * @ref CPP11 "C++11 standard".
844 *
845 * Range-based for loops are useful in that they require much less code than
846 * traditional loops (see <a href="http://en.wikipedia.org/wiki/C%2B%2B11
847 * #Range-based_for_loop">here</a> for a discussion of how they work). An
848 * example is that without range-based for loops, one often writes code such
849 * as the following:
850 * @code
851 * DoFHandler<dim> dof_handler;
852 * ...
853 * typename DoFHandler<dim>::active_cell_iterator
854 * cell = dof_handler.begin_active(),
855 * endc = dof_handler.end();
856 * for (; cell!=endc; ++cell)
857 * {
858 * fe_values.reinit (cell);
859 * ...do the local integration on 'cell'...;
860 * }
861 * @endcode
862 * Using C++11's range-based for loops, this is now entirely equivalent to
863 * the following:
864 * @code
865 * DoFHandler<dim> dof_handler;
866 * ...
867 * for (const auto &cell : dof_handler.active_cell_iterators())
868 * {
869 * fe_values.reinit (cell);
870 * ...do the local integration on 'cell'...;
871 * }
872 * @endcode
873 *
874 * @return The half open range <code>[this->begin_active(),
875 * this->end())</code>
876 *
877 * @ingroup CPP11
878 */
879 IteratorRange<active_cell_iterator>
880 active_cell_iterators() const;
881
882 /**
883 * Return an iterator range that contains all cells (active or not) that
884 * make up this DoFHandler in their level-cell form. Such a range is useful
885 * to initialize range-based for loops as supported by C++11. See the
886 * example in the documentation of active_cell_iterators().
887 *
888 * @return The half open range <code>[this->begin_mg(),
889 * this->end_mg())</code>
890 *
891 * @ingroup CPP11
892 */
893 IteratorRange<level_cell_iterator>
894 mg_cell_iterators() const;
895
896 /**
897 * Return an iterator range that contains all cells (active or not) that
898 * make up the given level of this DoFHandler. Such a range is useful to
899 * initialize range-based for loops as supported by C++11. See the example
900 * in the documentation of active_cell_iterators().
901 *
902 * @param[in] level A given level in the refinement hierarchy of this
903 * triangulation.
904 * @return The half open range <code>[this->begin(level),
905 * this->end(level))</code>
906 *
907 * @pre level must be less than this->n_levels().
908 *
909 * @ingroup CPP11
910 */
911 IteratorRange<cell_iterator>
912 cell_iterators_on_level(const unsigned int level) const;
913
914 /**
915 * Return an iterator range that contains all active cells that make up the
916 * given level of this DoFHandler. Such a range is useful to initialize
917 * range-based for loops as supported by C++11. See the example in the
918 * documentation of active_cell_iterators().
919 *
920 * @param[in] level A given level in the refinement hierarchy of this
921 * triangulation.
922 * @return The half open range <code>[this->begin_active(level),
923 * this->end(level))</code>
924 *
925 * @pre level must be less than this->n_levels().
926 *
927 * @ingroup CPP11
928 */
929 IteratorRange<active_cell_iterator>
930 active_cell_iterators_on_level(const unsigned int level) const;
931
932 /**
933 * Return an iterator range that contains all cells (active or not) that
934 * make up the given level of this DoFHandler in their level-cell form. Such
935 * a range is useful to initialize range-based for loops as supported by
936 * C++11. See the example in the documentation of active_cell_iterators().
937 *
938 * @param[in] level A given level in the refinement hierarchy of this
939 * triangulation.
940 * @return The half open range <code>[this->begin_mg(level),
941 * this->end_mg(level))</code>
942 *
943 * @pre level must be less than this->n_levels().
944 *
945 * @ingroup CPP11
946 */
947 IteratorRange<level_cell_iterator>
948 mg_cell_iterators_on_level(const unsigned int level) const;
949
950 /*
951 * @}
952 */
953
954
955 /*---------------------------------------*/
956
957
958 /**
959 * Return the global number of degrees of freedom. If the current object
960 * handles all degrees of freedom itself (even if you may intend to solve
961 * your linear system in parallel, such as in step-17 or step-18), then this
962 * number equals the number of locally owned degrees of freedom since this
963 * object doesn't know anything about what you want to do with it and
964 * believes that it owns every degree of freedom it knows about.
965 *
966 * On the other hand, if this object operates on a
967 * parallel::distributed::Triangulation object, then this function returns
968 * the global number of degrees of freedom, accumulated over all processors.
969 *
970 * In either case, included in the returned number are those DoFs which are
971 * constrained by hanging nodes, see
972 * @ref constraints.
973 *
974 * Mathematically speaking, the number returned by this function equals the
975 * dimension of the finite element space (without taking into account
976 * constraints) that corresponds to (i) the mesh on which it is defined,
977 * and (ii) the finite element that is used by the current object. It
978 * also, of course, equals the number of shape functions that span this
979 * space.
980 */
981 types::global_dof_index
982 n_dofs() const;
983
984 /**
985 * The (global) number of multilevel degrees of freedom on a given level.
986 *
987 * If no level degrees of freedom have been assigned to this level, returns
988 * numbers::invalid_dof_index. Else returns the number of degrees of freedom
989 * on this level.
990 */
991 types::global_dof_index
992 n_dofs(const unsigned int level) const;
993
994 /**
995 * Return the number of locally owned degrees of freedom located on the
996 * boundary.
997 */
998 types::global_dof_index
999 n_boundary_dofs() const;
1000
1001 /**
1002 * Return the number of locally owned degrees of freedom located on those
1003 * parts of the boundary which have a boundary indicator listed in the given
1004 * set.
1005 * The reason that a @p map rather than a @p set is used is the same as
1006 * described in the documentation of that variant of
1007 * DoFTools::make_boundary_sparsity_pattern() that takes a map.
1008 *
1009 * There is, however, another overload of this function that takes
1010 * a @p set argument (see below).
1011 */
1012 template <typename number>
1013 types::global_dof_index
1014 n_boundary_dofs(
1015 const std::map<types::boundary_id, const Function<spacedim, number> *>
1016 &boundary_ids) const;
1017
1018 /**
1019 * Return the number of degrees of freedom located on those parts of the
1020 * boundary which have a boundary indicator listed in the given set. The
1021 */
1022 types::global_dof_index
1023 n_boundary_dofs(const std::set<types::boundary_id> &boundary_ids) const;
1024
1025 /**
1026 * Access to an object informing of the block structure of the dof handler.
1027 *
1028 * If an FESystem is used in distribute_dofs(), degrees of freedom naturally
1029 * split into several
1030 * @ref GlossBlock "blocks".
1031 * For each base element as many blocks appear as its multiplicity.
1032 *
1033 * At the end of distribute_dofs(), the number of degrees of freedom in each
1034 * block is counted, and stored in a BlockInfo object, which can be accessed
1035 * here. If you have previously called distribute_mg_dofs(), the same is
1036 * done on each level of the multigrid hierarchy. Additionally, the block
1037 * structure on each cell can be generated in this object by calling
1038 * initialize_local_block_info().
1039 */
1040 const BlockInfo &
1041 block_info() const;
1042
1043 /**
1044 * Return the number of degrees of freedom that belong to this process.
1045 *
1046 * If this is a sequential DoFHandler, then the result equals that produced by
1047 * n_dofs(). (Here, "sequential" means that either
1048 * the whole program does not use MPI, or that it uses MPI
1049 * but only uses a single MPI process, or that there are multiple MPI
1050 * processes but the Triangulation on which this DoFHandler builds
1051 * works only on one MPI process.)
1052 * On the other hand, if we are operating on a
1053 * parallel::distributed::Triangulation or parallel::shared::Triangulation,
1054 * then it includes only the degrees
1055 * of freedom that the current processor owns. Note that in this case this
1056 * does not include all degrees of freedom that have been distributed on the
1057 * current processor's image of the mesh: in particular, some of the degrees
1058 * of freedom on the interface between the cells owned by this processor and
1059 * cells owned by other processors may be theirs, and degrees of freedom on
1060 * ghost cells are also not necessarily included.
1061 */
1062 types::global_dof_index
1063 n_locally_owned_dofs() const;
1064
1065 /**
1066 * Return an IndexSet describing the set of locally owned DoFs as a subset
1067 * of 0..n_dofs(). The number of elements of this set equals
1068 * n_locally_owned_dofs().
1069 */
1070 const IndexSet &
1071 locally_owned_dofs() const;
1072
1073 /**
1074 * Return an IndexSet describing the set of locally owned DoFs used for the
1075 * given multigrid level as a subset of 0..n_dofs(level).
1076 */
1077 const IndexSet &
1078 locally_owned_mg_dofs(const unsigned int level) const;
1079
1080 /**
1081 * Return a vector that stores the locally owned DoFs of each processor.
1082 *
1083 * @deprecated As of deal.II version 9.2, we do not populate a vector with
1084 * the index sets of all processors by default any more due to a possibly
1085 * large memory footprint on many processors. As a consequence, this
1086 * function needs to call `Utilities::all_gather(comm, locally_owned_dofs())`
1087 * upon the first invocation, including global communication. Use
1088 * `Utilities::all_gather(comm, dof_handler.locally_owned_dofs())` instead if
1089 * using up to a few thousands of MPI ranks or some variant involving local
1090 * communication with more processors.
1091 */
1092 DEAL_II_DEPRECATED const std::vector<IndexSet> &
1093 locally_owned_dofs_per_processor() const;
1094
1095 /**
1096 * Return a vector that stores the number of degrees of freedom each
1097 * processor that participates in this triangulation owns locally. The sum
1098 * of all these numbers equals the number of degrees of freedom that exist
1099 * globally, i.e. what n_dofs() returns.
1100 *
1101 * @deprecated As of deal.II version 9.2, we do not populate a vector with
1102 * the numbers of dofs of all processors by default any more due to a
1103 * possibly large memory footprint on many processors. As a consequence,
1104 * this function needs to call `Utilities::all_gather(comm,
1105 * n_locally_owned_dofs()` upon the first invocation, including global
1106 * communication. Use `Utilities::all_gather(comm,
1107 * dof_handler.n_locally_owned_dofs()` instead if using up to a few thousands
1108 * of MPI ranks or some variant involving local communication with more
1109 * processors.
1110 */
1111 DEAL_II_DEPRECATED const std::vector<types::global_dof_index> &
1112 n_locally_owned_dofs_per_processor() const;
1113
1114 /**
1115 * Return a vector that stores the locally owned DoFs of each processor on
1116 * the given level @p level.
1117 *
1118 * @deprecated As of deal.II version 9.2, we do not populate a vector with
1119 * the index sets of all processors by default any more due to a possibly
1120 * large memory footprint on many processors. As a consequence, this
1121 * function needs to call `Utilities::all_gather(comm,
1122 * locally_owned_dofs_mg())` upon the first invocation, including global
1123 * communication. Use `Utilities::all_gather(comm,
1124 * dof_handler.locally_owned_dofs_mg())` instead if using up to a few
1125 * thousands of MPI ranks or some variant involving local communication with
1126 * more processors.
1127 */
1128 DEAL_II_DEPRECATED const std::vector<IndexSet> &
1129 locally_owned_mg_dofs_per_processor(const unsigned int level) const;
1130
1131 /**
1132 * Return a constant reference to the selected finite element object.
1133 * Since there is only one FiniteElement @p index must be equal to zero
1134 * which is also the default value.
1135 */
1136 const FiniteElement<dim, spacedim> &
1137 get_fe(const unsigned int index = 0) const;
1138
1139 /**
1140 * Return a constant reference to the set of finite element objects that
1141 * are used by this @p DoFHandler. Since this object only contains one
1142 * FiniteElement, only this one object is returned wrapped in a
1143 * hp::FECollection.
1144 */
1145 const hp::FECollection<dim, spacedim> &
1146 get_fe_collection() const;
1147
1148 /**
1149 * Return a constant reference to the triangulation underlying this object.
1150 */
1151 const Triangulation<dim, spacedim> &
1152 get_triangulation() const;
1153
1154 /**
1155 * Whenever serialization with a parallel::distributed::Triangulation as the
1156 * underlying triangulation is considered, we also need to consider storing
1157 * the active_fe_indices on all active cells as well.
1158 *
1159 * This function registers that these indices are to be stored whenever the
1160 * parallel::distributed::Triangulation::save() function is called on the
1161 * underlying triangulation.
1162 *
1163 * @note Currently only implemented for triangulations of type
1164 * parallel::distributed::Triangulation. An assertion will be triggered if
1165 * a different type is registered.
1166 *
1167 * @see The documentation of parallel::distributed::SolutionTransfer has further
1168 * information on serialization.
1169 */
1170 void
1171 prepare_for_serialization_of_active_fe_indices();
1172
1173 /**
1174 * Whenever serialization with a parallel::distributed::Triangulation as the
1175 * underlying triangulation is considered, we also need to consider storing
1176 * the active_fe_indices on all active cells as well.
1177 *
1178 * This function deserializes and distributes the previously stored
1179 * active_fe_indices on all active cells.
1180 *
1181 * @note Currently only implemented for triangulations of type
1182 * parallel::distributed::Triangulation. An assertion will be triggered if
1183 * a different type is registered.
1184 *
1185 * @see The documentation of parallel::distributed::SolutionTransfer has further
1186 * information on serialization.
1187 */
1188 void
1189 deserialize_active_fe_indices();
1190
1191 /**
1192 * Determine an estimate for the memory consumption (in bytes) of this
1193 * object.
1194 *
1195 * This function is made virtual, since a dof handler object might be
1196 * accessed through a pointers to this base class, although the actual
1197 * object might be a derived class.
1198 */
1199 virtual std::size_t
1200 memory_consumption() const;
1201
1202 /**
1203 * Write the data of this object to a stream for the purpose of
1204 * serialization.
1205 */
1206 template <class Archive>
1207 void
1208 save(Archive &ar, const unsigned int version) const;
1209
1210 /**
1211 * Read the data of this object from a stream for the purpose of
1212 * serialization.
1213 */
1214 template <class Archive>
1215 void
1216 load(Archive &ar, const unsigned int version);
1217
1218 #ifdef DOXYGEN
1219 /**
1220 * Write and read the data of this object from a stream for the purpose
1221 * of serialization.
1222 */
1223 template <class Archive>
1224 void
1225 serialize(Archive &archive, const unsigned int version);
1226 #else
1227 // This macro defines the serialize() method that is compatible with
1228 // the templated save() and load() method that have been implemented.
1229 BOOST_SERIALIZATION_SPLIT_MEMBER()
1230 #endif
1231
1232 /**
1233 * Exception
1234 */
1235 DeclException0(ExcNoFESelected);
1236 /**
1237 * Exception
1238 * @ingroup Exceptions
1239 */
1240 DeclException0(ExcInvalidBoundaryIndicator);
1241 /**
1242 * Exception
1243 * @ingroup Exceptions
1244 */
1245 DeclException1(ExcInvalidLevel,
1246 int,
1247 << "The given level " << arg1
1248 << " is not in the valid range!");
1249 /**
1250 * Exception
1251 * @ingroup Exceptions
1252 */
1253 DeclException1(ExcNewNumbersNotConsecutive,
1254 types::global_dof_index,
1255 << "The given list of new dof indices is not consecutive: "
1256 << "the index " << arg1 << " does not exist.");
1257 /**
1258 * Exception
1259 */
1260 DeclException2(ExcInvalidFEIndex,
1261 int,
1262 int,
1263 << "The mesh contains a cell with an active_fe_index of "
1264 << arg1 << ", but the finite element collection only has "
1265 << arg2 << " elements");
1266
1267 /**
1268 * Exception used when a certain feature doesn't make sense when
1269 * DoFHandler does not hp capabilities.
1270 */
1271 DeclExceptionMsg(ExcNotAvailableWithoutHP,
1272 "The current function doesn't make sense when used with a "
1273 "DoFHandler without hp capabilities.");
1274
1275 /**
1276 * Exception used when a certain feature is not implemented when the
1277 * DoFHandler has hp capabilities.
1278 */
1279 DeclExceptionMsg(ExcNotImplementedWithHP,
1280 "The current function has not yet been implemented for a "
1281 "DoFHandler with hp capabilities.");
1282
1283 private:
1284 /**
1285 * A data structure that is used to store the DoF indices associated with a
1286 * particular vertex. Unlike cells, vertices live on several levels of a
1287 * multigrid hierarchy; consequently, we need to store DoF indices for each
1288 * vertex for each of the levels it lives on. This class does this.
1289 */
1290 class MGVertexDoFs
1291 {
1292 public:
1293 /**
1294 * Constructor.
1295 */
1296 MGVertexDoFs();
1297
1298 /**
1299 * A function that is called to allocate the necessary amount of memory to
1300 * store the indices of the DoFs that live on this vertex for the given
1301 * (inclusive) range of levels.
1302 */
1303 void
1304 init(const unsigned int coarsest_level,
1305 const unsigned int finest_level,
1306 const unsigned int dofs_per_vertex);
1307
1308 /**
1309 * Return the coarsest level for which this structure stores data.
1310 */
1311 unsigned int
1312 get_coarsest_level() const;
1313
1314 /**
1315 * Return the finest level for which this structure stores data.
1316 */
1317 unsigned int
1318 get_finest_level() const;
1319
1320 /**
1321 * Return the index of the <code>dof_number</code>th degree of freedom for
1322 * the given level stored for the current vertex.
1323 */
1324 types::global_dof_index
1325 get_index(const unsigned int level,
1326 const unsigned int dof_number,
1327 const unsigned int dofs_per_vertex) const;
1328
1329 /**
1330 * Set the index of the <code>dof_number</code>th degree of freedom for
1331 * the given level stored for the current vertex to <code>index</code>.
1332 */
1333 void
1334 set_index(const unsigned int level,
1335 const unsigned int dof_number,
1336 const unsigned int dofs_per_vertex,
1337 const types::global_dof_index index);
1338
1339 private:
1340 /**
1341 * Coarsest level for which this object stores DoF indices.
1342 */
1343 unsigned int coarsest_level;
1344
1345 /**
1346 * Finest level for which this object stores DoF indices.
1347 */
1348 unsigned int finest_level;
1349
1350 /**
1351 * A pointer to an array where we store the indices of the DoFs that live
1352 * on the various levels this vertex exists on.
1353 *
1354 * The starting offset of the DoFs that belong to a @p level are given by
1355 * <code>n_dofs_per_vertex() * (level-coarsest_level)</code>. @p n_dofs_per_vertex()
1356 * must therefore be passed as an argument to the functions that set or
1357 * read an index.
1358 */
1359 std::unique_ptr<types::global_dof_index[]> indices;
1360 };
1361
1362 /**
1363 * Whenever the underlying triangulation changes by either
1364 * h/p refinement/coarsening and serialization, the active_fe_index of cells
1365 * needs to be transferred. This structure stores all temporary information
1366 * required during that process.
1367 */
1368 struct ActiveFEIndexTransfer
1369 {
1370 /**
1371 * Container to temporarily store the iterator and future active FE index
1372 * of cells that persist.
1373 */
1374 std::map<const cell_iterator, const unsigned int> persisting_cells_fe_index;
1375
1376 /**
1377 * Container to temporarily store the iterator and future active FE index
1378 * of cells that will be refined.
1379 */
1380 std::map<const cell_iterator, const unsigned int> refined_cells_fe_index;
1381
1382 /**
1383 * Container to temporarily store the iterator and future active FE index
1384 * of parent cells that will remain after coarsening.
1385 */
1386 std::map<const cell_iterator, const unsigned int> coarsened_cells_fe_index;
1387
1388 /**
1389 * Container to temporarily store the active_fe_index of every locally
1390 * owned cell for transfer across parallel::distributed::Triangulation
1391 * objects.
1392 */
1393 std::vector<unsigned int> active_fe_indices;
1394
1395 /**
1396 * Helper object to transfer all active_fe_indices on
1397 * parallel::distributed::Triangulation objects during
1398 * refinement/coarsening and serialization.
1399 */
1400 std::unique_ptr<
1401 parallel::distributed::
1402 CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>
1403 cell_data_transfer;
1404 };
1405
1406 /**
1407 * An object containing information on the block structure.
1408 */
1409 BlockInfo block_info_object;
1410
1411 /**
1412 * Address of the triangulation to work on.
1413 */
1414 SmartPointer<const Triangulation<dim, spacedim>, DoFHandler<dim, spacedim>>
1415 tria;
1416
1417 /**
1418 * Store a hp::FECollection object. If only a single FiniteElement is
1419 * used during initialization of this object, it contains the (one)
1420 * FiniteElement.
1421 */
1422 hp::FECollection<dim, spacedim> fe_collection;
1423
1424 /**
1425 * An object that describes how degrees of freedom should be distributed and
1426 * renumbered.
1427 */
1428 std::unique_ptr<dealii::internal::DoFHandlerImplementation::Policy::
1429 PolicyBase<dim, spacedim>>
1430 policy;
1431
1432 /**
1433 * A structure that contains all sorts of numbers that characterize the
1434 * degrees of freedom this object works on.
1435 *
1436 * For most members of this structure, there is an accessor function in this
1437 * class that returns its value.
1438 */
1439 dealii::internal::DoFHandlerImplementation::NumberCache number_cache;
1440
1441 /**
1442 * Data structure like number_cache, but for each multigrid level.
1443 */
1444 std::vector<dealii::internal::DoFHandlerImplementation::NumberCache>
1445 mg_number_cache;
1446
1447 /**
1448 * Cached indices of the degrees of freedom of all active cell. Identification
1449 * of the appropriate position of a cell in the vectors is done via
1450 * cell_dof_cache_ptr (CRS scheme).
1451 */
1452 mutable std::vector<std::vector<types::global_dof_index>>
1453 cell_dof_cache_indices;
1454
1455 /**
1456 * Pointer to the first cached degree of freedom of an active cell
1457 * (identified by level and level index) within cell_dof_cache_indices.
1458 */
1459 mutable std::vector<std::vector<offset_type>> cell_dof_cache_ptr;
1460
1461 /**
1462 * Indices of degree of freedom of each d+1 geometric object (3D: vertex,
1463 * line, quad, hex) for all relevant active finite elements. Identification
1464 * of the appropriate position is done via object_dof_ptr (CRS scheme).
1465 */
1466 mutable std::vector<std::array<std::vector<types::global_dof_index>, dim + 1>>
1467 object_dof_indices;
1468
1469 /**
1470 * Pointer to the first cached degree of freedom of a geometric object for all
1471 * relevant active finite elements.
1472 *
1473 * @note In normal mode it is possible to access this data strucutre directly.
1474 * In hp mode, an indirection via hp_object_fe_indices/hp_object_fe_ptr is
1475 * necessary.
1476 */
1477 mutable std::vector<std::array<std::vector<offset_type>, dim + 1>>
1478 object_dof_ptr;
1479
1480 /**
1481 * Active fe indices of each geometric object. Identification
1482 * of the appropriate position of a cell in the vectors is done via
1483 * hp_object_fe_ptr (CRS scheme).
1484 */
1485 mutable std::array<std::vector<active_fe_index_type>, dim + 1>
1486 hp_object_fe_indices;
1487
1488 /**
1489 * Pointer to the first fe index of a geometric object.
1490 */
1491 mutable std::array<std::vector<offset_type>, dim + 1> hp_object_fe_ptr;
1492
1493 /**
1494 * Active fe index of an active cell (identified by level and level index).
1495 * This vector is only used in hp mode.
1496 */
1497 mutable std::vector<std::vector<active_fe_index_type>>
1498 hp_cell_active_fe_indices;
1499
1500 /**
1501 * Future fe index of an active cell (identified by level and level index).
1502 * This vector is only used in hp mode.
1503 */
1504 mutable std::vector<std::vector<active_fe_index_type>>
1505 hp_cell_future_fe_indices;
1506
1507 /**
1508 * An array to store the indices for level degrees of freedom located at
1509 * vertices.
1510 */
1511 std::vector<MGVertexDoFs> mg_vertex_dofs;
1512
1513 /**
1514 * Space to store the DoF numbers for the different multigrid levels.
1515 */
1516 std::vector<
1517 std::unique_ptr<dealii::internal::DoFHandlerImplementation::DoFLevel<dim>>>
1518 mg_levels;
1519
1520 /**
1521 * Space to store DoF numbers of faces in the multigrid context.
1522 */
1523 std::unique_ptr<dealii::internal::DoFHandlerImplementation::DoFFaces<dim>>
1524 mg_faces;
1525
1526 /**
1527 * We embed our data structure into a pointer to control that
1528 * all transfer related data only exists during the actual transfer process.
1529 */
1530 std::unique_ptr<ActiveFEIndexTransfer> active_fe_index_transfer;
1531
1532 /**
1533 * A list of connections with which this object connects to the
1534 * triangulation to get information about when the triangulation changes.
1535 */
1536 std::vector<boost::signals2::connection> tria_listeners;
1537
1538 /**
1539 * Free all memory used for non-multigrid data structures.
1540 */
1541 void
1542 clear_space();
1543
1544 /**
1545 * Free all memory used for multigrid data structures.
1546 */
1547 void
1548 clear_mg_space();
1549
1550 /**
1551 * Return dof index of specified object.
1552 */
1553 template <int structdim>
1554 types::global_dof_index
1555 get_dof_index(const unsigned int obj_level,
1556 const unsigned int obj_index,
1557 const unsigned int fe_index,
1558 const unsigned int local_index) const;
1559
1560 /**
1561 * Return dof index of specified object.
1562 */
1563 template <int structdim>
1564 void
1565 set_dof_index(const unsigned int obj_level,
1566 const unsigned int obj_index,
1567 const unsigned int fe_index,
1568 const unsigned int local_index,
1569 const types::global_dof_index global_index) const;
1570
1571 /**
1572 * Setup DoFHandler policy.
1573 */
1574 void
1575 setup_policy();
1576
1577 /**
1578 * Setup DoFHandler policy and listeners (in the hp-context).
1579 */
1580 void
1581 setup_policy_and_listeners();
1582
1583 /**
1584 * Create default tables for the active_fe_indices.
1585 * They are initialized with a zero
1586 * indicator, meaning that fe[0] is going to be used by default. This
1587 * method is called before refinement and while setting the finite elements
1588 * via set_fe(). It ensures each cell has a valid active_fe_index.
1589 */
1590 void
1591 create_active_fe_table();
1592
1593 /**
1594 * A function that will be triggered through a triangulation
1595 * signal just before the triangulation is modified.
1596 *
1597 * The function that stores the active_fe_flags of all cells that will
1598 * be refined or coarsened before the refinement happens, so that
1599 * they can be set again after refinement.
1600 */
1601 void
1602 pre_refinement_action();
1603
1604 /**
1605 * A function that will be triggered through a triangulation
1606 * signal just after the triangulation is modified.
1607 *
1608 * The function that restores the active_fe_flags of all cells that
1609 * were refined.
1610 */
1611 void
1612 post_refinement_action();
1613
1614 /**
1615 * A function that will be triggered through a triangulation
1616 * signal just before the associated Triangulation or
1617 * parallel::shared::Triangulation is modified.
1618 *
1619 * The function that stores the active_fe_indices of all cells that will
1620 * be refined or coarsened before the refinement happens, so that
1621 * they can be set again after refinement.
1622 */
1623 void
1624 pre_active_fe_index_transfer();
1625
1626 /**
1627 * A function that will be triggered through a triangulation
1628 * signal just before the associated parallel::distributed::Triangulation is
1629 * modified.
1630 *
1631 * The function that stores all active_fe_indices on locally owned cells for
1632 * distribution over all participating processors.
1633 */
1634 void
1635 pre_distributed_active_fe_index_transfer();
1636
1637 /**
1638 * A function that will be triggered through a triangulation
1639 * signal just after the associated Triangulation or
1640 * parallel::shared::Triangulation is modified.
1641 *
1642 * The function that restores the active_fe_indices of all cells that
1643 * were refined or coarsened.
1644 */
1645 void
1646 post_active_fe_index_transfer();
1647
1648 /**
1649 * A function that will be triggered through a triangulation
1650 * signal just after the associated parallel::distributed::Triangulation is
1651 * modified.
1652 *
1653 * The function that restores all active_fe_indices on locally owned cells
1654 * that have been communicated.
1655 */
1656 void
1657 post_distributed_active_fe_index_transfer();
1658
1659 /**
1660 * A function that will be triggered through a triangulation
1661 * signal just after the associated parallel::distributed::Triangulation has
1662 * been saved.
1663 *
1664 * The function frees all memory related to the transfer of
1665 * active_fe_indices.
1666 */
1667 void
1668 post_distributed_serialization_of_active_fe_indices();
1669
1670
1671 // Make accessor objects friends.
1672 template <int, int, int, bool>
1673 friend class dealii::DoFAccessor;
1674 template <int, int, bool>
1675 friend class dealii::DoFCellAccessor;
1676 friend struct dealii::internal::DoFAccessorImplementation::Implementation;
1677 friend struct dealii::internal::DoFCellAccessorImplementation::Implementation;
1678
1679 // Likewise for DoFLevel objects since they need to access the vertex dofs
1680 // in the functions that set and retrieve vertex dof indices.
1681 friend struct dealii::internal::DoFHandlerImplementation::Implementation;
1682 friend struct dealii::internal::hp::DoFHandlerImplementation::Implementation;
1683 friend struct dealii::internal::DoFHandlerImplementation::Policy::
1684 Implementation;
1685
1686 // explicitly check for sensible template arguments, but not on windows
1687 // because MSVC creates bogus warnings during normal compilation
1688 #ifndef DEAL_II_MSVC
1689 static_assert(dim <= spacedim,
1690 "The dimension <dim> of a DoFHandler must be less than or "
1691 "equal to the space dimension <spacedim> in which it lives.");
1692 #endif
1693 };
1694
1695 #ifndef DOXYGEN
1696
1697 /* ----------------------- Inline functions ----------------------------------
1698 */
1699
1700
1701 template <int dim, int spacedim>
1702 inline bool
has_level_dofs()1703 DoFHandler<dim, spacedim>::has_level_dofs() const
1704 {
1705 return mg_number_cache.size() > 0;
1706 }
1707
1708
1709
1710 template <int dim, int spacedim>
1711 inline bool
has_active_dofs()1712 DoFHandler<dim, spacedim>::has_active_dofs() const
1713 {
1714 return number_cache.n_global_dofs > 0;
1715 }
1716
1717
1718
1719 template <int dim, int spacedim>
1720 inline types::global_dof_index
n_dofs()1721 DoFHandler<dim, spacedim>::n_dofs() const
1722 {
1723 return number_cache.n_global_dofs;
1724 }
1725
1726
1727
1728 template <int dim, int spacedim>
1729 inline types::global_dof_index
n_dofs(const unsigned int level)1730 DoFHandler<dim, spacedim>::n_dofs(const unsigned int level) const
1731 {
1732 Assert(has_level_dofs(),
1733 ExcMessage(
1734 "n_dofs(level) can only be called after distribute_mg_dofs()"));
1735 Assert(level < mg_number_cache.size(), ExcInvalidLevel(level));
1736 return mg_number_cache[level].n_global_dofs;
1737 }
1738
1739
1740
1741 template <int dim, int spacedim>
1742 types::global_dof_index
n_locally_owned_dofs()1743 DoFHandler<dim, spacedim>::n_locally_owned_dofs() const
1744 {
1745 return number_cache.n_locally_owned_dofs;
1746 }
1747
1748
1749
1750 template <int dim, int spacedim>
1751 const IndexSet &
locally_owned_dofs()1752 DoFHandler<dim, spacedim>::locally_owned_dofs() const
1753 {
1754 return number_cache.locally_owned_dofs;
1755 }
1756
1757
1758
1759 template <int dim, int spacedim>
1760 const IndexSet &
locally_owned_mg_dofs(const unsigned int level)1761 DoFHandler<dim, spacedim>::locally_owned_mg_dofs(const unsigned int level) const
1762 {
1763 Assert(level < this->get_triangulation().n_global_levels(),
1764 ExcMessage("The given level index exceeds the number of levels "
1765 "present in the triangulation"));
1766 Assert(
1767 mg_number_cache.size() == this->get_triangulation().n_global_levels(),
1768 ExcMessage(
1769 "The level dofs are not set up properly! Did you call distribute_mg_dofs()?"));
1770 return mg_number_cache[level].locally_owned_dofs;
1771 }
1772
1773
1774
1775 template <int dim, int spacedim>
1776 const std::vector<types::global_dof_index> &
n_locally_owned_dofs_per_processor()1777 DoFHandler<dim, spacedim>::n_locally_owned_dofs_per_processor() const
1778 {
1779 if (number_cache.n_locally_owned_dofs_per_processor.empty() &&
1780 number_cache.n_global_dofs > 0)
1781 {
1782 MPI_Comm comm;
1783
1784 const parallel::TriangulationBase<dim, spacedim> *tr =
1785 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1786 &this->get_triangulation()));
1787 if (tr != nullptr)
1788 comm = tr->get_communicator();
1789 else
1790 comm = MPI_COMM_SELF;
1791
1792 const_cast<dealii::internal::DoFHandlerImplementation::NumberCache &>(
1793 number_cache)
1794 .n_locally_owned_dofs_per_processor =
1795 number_cache.get_n_locally_owned_dofs_per_processor(comm);
1796 }
1797 return number_cache.n_locally_owned_dofs_per_processor;
1798 }
1799
1800
1801
1802 template <int dim, int spacedim>
1803 const std::vector<IndexSet> &
locally_owned_dofs_per_processor()1804 DoFHandler<dim, spacedim>::locally_owned_dofs_per_processor() const
1805 {
1806 if (number_cache.locally_owned_dofs_per_processor.empty() &&
1807 number_cache.n_global_dofs > 0)
1808 {
1809 MPI_Comm comm;
1810
1811 const parallel::TriangulationBase<dim, spacedim> *tr =
1812 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1813 &this->get_triangulation()));
1814 if (tr != nullptr)
1815 comm = tr->get_communicator();
1816 else
1817 comm = MPI_COMM_SELF;
1818
1819 const_cast<dealii::internal::DoFHandlerImplementation::NumberCache &>(
1820 number_cache)
1821 .locally_owned_dofs_per_processor =
1822 number_cache.get_locally_owned_dofs_per_processor(comm);
1823 }
1824 return number_cache.locally_owned_dofs_per_processor;
1825 }
1826
1827
1828
1829 template <int dim, int spacedim>
1830 const std::vector<IndexSet> &
locally_owned_mg_dofs_per_processor(const unsigned int level)1831 DoFHandler<dim, spacedim>::locally_owned_mg_dofs_per_processor(
1832 const unsigned int level) const
1833 {
1834 Assert(level < this->get_triangulation().n_global_levels(),
1835 ExcMessage("The given level index exceeds the number of levels "
1836 "present in the triangulation"));
1837 Assert(
1838 mg_number_cache.size() == this->get_triangulation().n_global_levels(),
1839 ExcMessage(
1840 "The level dofs are not set up properly! Did you call distribute_mg_dofs()?"));
1841 if (mg_number_cache[level].locally_owned_dofs_per_processor.empty() &&
1842 mg_number_cache[level].n_global_dofs > 0)
1843 {
1844 MPI_Comm comm;
1845
1846 const parallel::TriangulationBase<dim, spacedim> *tr =
1847 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1848 &this->get_triangulation()));
1849 if (tr != nullptr)
1850 comm = tr->get_communicator();
1851 else
1852 comm = MPI_COMM_SELF;
1853
1854 const_cast<dealii::internal::DoFHandlerImplementation::NumberCache &>(
1855 mg_number_cache[level])
1856 .locally_owned_dofs_per_processor =
1857 mg_number_cache[level].get_locally_owned_dofs_per_processor(comm);
1858 }
1859 return mg_number_cache[level].locally_owned_dofs_per_processor;
1860 }
1861
1862
1863
1864 template <int dim, int spacedim>
1865 inline const FiniteElement<dim, spacedim> &
get_fe(const unsigned int number)1866 DoFHandler<dim, spacedim>::get_fe(const unsigned int number) const
1867 {
1868 Assert(fe_collection.size() > 0,
1869 ExcMessage("No finite element collection is associated with "
1870 "this DoFHandler"));
1871 return fe_collection[number];
1872 }
1873
1874
1875
1876 template <int dim, int spacedim>
1877 inline const hp::FECollection<dim, spacedim> &
get_fe_collection()1878 DoFHandler<dim, spacedim>::get_fe_collection() const
1879 {
1880 Assert(fe_collection.size() > 0,
1881 ExcMessage("No finite element collection is associated with "
1882 "this DoFHandler"));
1883 return fe_collection;
1884 }
1885
1886
1887
1888 template <int dim, int spacedim>
1889 inline const Triangulation<dim, spacedim> &
get_triangulation()1890 DoFHandler<dim, spacedim>::get_triangulation() const
1891 {
1892 Assert(tria != nullptr,
1893 ExcMessage("This DoFHandler object has not been associated "
1894 "with a triangulation."));
1895 return *tria;
1896 }
1897
1898
1899
1900 template <int dim, int spacedim>
1901 inline const BlockInfo &
block_info()1902 DoFHandler<dim, spacedim>::block_info() const
1903 {
1904 Assert(this->hp_capability_enabled == false, ExcNotImplementedWithHP());
1905
1906 return block_info_object;
1907 }
1908
1909
1910
1911 template <int dim, int spacedim>
1912 template <typename number>
1913 types::global_dof_index
n_boundary_dofs(const std::map<types::boundary_id,const Function<spacedim,number> * > & boundary_ids)1914 DoFHandler<dim, spacedim>::n_boundary_dofs(
1915 const std::map<types::boundary_id, const Function<spacedim, number> *>
1916 &boundary_ids) const
1917 {
1918 Assert(!(dim == 2 && spacedim == 3) || this->hp_capability_enabled == false,
1919 ExcNotImplementedWithHP());
1920
1921 // extract the set of boundary ids and forget about the function object
1922 // pointers
1923 std::set<types::boundary_id> boundary_ids_only;
1924 for (typename std::map<types::boundary_id,
1925 const Function<spacedim, number> *>::const_iterator p =
1926 boundary_ids.begin();
1927 p != boundary_ids.end();
1928 ++p)
1929 boundary_ids_only.insert(p->first);
1930
1931 // then just hand everything over to the other function that does the work
1932 return n_boundary_dofs(boundary_ids_only);
1933 }
1934
1935
1936
1937 namespace internal
1938 {
1939 /**
1940 * Return a string representing the dynamic type of the given argument.
1941 * This is basically the same what typeid(...).name() does, but it turns out
1942 * this is broken on Intel 13+.
1943 *
1944 * Defined in dof_handler.cc.
1945 */
1946 template <int dim, int spacedim>
1947 std::string
1948 policy_to_string(const dealii::internal::DoFHandlerImplementation::Policy::
1949 PolicyBase<dim, spacedim> &policy);
1950 } // namespace internal
1951
1952
1953
1954 template <int dim, int spacedim>
1955 template <class Archive>
1956 void
save(Archive & ar,const unsigned int)1957 DoFHandler<dim, spacedim>::save(Archive &ar, const unsigned int) const
1958 {
1959 if (this->hp_capability_enabled)
1960 {
1961 ar & this->object_dof_indices;
1962 ar & this->object_dof_ptr;
1963
1964 ar & this->cell_dof_cache_indices;
1965 ar & this->cell_dof_cache_ptr;
1966
1967 ar & this->hp_cell_active_fe_indices;
1968 ar & this->hp_cell_future_fe_indices;
1969
1970 ar &hp_object_fe_ptr;
1971 ar &hp_object_fe_indices;
1972
1973 ar &number_cache;
1974
1975 ar &mg_number_cache;
1976
1977 // write out the number of triangulation cells and later check during
1978 // loading that this number is indeed correct; same with something that
1979 // identifies the policy
1980 const unsigned int n_cells = this->tria->n_cells();
1981 std::string policy_name =
1982 dealii::internal::policy_to_string(*this->policy);
1983
1984 ar &n_cells &policy_name;
1985 }
1986 else
1987 {
1988 ar & this->block_info_object;
1989 ar &number_cache;
1990
1991 ar & this->object_dof_indices;
1992 ar & this->object_dof_ptr;
1993
1994 ar & this->cell_dof_cache_indices;
1995 ar & this->cell_dof_cache_ptr;
1996
1997 // write out the number of triangulation cells and later check during
1998 // loading that this number is indeed correct; same with something that
1999 // identifies the FE and the policy
2000 unsigned int n_cells = this->tria->n_cells();
2001 std::string fe_name = this->get_fe(0).get_name();
2002 std::string policy_name = internal::policy_to_string(*this->policy);
2003
2004 ar &n_cells &fe_name &policy_name;
2005 }
2006 }
2007
2008
2009
2010 template <int dim, int spacedim>
2011 template <class Archive>
2012 void
load(Archive & ar,const unsigned int)2013 DoFHandler<dim, spacedim>::load(Archive &ar, const unsigned int)
2014 {
2015 if (this->hp_capability_enabled)
2016 {
2017 ar & this->object_dof_indices;
2018 ar & this->object_dof_ptr;
2019
2020 ar & this->cell_dof_cache_indices;
2021 ar & this->cell_dof_cache_ptr;
2022
2023 ar & this->hp_cell_active_fe_indices;
2024 ar & this->hp_cell_future_fe_indices;
2025
2026 ar &hp_object_fe_ptr;
2027 ar &hp_object_fe_indices;
2028
2029 ar &number_cache;
2030
2031 ar &mg_number_cache;
2032
2033 // these are the checks that correspond to the last block in the save()
2034 // function
2035 unsigned int n_cells;
2036 std::string policy_name;
2037
2038 ar &n_cells &policy_name;
2039
2040 AssertThrow(
2041 n_cells == this->tria->n_cells(),
2042 ExcMessage(
2043 "The object being loaded into does not match the triangulation "
2044 "that has been stored previously."));
2045 AssertThrow(
2046 policy_name == dealii::internal::policy_to_string(*this->policy),
2047 ExcMessage("The policy currently associated with this DoFHandler (" +
2048 dealii::internal::policy_to_string(*this->policy) +
2049 ") does not match the one that was associated with the "
2050 "DoFHandler previously stored (" +
2051 policy_name + ")."));
2052 }
2053 else
2054 {
2055 ar & this->block_info_object;
2056 ar &number_cache;
2057
2058 object_dof_indices.clear();
2059
2060 object_dof_ptr.clear();
2061
2062 ar & this->object_dof_indices;
2063 ar & this->object_dof_ptr;
2064
2065 ar & this->cell_dof_cache_indices;
2066 ar & this->cell_dof_cache_ptr;
2067
2068 // these are the checks that correspond to the last block in the save()
2069 // function
2070 unsigned int n_cells;
2071 std::string fe_name;
2072 std::string policy_name;
2073
2074 ar &n_cells &fe_name &policy_name;
2075
2076 AssertThrow(
2077 n_cells == this->tria->n_cells(),
2078 ExcMessage(
2079 "The object being loaded into does not match the triangulation "
2080 "that has been stored previously."));
2081 AssertThrow(
2082 fe_name == this->get_fe(0).get_name(),
2083 ExcMessage(
2084 "The finite element associated with this DoFHandler does not match "
2085 "the one that was associated with the DoFHandler previously stored."));
2086 AssertThrow(policy_name == internal::policy_to_string(*this->policy),
2087 ExcMessage(
2088 "The policy currently associated with this DoFHandler (" +
2089 internal::policy_to_string(*this->policy) +
2090 ") does not match the one that was associated with the "
2091 "DoFHandler previously stored (" +
2092 policy_name + ")."));
2093 }
2094 }
2095
2096
2097
2098 template <int dim, int spacedim>
2099 inline types::global_dof_index
get_index(const unsigned int level,const unsigned int dof_number,const unsigned int dofs_per_vertex)2100 DoFHandler<dim, spacedim>::MGVertexDoFs::get_index(
2101 const unsigned int level,
2102 const unsigned int dof_number,
2103 const unsigned int dofs_per_vertex) const
2104 {
2105 Assert((level >= coarsest_level) && (level <= finest_level),
2106 ExcInvalidLevel(level));
2107 return indices[dofs_per_vertex * (level - coarsest_level) + dof_number];
2108 }
2109
2110
2111
2112 template <int dim, int spacedim>
2113 inline void
set_index(const unsigned int level,const unsigned int dof_number,const unsigned int dofs_per_vertex,const types::global_dof_index index)2114 DoFHandler<dim, spacedim>::MGVertexDoFs::set_index(
2115 const unsigned int level,
2116 const unsigned int dof_number,
2117 const unsigned int dofs_per_vertex,
2118 const types::global_dof_index index)
2119 {
2120 Assert((level >= coarsest_level) && (level <= finest_level),
2121 ExcInvalidLevel(level));
2122 indices[dofs_per_vertex * (level - coarsest_level) + dof_number] = index;
2123 }
2124
2125
2126 #endif // DOXYGEN
2127
2128 DEAL_II_NAMESPACE_CLOSE
2129
2130 #endif
2131