Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
dof_tools.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#include <deal.II/base/mpi.h>
18#include <deal.II/base/table.h>
20
23
27
28#include <deal.II/fe/fe.h>
29#include <deal.II/fe/fe_tools.h>
31
34#include <deal.II/grid/tria.h>
36
41
45#include <deal.II/lac/vector.h>
46
47#include <algorithm>
48#include <numeric>
49
51
52
53
54namespace DoFTools
55{
56 namespace internal
57 {
72 template <int dim, typename Number = double>
74 {
80 bool
82 const Point<dim, Number> &rhs) const
83 {
84 double downstream_size = 0;
85 double weight = 1.;
86 for (unsigned int d = 0; d < dim; ++d)
87 {
88 downstream_size += (rhs[d] - lhs[d]) * weight;
89 weight *= 1e-5;
90 }
91 if (downstream_size < 0)
92 return false;
93 else if (downstream_size > 0)
94 return true;
95 else
96 {
97 for (unsigned int d = 0; d < dim; ++d)
98 {
99 if (lhs[d] == rhs[d])
100 continue;
101 return lhs[d] < rhs[d];
102 }
103 return false;
104 }
105 }
106 };
107
108
109
110 // return an array that for each dof on the reference cell
111 // lists the corresponding vector component.
112 //
113 // if an element is non-primitive then we assign to each degree of freedom
114 // the following component:
115 // - if the nonzero components that belong to a shape function are not
116 // selected in the component_mask, then the shape function is assigned
117 // to the first nonzero vector component that corresponds to this
118 // shape function
119 // - otherwise, the shape function is assigned the first component selected
120 // in the component_mask that corresponds to this shape function
121 template <int dim, int spacedim>
122 std::vector<unsigned char>
124 const ComponentMask &component_mask)
125 {
126 std::vector<unsigned char> local_component_association(
127 fe.n_dofs_per_cell(), static_cast<unsigned char>(-1));
128
129 // compute the component each local dof belongs to.
130 // if the shape function is primitive, then this
131 // is simple and we can just associate it with
132 // what system_to_component_index gives us
133 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
134 if (fe.is_primitive(i))
135 local_component_association[i] =
136 fe.system_to_component_index(i).first;
137 else
138 // if the shape function is not primitive, then either use the
139 // component of the first nonzero component corresponding
140 // to this shape function (if the component is not specified
141 // in the component_mask), or the first component of this block
142 // that is listed in the component_mask (if the block this
143 // component corresponds to is indeed specified in the component
144 // mask)
145 {
146 const unsigned int first_comp =
148
149 if ((fe.get_nonzero_components(i) & component_mask)
150 .n_selected_components(fe.n_components()) == 0)
151 local_component_association[i] = first_comp;
152 else
153 // pick the component selected. we know from the previous 'if'
154 // that within the components that are nonzero for this
155 // shape function there must be at least one for which the
156 // mask is true, so we will for sure run into the break()
157 // at one point
158 for (unsigned int c = first_comp; c < fe.n_components(); ++c)
159 if (component_mask[c] == true)
160 {
161 local_component_association[i] = c;
162 break;
163 }
164 }
165
166 Assert(std::find(local_component_association.begin(),
167 local_component_association.end(),
168 static_cast<unsigned char>(-1)) ==
169 local_component_association.end(),
171
172 return local_component_association;
173 }
174
175
176 // this internal function assigns to each dof the respective component
177 // of the vector system.
178 //
179 // the output array dofs_by_component lists for each dof the
180 // corresponding vector component. if the DoFHandler is based on a
181 // parallel distributed triangulation then the output array is index by
182 // dof.locally_owned_dofs().index_within_set(indices[i])
183 //
184 // if an element is non-primitive then we assign to each degree of
185 // freedom the following component:
186 // - if the nonzero components that belong to a shape function are not
187 // selected in the component_mask, then the shape function is assigned
188 // to the first nonzero vector component that corresponds to this
189 // shape function
190 // - otherwise, the shape function is assigned the first component selected
191 // in the component_mask that corresponds to this shape function
192 template <int dim, int spacedim>
193 void
195 const ComponentMask & component_mask,
196 std::vector<unsigned char> &dofs_by_component)
197 {
198 const ::hp::FECollection<dim, spacedim> &fe_collection =
199 dof.get_fe_collection();
200 Assert(fe_collection.n_components() < 256, ExcNotImplemented());
201 Assert(dofs_by_component.size() == dof.n_locally_owned_dofs(),
202 ExcDimensionMismatch(dofs_by_component.size(),
203 dof.n_locally_owned_dofs()));
204
205 // next set up a table for the degrees of freedom on each of the
206 // cells (regardless of the fact whether it is listed in the
207 // component_select argument or not)
208 //
209 // for each element 'f' of the FECollection,
210 // local_component_association[f][d] then returns the vector
211 // component that degree of freedom 'd' belongs to
212 std::vector<std::vector<unsigned char>> local_component_association(
213 fe_collection.size());
214 for (unsigned int f = 0; f < fe_collection.size(); ++f)
215 {
216 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
217 local_component_association[f] =
218 get_local_component_association(fe, component_mask);
219 }
220
221 // then loop over all cells and do the work
222 std::vector<types::global_dof_index> indices;
223 for (const auto &c :
225 {
226 const types::fe_index fe_index = c->active_fe_index();
227 const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell();
228 indices.resize(dofs_per_cell);
229 c->get_dof_indices(indices);
230 for (unsigned int i = 0; i < dofs_per_cell; ++i)
231 if (dof.locally_owned_dofs().is_element(indices[i]))
232 dofs_by_component[dof.locally_owned_dofs().index_within_set(
233 indices[i])] = local_component_association[fe_index][i];
234 }
235 }
236
237
238 // this is the function corresponding to the one above but working on
239 // blocks instead of components.
240 //
241 // the output array dofs_by_block lists for each dof the corresponding
242 // vector block. if the DoFHandler is based on a parallel distributed
243 // triangulation then the output array is index by
244 // dof.locally_owned_dofs().index_within_set(indices[i])
245 template <int dim, int spacedim>
246 inline void
248 std::vector<unsigned char> & dofs_by_block)
249 {
250 const ::hp::FECollection<dim, spacedim> &fe_collection =
251 dof.get_fe_collection();
252 Assert(fe_collection.n_components() < 256, ExcNotImplemented());
253 Assert(dofs_by_block.size() == dof.n_locally_owned_dofs(),
254 ExcDimensionMismatch(dofs_by_block.size(),
255 dof.n_locally_owned_dofs()));
256
257 // next set up a table for the degrees of freedom on each of the
258 // cells (regardless of the fact whether it is listed in the
259 // component_select argument or not)
260 //
261 // for each element 'f' of the FECollection,
262 // local_block_association[f][d] then returns the vector block that
263 // degree of freedom 'd' belongs to
264 std::vector<std::vector<unsigned char>> local_block_association(
265 fe_collection.size());
266 for (unsigned int f = 0; f < fe_collection.size(); ++f)
267 {
268 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
269 local_block_association[f].resize(fe.n_dofs_per_cell(),
270 static_cast<unsigned char>(-1));
271 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
272 local_block_association[f][i] = fe.system_to_block_index(i).first;
273
274 Assert(std::find(local_block_association[f].begin(),
275 local_block_association[f].end(),
276 static_cast<unsigned char>(-1)) ==
277 local_block_association[f].end(),
279 }
280
281 // then loop over all cells and do the work
282 std::vector<types::global_dof_index> indices;
283 for (const auto &cell : dof.active_cell_iterators())
284 if (cell->is_locally_owned())
285 {
286 const types::fe_index fe_index = cell->active_fe_index();
287 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
288 indices.resize(dofs_per_cell);
289 cell->get_dof_indices(indices);
290 for (unsigned int i = 0; i < dofs_per_cell; ++i)
291 if (dof.locally_owned_dofs().is_element(indices[i]))
292 dofs_by_block[dof.locally_owned_dofs().index_within_set(
293 indices[i])] = local_block_association[fe_index][i];
294 }
295 }
296 } // namespace internal
297
298
299
300 template <int dim, int spacedim, typename Number>
301 void
303 const Vector<Number> & cell_data,
304 Vector<double> & dof_data,
305 const unsigned int component)
306 {
308 (void)tria;
309
310 AssertDimension(cell_data.size(), tria.n_active_cells());
311 AssertDimension(dof_data.size(), dof_handler.n_dofs());
312 const auto &fe_collection = dof_handler.get_fe_collection();
313 AssertIndexRange(component, fe_collection.n_components());
314 for (unsigned int i = 0; i < fe_collection.size(); ++i)
315 {
316 Assert(fe_collection[i].is_primitive() == true,
318 }
319
320 // store a flag whether we should care about different components. this
321 // is just a simplification, we could ask for this at every single
322 // place equally well
323 const bool consider_components =
324 (dof_handler.get_fe_collection().n_components() != 1);
325
326 // zero out the components that we will touch
327 if (consider_components == false)
328 dof_data = 0;
329 else
330 {
331 std::vector<unsigned char> component_dofs(
332 dof_handler.n_locally_owned_dofs());
334 dof_handler,
335 fe_collection.component_mask(FEValuesExtractors::Scalar(component)),
336 component_dofs);
337
338 for (unsigned int i = 0; i < dof_data.size(); ++i)
339 if (component_dofs[i] == static_cast<unsigned char>(component))
340 dof_data(i) = 0;
341 }
342
343 // count how often we have added a value in the sum for each dof
344 std::vector<unsigned char> touch_count(dof_handler.n_dofs(), 0);
345
346 std::vector<types::global_dof_index> dof_indices;
347 dof_indices.reserve(fe_collection.max_dofs_per_cell());
348
349 for (const auto &cell : dof_handler.active_cell_iterators())
350 {
351 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
352 dof_indices.resize(dofs_per_cell);
353 cell->get_dof_indices(dof_indices);
354
355 for (unsigned int i = 0; i < dofs_per_cell; ++i)
356 // consider this dof only if it is the right component. if there
357 // is only one component, short cut the test
358 if (!consider_components ||
359 (cell->get_fe().system_to_component_index(i).first == component))
360 {
361 // sum up contribution of the present_cell to this dof
362 dof_data(dof_indices[i]) += cell_data(cell->active_cell_index());
363
364 // note that we added another summand
365 ++touch_count[dof_indices[i]];
366 }
367 }
368
369 // compute the mean value on all the dofs by dividing with the number
370 // of summands.
371 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
372 {
373 // assert that each dof was used at least once. this needs not be
374 // the case if the vector has more than one component
375 Assert(consider_components || (touch_count[i] != 0),
377 if (touch_count[i] != 0)
378 dof_data(i) /= touch_count[i];
379 }
380 }
381
382
383
384 template <int dim, int spacedim>
387 const ComponentMask & component_mask)
388 {
389 Assert(component_mask.represents_n_components(
392 "The given component mask is not sized correctly to represent the "
393 "components of the given finite element."));
394
395 // Two special cases: no component is selected, and all components are
396 // selected; both rather stupid, but easy to catch
397 if (component_mask.n_selected_components(
398 dof.get_fe_collection().n_components()) == 0)
399 return IndexSet(dof.n_dofs());
400 else if (component_mask.n_selected_components(
403 return dof.locally_owned_dofs();
404
405 // get the component association of each DoF and then select the ones
406 // that match the given set of components
407 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
408 internal::get_component_association(dof, component_mask, dofs_by_component);
409
410 // fill the selected components in a vector
411 std::vector<types::global_dof_index> selected_dofs;
412 selected_dofs.reserve(dof.n_locally_owned_dofs());
413 for (types::global_dof_index i = 0; i < dofs_by_component.size(); ++i)
414 if (component_mask[dofs_by_component[i]] == true)
415 selected_dofs.push_back(dof.locally_owned_dofs().nth_index_in_set(i));
416
417 // fill vector of indices to return argument
418 IndexSet result(dof.n_dofs());
419 result.add_indices(selected_dofs.begin(), selected_dofs.end());
420 return result;
421 }
422
423
424
425 template <int dim, int spacedim>
428 const BlockMask & block_mask)
429 {
430 // simply forward to the function that works based on a component mask
431 return extract_dofs<dim, spacedim>(
432 dof, dof.get_fe_collection().component_mask(block_mask));
433 }
434
435
436
437 template <int dim, int spacedim>
438 std::vector<IndexSet>
440 const ComponentMask &component_mask)
441 {
442 const auto n_comps = dof.get_fe_collection().n_components();
443 Assert(component_mask.represents_n_components(n_comps),
445 "The given component mask is not sized correctly to represent the "
446 "components of the given finite element."));
447
448 const auto &locally_owned_dofs = dof.locally_owned_dofs();
449
450 // get the component association of each DoF and then select the ones
451 // that match the given set of components
452 std::vector<unsigned char> dofs_by_component(dof.n_locally_owned_dofs());
453 internal::get_component_association(dof, component_mask, dofs_by_component);
454
455 std::vector<IndexSet> index_per_comp(n_comps, IndexSet(dof.n_dofs()));
456
457 for (types::global_dof_index i = 0; i < dof.n_locally_owned_dofs(); ++i)
458 {
459 const auto &comp_i = dofs_by_component[i];
460 if (component_mask[comp_i])
461 index_per_comp[comp_i].add_index(
462 locally_owned_dofs.nth_index_in_set(i));
463 }
464 for (auto &c : index_per_comp)
465 c.compress();
466 return index_per_comp;
467 }
468
469
470
471 template <int dim, int spacedim>
472 void
473 extract_level_dofs(const unsigned int level,
474 const DoFHandler<dim, spacedim> &dof,
475 const ComponentMask & component_mask,
476 std::vector<bool> & selected_dofs)
477 {
478 const FiniteElement<dim, spacedim> &fe = dof.get_fe();
479
480 Assert(component_mask.represents_n_components(
483 "The given component mask is not sized correctly to represent the "
484 "components of the given finite element."));
485 Assert(selected_dofs.size() == dof.n_dofs(level),
486 ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
487
488 // two special cases: no component is selected, and all components are
489 // selected, both rather stupid, but easy to catch
490 if (component_mask.n_selected_components(
491 dof.get_fe_collection().n_components()) == 0)
492 {
493 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
494 return;
495 }
496 else if (component_mask.n_selected_components(
499 {
500 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), true);
501 return;
502 }
503
504 // preset all values by false
505 std::fill_n(selected_dofs.begin(), dof.n_dofs(level), false);
506
507 // next set up a table for the degrees of freedom on each of the cells
508 // whether it is something interesting or not
509 std::vector<unsigned char> local_component_association =
511 std::vector<bool> local_selected_dofs(fe.n_dofs_per_cell());
512 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
513 local_selected_dofs[i] = component_mask[local_component_association[i]];
514
515 // then loop over all cells and do work
516 std::vector<types::global_dof_index> indices(fe.n_dofs_per_cell());
517 for (const auto &cell : dof.cell_iterators_on_level(level))
518 {
519 cell->get_mg_dof_indices(indices);
520 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
521 selected_dofs[indices[i]] = local_selected_dofs[i];
522 }
523 }
524
525
526
527 template <int dim, int spacedim>
528 void
529 extract_level_dofs(const unsigned int level,
530 const DoFHandler<dim, spacedim> &dof,
531 const BlockMask & block_mask,
532 std::vector<bool> & selected_dofs)
533 {
534 // simply defer to the other extract_level_dofs() function
536 dof,
537 dof.get_fe().component_mask(block_mask),
538 selected_dofs);
539 }
540
541
542
543 template <int dim, int spacedim>
544 void
546 const ComponentMask & component_mask,
547 std::vector<bool> & selected_dofs,
548 const std::set<types::boundary_id> &boundary_ids)
549 {
550 Assert((dynamic_cast<
552 &dof_handler.get_triangulation()) == nullptr),
554 "This function can not be used with distributed triangulations. "
555 "See the documentation for more information."));
556
557 IndexSet indices;
558 extract_boundary_dofs(dof_handler, component_mask, indices, boundary_ids);
559
560 // clear and reset array by default values
561 selected_dofs.clear();
562 selected_dofs.resize(dof_handler.n_dofs(), false);
563
564 // then convert the values computed above to the binary vector
565 indices.fill_binary_vector(selected_dofs);
566 }
567
568
569
570 template <int dim, int spacedim>
571 void
573 const ComponentMask & component_mask,
574 IndexSet & selected_dofs,
575 const std::set<types::boundary_id> &boundary_ids)
576 {
577 // Simply forward to the other function
578 selected_dofs =
579 extract_boundary_dofs(dof_handler, component_mask, boundary_ids);
580 }
581
582
583
584 template <int dim, int spacedim>
587 const ComponentMask & component_mask,
588 const std::set<types::boundary_id> &boundary_ids)
589 {
590 Assert(component_mask.represents_n_components(
591 dof_handler.get_fe_collection().n_components()),
592 ExcMessage("Component mask has invalid size."));
593 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
594 boundary_ids.end(),
596
597 IndexSet selected_dofs(dof_handler.n_dofs());
598
599 // let's see whether we have to check for certain boundary indicators
600 // or whether we can accept all
601 const bool check_boundary_id = (boundary_ids.size() != 0);
602
603 // also see whether we have to check whether a certain vector component
604 // is selected, or all
605 const bool check_vector_component =
606 ((component_mask.represents_the_all_selected_mask() == false) ||
607 (component_mask.n_selected_components(
608 dof_handler.get_fe_collection().n_components()) !=
609 dof_handler.get_fe_collection().n_components()));
610
611 std::vector<types::global_dof_index> face_dof_indices;
612 face_dof_indices.reserve(
613 dof_handler.get_fe_collection().max_dofs_per_face());
614
615 // now loop over all cells and check whether their faces are at the
616 // boundary. note that we need not take special care of single lines
617 // being at the boundary (using @p{cell->has_boundary_lines}), since we
618 // do not support boundaries of dimension dim-2, and so every isolated
619 // boundary line is also part of a boundary face which we will be
620 // visiting sooner or later
621 for (const auto &cell : dof_handler.active_cell_iterators())
622 // only work on cells that are either locally owned or at least ghost
623 // cells
624 if (cell->is_artificial() == false)
625 for (const unsigned int face : cell->face_indices())
626 if (cell->at_boundary(face))
627 if (!check_boundary_id ||
628 (boundary_ids.find(cell->face(face)->boundary_id()) !=
629 boundary_ids.end()))
630 {
631 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
632
633 const auto reference_cell = cell->reference_cell();
634
635 const unsigned int n_vertices_per_cell =
636 reference_cell.n_vertices();
637 const unsigned int n_lines_per_cell = reference_cell.n_lines();
638 const unsigned int n_vertices_per_face =
639 reference_cell.face_reference_cell(face).n_vertices();
640 const unsigned int n_lines_per_face =
641 reference_cell.face_reference_cell(face).n_lines();
642
643 const unsigned int dofs_per_face = fe.n_dofs_per_face(face);
644 face_dof_indices.resize(dofs_per_face);
645 cell->face(face)->get_dof_indices(face_dof_indices,
646 cell->active_fe_index());
647
648 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
649 if (!check_vector_component)
650 selected_dofs.add_index(face_dof_indices[i]);
651 else
652 // check for component is required. somewhat tricky as
653 // usual for the case that the shape function is
654 // non-primitive, but use usual convention (see docs)
655 {
656 // first get at the cell-global number of a face dof,
657 // to ask the FE certain questions
658 const unsigned int cell_index =
659 (dim == 1 ?
660 i :
661 (dim == 2 ?
662 (i < 2 * fe.n_dofs_per_vertex() ?
663 i :
664 i + 2 * fe.n_dofs_per_vertex()) :
665 (dim == 3 ? (i < n_vertices_per_face *
666 fe.n_dofs_per_vertex() ?
667 i :
668 (i < n_vertices_per_face *
669 fe.n_dofs_per_vertex() +
670 n_lines_per_face *
671 fe.n_dofs_per_line() ?
672 (i - n_vertices_per_face *
673 fe.n_dofs_per_vertex()) +
674 n_vertices_per_cell *
675 fe.n_dofs_per_vertex() :
676 (i -
677 n_vertices_per_face *
678 fe.n_dofs_per_vertex() -
679 n_lines_per_face *
680 fe.n_dofs_per_line()) +
681 n_vertices_per_cell *
682 fe.n_dofs_per_vertex() +
683 n_lines_per_cell *
684 fe.n_dofs_per_line())) :
686 if (fe.is_primitive(cell_index))
687 {
688 if (component_mask[fe.face_system_to_component_index(
689 i, face)
690 .first] == true)
691 selected_dofs.add_index(face_dof_indices[i]);
692 }
693 else // not primitive
694 {
695 const unsigned int first_nonzero_comp =
698 Assert(first_nonzero_comp < fe.n_components(),
700
701 if (component_mask[first_nonzero_comp] == true)
702 selected_dofs.add_index(face_dof_indices[i]);
703 }
704 }
705 }
706
707 return selected_dofs;
708 }
709
710
711
712 template <int dim, int spacedim>
713 void
715 const DoFHandler<dim, spacedim> & dof_handler,
716 const ComponentMask & component_mask,
717 std::vector<bool> & selected_dofs,
718 const std::set<types::boundary_id> &boundary_ids)
719 {
720 Assert(component_mask.represents_n_components(
721 dof_handler.get_fe_collection().n_components()),
722 ExcMessage("This component mask has the wrong size."));
723 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
724 boundary_ids.end(),
726
727 // let's see whether we have to check for certain boundary indicators
728 // or whether we can accept all
729 const bool check_boundary_id = (boundary_ids.size() != 0);
730
731 // also see whether we have to check whether a certain vector component
732 // is selected, or all
733 const bool check_vector_component =
734 (component_mask.represents_the_all_selected_mask() == false);
735
736 // clear and reset array by default values
737 selected_dofs.clear();
738 selected_dofs.resize(dof_handler.n_dofs(), false);
739 std::vector<types::global_dof_index> cell_dof_indices;
740 cell_dof_indices.reserve(
741 dof_handler.get_fe_collection().max_dofs_per_cell());
742
743 // now loop over all cells and check whether their faces are at the
744 // boundary. note that we need not take special care of single lines
745 // being at the boundary (using @p{cell->has_boundary_lines}), since we
746 // do not support boundaries of dimension dim-2, and so every isolated
747 // boundary line is also part of a boundary face which we will be
748 // visiting sooner or later
749 for (const auto &cell : dof_handler.active_cell_iterators())
750 for (const unsigned int face : cell->face_indices())
751 if (cell->at_boundary(face))
752 if (!check_boundary_id ||
753 (boundary_ids.find(cell->face(face)->boundary_id()) !=
754 boundary_ids.end()))
755 {
756 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
757
758 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
759 cell_dof_indices.resize(dofs_per_cell);
760 cell->get_dof_indices(cell_dof_indices);
761
762 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
763 if (fe.has_support_on_face(i, face))
764 {
765 if (!check_vector_component)
766 selected_dofs[cell_dof_indices[i]] = true;
767 else
768 // check for component is required. somewhat tricky
769 // as usual for the case that the shape function is
770 // non-primitive, but use usual convention (see docs)
771 {
772 if (fe.is_primitive(i))
773 selected_dofs[cell_dof_indices[i]] =
774 (component_mask[fe.system_to_component_index(i)
775 .first] == true);
776 else // not primitive
777 {
778 const unsigned int first_nonzero_comp =
781 Assert(first_nonzero_comp < fe.n_components(),
783
784 selected_dofs[cell_dof_indices[i]] =
785 (component_mask[first_nonzero_comp] == true);
786 }
787 }
788 }
789 }
790 }
791
792
793
794 template <int dim, int spacedim, typename number>
797 const DoFHandler<dim, spacedim> &dof_handler,
798 const std::function<
800 & predicate,
802 {
803 const std::function<bool(
805 predicate_local =
806 [=](
808 -> bool { return cell->is_locally_owned() && predicate(cell); };
809
810 std::vector<types::global_dof_index> local_dof_indices;
811 local_dof_indices.reserve(
812 dof_handler.get_fe_collection().max_dofs_per_cell());
813
814 // Get all the dofs that live on the subdomain:
815 std::set<types::global_dof_index> predicate_dofs;
816
817 for (const auto &cell : dof_handler.active_cell_iterators())
818 if (!cell->is_artificial() && predicate(cell))
819 {
820 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
821 cell->get_dof_indices(local_dof_indices);
822 predicate_dofs.insert(local_dof_indices.begin(),
823 local_dof_indices.end());
824 }
825
826 // Get halo layer and accumulate its DoFs
827 std::set<types::global_dof_index> dofs_with_support_on_halo_cells;
828
829 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
830 halo_cells =
831 GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local);
832 for (typename std::vector<typename DoFHandler<dim, spacedim>::
833 active_cell_iterator>::const_iterator it =
834 halo_cells.begin();
835 it != halo_cells.end();
836 ++it)
837 {
838 // skip halo cells that satisfy the given predicate and thereby
839 // shall be a part of the index set on another MPI core.
840 // Those could only be ghost cells with p::d::Tria.
841 if (predicate(*it))
842 {
843 Assert((*it)->is_ghost(), ExcInternalError());
844 continue;
845 }
846
847 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell();
848 local_dof_indices.resize(dofs_per_cell);
849 (*it)->get_dof_indices(local_dof_indices);
850 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(),
851 local_dof_indices.end());
852 }
853
854 // A complication coming from the fact that DoFs living on locally
855 // owned cells which satisfy predicate may participate in constraints for
856 // DoFs outside of this region.
857 if (cm.n_constraints() > 0)
858 {
859 // Remove DoFs that are part of constraints for DoFs outside
860 // of predicate. Since we will subtract halo_dofs from predicate_dofs,
861 // achieve this by extending halo_dofs with DoFs to which
862 // halo_dofs are constrained.
863 std::set<types::global_dof_index> extra_halo;
864 for (const auto dof : dofs_with_support_on_halo_cells)
865 // if halo DoF is constrained, add all DoFs to which it's constrained
866 // because after resolving constraints, the support of the DoFs that
867 // constrain the current DoF will extend to the halo cells.
868 if (const auto *line_ptr = cm.get_constraint_entries(dof))
869 {
870 const unsigned int line_size = line_ptr->size();
871 for (unsigned int j = 0; j < line_size; ++j)
872 extra_halo.insert((*line_ptr)[j].first);
873 }
874
875 dofs_with_support_on_halo_cells.insert(extra_halo.begin(),
876 extra_halo.end());
877 }
878
879 // Rework our std::set's into IndexSet and subtract halo DoFs from the
880 // predicate_dofs:
881 IndexSet support_set(dof_handler.n_dofs());
882 support_set.add_indices(predicate_dofs.begin(), predicate_dofs.end());
883 support_set.compress();
884
885 IndexSet halo_set(dof_handler.n_dofs());
886 halo_set.add_indices(dofs_with_support_on_halo_cells.begin(),
887 dofs_with_support_on_halo_cells.end());
888 halo_set.compress();
889
890 support_set.subtract_set(halo_set);
891
892 // we intentionally do not want to limit the output to locally owned DoFs.
893 return support_set;
894 }
895
896
897
898 namespace internal
899 {
900 namespace
901 {
902 template <int spacedim>
905 {
906 // there are no hanging nodes in 1d
907 return IndexSet(dof_handler.n_dofs());
908 }
909
910
911 template <int spacedim>
914 {
915 const unsigned int dim = 2;
916
917 IndexSet selected_dofs(dof_handler.n_dofs());
918
919 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
920
921 // this function is similar to the make_sparsity_pattern function,
922 // see there for more information
923 for (const auto &cell : dof_handler.active_cell_iterators())
924 if (!cell->is_artificial())
925 {
926 for (const unsigned int face : cell->face_indices())
927 if (cell->face(face)->has_children())
928 {
930 line = cell->face(face);
931
932 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
933 ++dof)
934 selected_dofs.add_index(
935 line->child(0)->vertex_dof_index(1, dof));
936
937 for (unsigned int child = 0; child < 2; ++child)
938 {
939 if (cell->neighbor_child_on_subface(face, child)
940 ->is_artificial())
941 continue;
942 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
943 ++dof)
944 selected_dofs.add_index(
945 line->child(child)->dof_index(dof));
946 }
947 }
948 }
949
950 selected_dofs.compress();
951 return selected_dofs;
952 }
953
954
955 template <int spacedim>
958 {
959 const unsigned int dim = 3;
960
961 IndexSet selected_dofs(dof_handler.n_dofs());
962 IndexSet unconstrained_dofs(dof_handler.n_dofs());
963
964 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
965
966 for (const auto &cell : dof_handler.active_cell_iterators())
967 if (!cell->is_artificial())
968 for (auto f : cell->face_indices())
969 {
970 const typename DoFHandler<dim, spacedim>::face_iterator face =
971 cell->face(f);
972 if (cell->face(f)->has_children())
973 {
974 for (unsigned int child = 0; child < 4; ++child)
975 if (!cell->neighbor_child_on_subface(f, child)
976 ->is_artificial())
977 {
978 // simply take all DoFs that live on this subface
979 std::vector<types::global_dof_index> ldi(
980 fe.n_dofs_per_face(f, child));
981 face->child(child)->get_dof_indices(ldi);
982 selected_dofs.add_indices(ldi.begin(), ldi.end());
983 }
984
985 // and subtract (in the end) all the indices which a shared
986 // between this face and its subfaces
987 for (unsigned int vertex = 0; vertex < 4; ++vertex)
988 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
989 ++dof)
990 unconstrained_dofs.add_index(
991 face->vertex_dof_index(vertex, dof));
992 }
993 }
994 selected_dofs.subtract_set(unconstrained_dofs);
995 return selected_dofs;
996 }
997 } // namespace
998 } // namespace internal
999
1000
1001
1002 template <int dim, int spacedim>
1003 IndexSet
1005 {
1006 return internal::extract_hanging_node_dofs(dof_handler);
1007 }
1008
1009
1010
1011 template <int dim, int spacedim>
1012 void
1014 const types::subdomain_id subdomain_id,
1015 std::vector<bool> & selected_dofs)
1016 {
1017 Assert(selected_dofs.size() == dof_handler.n_dofs(),
1018 ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
1019
1020 // preset all values by false
1021 std::fill_n(selected_dofs.begin(), dof_handler.n_dofs(), false);
1022
1023 std::vector<types::global_dof_index> local_dof_indices;
1024 local_dof_indices.reserve(
1025 dof_handler.get_fe_collection().max_dofs_per_cell());
1026
1027 // this function is similar to the make_sparsity_pattern function, see
1028 // there for more information
1029 for (const auto &cell : dof_handler.active_cell_iterators())
1030 if (cell->subdomain_id() == subdomain_id)
1031 {
1032 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1033 local_dof_indices.resize(dofs_per_cell);
1034 cell->get_dof_indices(local_dof_indices);
1035 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1036 selected_dofs[local_dof_indices[i]] = true;
1037 }
1038 }
1039
1040
1041
1042 template <int dim, int spacedim>
1043 IndexSet
1045 {
1046 // collect all the locally owned dofs
1047 IndexSet dof_set = dof_handler.locally_owned_dofs();
1048
1049 // add the DoF on the adjacent ghost cells to the IndexSet, cache them
1050 // in a set. need to check each dof manually because we can't be sure
1051 // that the dof range of locally_owned_dofs is really contiguous.
1052 std::vector<types::global_dof_index> dof_indices;
1053 std::set<types::global_dof_index> global_dof_indices;
1054
1055 for (const auto &cell : dof_handler.active_cell_iterators() |
1057 {
1058 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1059 cell->get_dof_indices(dof_indices);
1060
1061 for (const types::global_dof_index dof_index : dof_indices)
1062 if (!dof_set.is_element(dof_index))
1063 global_dof_indices.insert(dof_index);
1064 }
1065
1066 dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
1067
1068 dof_set.compress();
1069
1070 return dof_set;
1071 }
1072
1073
1074
1075 template <int dim, int spacedim>
1076 void
1078 IndexSet & dof_set)
1079 {
1080 dof_set = extract_locally_active_dofs(dof_handler);
1081 }
1082
1083
1084
1085 template <int dim, int spacedim>
1086 IndexSet
1088 const DoFHandler<dim, spacedim> &dof_handler,
1089 const unsigned int level)
1090 {
1091 // collect all the locally owned dofs
1092 IndexSet dof_set = dof_handler.locally_owned_mg_dofs(level);
1093
1094 // add the DoF on the adjacent ghost cells to the IndexSet, cache them
1095 // in a set. need to check each dof manually because we can't be sure
1096 // that the dof range of locally_owned_dofs is really contiguous.
1097 std::vector<types::global_dof_index> dof_indices;
1098 std::set<types::global_dof_index> global_dof_indices;
1099
1100 const auto filtered_iterators_range =
1103 for (const auto &cell : filtered_iterators_range)
1104 {
1105 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1106 cell->get_mg_dof_indices(dof_indices);
1107
1108 for (const types::global_dof_index dof_index : dof_indices)
1109 if (!dof_set.is_element(dof_index))
1110 global_dof_indices.insert(dof_index);
1111 }
1112
1113 dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
1114
1115 dof_set.compress();
1116
1117 return dof_set;
1118 }
1119
1120
1121
1122 template <int dim, int spacedim>
1123 void
1125 const DoFHandler<dim, spacedim> &dof_handler,
1126 IndexSet & dof_set,
1127 const unsigned int level)
1128 {
1129 dof_set = extract_locally_active_level_dofs(dof_handler, level);
1130 }
1131
1132
1133
1134 template <int dim, int spacedim>
1135 IndexSet
1137 {
1138 // collect all the locally owned dofs
1139 IndexSet dof_set = dof_handler.locally_owned_dofs();
1140
1141 // now add the DoF on the adjacent ghost cells to the IndexSet
1142
1143 // Note: For certain meshes (in particular in 3d and with many
1144 // processors), it is really necessary to cache intermediate data. After
1145 // trying several objects such as std::set, a vector that is always kept
1146 // sorted, and a vector that is initially unsorted and sorted once at the
1147 // end, the latter has been identified to provide the best performance.
1148 // Martin Kronbichler
1149 std::vector<types::global_dof_index> dof_indices;
1150 std::vector<types::global_dof_index> dofs_on_ghosts;
1151
1152 for (const auto &cell : dof_handler.active_cell_iterators())
1153 if (cell->is_ghost())
1154 {
1155 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1156 cell->get_dof_indices(dof_indices);
1157 for (const auto dof_index : dof_indices)
1158 if (!dof_set.is_element(dof_index))
1159 dofs_on_ghosts.push_back(dof_index);
1160 }
1161
1162 // sort, compress out duplicates, fill into index set
1163 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1164 dof_set.add_indices(dofs_on_ghosts.begin(),
1165 std::unique(dofs_on_ghosts.begin(),
1166 dofs_on_ghosts.end()));
1167 dof_set.compress();
1168
1169 return dof_set;
1170 }
1171
1172
1173
1174 template <int dim, int spacedim>
1175 void
1177 IndexSet & dof_set)
1178 {
1179 dof_set = extract_locally_relevant_dofs(dof_handler);
1180 }
1181
1182
1183
1184 template <int dim, int spacedim>
1185 IndexSet
1187 const DoFHandler<dim, spacedim> &dof_handler,
1188 const unsigned int level)
1189 {
1190 // collect all the locally owned dofs
1191 IndexSet dof_set = dof_handler.locally_owned_mg_dofs(level);
1192
1193 // add the DoF on the adjacent ghost cells to the IndexSet
1194
1195 // Note: For certain meshes (in particular in 3d and with many
1196 // processors), it is really necessary to cache intermediate data. After
1197 // trying several objects such as std::set, a vector that is always kept
1198 // sorted, and a vector that is initially unsorted and sorted once at the
1199 // end, the latter has been identified to provide the best performance.
1200 // Martin Kronbichler
1201 std::vector<types::global_dof_index> dof_indices;
1202 std::vector<types::global_dof_index> dofs_on_ghosts;
1203
1204 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
1205 {
1206 const types::subdomain_id id = cell->level_subdomain_id();
1207
1208 // skip artificial and own cells (only look at ghost cells)
1209 if (id == dof_handler.get_triangulation().locally_owned_subdomain() ||
1211 continue;
1212
1213 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1214 cell->get_mg_dof_indices(dof_indices);
1215 for (const auto dof_index : dof_indices)
1216 if (!dof_set.is_element(dof_index))
1217 dofs_on_ghosts.push_back(dof_index);
1218 }
1219
1220 // sort, compress out duplicates, fill into index set
1221 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1222 dof_set.add_indices(dofs_on_ghosts.begin(),
1223 std::unique(dofs_on_ghosts.begin(),
1224 dofs_on_ghosts.end()));
1225
1226 dof_set.compress();
1227
1228 return dof_set;
1229 }
1230
1231
1232
1233 template <int dim, int spacedim>
1234 void
1236 const DoFHandler<dim, spacedim> &dof_handler,
1237 const unsigned int level,
1238 IndexSet & dof_set)
1239 {
1240 dof_set = extract_locally_relevant_level_dofs(dof_handler, level);
1241 }
1242
1243
1244
1245 template <int dim, int spacedim>
1246 void
1248 const ComponentMask & component_mask,
1249 std::vector<std::vector<bool>> & constant_modes)
1250 {
1251 // If there are no locally owned DoFs, return with an empty
1252 // constant_modes object:
1253 if (dof_handler.n_locally_owned_dofs() == 0)
1254 {
1255 constant_modes = std::vector<std::vector<bool>>(0);
1256 return;
1257 }
1258
1259 const unsigned int n_components = dof_handler.get_fe(0).n_components();
1260 Assert(component_mask.represents_n_components(n_components),
1261 ExcDimensionMismatch(n_components, component_mask.size()));
1262
1263 std::vector<unsigned char> dofs_by_component(
1264 dof_handler.n_locally_owned_dofs());
1266 component_mask,
1267 dofs_by_component);
1268 unsigned int n_selected_dofs = 0;
1269 for (unsigned int i = 0; i < n_components; ++i)
1270 if (component_mask[i] == true)
1271 n_selected_dofs +=
1272 std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1273
1274 // Find local numbering within the selected components
1275 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
1276 std::vector<unsigned int> component_numbering(
1277 locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int);
1278 for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements();
1279 ++i)
1280 if (component_mask[dofs_by_component[i]])
1281 component_numbering[i] = count++;
1282
1283 // get the element constant modes and find a translation table between
1284 // index in the constant modes and the components.
1285 //
1286 // TODO: We might be able to extend this also for elements which do not
1287 // have the same constant modes, but that is messy...
1288 const ::hp::FECollection<dim, spacedim> &fe_collection =
1289 dof_handler.get_fe_collection();
1290 std::vector<Table<2, bool>> element_constant_modes;
1291 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1292 constant_mode_to_component_translation(n_components);
1293 {
1294 unsigned int n_constant_modes = 0;
1295 int first_non_empty_constant_mode = -1;
1296 for (unsigned int f = 0; f < fe_collection.size(); ++f)
1297 {
1298 std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1299 fe_collection[f].get_constant_modes();
1300
1301 // Store the index of the current element if it is the first that has
1302 // non-empty constant modes.
1303 if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
1304 {
1305 first_non_empty_constant_mode = f;
1306 // This is the first non-empty constant mode, so we figure out the
1307 // translation between index in the constant modes and the
1308 // components
1309 for (unsigned int i = 0; i < data.second.size(); ++i)
1310 if (component_mask[data.second[i]])
1311 constant_mode_to_component_translation[data.second[i]]
1312 .emplace_back(n_constant_modes++, i);
1313 }
1314
1315 // Add the constant modes of this element to the list and assert that
1316 // there are as many constant modes as for the other elements (or zero
1317 // constant modes).
1318 element_constant_modes.push_back(data.first);
1319 Assert(
1320 element_constant_modes.back().n_rows() == 0 ||
1321 element_constant_modes.back().n_rows() ==
1322 element_constant_modes[first_non_empty_constant_mode].n_rows(),
1324 }
1325 AssertIndexRange(first_non_empty_constant_mode, fe_collection.size());
1326
1327 // Now we know the number of constant modes and resize the return vector
1328 // accordingly
1329 constant_modes.clear();
1330 constant_modes.resize(n_constant_modes,
1331 std::vector<bool>(n_selected_dofs, false));
1332 }
1333
1334 // Loop over all owned cells and ask the element for the constant modes
1335 std::vector<types::global_dof_index> dof_indices;
1336 for (const auto &cell : dof_handler.active_cell_iterators() |
1338 {
1339 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1340 cell->get_dof_indices(dof_indices);
1341
1342 for (unsigned int i = 0; i < dof_indices.size(); ++i)
1343 if (locally_owned_dofs.is_element(dof_indices[i]))
1344 {
1345 const unsigned int loc_index =
1346 locally_owned_dofs.index_within_set(dof_indices[i]);
1347 const unsigned int comp = dofs_by_component[loc_index];
1348 if (component_mask[comp])
1349 for (auto &indices :
1350 constant_mode_to_component_translation[comp])
1351 constant_modes[indices
1352 .first][component_numbering[loc_index]] =
1353 element_constant_modes[cell->active_fe_index()](
1354 indices.second, i);
1355 }
1356 }
1357 }
1358
1359
1360
1361 template <int dim, int spacedim>
1362 void
1364 std::vector<unsigned int> & active_fe_indices)
1365 {
1366 AssertDimension(active_fe_indices.size(),
1367 dof_handler.get_triangulation().n_active_cells());
1368
1369 std::vector<types::fe_index> indices = dof_handler.get_active_fe_indices();
1370
1371 active_fe_indices.assign(indices.begin(), indices.end());
1372 }
1373
1374 template <int dim, int spacedim>
1375 std::vector<IndexSet>
1377 {
1378 Assert(dof_handler.n_dofs() > 0,
1379 ExcMessage("The given DoFHandler has no DoFs."));
1380
1381 // If the Triangulation is distributed, the only thing we can usefully
1382 // ask is for its locally owned subdomain
1383 Assert((dynamic_cast<
1385 &dof_handler.get_triangulation()) == nullptr),
1386 ExcMessage(
1387 "For parallel::distributed::Triangulation objects and "
1388 "associated DoF handler objects, asking for any information "
1389 "related to a subdomain other than the locally owned one does "
1390 "not make sense."));
1391
1392 // The following is a random process (flip of a coin), thus should be called
1393 // once only.
1394 std::vector<::types::subdomain_id> subdomain_association(
1395 dof_handler.n_dofs());
1397 subdomain_association);
1398
1399 // Figure out how many subdomain ids there are.
1400 //
1401 // if this is a parallel triangulation, then we can just ask the
1402 // triangulation for this. if this is a sequential triangulation, we loop
1403 // over all cells and take the largest subdomain_id value we find; the
1404 // number of subdomains is then the largest found value plus one. (we here
1405 // assume that all subdomain ids up to the largest are actually used; this
1406 // may not be true for a sequential triangulation where these values have
1407 // been set by hand and not in accordance with some MPI communicator; but
1408 // the function returns an array indexed starting at zero, so we need to
1409 // collect information for each subdomain index anyway, not just for the
1410 // used one.)
1411 const unsigned int n_subdomains =
1412 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1413 &dof_handler.get_triangulation()) == nullptr ?
1414 [&dof_handler]() {
1415 unsigned int max_subdomain_id = 0;
1416 for (const auto &cell : dof_handler.active_cell_iterators())
1417 max_subdomain_id =
1418 std::max(max_subdomain_id, cell->subdomain_id());
1419 return max_subdomain_id + 1;
1420 }() :
1422 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1423 &dof_handler.get_triangulation())
1424 ->get_communicator()));
1425 Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1426 subdomain_association.end()),
1428
1429 std::vector<::IndexSet> index_sets(
1430 n_subdomains, ::IndexSet(dof_handler.n_dofs()));
1431
1432 // loop over subdomain_association and populate IndexSet when a
1433 // change in subdomain ID is found
1434 ::types::global_dof_index i_min = 0;
1435 ::types::global_dof_index this_subdomain = subdomain_association[0];
1436
1437 for (::types::global_dof_index index = 1;
1438 index < subdomain_association.size();
1439 ++index)
1440 {
1441 // found index different from the current one
1442 if (subdomain_association[index] != this_subdomain)
1443 {
1444 index_sets[this_subdomain].add_range(i_min, index);
1445 i_min = index;
1446 this_subdomain = subdomain_association[index];
1447 }
1448 }
1449
1450 // the very last element is of different index
1451 if (i_min == subdomain_association.size() - 1)
1452 {
1453 index_sets[this_subdomain].add_index(i_min);
1454 }
1455
1456 // otherwise there are at least two different indices
1457 else
1458 {
1459 index_sets[this_subdomain].add_range(i_min,
1460 subdomain_association.size());
1461 }
1462
1463 for (unsigned int i = 0; i < n_subdomains; ++i)
1464 index_sets[i].compress();
1465
1466 return index_sets;
1467 }
1468
1469 template <int dim, int spacedim>
1470 std::vector<IndexSet>
1472 const DoFHandler<dim, spacedim> &dof_handler)
1473 {
1474 // If the Triangulation is distributed, the only thing we can usefully
1475 // ask is for its locally owned subdomain
1476 Assert((dynamic_cast<
1478 &dof_handler.get_triangulation()) == nullptr),
1479 ExcMessage(
1480 "For parallel::distributed::Triangulation objects and "
1481 "associated DoF handler objects, asking for any information "
1482 "related to a subdomain other than the locally owned one does "
1483 "not make sense."));
1484
1485 // Collect all the locally owned DoFs
1486 // Note: Even though the distribution of DoFs by the
1487 // locally_owned_dofs_per_subdomain function is pseudo-random, we will
1488 // collect all the DoFs on the subdomain and its layer cell. Therefore, the
1489 // random nature of this function does not play a role in the extraction of
1490 // the locally relevant DoFs
1491 std::vector<IndexSet> dof_set =
1493 const ::types::subdomain_id n_subdomains = dof_set.size();
1494
1495 // Add the DoFs on the adjacent (equivalent ghost) cells to the IndexSet,
1496 // cache them in a set. Need to check each DoF manually because we can't
1497 // be sure that the DoF range of locally_owned_dofs is really contiguous.
1498 for (::types::subdomain_id subdomain_id = 0;
1499 subdomain_id < n_subdomains;
1500 ++subdomain_id)
1501 {
1502 // Extract the layer of cells around this subdomain
1503 std::function<bool(
1505 predicate = IteratorFilters::SubdomainEqualTo(subdomain_id);
1506 const std::vector<
1508 active_halo_layer =
1509 GridTools::compute_active_cell_halo_layer(dof_handler, predicate);
1510
1511 // Extract DoFs associated with halo layer
1512 std::vector<types::global_dof_index> local_dof_indices;
1513 std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1514 for (typename std::vector<
1516 const_iterator it_cell = active_halo_layer.begin();
1517 it_cell != active_halo_layer.end();
1518 ++it_cell)
1519 {
1521 &cell = *it_cell;
1522 Assert(
1523 cell->subdomain_id() != subdomain_id,
1524 ExcMessage(
1525 "The subdomain ID of the halo cell should not match that of the vector entry."));
1526
1527 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1528 cell->get_dof_indices(local_dof_indices);
1529
1530 for (const types::global_dof_index local_dof_index :
1531 local_dof_indices)
1532 subdomain_halo_global_dof_indices.insert(local_dof_index);
1533 }
1534
1535 dof_set[subdomain_id].add_indices(
1536 subdomain_halo_global_dof_indices.begin(),
1537 subdomain_halo_global_dof_indices.end());
1538
1539 dof_set[subdomain_id].compress();
1540 }
1541
1542 return dof_set;
1543 }
1544
1545 template <int dim, int spacedim>
1546 void
1548 const DoFHandler<dim, spacedim> & dof_handler,
1549 std::vector<types::subdomain_id> &subdomain_association)
1550 {
1551 // if the Triangulation is distributed, the only thing we can usefully
1552 // ask is for its locally owned subdomain
1553 Assert((dynamic_cast<
1555 &dof_handler.get_triangulation()) == nullptr),
1556 ExcMessage(
1557 "For parallel::distributed::Triangulation objects and "
1558 "associated DoF handler objects, asking for any subdomain other "
1559 "than the locally owned one does not make sense."));
1560
1561 Assert(subdomain_association.size() == dof_handler.n_dofs(),
1562 ExcDimensionMismatch(subdomain_association.size(),
1563 dof_handler.n_dofs()));
1564
1565 // catch an error that happened in some versions of the shared tria
1566 // distribute_dofs() function where we were trying to call this
1567 // function at a point in time when not all internal DoFHandler
1568 // structures were quite set up yet.
1569 Assert(dof_handler.n_dofs() > 0, ExcInternalError());
1570
1571 // In case this function is executed with parallel::shared::Triangulation
1572 // with possibly artificial cells, we need to take "true" subdomain IDs
1573 // (i.e. without artificial cells). Otherwise we are good to use
1574 // subdomain_id as stored in cell->subdomain_id().
1575 std::vector<types::subdomain_id> cell_owners(
1576 dof_handler.get_triangulation().n_active_cells());
1579 &dof_handler.get_triangulation())))
1580 {
1581 cell_owners = tr->get_true_subdomain_ids_of_cells();
1582 Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1583 tr->n_active_cells(),
1585 }
1586 else
1587 {
1588 for (const auto &cell : dof_handler.active_cell_iterators() |
1590 cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1591 }
1592
1593 // preset all values by an invalid value
1594 std::fill_n(subdomain_association.begin(),
1595 dof_handler.n_dofs(),
1597
1598 std::vector<types::global_dof_index> local_dof_indices;
1599 local_dof_indices.reserve(
1600 dof_handler.get_fe_collection().max_dofs_per_cell());
1601
1602 // loop over all cells and record which subdomain a DoF belongs to.
1603 // give to the smaller subdomain_id in case it is on an interface
1604 for (const auto &cell : dof_handler.active_cell_iterators())
1605 {
1606 const types::subdomain_id subdomain_id =
1607 cell_owners[cell->active_cell_index()];
1608 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1609 local_dof_indices.resize(dofs_per_cell);
1610 cell->get_dof_indices(local_dof_indices);
1611
1612 // set subdomain ids. if dofs already have their values set then
1613 // they must be on partition interfaces. in that case assign them
1614 // to either the previous association or the current processor
1615 // with the smaller subdomain id.
1616 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1617 if (subdomain_association[local_dof_indices[i]] ==
1619 subdomain_association[local_dof_indices[i]] = subdomain_id;
1620 else if (subdomain_association[local_dof_indices[i]] > subdomain_id)
1621 {
1622 subdomain_association[local_dof_indices[i]] = subdomain_id;
1623 }
1624 }
1625
1626 Assert(std::find(subdomain_association.begin(),
1627 subdomain_association.end(),
1629 subdomain_association.end(),
1631 }
1632
1633
1634
1635 template <int dim, int spacedim>
1636 unsigned int
1638 const DoFHandler<dim, spacedim> &dof_handler,
1639 const types::subdomain_id subdomain)
1640 {
1641 std::vector<types::subdomain_id> subdomain_association(
1642 dof_handler.n_dofs());
1643 get_subdomain_association(dof_handler, subdomain_association);
1644
1645 return std::count(subdomain_association.begin(),
1646 subdomain_association.end(),
1647 subdomain);
1648 }
1649
1650
1651
1652 template <int dim, int spacedim>
1653 IndexSet
1655 const DoFHandler<dim, spacedim> &dof_handler,
1656 const types::subdomain_id subdomain)
1657 {
1658 // If we have a distributed::Triangulation only allow locally_owned
1659 // subdomain.
1662 (subdomain ==
1664 ExcMessage(
1665 "For parallel::distributed::Triangulation objects and "
1666 "associated DoF handler objects, asking for any subdomain other "
1667 "than the locally owned one does not make sense."));
1668
1669 IndexSet index_set(dof_handler.n_dofs());
1670
1671 std::vector<types::global_dof_index> local_dof_indices;
1672 local_dof_indices.reserve(
1673 dof_handler.get_fe_collection().max_dofs_per_cell());
1674
1675 // first generate an unsorted list of all indices which we fill from
1676 // the back. could also insert them directly into the IndexSet, but
1677 // that inserts indices in the middle, which is an O(n^2) algorithm and
1678 // hence too expensive. Could also use std::set, but that is in general
1679 // more expensive than a vector
1680 std::vector<types::global_dof_index> subdomain_indices;
1681
1682 for (const auto &cell : dof_handler.active_cell_iterators())
1683 if ((cell->is_artificial() == false) &&
1684 (cell->subdomain_id() == subdomain))
1685 {
1686 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1687 local_dof_indices.resize(dofs_per_cell);
1688 cell->get_dof_indices(local_dof_indices);
1689 subdomain_indices.insert(subdomain_indices.end(),
1690 local_dof_indices.begin(),
1691 local_dof_indices.end());
1692 }
1693 // sort indices and remove duplicates
1694 std::sort(subdomain_indices.begin(), subdomain_indices.end());
1695 subdomain_indices.erase(std::unique(subdomain_indices.begin(),
1696 subdomain_indices.end()),
1697 subdomain_indices.end());
1698
1699 // insert into IndexSet
1700 index_set.add_indices(subdomain_indices.begin(), subdomain_indices.end());
1701 index_set.compress();
1702
1703 return index_set;
1704 }
1705
1706
1707
1708 template <int dim, int spacedim>
1709 void
1711 const DoFHandler<dim, spacedim> &dof_handler,
1712 const types::subdomain_id subdomain,
1713 std::vector<unsigned int> & n_dofs_on_subdomain)
1714 {
1715 Assert(n_dofs_on_subdomain.size() == dof_handler.get_fe(0).n_components(),
1716 ExcDimensionMismatch(n_dofs_on_subdomain.size(),
1717 dof_handler.get_fe(0).n_components()));
1718 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1719
1720 // Make sure there are at least some cells with this subdomain id
1721 Assert(std::any_of(
1722 dof_handler.begin_active(),
1724 dof_handler.end()},
1725 [subdomain](
1726 const typename DoFHandler<dim, spacedim>::cell_accessor &cell) {
1727 return cell.subdomain_id() == subdomain;
1728 }),
1729 ExcMessage("There are no cells for the given subdomain!"));
1730
1731 std::vector<types::subdomain_id> subdomain_association(
1732 dof_handler.n_dofs());
1733 get_subdomain_association(dof_handler, subdomain_association);
1734
1735 std::vector<unsigned char> component_association(dof_handler.n_dofs());
1737 std::vector<bool>(),
1738 component_association);
1739
1740 for (unsigned int c = 0; c < dof_handler.get_fe(0).n_components(); ++c)
1741 {
1742 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
1743 if ((subdomain_association[i] == subdomain) &&
1744 (component_association[i] == static_cast<unsigned char>(c)))
1745 ++n_dofs_on_subdomain[c];
1746 }
1747 }
1748
1749
1750
1751 namespace internal
1752 {
1753 // TODO: why is this function so complicated? It would be nice to have
1754 // comments that explain why we can't just loop over all components and
1755 // count the entries in dofs_by_component that have this component's
1756 // index
1757 template <int dim, int spacedim>
1758 void
1760 const std::vector<unsigned char> & dofs_by_component,
1761 const std::vector<unsigned int> & target_component,
1762 const bool only_once,
1763 std::vector<types::global_dof_index> &dofs_per_component,
1764 unsigned int & component)
1765 {
1766 for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1767 {
1768 const FiniteElement<dim, spacedim> &base = fe.base_element(b);
1769 // Dimension of base element
1770 unsigned int d = base.n_components();
1771
1772 for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1773 {
1774 if (base.n_base_elements() > 1)
1775 resolve_components(base,
1776 dofs_by_component,
1777 target_component,
1778 only_once,
1779 dofs_per_component,
1780 component);
1781 else
1782 {
1783 for (unsigned int dd = 0; dd < d; ++dd, ++component)
1784 dofs_per_component[target_component[component]] +=
1785 std::count(dofs_by_component.begin(),
1786 dofs_by_component.end(),
1787 component);
1788
1789 // if we have non-primitive FEs and want all components
1790 // to show the number of dofs, need to copy the result to
1791 // those components
1792 if (!base.is_primitive() && !only_once)
1793 for (unsigned int dd = 1; dd < d; ++dd)
1794 dofs_per_component[target_component[component - d + dd]] =
1795 dofs_per_component[target_component[component - d]];
1796 }
1797 }
1798 }
1799 }
1800
1801
1802 template <int dim, int spacedim>
1803 void
1805 const std::vector<unsigned char> & dofs_by_component,
1806 const std::vector<unsigned int> & target_component,
1807 const bool only_once,
1808 std::vector<types::global_dof_index> &dofs_per_component,
1809 unsigned int & component)
1810 {
1811 // assert that all elements in the collection have the same structure
1812 // (base elements and multiplicity, components per base element) and
1813 // then simply call the function above
1814 for (unsigned int fe = 1; fe < fe_collection.size(); ++fe)
1815 {
1816 Assert(fe_collection[fe].n_components() ==
1817 fe_collection[0].n_components(),
1819 Assert(fe_collection[fe].n_base_elements() ==
1820 fe_collection[0].n_base_elements(),
1822 for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b)
1823 {
1824 Assert(fe_collection[fe].base_element(b).n_components() ==
1825 fe_collection[0].base_element(b).n_components(),
1827 Assert(fe_collection[fe].base_element(b).n_base_elements() ==
1828 fe_collection[0].base_element(b).n_base_elements(),
1830 }
1831 }
1832
1833 resolve_components(fe_collection[0],
1834 dofs_by_component,
1835 target_component,
1836 only_once,
1837 dofs_per_component,
1838 component);
1839 }
1840 } // namespace internal
1841
1842
1843
1844 namespace internal
1845 {
1846 namespace
1847 {
1851 template <int dim, int spacedim>
1852 bool
1853 all_elements_are_primitive(const FiniteElement<dim, spacedim> &fe)
1854 {
1855 return fe.is_primitive();
1856 }
1857
1858
1863 template <int dim, int spacedim>
1864 bool
1865 all_elements_are_primitive(
1866 const ::hp::FECollection<dim, spacedim> &fe_collection)
1867 {
1868 for (unsigned int i = 0; i < fe_collection.size(); ++i)
1869 if (fe_collection[i].is_primitive() == false)
1870 return false;
1871
1872 return true;
1873 }
1874 } // namespace
1875 } // namespace internal
1876
1877
1878
1879 template <int dim, int spacedim>
1880 std::vector<types::global_dof_index>
1882 const DoFHandler<dim, spacedim> &dof_handler,
1883 const bool only_once,
1884 const std::vector<unsigned int> &target_component_)
1885 {
1886 const unsigned int n_components = dof_handler.get_fe(0).n_components();
1887
1888 // If the empty vector was given as default argument, set up this
1889 // vector as identity.
1890 std::vector<unsigned int> target_component = target_component_;
1891 if (target_component.size() == 0)
1892 {
1893 target_component.resize(n_components);
1894 for (unsigned int i = 0; i < n_components; ++i)
1895 target_component[i] = i;
1896 }
1897 else
1898 Assert(target_component.size() == n_components,
1899 ExcDimensionMismatch(target_component.size(), n_components));
1900
1901
1902 const unsigned int max_component =
1903 *std::max_element(target_component.begin(), target_component.end());
1904 const unsigned int n_target_components = max_component + 1;
1905
1906 std::vector<types::global_dof_index> dofs_per_component(
1907 n_target_components, types::global_dof_index(0));
1908
1909 // special case for only one component. treat this first since it does
1910 // not require any computations
1911 if (n_components == 1)
1912 {
1913 dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
1914 return dofs_per_component;
1915 }
1916
1917
1918 // otherwise determine the number of dofs in each component separately.
1919 // do so in parallel
1920 std::vector<unsigned char> dofs_by_component(
1921 dof_handler.n_locally_owned_dofs());
1923 ComponentMask(),
1924 dofs_by_component);
1925
1926 // next count what we got
1927 unsigned int component = 0;
1929 dofs_by_component,
1930 target_component,
1931 only_once,
1932 dofs_per_component,
1933 component);
1934 Assert(n_components == component, ExcInternalError());
1935
1936 // finally sanity check. this is only valid if the finite element is
1937 // actually primitive, so exclude other elements from this
1938 Assert((internal::all_elements_are_primitive(
1939 dof_handler.get_fe_collection()) == false) ||
1940 (std::accumulate(dofs_per_component.begin(),
1941 dofs_per_component.end(),
1943 dof_handler.n_locally_owned_dofs()),
1945
1946 // reduce information from all CPUs
1947#ifdef DEAL_II_WITH_MPI
1948
1950 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1951 &dof_handler.get_triangulation())))
1952 {
1953 std::vector<types::global_dof_index> local_dof_count =
1954 dofs_per_component;
1955
1956 const int ierr = MPI_Allreduce(local_dof_count.data(),
1957 dofs_per_component.data(),
1958 n_target_components,
1960 MPI_SUM,
1962 AssertThrowMPI(ierr);
1963 }
1964#endif
1965
1966 return dofs_per_component;
1967 }
1968
1969
1970
1971 template <int dim, int spacedim>
1972 std::vector<types::global_dof_index>
1974 const std::vector<unsigned int> &target_block_)
1975 {
1976 const ::hp::FECollection<dim, spacedim> &fe_collection =
1977 dof_handler.get_fe_collection();
1978 Assert(fe_collection.size() < 256, ExcNotImplemented());
1979
1980 // If the empty vector for target_block(e.g., as default argument), then
1981 // set up this vector as identity. We do this set up with the first
1982 // element of the collection, but the whole thing can only work if
1983 // all elements have the same number of blocks anyway -- so check
1984 // that right after
1985 const unsigned int n_blocks = fe_collection[0].n_blocks();
1986
1987 std::vector<unsigned int> target_block = target_block_;
1988 if (target_block.size() == 0)
1989 {
1990 target_block.resize(fe_collection[0].n_blocks());
1991 for (unsigned int i = 0; i < n_blocks; ++i)
1992 target_block[i] = i;
1993 }
1994 else
1995 Assert(target_block.size() == n_blocks,
1996 ExcDimensionMismatch(target_block.size(), n_blocks));
1997 for (unsigned int f = 1; f < fe_collection.size(); ++f)
1998 Assert(fe_collection[0].n_blocks() == fe_collection[f].n_blocks(),
1999 ExcMessage("This function can only work if all elements in a "
2000 "collection have the same number of blocks."));
2001
2002 // special case for only one block. treat this first since it does
2003 // not require any computations
2004 if (n_blocks == 1)
2005 {
2006 std::vector<types::global_dof_index> dofs_per_block(1);
2007 dofs_per_block[0] = dof_handler.n_dofs();
2008 return dofs_per_block;
2009 }
2010
2011 // Otherwise set up the right-sized object and start working
2012 const unsigned int max_block =
2013 *std::max_element(target_block.begin(), target_block.end());
2014 const unsigned int n_target_blocks = max_block + 1;
2015
2016 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2017
2018 // Loop over the elements of the collection, but really only consider
2019 // the last element (see #9271)
2020 for (unsigned int this_fe = fe_collection.size() - 1;
2021 this_fe < fe_collection.size();
2022 ++this_fe)
2023 {
2024 const FiniteElement<dim, spacedim> &fe = fe_collection[this_fe];
2025
2026 std::vector<unsigned char> dofs_by_block(
2027 dof_handler.n_locally_owned_dofs());
2028 internal::get_block_association(dof_handler, dofs_by_block);
2029
2030 // next count what we got
2031 for (unsigned int block = 0; block < fe.n_blocks(); ++block)
2032 dofs_per_block[target_block[block]] +=
2033 std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2034
2035#ifdef DEAL_II_WITH_MPI
2036 // if we are working on a parallel mesh, we now need to collect
2037 // this information from all processors
2039 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2040 &dof_handler.get_triangulation())))
2041 {
2042 std::vector<types::global_dof_index> local_dof_count =
2043 dofs_per_block;
2044 const int ierr = MPI_Allreduce(local_dof_count.data(),
2045 dofs_per_block.data(),
2046 n_target_blocks,
2048 MPI_SUM,
2050 AssertThrowMPI(ierr);
2051 }
2052#endif
2053 }
2054
2055 return dofs_per_block;
2056 }
2057
2058
2059
2060 template <int dim, int spacedim>
2061 void
2063 std::vector<types::global_dof_index> &mapping)
2064 {
2065 mapping.clear();
2066 mapping.insert(mapping.end(),
2067 dof_handler.n_dofs(),
2069
2070 std::vector<types::global_dof_index> dofs_on_face;
2071 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2072 types::global_dof_index next_boundary_index = 0;
2073
2074 // now loop over all cells and check whether their faces are at the
2075 // boundary. note that we need not take special care of single lines
2076 // being at the boundary (using @p{cell->has_boundary_lines}), since we
2077 // do not support boundaries of dimension dim-2, and so every isolated
2078 // boundary line is also part of a boundary face which we will be
2079 // visiting sooner or later
2080 for (const auto &cell : dof_handler.active_cell_iterators())
2081 for (const unsigned int f : cell->face_indices())
2082 if (cell->at_boundary(f))
2083 {
2084 const unsigned int dofs_per_face =
2085 cell->get_fe().n_dofs_per_face(f);
2086 dofs_on_face.resize(dofs_per_face);
2087 cell->face(f)->get_dof_indices(dofs_on_face,
2088 cell->active_fe_index());
2089 for (unsigned int i = 0; i < dofs_per_face; ++i)
2090 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2091 mapping[dofs_on_face[i]] = next_boundary_index++;
2092 }
2093
2094 AssertDimension(next_boundary_index, dof_handler.n_boundary_dofs());
2095 }
2096
2097
2098
2099 template <int dim, int spacedim>
2100 void
2102 const std::set<types::boundary_id> &boundary_ids,
2103 std::vector<types::global_dof_index> &mapping)
2104 {
2105 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2106 boundary_ids.end(),
2108
2109 mapping.clear();
2110 mapping.insert(mapping.end(),
2111 dof_handler.n_dofs(),
2113
2114 // return if there is nothing to do
2115 if (boundary_ids.size() == 0)
2116 return;
2117
2118 std::vector<types::global_dof_index> dofs_on_face;
2119 dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2120 types::global_dof_index next_boundary_index = 0;
2121
2122 for (const auto &cell : dof_handler.active_cell_iterators())
2123 for (const unsigned int f : cell->face_indices())
2124 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2125 boundary_ids.end())
2126 {
2127 const unsigned int dofs_per_face =
2128 cell->get_fe().n_dofs_per_face(f);
2129 dofs_on_face.resize(dofs_per_face);
2130 cell->face(f)->get_dof_indices(dofs_on_face,
2131 cell->active_fe_index());
2132 for (unsigned int i = 0; i < dofs_per_face; ++i)
2133 if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2134 mapping[dofs_on_face[i]] = next_boundary_index++;
2135 }
2136
2137 AssertDimension(next_boundary_index,
2138 dof_handler.n_boundary_dofs(boundary_ids));
2139 }
2140
2141 namespace internal
2142 {
2143 namespace
2144 {
2145 template <int dim, int spacedim>
2146 std::map<types::global_dof_index, Point<spacedim>>
2149 const DoFHandler<dim, spacedim> & dof_handler,
2150 const ComponentMask & in_mask)
2151 {
2152 std::map<types::global_dof_index, Point<spacedim>> support_points;
2153
2154 const hp::FECollection<dim, spacedim> &fe_collection =
2155 dof_handler.get_fe_collection();
2156 hp::QCollection<dim> q_coll_dummy;
2157
2158 for (unsigned int fe_index = 0; fe_index < fe_collection.size();
2159 ++fe_index)
2160 {
2161 // check whether every FE in the collection has support points
2162 Assert((fe_collection[fe_index].n_dofs_per_cell() == 0) ||
2163 (fe_collection[fe_index].has_support_points()),
2165 q_coll_dummy.push_back(Quadrature<dim>(
2166 fe_collection[fe_index].get_unit_support_points()));
2167 }
2168
2169 // Take care of components
2170 const ComponentMask mask =
2171 (in_mask.size() == 0 ?
2172 ComponentMask(fe_collection.n_components(), true) :
2173 in_mask);
2174
2175 // Now loop over all cells and enquire the support points on each
2176 // of these. we use dummy quadrature formulas where the quadrature
2177 // points are located at the unit support points to enquire the
2178 // location of the support points in real space.
2179 //
2180 // The weights of the quadrature rule have been set to invalid
2181 // values by the used constructor.
2182 hp::FEValues<dim, spacedim> hp_fe_values(mapping,
2183 fe_collection,
2184 q_coll_dummy,
2186
2187 std::vector<types::global_dof_index> local_dof_indices;
2188 for (const auto &cell : dof_handler.active_cell_iterators())
2189 // only work on locally relevant cells
2190 if (cell->is_artificial() == false)
2191 {
2192 hp_fe_values.reinit(cell);
2193 const FEValues<dim, spacedim> &fe_values =
2194 hp_fe_values.get_present_fe_values();
2195
2196 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2197 cell->get_dof_indices(local_dof_indices);
2198
2199 const std::vector<Point<spacedim>> &points =
2200 fe_values.get_quadrature_points();
2201 for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2202 ++i)
2203 {
2204 const unsigned int dof_comp =
2205 cell->get_fe().system_to_component_index(i).first;
2206
2207 // insert the values into the map if it is a valid component
2208 if (mask[dof_comp])
2209 support_points[local_dof_indices[i]] = points[i];
2210 }
2211 }
2212
2213 return support_points;
2214 }
2215
2216
2217 template <int dim, int spacedim>
2218 std::vector<Point<spacedim>>
2219 map_dofs_to_support_points_vector(
2221 const DoFHandler<dim, spacedim> & dof_handler,
2222 const ComponentMask & mask)
2223 {
2224 std::vector<Point<spacedim>> support_points(dof_handler.n_dofs());
2225
2226 // get the data in the form of the map as above
2227 const std::map<types::global_dof_index, Point<spacedim>>
2228 x_support_points =
2229 map_dofs_to_support_points(mapping, dof_handler, mask);
2230
2231 // now convert from the map to the linear vector. make sure every
2232 // entry really appeared in the map
2233 for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
2234 {
2235 Assert(x_support_points.find(i) != x_support_points.end(),
2237
2238 support_points[i] = x_support_points.find(i)->second;
2239 }
2240
2241 return support_points;
2242 }
2243 } // namespace
2244 } // namespace internal
2245
2246
2247 template <int dim, int spacedim>
2248 void
2250 const DoFHandler<dim, spacedim> &dof_handler,
2251 std::vector<Point<spacedim>> & support_points,
2252 const ComponentMask & mask)
2253 {
2254 AssertDimension(support_points.size(), dof_handler.n_dofs());
2255 Assert((dynamic_cast<
2257 &dof_handler.get_triangulation()) == nullptr),
2258 ExcMessage(
2259 "This function can not be used with distributed triangulations. "
2260 "See the documentation for more information."));
2261
2262 // Let the internal function do all the work, just make sure that it
2263 // gets a MappingCollection
2264 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2265
2266 support_points =
2267 internal::map_dofs_to_support_points_vector(mapping_collection,
2268 dof_handler,
2269 mask);
2270 }
2271
2272
2273 template <int dim, int spacedim>
2274 void
2277 const DoFHandler<dim, spacedim> & dof_handler,
2278 std::vector<Point<spacedim>> & support_points,
2279 const ComponentMask & mask)
2280 {
2281 AssertDimension(support_points.size(), dof_handler.n_dofs());
2282 Assert((dynamic_cast<
2284 &dof_handler.get_triangulation()) == nullptr),
2285 ExcMessage(
2286 "This function can not be used with distributed triangulations. "
2287 "See the documentation for more information."));
2288
2289 // Let the internal function do all the work, just make sure that it
2290 // gets a MappingCollection
2291 support_points =
2292 internal::map_dofs_to_support_points_vector(mapping, dof_handler, mask);
2293 }
2294
2295
2296 // This function is deprecated:
2297 template <int dim, int spacedim>
2298 DEAL_II_DEPRECATED_EARLY void
2300 const Mapping<dim, spacedim> & mapping,
2301 const DoFHandler<dim, spacedim> & dof_handler,
2302 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2303 const ComponentMask & mask)
2304 {
2305 support_points.clear();
2306
2307 // Let the internal function do all the work, just make sure that it
2308 // gets a MappingCollection
2309 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2310
2311 support_points = internal::map_dofs_to_support_points(mapping_collection,
2312 dof_handler,
2313 mask);
2314 }
2315
2316
2317 // This function is deprecated:
2318 template <int dim, int spacedim>
2319 DEAL_II_DEPRECATED_EARLY void
2322 const DoFHandler<dim, spacedim> & dof_handler,
2323 std::map<types::global_dof_index, Point<spacedim>> &support_points,
2324 const ComponentMask & mask)
2325 {
2326 support_points.clear();
2327
2328 // Let the internal function do all the work, just make sure that it
2329 // gets a MappingCollection
2330 support_points =
2331 internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2332 }
2333
2334
2335 template <int dim, int spacedim>
2336 std::map<types::global_dof_index, Point<spacedim>>
2338 const DoFHandler<dim, spacedim> &dof_handler,
2339 const ComponentMask & mask)
2340 {
2341 // Let the internal function do all the work, just make sure that it
2342 // gets a MappingCollection
2343 const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2344
2345 return internal::map_dofs_to_support_points(mapping_collection,
2346 dof_handler,
2347 mask);
2348 }
2349
2350
2351 template <int dim, int spacedim>
2352 std::map<types::global_dof_index, Point<spacedim>>
2355 const DoFHandler<dim, spacedim> & dof_handler,
2356 const ComponentMask & mask)
2357 {
2358 return internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2359 }
2360
2361
2362 template <int spacedim>
2363 void
2365 std::ostream & out,
2366 const std::map<types::global_dof_index, Point<spacedim>> &support_points)
2367 {
2368 AssertThrow(out.fail() == false, ExcIO());
2369
2370 std::map<Point<spacedim>,
2371 std::vector<types::global_dof_index>,
2373 point_map;
2374
2375 // convert to map point -> list of DoFs
2376 for (const auto &it : support_points)
2377 {
2378 std::vector<types::global_dof_index> &v = point_map[it.second];
2379 v.push_back(it.first);
2380 }
2381
2382 // print the newly created map:
2383 for (const auto &it : point_map)
2384 {
2385 out << it.first << " \"";
2386 const std::vector<types::global_dof_index> &v = it.second;
2387 for (unsigned int i = 0; i < v.size(); ++i)
2388 {
2389 if (i > 0)
2390 out << ", ";
2391 out << v[i];
2392 }
2393
2394 out << "\"\n";
2395 }
2396
2397 out << std::flush;
2398 }
2399
2400
2401 template <int dim, int spacedim>
2402 void
2404 const Table<2, Coupling> & table,
2405 std::vector<Table<2, Coupling>> &tables_by_block)
2406 {
2407 if (dof_handler.has_hp_capabilities() == false)
2408 {
2409 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
2410 const unsigned int nb = fe.n_blocks();
2411
2412 tables_by_block.resize(1);
2413 tables_by_block[0].reinit(nb, nb);
2414 tables_by_block[0].fill(none);
2415
2416 for (unsigned int i = 0; i < fe.n_components(); ++i)
2417 {
2418 const unsigned int ib = fe.component_to_block_index(i);
2419 for (unsigned int j = 0; j < fe.n_components(); ++j)
2420 {
2421 const unsigned int jb = fe.component_to_block_index(j);
2422 tables_by_block[0](ib, jb) |= table(i, j);
2423 }
2424 }
2425 }
2426 else
2427 {
2428 const hp::FECollection<dim> &fe_collection =
2429 dof_handler.get_fe_collection();
2430 tables_by_block.resize(fe_collection.size());
2431
2432 for (unsigned int f = 0; f < fe_collection.size(); ++f)
2433 {
2434 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
2435
2436 const unsigned int nb = fe.n_blocks();
2437 tables_by_block[f].reinit(nb, nb);
2438 tables_by_block[f].fill(none);
2439 for (unsigned int i = 0; i < fe.n_components(); ++i)
2440 {
2441 const unsigned int ib = fe.component_to_block_index(i);
2442 for (unsigned int j = 0; j < fe.n_components(); ++j)
2443 {
2444 const unsigned int jb = fe.component_to_block_index(j);
2445 tables_by_block[f](ib, jb) |= table(i, j);
2446 }
2447 }
2448 }
2449 }
2450 }
2451
2452
2453
2454 template <int dim, int spacedim>
2455 void
2457 const DoFHandler<dim, spacedim> &dof_handler,
2458 const unsigned int level,
2459 const std::vector<bool> & selected_dofs,
2460 const types::global_dof_index offset)
2461 {
2462 std::vector<types::global_dof_index> indices;
2463
2464 unsigned int i = 0;
2465
2466 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2467 if (cell->is_locally_owned_on_level())
2468 ++i;
2469 block_list.reinit(i,
2470 dof_handler.n_dofs(),
2471 dof_handler.get_fe().n_dofs_per_cell());
2472 i = 0;
2473 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2474 if (cell->is_locally_owned_on_level())
2475 {
2476 indices.resize(cell->get_fe().n_dofs_per_cell());
2477 cell->get_mg_dof_indices(indices);
2478
2479 if (selected_dofs.size() != 0)
2480 AssertDimension(indices.size(), selected_dofs.size());
2481
2482 for (types::global_dof_index j = 0; j < indices.size(); ++j)
2483 {
2484 if (selected_dofs.size() == 0)
2485 block_list.add(i, indices[j] - offset);
2486 else
2487 {
2488 if (selected_dofs[j])
2489 block_list.add(i, indices[j] - offset);
2490 }
2491 }
2492 ++i;
2493 }
2494 }
2495
2496
2497 template <int dim, int spacedim>
2498 void
2500 const DoFHandler<dim, spacedim> &dof_handler,
2501 const unsigned int level,
2502 const bool interior_only)
2503 {
2504 const FiniteElement<dim> &fe = dof_handler.get_fe();
2505 block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level));
2506
2507 std::vector<types::global_dof_index> indices;
2508 std::vector<bool> exclude;
2509
2510 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2511 {
2512 indices.resize(cell->get_fe().n_dofs_per_cell());
2513 cell->get_mg_dof_indices(indices);
2514
2515 if (interior_only)
2516 {
2517 // Exclude degrees of freedom on faces opposite to the vertex
2518 exclude.resize(fe.n_dofs_per_cell());
2519 std::fill(exclude.begin(), exclude.end(), false);
2520
2521 for (const unsigned int face : cell->face_indices())
2522 if (cell->at_boundary(face) ||
2523 cell->neighbor(face)->level() != cell->level())
2524 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2525 exclude[fe.face_to_cell_index(i, face)] = true;
2526 for (types::global_dof_index j = 0; j < indices.size(); ++j)
2527 if (!exclude[j])
2528 block_list.add(0, indices[j]);
2529 }
2530 else
2531 {
2532 for (const auto index : indices)
2533 block_list.add(0, index);
2534 }
2535 }
2536 }
2537
2538
2539 template <int dim, int spacedim>
2540 void
2542 const DoFHandler<dim, spacedim> &dof_handler,
2543 const unsigned int level,
2544 const bool interior_dofs_only,
2545 const bool boundary_dofs)
2546 {
2547 Assert(level > 0 && level < dof_handler.get_triangulation().n_levels(),
2548 ExcIndexRange(level, 1, dof_handler.get_triangulation().n_levels()));
2549
2550 std::vector<types::global_dof_index> indices;
2551 std::vector<bool> exclude;
2552
2553 unsigned int block = 0;
2554 for (const auto &pcell : dof_handler.cell_iterators_on_level(level - 1))
2555 {
2556 if (pcell->is_active())
2557 continue;
2558
2559 for (unsigned int child = 0; child < pcell->n_children(); ++child)
2560 {
2561 const auto cell = pcell->child(child);
2562
2563 // For hp, only this line here would have to be replaced.
2564 const FiniteElement<dim> &fe = dof_handler.get_fe();
2565 const unsigned int n_dofs = fe.n_dofs_per_cell();
2566 indices.resize(n_dofs);
2567 exclude.resize(n_dofs);
2568 std::fill(exclude.begin(), exclude.end(), false);
2569 cell->get_mg_dof_indices(indices);
2570
2571 if (interior_dofs_only)
2572 {
2573 // Eliminate dofs on faces of the child which are on faces
2574 // of the parent
2575 for (unsigned int d = 0; d < dim; ++d)
2576 {
2577 const unsigned int face =
2579 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2580 exclude[fe.face_to_cell_index(i, face)] = true;
2581 }
2582
2583 // Now remove all degrees of freedom on the domain boundary
2584 // from the exclusion list
2585 if (boundary_dofs)
2586 for (const unsigned int face :
2588 if (cell->at_boundary(face))
2589 for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
2590 ++i)
2591 exclude[fe.face_to_cell_index(i, face)] = false;
2592 }
2593
2594 for (unsigned int i = 0; i < n_dofs; ++i)
2595 if (!exclude[i])
2596 block_list.add(block, indices[i]);
2597 }
2598 ++block;
2599 }
2600 }
2601
2602 template <int dim, int spacedim>
2603 std::vector<unsigned int>
2605 const DoFHandler<dim, spacedim> &dof_handler,
2606 const unsigned int level,
2607 const bool interior_only,
2608 const bool boundary_patches,
2609 const bool level_boundary_patches,
2610 const bool single_cell_patches,
2611 const bool invert_vertex_mapping)
2612 {
2613 const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
2614 BlockMask exclude_boundary_dofs = BlockMask(n_blocks, interior_only);
2615 return make_vertex_patches(block_list,
2616 dof_handler,
2617 level,
2618 exclude_boundary_dofs,
2619 boundary_patches,
2620 level_boundary_patches,
2621 single_cell_patches,
2622 invert_vertex_mapping);
2623 }
2624
2625 template <int dim, int spacedim>
2626 std::vector<unsigned int>
2628 const DoFHandler<dim, spacedim> &dof_handler,
2629 const unsigned int level,
2630 const BlockMask & exclude_boundary_dofs,
2631 const bool boundary_patches,
2632 const bool level_boundary_patches,
2633 const bool single_cell_patches,
2634 const bool invert_vertex_mapping)
2635 {
2636 // Vector mapping from vertex index in the triangulation to consecutive
2637 // block indices on this level The number of cells at a vertex
2638 std::vector<unsigned int> vertex_cell_count(
2639 dof_handler.get_triangulation().n_vertices(), 0);
2640
2641 // Is a vertex at the boundary?
2642 std::vector<bool> vertex_boundary(
2643 dof_handler.get_triangulation().n_vertices(), false);
2644
2645 std::vector<unsigned int> vertex_mapping(
2646 dof_handler.get_triangulation().n_vertices(),
2648
2649 // Estimate for the number of dofs at this point
2650 std::vector<unsigned int> vertex_dof_count(
2651 dof_handler.get_triangulation().n_vertices(), 0);
2652
2653 // Identify all vertices active on this level and remember some data
2654 // about them
2655 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2656 for (const unsigned int v : cell->vertex_indices())
2657 {
2658 const unsigned int vg = cell->vertex_index(v);
2659 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2660 ++vertex_cell_count[vg];
2661 for (unsigned int d = 0; d < dim; ++d)
2662 {
2663 const unsigned int face = GeometryInfo<dim>::vertex_to_face[v][d];
2664 if (cell->at_boundary(face))
2665 vertex_boundary[vg] = true;
2666 else if ((!level_boundary_patches) &&
2667 (cell->neighbor(face)->level() !=
2668 static_cast<int>(level)))
2669 vertex_boundary[vg] = true;
2670 }
2671 }
2672 // From now on, only vertices with positive dof count are "in".
2673
2674 // Remove vertices at boundaries or in corners
2675 for (unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2676 if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2677 (!boundary_patches && vertex_boundary[vg]))
2678 vertex_dof_count[vg] = 0;
2679
2680 // Create a mapping from all vertices to the ones used here
2681 unsigned int n_vertex_count = 0;
2682 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2683 if (vertex_dof_count[vg] != 0)
2684 vertex_mapping[vg] = n_vertex_count++;
2685
2686 // Compactify dof count
2687 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2688 if (vertex_dof_count[vg] != 0)
2689 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2690
2691 // Now that we have all the data, we reduce it to the part we actually
2692 // want
2693 vertex_dof_count.resize(n_vertex_count);
2694
2695 // At this point, the list of patches is ready. Now we enter the dofs
2696 // into the sparsity pattern.
2697 block_list.reinit(vertex_dof_count.size(),
2698 dof_handler.n_dofs(level),
2699 vertex_dof_count);
2700
2701 std::vector<types::global_dof_index> indices;
2702 std::vector<bool> exclude;
2703
2704 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2705 {
2706 const FiniteElement<dim> &fe = cell->get_fe();
2707 indices.resize(fe.n_dofs_per_cell());
2708 cell->get_mg_dof_indices(indices);
2709
2710 for (const unsigned int v : cell->vertex_indices())
2711 {
2712 const unsigned int vg = cell->vertex_index(v);
2713 const unsigned int block = vertex_mapping[vg];
2714 if (block == numbers::invalid_unsigned_int)
2715 continue;
2716
2717 // Collect excluded dofs for some block(s) if boundary dofs
2718 // for a block are decided to be excluded
2719 if (exclude_boundary_dofs.size() == 0 ||
2720 exclude_boundary_dofs.n_selected_blocks() != 0)
2721 {
2722 // Exclude degrees of freedom on faces opposite to the
2723 // vertex
2724 exclude.resize(fe.n_dofs_per_cell());
2725 std::fill(exclude.begin(), exclude.end(), false);
2726
2727 for (unsigned int d = 0; d < dim; ++d)
2728 {
2729 const unsigned int a_face =
2731 const unsigned int face =
2733 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2734 {
2735 // For each dof, get the block it is in and decide to
2736 // exclude it or not
2737 if (exclude_boundary_dofs[fe.system_to_block_index(
2739 i, face))
2740 .first] == true)
2741 exclude[fe.face_to_cell_index(i, face)] = true;
2742 }
2743 }
2744 for (unsigned int j = 0; j < indices.size(); ++j)
2745 if (!exclude[j])
2746 block_list.add(block, indices[j]);
2747 }
2748 else
2749 {
2750 for (const auto index : indices)
2751 block_list.add(block, index);
2752 }
2753 }
2754 }
2755
2756 if (invert_vertex_mapping)
2757 {
2758 // Compress vertex mapping
2759 unsigned int n_vertex_count = 0;
2760 for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2761 if (vertex_mapping[vg] != numbers::invalid_unsigned_int)
2762 vertex_mapping[n_vertex_count++] = vg;
2763
2764 // Now we reduce it to the part we actually want
2765 vertex_mapping.resize(n_vertex_count);
2766 }
2767
2768 return vertex_mapping;
2769 }
2770
2771
2772 template <int dim, int spacedim>
2773 unsigned int
2775 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2776 &patch)
2777 {
2778 std::set<types::global_dof_index> dofs_on_patch;
2779 std::vector<types::global_dof_index> local_dof_indices;
2780
2781 // loop over the cells in the patch and get the DoFs on each.
2782 // add all of them to a std::set which automatically makes sure
2783 // all duplicates are ignored
2784 for (unsigned int i = 0; i < patch.size(); ++i)
2785 {
2787 patch[i];
2788 Assert(cell->is_artificial() == false,
2789 ExcMessage("This function can not be called with cells that are "
2790 "not either locally owned or ghost cells."));
2791 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2792 cell->get_dof_indices(local_dof_indices);
2793 dofs_on_patch.insert(local_dof_indices.begin(),
2794 local_dof_indices.end());
2795 }
2796
2797 // now return the number of DoFs (duplicates were ignored)
2798 return dofs_on_patch.size();
2799 }
2800
2801
2802
2803 template <int dim, int spacedim>
2804 std::vector<types::global_dof_index>
2806 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2807 &patch)
2808 {
2809 std::set<types::global_dof_index> dofs_on_patch;
2810 std::vector<types::global_dof_index> local_dof_indices;
2811
2812 // loop over the cells in the patch and get the DoFs on each.
2813 // add all of them to a std::set which automatically makes sure
2814 // all duplicates are ignored
2815 for (unsigned int i = 0; i < patch.size(); ++i)
2816 {
2818 patch[i];
2819 Assert(cell->is_artificial() == false,
2820 ExcMessage("This function can not be called with cells that are "
2821 "not either locally owned or ghost cells."));
2822 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2823 cell->get_dof_indices(local_dof_indices);
2824 dofs_on_patch.insert(local_dof_indices.begin(),
2825 local_dof_indices.end());
2826 }
2827
2828 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2830
2831 // return a vector with the content of the set above. copying
2832 // also ensures that we retain sortedness as promised in the
2833 // documentation and as necessary to retain the block structure
2834 // also on the local system
2835 return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2836 dofs_on_patch.end());
2837 }
2838
2839
2840} // end of namespace DoFTools
2841
2842
2843
2844// explicit instantiations
2845
2846#include "dof_tools.inst"
2847
2848
2849
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type n_constraints() const
unsigned int size() const
unsigned int n_selected_blocks(const unsigned int overall_number_of_blocks=numbers::invalid_unsigned_int) const
bool represents_n_components(const unsigned int n) const
bool represents_the_all_selected_mask() const
unsigned int size() const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_dofs() const
active_cell_iterator begin_active(const unsigned int level=0) const
std::vector< types::fe_index > get_active_fe_indices() const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
types::global_dof_index n_locally_owned_dofs() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
const FEValues< dim, spacedim > & get_present_fe_values() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_blocks() const
unsigned int n_components() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
bool is_primitive() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
unsigned int component_to_block_index(const unsigned int component) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
unsigned int element_multiplicity(const unsigned int index) const
unsigned int n_base_elements() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
size_type n_elements() const
Definition index_set.h:1816
bool is_element(const size_type index) const
Definition index_set.h:1776
void add_index(const size_type index)
Definition index_set.h:1680
void fill_binary_vector(VectorType &vector) const
void subtract_set(const IndexSet &other)
Definition index_set.cc:268
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1864
void compress() const
Definition index_set.h:1669
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1727
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
void add(const size_type i, const size_type j)
virtual MPI_Comm get_communicator() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
unsigned int n_levels() const
unsigned int n_vertices() const
size_type size() const
unsigned int size() const
Definition collection.h:265
unsigned int max_dofs_per_face() const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
unsigned int n_components() const
unsigned int max_dofs_per_cell() const
void push_back(const Quadrature< dim_in > &new_quadrature)
virtual MPI_Comm get_communicator() const override
Definition tria_base.cc:160
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
unsigned int cell_index
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< FilteredIterator< BaseIterator > > filter_iterators(IteratorRange< BaseIterator > i, const Predicate &p)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::line_iterator line_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
typename ActiveSelector::CellAccessor cell_accessor
@ update_quadrature_points
Transformed quadrature points.
void get_block_association(const DoFHandler< dim, spacedim > &dof, std::vector< unsigned char > &dofs_by_block)
Definition dof_tools.cc:247
std::vector< unsigned char > get_local_component_association(const FiniteElement< dim, spacedim > &fe, const ComponentMask &component_mask)
Definition dof_tools.cc:123
void resolve_components(const FiniteElement< dim, spacedim > &fe, const std::vector< unsigned char > &dofs_by_component, const std::vector< unsigned int > &target_component, const bool only_once, std::vector< types::global_dof_index > &dofs_per_component, unsigned int &component)
void get_component_association(const DoFHandler< dim, spacedim > &dof, const ComponentMask &component_mask, std::vector< unsigned char > &dofs_by_component)
Definition dof_tools.cc:194
void get_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::subdomain_id > &subdomain)
void extract_boundary_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs, const std::set< types::boundary_id > &boundary_ids={})
Definition dof_tools.cc:545
IndexSet extract_dofs_with_support_contained_within(const DoFHandler< dim, spacedim > &dof_handler, const std::function< bool(const typename DoFHandler< dim, spacedim >::active_cell_iterator &)> &predicate, const AffineConstraints< number > &constraints=AffineConstraints< number >())
Definition dof_tools.cc:796
IndexSet dof_indices_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
void map_dofs_to_support_points(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, std::vector< Point< spacedim > > &support_points, const ComponentMask &mask=ComponentMask())
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
void make_cell_patches(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const std::vector< bool > &selected_dofs={}, const types::global_dof_index offset=0)
IndexSet extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask)
Definition dof_tools.cc:386
IndexSet extract_locally_active_dofs(const DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
std::vector< IndexSet > locally_owned_dofs_per_component(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &components=ComponentMask())
Definition dof_tools.cc:439
void write_gnuplot_dof_support_point_info(std::ostream &out, const std::map< types::global_dof_index, Point< spacedim > > &support_points)
void extract_subdomain_dofs(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain_id, std::vector< bool > &selected_dofs)
IndexSet extract_locally_active_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
std::vector< unsigned int > make_vertex_patches(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only, const bool boundary_patches=false, const bool level_boundary_patches=false, const bool single_cell_patches=false, const bool invert_vertex_mapping=false)
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
void extract_level_dofs(const unsigned int level, const DoFHandler< dim, spacedim > &dof, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition dof_tools.cc:473
void extract_dofs_with_support_on_boundary(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
Definition dof_tools.cc:714
void distribute_cell_to_dof_vector(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &cell_data, Vector< double > &dof_data, const unsigned int component=0)
Definition dof_tools.cc:302
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
void make_child_patches(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only, const bool boundary_dofs=false)
void map_dof_to_boundary_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &mapping)
std::vector< IndexSet > locally_relevant_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
void make_single_patch(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only=false)
std::vector< types::global_dof_index > get_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
unsigned int count_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
void convert_couplings_to_blocks(const DoFHandler< dim, spacedim > &dof_handler, const Table< 2, Coupling > &table_by_component, std::vector< Table< 2, Coupling > > &tables_by_block)
void extract_constant_modes(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool > > &constant_modes)
IndexSet extract_hanging_node_dofs(const DoFHandler< dim, spacedim > &dof_handler)
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:150
const types::boundary_id internal_face_boundary_id
Definition types.h:277
const types::subdomain_id artificial_subdomain_id
Definition types.h:315
const types::subdomain_id invalid_subdomain_id
Definition types.h:298
static const unsigned int invalid_unsigned_int
Definition types.h:213
const types::global_dof_index invalid_dof_index
Definition types.h:233
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:60
bool operator()(const Point< dim, Number > &lhs, const Point< dim, Number > &rhs) const
Definition dof_tools.cc:81
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
const ::Triangulation< dim, spacedim > & tria
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition types.h:99