Reference documentation for deal.II version GIT 9d46a3bd8c 2023-10-03 12:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 
34 #include <deal.II/grid/tria.h>
36 
38 #include <deal.II/hp/fe_values.h>
41 
45 #include <deal.II/lac/vector.h>
46 
47 #include <algorithm>
48 #include <numeric>
49 
51 
52 
53 
54 namespace 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(),
170  ExcInternalError());
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(),
278  ExcInternalError());
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  {
307  const Triangulation<dim, spacedim> &tria = dof_handler.get_triangulation();
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),
376  ExcInternalError());
377  if (touch_count[i] != 0)
378  dof_data(i) /= touch_count[i];
379  }
380  }
381 
382 
383 
384  template <int dim, int spacedim>
385  IndexSet
387  const ComponentMask &component_mask)
388  {
389  Assert(component_mask.represents_n_components(
391  ExcMessage(
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(
401  dof.get_fe_collection().n_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>
426  IndexSet
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),
444  ExcMessage(
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(
482  ExcMessage(
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(
497  dof.get_fe_collection().n_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 =
510  internal::get_local_component_association(fe, component_mask);
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),
553  ExcMessage(
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, 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>
585  IndexSet
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(),
699  ExcInternalError());
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(),
782  ExcInternalError());
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>
795  IndexSet
797  const DoFHandler<dim, spacedim> &dof_handler,
798  const std::function<
800  &predicate,
801  const AffineConstraints<number> &cm)
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>
903  IndexSet
905  {
906  // there are no hanging nodes in 1d
907  return IndexSet(dof_handler.n_dofs());
908  }
909 
910 
911  template <int spacedim>
912  IndexSet
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>
956  IndexSet
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
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 and put into an index set
1163  std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1164  dof_set.add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1165  dof_set.compress();
1166 
1167  return dof_set;
1168  }
1169 
1170 
1171 
1172  template <int dim, int spacedim>
1173  void
1175  IndexSet &dof_set)
1176  {
1177  dof_set = extract_locally_relevant_dofs(dof_handler);
1178  }
1179 
1180 
1181 
1182  template <int dim, int spacedim>
1183  IndexSet
1185  const DoFHandler<dim, spacedim> &dof_handler,
1186  const unsigned int level)
1187  {
1188  // collect all the locally owned dofs
1189  IndexSet dof_set = dof_handler.locally_owned_mg_dofs(level);
1190 
1191  // add the DoF on the adjacent ghost cells to the IndexSet
1192 
1193  // Note: For certain meshes (in particular in 3d and with many
1194  // processors), it is really necessary to cache intermediate data. After
1195  // trying several objects such as std::set, a vector that is always kept
1196  // sorted, and a vector that is initially unsorted and sorted once at the
1197  // end, the latter has been identified to provide the best performance.
1198  // Martin Kronbichler
1199  std::vector<types::global_dof_index> dof_indices;
1200  std::vector<types::global_dof_index> dofs_on_ghosts;
1201 
1202  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
1203  {
1204  const types::subdomain_id id = cell->level_subdomain_id();
1205 
1206  // skip artificial and own cells (only look at ghost cells)
1207  if (id == dof_handler.get_triangulation().locally_owned_subdomain() ||
1209  continue;
1210 
1211  dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1212  cell->get_mg_dof_indices(dof_indices);
1213  for (const auto dof_index : dof_indices)
1214  if (!dof_set.is_element(dof_index))
1215  dofs_on_ghosts.push_back(dof_index);
1216  }
1217 
1218  // sort and fill into an index set
1219  std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1220  dof_set.add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1221  dof_set.compress();
1222 
1223  return dof_set;
1224  }
1225 
1226 
1227 
1228  template <int dim, int spacedim>
1229  void
1231  const DoFHandler<dim, spacedim> &dof_handler,
1232  const unsigned int level,
1233  IndexSet &dof_set)
1234  {
1235  dof_set = extract_locally_relevant_level_dofs(dof_handler, level);
1236  }
1237 
1238 
1239 
1240  template <int dim, int spacedim>
1241  void
1243  const ComponentMask &component_mask,
1244  std::vector<std::vector<bool>> &constant_modes)
1245  {
1246  // If there are no locally owned DoFs, return with an empty
1247  // constant_modes object:
1248  if (dof_handler.n_locally_owned_dofs() == 0)
1249  {
1250  constant_modes = std::vector<std::vector<bool>>(0);
1251  return;
1252  }
1253 
1254  const unsigned int n_components = dof_handler.get_fe(0).n_components();
1255  Assert(component_mask.represents_n_components(n_components),
1256  ExcDimensionMismatch(n_components, component_mask.size()));
1257 
1258  std::vector<unsigned char> dofs_by_component(
1259  dof_handler.n_locally_owned_dofs());
1261  component_mask,
1262  dofs_by_component);
1263  unsigned int n_selected_dofs = 0;
1264  for (unsigned int i = 0; i < n_components; ++i)
1265  if (component_mask[i] == true)
1266  n_selected_dofs +=
1267  std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1268 
1269  // Find local numbering within the selected components
1270  const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
1271  std::vector<unsigned int> component_numbering(
1272  locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int);
1273  for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements();
1274  ++i)
1275  if (component_mask[dofs_by_component[i]])
1276  component_numbering[i] = count++;
1277 
1278  // get the element constant modes and find a translation table between
1279  // index in the constant modes and the components.
1280  //
1281  // TODO: We might be able to extend this also for elements which do not
1282  // have the same constant modes, but that is messy...
1283  const ::hp::FECollection<dim, spacedim> &fe_collection =
1284  dof_handler.get_fe_collection();
1285  std::vector<Table<2, bool>> element_constant_modes;
1286  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1287  constant_mode_to_component_translation(n_components);
1288  {
1289  unsigned int n_constant_modes = 0;
1290  int first_non_empty_constant_mode = -1;
1291  for (unsigned int f = 0; f < fe_collection.size(); ++f)
1292  {
1293  std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1294  fe_collection[f].get_constant_modes();
1295 
1296  // Store the index of the current element if it is the first that has
1297  // non-empty constant modes.
1298  if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
1299  {
1300  first_non_empty_constant_mode = f;
1301  // This is the first non-empty constant mode, so we figure out the
1302  // translation between index in the constant modes and the
1303  // components
1304  for (unsigned int i = 0; i < data.second.size(); ++i)
1305  if (component_mask[data.second[i]])
1306  constant_mode_to_component_translation[data.second[i]]
1307  .emplace_back(n_constant_modes++, i);
1308  }
1309 
1310  // Add the constant modes of this element to the list and assert that
1311  // there are as many constant modes as for the other elements (or zero
1312  // constant modes).
1313  element_constant_modes.push_back(data.first);
1314  Assert(
1315  element_constant_modes.back().n_rows() == 0 ||
1316  element_constant_modes.back().n_rows() ==
1317  element_constant_modes[first_non_empty_constant_mode].n_rows(),
1318  ExcInternalError());
1319  }
1320  AssertIndexRange(first_non_empty_constant_mode, fe_collection.size());
1321 
1322  // Now we know the number of constant modes and resize the return vector
1323  // accordingly
1324  constant_modes.clear();
1325  constant_modes.resize(n_constant_modes,
1326  std::vector<bool>(n_selected_dofs, false));
1327  }
1328 
1329  // Loop over all owned cells and ask the element for the constant modes
1330  std::vector<types::global_dof_index> dof_indices;
1331  for (const auto &cell : dof_handler.active_cell_iterators() |
1333  {
1334  dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1335  cell->get_dof_indices(dof_indices);
1336 
1337  for (unsigned int i = 0; i < dof_indices.size(); ++i)
1338  if (locally_owned_dofs.is_element(dof_indices[i]))
1339  {
1340  const unsigned int loc_index =
1341  locally_owned_dofs.index_within_set(dof_indices[i]);
1342  const unsigned int comp = dofs_by_component[loc_index];
1343  if (component_mask[comp])
1344  for (auto &indices :
1345  constant_mode_to_component_translation[comp])
1346  constant_modes[indices
1347  .first][component_numbering[loc_index]] =
1348  element_constant_modes[cell->active_fe_index()](
1349  indices.second, i);
1350  }
1351  }
1352  }
1353 
1354 
1355 
1356  template <int dim, int spacedim>
1357  std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1358  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1359  unsigned int>>
1361  const std::map<typename Triangulation<dim - 1, spacedim>::cell_iterator,
1363  &c1_to_c0_face,
1364  const DoFHandler<dim, spacedim> &c0_dh,
1365  const DoFHandler<dim - 1, spacedim> &c1_dh)
1366  {
1367  // This is the returned object: a map of codimension-1 active dof cell
1368  // iterators to codimension-0 cells and face indices
1369  std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1370  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1371  unsigned int>>
1372  c1_to_c0_cell_and_face;
1373 
1374  // Shortcut if there are no faces to check
1375  if (c1_to_c0_face.empty())
1376  return c1_to_c0_cell_and_face;
1377 
1378  // This is the partial inverse of the map passed as input, for dh
1379  std::map<typename Triangulation<dim, spacedim>::face_iterator,
1380  typename DoFHandler<dim - 1, spacedim>::active_cell_iterator>
1381  c0_to_c1;
1382 
1383  // map volume mesh face -> codimension 1 dof cell
1384  for (const auto &[c1_cell, c0_cell] : c1_to_c0_face)
1385  if (!c1_cell->has_children())
1386  c0_to_c1[c0_cell] = c1_cell->as_dof_handler_iterator(c1_dh);
1387 
1388  // generate a mapping that maps codimension-1 cells
1389  // to codimension-0 cells and faces
1390  for (const auto &cell :
1391  c0_dh.active_cell_iterators()) // disp_dof.active_cell_iterators())
1392  for (const auto f : cell->face_indices())
1393  if (cell->face(f)->at_boundary())
1394  {
1395  const auto &it = c0_to_c1.find(cell->face(f));
1396  if (it != c0_to_c1.end())
1397  {
1398  const auto &c1_cell = it->second;
1399  c1_to_c0_cell_and_face[c1_cell] = {cell, f};
1400  c0_to_c1.erase(it);
1401  }
1402  }
1403  // Check the dimensions: make sure all active cells we had have been mapped.
1404  AssertDimension(c0_to_c1.size(), 0);
1405  return c1_to_c0_cell_and_face;
1406  }
1407 
1408 
1409 
1410  template <int dim, int spacedim>
1411  void
1413  std::vector<unsigned int> &active_fe_indices)
1414  {
1415  AssertDimension(active_fe_indices.size(),
1416  dof_handler.get_triangulation().n_active_cells());
1417 
1418  std::vector<types::fe_index> indices = dof_handler.get_active_fe_indices();
1419 
1420  active_fe_indices.assign(indices.begin(), indices.end());
1421  }
1422 
1423  template <int dim, int spacedim>
1424  std::vector<IndexSet>
1426  {
1427  Assert(dof_handler.n_dofs() > 0,
1428  ExcMessage("The given DoFHandler has no DoFs."));
1429 
1430  // If the Triangulation is distributed, the only thing we can usefully
1431  // ask is for its locally owned subdomain
1432  Assert((dynamic_cast<
1434  &dof_handler.get_triangulation()) == nullptr),
1435  ExcMessage(
1436  "For parallel::distributed::Triangulation objects and "
1437  "associated DoF handler objects, asking for any information "
1438  "related to a subdomain other than the locally owned one does "
1439  "not make sense."));
1440 
1441  // The following is a random process (flip of a coin), thus should be called
1442  // once only.
1443  std::vector<::types::subdomain_id> subdomain_association(
1444  dof_handler.n_dofs());
1446  subdomain_association);
1447 
1448  // Figure out how many subdomain ids there are.
1449  //
1450  // if this is a parallel triangulation, then we can just ask the
1451  // triangulation for this. if this is a sequential triangulation, we loop
1452  // over all cells and take the largest subdomain_id value we find; the
1453  // number of subdomains is then the largest found value plus one. (we here
1454  // assume that all subdomain ids up to the largest are actually used; this
1455  // may not be true for a sequential triangulation where these values have
1456  // been set by hand and not in accordance with some MPI communicator; but
1457  // the function returns an array indexed starting at zero, so we need to
1458  // collect information for each subdomain index anyway, not just for the
1459  // used one.)
1460  const unsigned int n_subdomains =
1461  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1462  &dof_handler.get_triangulation()) == nullptr ?
1463  [&dof_handler]() {
1464  unsigned int max_subdomain_id = 0;
1465  for (const auto &cell : dof_handler.active_cell_iterators())
1466  max_subdomain_id =
1467  std::max(max_subdomain_id, cell->subdomain_id());
1468  return max_subdomain_id + 1;
1469  }() :
1471  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1472  &dof_handler.get_triangulation())
1473  ->get_communicator()));
1474  Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1475  subdomain_association.end()),
1476  ExcInternalError());
1477 
1478  std::vector<::IndexSet> index_sets(
1479  n_subdomains, ::IndexSet(dof_handler.n_dofs()));
1480 
1481  // loop over subdomain_association and populate IndexSet when a
1482  // change in subdomain ID is found
1483  ::types::global_dof_index i_min = 0;
1484  ::types::global_dof_index this_subdomain = subdomain_association[0];
1485 
1486  for (::types::global_dof_index index = 1;
1487  index < subdomain_association.size();
1488  ++index)
1489  {
1490  // found index different from the current one
1491  if (subdomain_association[index] != this_subdomain)
1492  {
1493  index_sets[this_subdomain].add_range(i_min, index);
1494  i_min = index;
1495  this_subdomain = subdomain_association[index];
1496  }
1497  }
1498 
1499  // the very last element is of different index
1500  if (i_min == subdomain_association.size() - 1)
1501  {
1502  index_sets[this_subdomain].add_index(i_min);
1503  }
1504 
1505  // otherwise there are at least two different indices
1506  else
1507  {
1508  index_sets[this_subdomain].add_range(i_min,
1509  subdomain_association.size());
1510  }
1511 
1512  for (unsigned int i = 0; i < n_subdomains; ++i)
1513  index_sets[i].compress();
1514 
1515  return index_sets;
1516  }
1517 
1518  template <int dim, int spacedim>
1519  std::vector<IndexSet>
1521  const DoFHandler<dim, spacedim> &dof_handler)
1522  {
1523  // If the Triangulation is distributed, the only thing we can usefully
1524  // ask is for its locally owned subdomain
1525  Assert((dynamic_cast<
1527  &dof_handler.get_triangulation()) == nullptr),
1528  ExcMessage(
1529  "For parallel::distributed::Triangulation objects and "
1530  "associated DoF handler objects, asking for any information "
1531  "related to a subdomain other than the locally owned one does "
1532  "not make sense."));
1533 
1534  // Collect all the locally owned DoFs
1535  // Note: Even though the distribution of DoFs by the
1536  // locally_owned_dofs_per_subdomain function is pseudo-random, we will
1537  // collect all the DoFs on the subdomain and its layer cell. Therefore, the
1538  // random nature of this function does not play a role in the extraction of
1539  // the locally relevant DoFs
1540  std::vector<IndexSet> dof_set =
1541  locally_owned_dofs_per_subdomain(dof_handler);
1542  const ::types::subdomain_id n_subdomains = dof_set.size();
1543 
1544  // Add the DoFs on the adjacent (equivalent ghost) cells to the IndexSet,
1545  // cache them in a set. Need to check each DoF manually because we can't
1546  // be sure that the DoF range of locally_owned_dofs is really contiguous.
1547  for (::types::subdomain_id subdomain_id = 0;
1548  subdomain_id < n_subdomains;
1549  ++subdomain_id)
1550  {
1551  // Extract the layer of cells around this subdomain
1552  std::function<bool(
1555  const std::vector<
1557  active_halo_layer =
1558  GridTools::compute_active_cell_halo_layer(dof_handler, predicate);
1559 
1560  // Extract DoFs associated with halo layer
1561  std::vector<types::global_dof_index> local_dof_indices;
1562  std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1563  for (typename std::vector<
1565  const_iterator it_cell = active_halo_layer.begin();
1566  it_cell != active_halo_layer.end();
1567  ++it_cell)
1568  {
1570  &cell = *it_cell;
1571  Assert(
1572  cell->subdomain_id() != subdomain_id,
1573  ExcMessage(
1574  "The subdomain ID of the halo cell should not match that of the vector entry."));
1575 
1576  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1577  cell->get_dof_indices(local_dof_indices);
1578 
1579  for (const types::global_dof_index local_dof_index :
1580  local_dof_indices)
1581  subdomain_halo_global_dof_indices.insert(local_dof_index);
1582  }
1583 
1584  dof_set[subdomain_id].add_indices(
1585  subdomain_halo_global_dof_indices.begin(),
1586  subdomain_halo_global_dof_indices.end());
1587 
1588  dof_set[subdomain_id].compress();
1589  }
1590 
1591  return dof_set;
1592  }
1593 
1594  template <int dim, int spacedim>
1595  void
1597  const DoFHandler<dim, spacedim> &dof_handler,
1598  std::vector<types::subdomain_id> &subdomain_association)
1599  {
1600  // if the Triangulation is distributed, the only thing we can usefully
1601  // ask is for its locally owned subdomain
1602  Assert((dynamic_cast<
1604  &dof_handler.get_triangulation()) == nullptr),
1605  ExcMessage(
1606  "For parallel::distributed::Triangulation objects and "
1607  "associated DoF handler objects, asking for any subdomain other "
1608  "than the locally owned one does not make sense."));
1609 
1610  Assert(subdomain_association.size() == dof_handler.n_dofs(),
1611  ExcDimensionMismatch(subdomain_association.size(),
1612  dof_handler.n_dofs()));
1613 
1614  // catch an error that happened in some versions of the shared tria
1615  // distribute_dofs() function where we were trying to call this
1616  // function at a point in time when not all internal DoFHandler
1617  // structures were quite set up yet.
1618  Assert(dof_handler.n_dofs() > 0, ExcInternalError());
1619 
1620  // In case this function is executed with parallel::shared::Triangulation
1621  // with possibly artificial cells, we need to take "true" subdomain IDs
1622  // (i.e. without artificial cells). Otherwise we are good to use
1623  // subdomain_id as stored in cell->subdomain_id().
1624  std::vector<types::subdomain_id> cell_owners(
1625  dof_handler.get_triangulation().n_active_cells());
1627  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
1628  &dof_handler.get_triangulation())))
1629  {
1630  cell_owners = tr->get_true_subdomain_ids_of_cells();
1631  Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1632  tr->n_active_cells(),
1633  ExcInternalError());
1634  }
1635  else
1636  {
1637  for (const auto &cell : dof_handler.active_cell_iterators() |
1639  cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1640  }
1641 
1642  // preset all values by an invalid value
1643  std::fill_n(subdomain_association.begin(),
1644  dof_handler.n_dofs(),
1646 
1647  std::vector<types::global_dof_index> local_dof_indices;
1648  local_dof_indices.reserve(
1649  dof_handler.get_fe_collection().max_dofs_per_cell());
1650 
1651  // loop over all cells and record which subdomain a DoF belongs to.
1652  // give to the smaller subdomain_id in case it is on an interface
1653  for (const auto &cell : dof_handler.active_cell_iterators())
1654  {
1656  cell_owners[cell->active_cell_index()];
1657  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1658  local_dof_indices.resize(dofs_per_cell);
1659  cell->get_dof_indices(local_dof_indices);
1660 
1661  // set subdomain ids. if dofs already have their values set then
1662  // they must be on partition interfaces. in that case assign them
1663  // to either the previous association or the current processor
1664  // with the smaller subdomain id.
1665  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1666  if (subdomain_association[local_dof_indices[i]] ==
1668  subdomain_association[local_dof_indices[i]] = subdomain_id;
1669  else if (subdomain_association[local_dof_indices[i]] > subdomain_id)
1670  {
1671  subdomain_association[local_dof_indices[i]] = subdomain_id;
1672  }
1673  }
1674 
1675  Assert(std::find(subdomain_association.begin(),
1676  subdomain_association.end(),
1678  subdomain_association.end(),
1679  ExcInternalError());
1680  }
1681 
1682 
1683 
1684  template <int dim, int spacedim>
1685  unsigned int
1687  const DoFHandler<dim, spacedim> &dof_handler,
1688  const types::subdomain_id subdomain)
1689  {
1690  std::vector<types::subdomain_id> subdomain_association(
1691  dof_handler.n_dofs());
1692  get_subdomain_association(dof_handler, subdomain_association);
1693 
1694  return std::count(subdomain_association.begin(),
1695  subdomain_association.end(),
1696  subdomain);
1697  }
1698 
1699 
1700 
1701  template <int dim, int spacedim>
1702  IndexSet
1704  const DoFHandler<dim, spacedim> &dof_handler,
1705  const types::subdomain_id subdomain)
1706  {
1707  // If we have a distributed::Triangulation only allow locally_owned
1708  // subdomain.
1711  (subdomain ==
1712  dof_handler.get_triangulation().locally_owned_subdomain()),
1713  ExcMessage(
1714  "For parallel::distributed::Triangulation objects and "
1715  "associated DoF handler objects, asking for any subdomain other "
1716  "than the locally owned one does not make sense."));
1717 
1718  IndexSet index_set(dof_handler.n_dofs());
1719 
1720  std::vector<types::global_dof_index> local_dof_indices;
1721  local_dof_indices.reserve(
1722  dof_handler.get_fe_collection().max_dofs_per_cell());
1723 
1724  // first generate an unsorted list of all indices which we fill from
1725  // the back. could also insert them directly into the IndexSet, but
1726  // that inserts indices in the middle, which is an O(n^2) algorithm and
1727  // hence too expensive. Could also use std::set, but that is in general
1728  // more expensive than a vector
1729  std::vector<types::global_dof_index> subdomain_indices;
1730 
1731  for (const auto &cell : dof_handler.active_cell_iterators())
1732  if ((cell->is_artificial() == false) &&
1733  (cell->subdomain_id() == subdomain))
1734  {
1735  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1736  local_dof_indices.resize(dofs_per_cell);
1737  cell->get_dof_indices(local_dof_indices);
1738  subdomain_indices.insert(subdomain_indices.end(),
1739  local_dof_indices.begin(),
1740  local_dof_indices.end());
1741  }
1742  // sort indices and put into an index set:
1743  std::sort(subdomain_indices.begin(), subdomain_indices.end());
1744  index_set.add_indices(subdomain_indices.begin(), subdomain_indices.end());
1745  index_set.compress();
1746 
1747  return index_set;
1748  }
1749 
1750 
1751 
1752  template <int dim, int spacedim>
1753  void
1755  const DoFHandler<dim, spacedim> &dof_handler,
1756  const types::subdomain_id subdomain,
1757  std::vector<unsigned int> &n_dofs_on_subdomain)
1758  {
1759  Assert(n_dofs_on_subdomain.size() == dof_handler.get_fe(0).n_components(),
1760  ExcDimensionMismatch(n_dofs_on_subdomain.size(),
1761  dof_handler.get_fe(0).n_components()));
1762  std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1763 
1764  // Make sure there are at least some cells with this subdomain id
1765  Assert(std::any_of(
1766  dof_handler.begin_active(),
1768  dof_handler.end()},
1769  [subdomain](
1770  const typename DoFHandler<dim, spacedim>::cell_accessor &cell) {
1771  return cell.subdomain_id() == subdomain;
1772  }),
1773  ExcMessage("There are no cells for the given subdomain!"));
1774 
1775  std::vector<types::subdomain_id> subdomain_association(
1776  dof_handler.n_dofs());
1777  get_subdomain_association(dof_handler, subdomain_association);
1778 
1779  std::vector<unsigned char> component_association(dof_handler.n_dofs());
1781  ComponentMask(std::vector<bool>()),
1782  component_association);
1783 
1784  for (unsigned int c = 0; c < dof_handler.get_fe(0).n_components(); ++c)
1785  {
1786  for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
1787  if ((subdomain_association[i] == subdomain) &&
1788  (component_association[i] == static_cast<unsigned char>(c)))
1789  ++n_dofs_on_subdomain[c];
1790  }
1791  }
1792 
1793 
1794 
1795  namespace internal
1796  {
1797  // TODO: why is this function so complicated? It would be nice to have
1798  // comments that explain why we can't just loop over all components and
1799  // count the entries in dofs_by_component that have this component's
1800  // index
1801  template <int dim, int spacedim>
1802  void
1804  const std::vector<unsigned char> &dofs_by_component,
1805  const std::vector<unsigned int> &target_component,
1806  const bool only_once,
1807  std::vector<types::global_dof_index> &dofs_per_component,
1808  unsigned int &component)
1809  {
1810  for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1811  {
1812  const FiniteElement<dim, spacedim> &base = fe.base_element(b);
1813  // Dimension of base element
1814  unsigned int d = base.n_components();
1815 
1816  for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1817  {
1818  if (base.n_base_elements() > 1)
1819  resolve_components(base,
1820  dofs_by_component,
1821  target_component,
1822  only_once,
1823  dofs_per_component,
1824  component);
1825  else
1826  {
1827  for (unsigned int dd = 0; dd < d; ++dd, ++component)
1828  dofs_per_component[target_component[component]] +=
1829  std::count(dofs_by_component.begin(),
1830  dofs_by_component.end(),
1831  component);
1832 
1833  // if we have non-primitive FEs and want all components
1834  // to show the number of dofs, need to copy the result to
1835  // those components
1836  if (!base.is_primitive() && !only_once)
1837  for (unsigned int dd = 1; dd < d; ++dd)
1838  dofs_per_component[target_component[component - d + dd]] =
1839  dofs_per_component[target_component[component - d]];
1840  }
1841  }
1842  }
1843  }
1844 
1845 
1846  template <int dim, int spacedim>
1847  void
1849  const std::vector<unsigned char> &dofs_by_component,
1850  const std::vector<unsigned int> &target_component,
1851  const bool only_once,
1852  std::vector<types::global_dof_index> &dofs_per_component,
1853  unsigned int &component)
1854  {
1855  // assert that all elements in the collection have the same structure
1856  // (base elements and multiplicity, components per base element) and
1857  // then simply call the function above
1858  for (unsigned int fe = 1; fe < fe_collection.size(); ++fe)
1859  {
1860  Assert(fe_collection[fe].n_components() ==
1861  fe_collection[0].n_components(),
1862  ExcNotImplemented());
1863  Assert(fe_collection[fe].n_base_elements() ==
1864  fe_collection[0].n_base_elements(),
1865  ExcNotImplemented());
1866  for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b)
1867  {
1868  Assert(fe_collection[fe].base_element(b).n_components() ==
1869  fe_collection[0].base_element(b).n_components(),
1870  ExcNotImplemented());
1871  Assert(fe_collection[fe].base_element(b).n_base_elements() ==
1872  fe_collection[0].base_element(b).n_base_elements(),
1873  ExcNotImplemented());
1874  }
1875  }
1876 
1877  resolve_components(fe_collection[0],
1878  dofs_by_component,
1879  target_component,
1880  only_once,
1881  dofs_per_component,
1882  component);
1883  }
1884  } // namespace internal
1885 
1886 
1887 
1888  namespace internal
1889  {
1890  namespace
1891  {
1895  template <int dim, int spacedim>
1896  bool
1897  all_elements_are_primitive(const FiniteElement<dim, spacedim> &fe)
1898  {
1899  return fe.is_primitive();
1900  }
1901 
1902 
1907  template <int dim, int spacedim>
1908  bool
1909  all_elements_are_primitive(
1910  const ::hp::FECollection<dim, spacedim> &fe_collection)
1911  {
1912  for (unsigned int i = 0; i < fe_collection.size(); ++i)
1913  if (fe_collection[i].is_primitive() == false)
1914  return false;
1915 
1916  return true;
1917  }
1918  } // namespace
1919  } // namespace internal
1920 
1921 
1922 
1923  template <int dim, int spacedim>
1924  std::vector<types::global_dof_index>
1926  const DoFHandler<dim, spacedim> &dof_handler,
1927  const bool only_once,
1928  const std::vector<unsigned int> &target_component_)
1929  {
1930  const unsigned int n_components = dof_handler.get_fe(0).n_components();
1931 
1932  // If the empty vector was given as default argument, set up this
1933  // vector as identity.
1934  std::vector<unsigned int> target_component = target_component_;
1935  if (target_component.empty())
1936  {
1937  target_component.resize(n_components);
1938  for (unsigned int i = 0; i < n_components; ++i)
1939  target_component[i] = i;
1940  }
1941  else
1942  Assert(target_component.size() == n_components,
1943  ExcDimensionMismatch(target_component.size(), n_components));
1944 
1945 
1946  const unsigned int max_component =
1947  *std::max_element(target_component.begin(), target_component.end());
1948  const unsigned int n_target_components = max_component + 1;
1949 
1950  std::vector<types::global_dof_index> dofs_per_component(
1951  n_target_components, types::global_dof_index(0));
1952 
1953  // special case for only one component. treat this first since it does
1954  // not require any computations
1955  if (n_components == 1)
1956  {
1957  dofs_per_component[0] = dof_handler.n_locally_owned_dofs();
1958  return dofs_per_component;
1959  }
1960 
1961 
1962  // otherwise determine the number of dofs in each component separately.
1963  // do so in parallel
1964  std::vector<unsigned char> dofs_by_component(
1965  dof_handler.n_locally_owned_dofs());
1967  ComponentMask(),
1968  dofs_by_component);
1969 
1970  // next count what we got
1971  unsigned int component = 0;
1973  dofs_by_component,
1974  target_component,
1975  only_once,
1976  dofs_per_component,
1977  component);
1978  Assert(n_components == component, ExcInternalError());
1979 
1980  // finally sanity check. this is only valid if the finite element is
1981  // actually primitive, so exclude other elements from this
1982  Assert((internal::all_elements_are_primitive(
1983  dof_handler.get_fe_collection()) == false) ||
1984  (std::accumulate(dofs_per_component.begin(),
1985  dofs_per_component.end(),
1987  dof_handler.n_locally_owned_dofs()),
1988  ExcInternalError());
1989 
1990  // reduce information from all CPUs
1991 #ifdef DEAL_II_WITH_MPI
1992 
1994  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1995  &dof_handler.get_triangulation())))
1996  {
1997  std::vector<types::global_dof_index> local_dof_count =
1998  dofs_per_component;
1999 
2000  const int ierr = MPI_Allreduce(local_dof_count.data(),
2001  dofs_per_component.data(),
2002  n_target_components,
2004  MPI_SUM,
2005  tria->get_communicator());
2006  AssertThrowMPI(ierr);
2007  }
2008 #endif
2009 
2010  return dofs_per_component;
2011  }
2012 
2013 
2014 
2015  template <int dim, int spacedim>
2016  std::vector<types::global_dof_index>
2018  const std::vector<unsigned int> &target_block_)
2019  {
2020  const ::hp::FECollection<dim, spacedim> &fe_collection =
2021  dof_handler.get_fe_collection();
2022  Assert(fe_collection.size() < 256, ExcNotImplemented());
2023 
2024  // If the empty vector for target_block(e.g., as default argument), then
2025  // set up this vector as identity. We do this set up with the first
2026  // element of the collection, but the whole thing can only work if
2027  // all elements have the same number of blocks anyway -- so check
2028  // that right after
2029  const unsigned int n_blocks = fe_collection[0].n_blocks();
2030 
2031  std::vector<unsigned int> target_block = target_block_;
2032  if (target_block.empty())
2033  {
2034  target_block.resize(fe_collection[0].n_blocks());
2035  for (unsigned int i = 0; i < n_blocks; ++i)
2036  target_block[i] = i;
2037  }
2038  else
2039  Assert(target_block.size() == n_blocks,
2040  ExcDimensionMismatch(target_block.size(), n_blocks));
2041  for (unsigned int f = 1; f < fe_collection.size(); ++f)
2042  Assert(fe_collection[0].n_blocks() == fe_collection[f].n_blocks(),
2043  ExcMessage("This function can only work if all elements in a "
2044  "collection have the same number of blocks."));
2045 
2046  // special case for only one block. treat this first since it does
2047  // not require any computations
2048  if (n_blocks == 1)
2049  {
2050  std::vector<types::global_dof_index> dofs_per_block(1);
2051  dofs_per_block[0] = dof_handler.n_dofs();
2052  return dofs_per_block;
2053  }
2054 
2055  // Otherwise set up the right-sized object and start working
2056  const unsigned int max_block =
2057  *std::max_element(target_block.begin(), target_block.end());
2058  const unsigned int n_target_blocks = max_block + 1;
2059 
2060  std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2061 
2062  // Loop over the elements of the collection, but really only consider
2063  // the last element (see #9271)
2064  for (unsigned int this_fe = fe_collection.size() - 1;
2065  this_fe < fe_collection.size();
2066  ++this_fe)
2067  {
2068  const FiniteElement<dim, spacedim> &fe = fe_collection[this_fe];
2069 
2070  std::vector<unsigned char> dofs_by_block(
2071  dof_handler.n_locally_owned_dofs());
2072  internal::get_block_association(dof_handler, dofs_by_block);
2073 
2074  // next count what we got
2075  for (unsigned int block = 0; block < fe.n_blocks(); ++block)
2076  dofs_per_block[target_block[block]] +=
2077  std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2078 
2079 #ifdef DEAL_II_WITH_MPI
2080  // if we are working on a parallel mesh, we now need to collect
2081  // this information from all processors
2083  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2084  &dof_handler.get_triangulation())))
2085  {
2086  std::vector<types::global_dof_index> local_dof_count =
2087  dofs_per_block;
2088  const int ierr = MPI_Allreduce(local_dof_count.data(),
2089  dofs_per_block.data(),
2090  n_target_blocks,
2092  MPI_SUM,
2093  tria->get_communicator());
2094  AssertThrowMPI(ierr);
2095  }
2096 #endif
2097  }
2098 
2099  return dofs_per_block;
2100  }
2101 
2102 
2103 
2104  template <int dim, int spacedim>
2105  void
2107  std::vector<types::global_dof_index> &mapping)
2108  {
2109  mapping.clear();
2110  mapping.insert(mapping.end(),
2111  dof_handler.n_dofs(),
2113 
2114  std::vector<types::global_dof_index> dofs_on_face;
2115  dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2116  types::global_dof_index next_boundary_index = 0;
2117 
2118  // now loop over all cells and check whether their faces are at the
2119  // boundary. note that we need not take special care of single lines
2120  // being at the boundary (using @p{cell->has_boundary_lines}), since we
2121  // do not support boundaries of dimension dim-2, and so every isolated
2122  // boundary line is also part of a boundary face which we will be
2123  // visiting sooner or later
2124  for (const auto &cell : dof_handler.active_cell_iterators())
2125  for (const unsigned int f : cell->face_indices())
2126  if (cell->at_boundary(f))
2127  {
2128  const unsigned int dofs_per_face =
2129  cell->get_fe().n_dofs_per_face(f);
2130  dofs_on_face.resize(dofs_per_face);
2131  cell->face(f)->get_dof_indices(dofs_on_face,
2132  cell->active_fe_index());
2133  for (unsigned int i = 0; i < dofs_per_face; ++i)
2134  if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2135  mapping[dofs_on_face[i]] = next_boundary_index++;
2136  }
2137 
2138  AssertDimension(next_boundary_index, dof_handler.n_boundary_dofs());
2139  }
2140 
2141 
2142 
2143  template <int dim, int spacedim>
2144  void
2146  const std::set<types::boundary_id> &boundary_ids,
2147  std::vector<types::global_dof_index> &mapping)
2148  {
2149  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2150  boundary_ids.end(),
2152 
2153  mapping.clear();
2154  mapping.insert(mapping.end(),
2155  dof_handler.n_dofs(),
2157 
2158  // return if there is nothing to do
2159  if (boundary_ids.empty())
2160  return;
2161 
2162  std::vector<types::global_dof_index> dofs_on_face;
2163  dofs_on_face.reserve(dof_handler.get_fe_collection().max_dofs_per_face());
2164  types::global_dof_index next_boundary_index = 0;
2165 
2166  for (const auto &cell : dof_handler.active_cell_iterators())
2167  for (const unsigned int f : cell->face_indices())
2168  if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2169  boundary_ids.end())
2170  {
2171  const unsigned int dofs_per_face =
2172  cell->get_fe().n_dofs_per_face(f);
2173  dofs_on_face.resize(dofs_per_face);
2174  cell->face(f)->get_dof_indices(dofs_on_face,
2175  cell->active_fe_index());
2176  for (unsigned int i = 0; i < dofs_per_face; ++i)
2177  if (mapping[dofs_on_face[i]] == numbers::invalid_dof_index)
2178  mapping[dofs_on_face[i]] = next_boundary_index++;
2179  }
2180 
2181  AssertDimension(next_boundary_index,
2182  dof_handler.n_boundary_dofs(boundary_ids));
2183  }
2184 
2185  namespace internal
2186  {
2187  namespace
2188  {
2189  template <int dim, int spacedim>
2190  std::map<types::global_dof_index, Point<spacedim>>
2192  const hp::MappingCollection<dim, spacedim> &mapping,
2193  const DoFHandler<dim, spacedim> &dof_handler,
2194  const ComponentMask &in_mask)
2195  {
2196  std::map<types::global_dof_index, Point<spacedim>> support_points;
2197 
2198  const hp::FECollection<dim, spacedim> &fe_collection =
2199  dof_handler.get_fe_collection();
2200  hp::QCollection<dim> q_coll_dummy;
2201 
2202  for (unsigned int fe_index = 0; fe_index < fe_collection.size();
2203  ++fe_index)
2204  {
2205  // check whether every FE in the collection has support points
2206  Assert((fe_collection[fe_index].n_dofs_per_cell() == 0) ||
2207  (fe_collection[fe_index].has_support_points()),
2209  q_coll_dummy.push_back(Quadrature<dim>(
2210  fe_collection[fe_index].get_unit_support_points()));
2211  }
2212 
2213  // Take care of components
2214  const ComponentMask mask =
2215  (in_mask.size() == 0 ?
2216  ComponentMask(fe_collection.n_components(), true) :
2217  in_mask);
2218 
2219  // Now loop over all cells and enquire the support points on each
2220  // of these. we use dummy quadrature formulas where the quadrature
2221  // points are located at the unit support points to enquire the
2222  // location of the support points in real space.
2223  //
2224  // The weights of the quadrature rule have been set to invalid
2225  // values by the used constructor.
2226  hp::FEValues<dim, spacedim> hp_fe_values(mapping,
2227  fe_collection,
2228  q_coll_dummy,
2230 
2231  std::vector<types::global_dof_index> local_dof_indices;
2232  for (const auto &cell : dof_handler.active_cell_iterators())
2233  // only work on locally relevant cells
2234  if (cell->is_artificial() == false)
2235  {
2236  hp_fe_values.reinit(cell);
2237  const FEValues<dim, spacedim> &fe_values =
2238  hp_fe_values.get_present_fe_values();
2239 
2240  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2241  cell->get_dof_indices(local_dof_indices);
2242 
2243  const std::vector<Point<spacedim>> &points =
2244  fe_values.get_quadrature_points();
2245  for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2246  ++i)
2247  {
2248  const unsigned int dof_comp =
2249  cell->get_fe().system_to_component_index(i).first;
2250 
2251  // insert the values into the map if it is a valid component
2252  if (mask[dof_comp])
2253  support_points[local_dof_indices[i]] = points[i];
2254  }
2255  }
2256 
2257  return support_points;
2258  }
2259 
2260 
2261  template <int dim, int spacedim>
2262  std::vector<Point<spacedim>>
2263  map_dofs_to_support_points_vector(
2264  const hp::MappingCollection<dim, spacedim> &mapping,
2265  const DoFHandler<dim, spacedim> &dof_handler,
2266  const ComponentMask &mask)
2267  {
2268  std::vector<Point<spacedim>> support_points(dof_handler.n_dofs());
2269 
2270  // get the data in the form of the map as above
2271  const std::map<types::global_dof_index, Point<spacedim>>
2272  x_support_points =
2273  map_dofs_to_support_points(mapping, dof_handler, mask);
2274 
2275  // now convert from the map to the linear vector. make sure every
2276  // entry really appeared in the map
2277  for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
2278  {
2279  Assert(x_support_points.find(i) != x_support_points.end(),
2280  ExcInternalError());
2281 
2282  support_points[i] = x_support_points.find(i)->second;
2283  }
2284 
2285  return support_points;
2286  }
2287  } // namespace
2288  } // namespace internal
2289 
2290 
2291  template <int dim, int spacedim>
2292  void
2294  const DoFHandler<dim, spacedim> &dof_handler,
2295  std::vector<Point<spacedim>> &support_points,
2296  const ComponentMask &mask)
2297  {
2298  AssertDimension(support_points.size(), dof_handler.n_dofs());
2299  Assert((dynamic_cast<
2301  &dof_handler.get_triangulation()) == nullptr),
2302  ExcMessage(
2303  "This function can not be used with distributed triangulations. "
2304  "See the documentation for more information."));
2305 
2306  // Let the internal function do all the work, just make sure that it
2307  // gets a MappingCollection
2308  const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2309 
2310  support_points =
2311  internal::map_dofs_to_support_points_vector(mapping_collection,
2312  dof_handler,
2313  mask);
2314  }
2315 
2316 
2317  template <int dim, int spacedim>
2318  void
2320  const hp::MappingCollection<dim, spacedim> &mapping,
2321  const DoFHandler<dim, spacedim> &dof_handler,
2322  std::vector<Point<spacedim>> &support_points,
2323  const ComponentMask &mask)
2324  {
2325  AssertDimension(support_points.size(), dof_handler.n_dofs());
2326  Assert((dynamic_cast<
2328  &dof_handler.get_triangulation()) == nullptr),
2329  ExcMessage(
2330  "This function can not be used with distributed triangulations. "
2331  "See the documentation for more information."));
2332 
2333  // Let the internal function do all the work, just make sure that it
2334  // gets a MappingCollection
2335  support_points =
2336  internal::map_dofs_to_support_points_vector(mapping, dof_handler, mask);
2337  }
2338 
2339 
2340  // This function is deprecated:
2341  template <int dim, int spacedim>
2342  DEAL_II_DEPRECATED void
2344  const Mapping<dim, spacedim> &mapping,
2345  const DoFHandler<dim, spacedim> &dof_handler,
2346  std::map<types::global_dof_index, Point<spacedim>> &support_points,
2347  const ComponentMask &mask)
2348  {
2349  support_points.clear();
2350 
2351  // Let the internal function do all the work, just make sure that it
2352  // gets a MappingCollection
2353  const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2354 
2355  support_points = internal::map_dofs_to_support_points(mapping_collection,
2356  dof_handler,
2357  mask);
2358  }
2359 
2360 
2361  // This function is deprecated:
2362  template <int dim, int spacedim>
2363  DEAL_II_DEPRECATED void
2365  const hp::MappingCollection<dim, spacedim> &mapping,
2366  const DoFHandler<dim, spacedim> &dof_handler,
2367  std::map<types::global_dof_index, Point<spacedim>> &support_points,
2368  const ComponentMask &mask)
2369  {
2370  support_points.clear();
2371 
2372  // Let the internal function do all the work, just make sure that it
2373  // gets a MappingCollection
2374  support_points =
2375  internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2376  }
2377 
2378 
2379  template <int dim, int spacedim>
2380  std::map<types::global_dof_index, Point<spacedim>>
2382  const DoFHandler<dim, spacedim> &dof_handler,
2383  const ComponentMask &mask)
2384  {
2385  // Let the internal function do all the work, just make sure that it
2386  // gets a MappingCollection
2387  const hp::MappingCollection<dim, spacedim> mapping_collection(mapping);
2388 
2389  return internal::map_dofs_to_support_points(mapping_collection,
2390  dof_handler,
2391  mask);
2392  }
2393 
2394 
2395  template <int dim, int spacedim>
2396  std::map<types::global_dof_index, Point<spacedim>>
2398  const hp::MappingCollection<dim, spacedim> &mapping,
2399  const DoFHandler<dim, spacedim> &dof_handler,
2400  const ComponentMask &mask)
2401  {
2402  return internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2403  }
2404 
2405 
2406  template <int spacedim>
2407  void
2409  std::ostream &out,
2410  const std::map<types::global_dof_index, Point<spacedim>> &support_points)
2411  {
2412  AssertThrow(out.fail() == false, ExcIO());
2413 
2414  std::map<Point<spacedim>,
2415  std::vector<types::global_dof_index>,
2417  point_map;
2418 
2419  // convert to map point -> list of DoFs
2420  for (const auto &it : support_points)
2421  {
2422  std::vector<types::global_dof_index> &v = point_map[it.second];
2423  v.push_back(it.first);
2424  }
2425 
2426  // print the newly created map:
2427  for (const auto &it : point_map)
2428  {
2429  out << it.first << " \"";
2430  const std::vector<types::global_dof_index> &v = it.second;
2431  for (unsigned int i = 0; i < v.size(); ++i)
2432  {
2433  if (i > 0)
2434  out << ", ";
2435  out << v[i];
2436  }
2437 
2438  out << "\"\n";
2439  }
2440 
2441  out << std::flush;
2442  }
2443 
2444 
2445  template <int dim, int spacedim>
2446  void
2448  const Table<2, Coupling> &table,
2449  std::vector<Table<2, Coupling>> &tables_by_block)
2450  {
2451  if (dof_handler.has_hp_capabilities() == false)
2452  {
2453  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
2454  const unsigned int nb = fe.n_blocks();
2455 
2456  tables_by_block.resize(1);
2457  tables_by_block[0].reinit(nb, nb);
2458  tables_by_block[0].fill(none);
2459 
2460  for (unsigned int i = 0; i < fe.n_components(); ++i)
2461  {
2462  const unsigned int ib = fe.component_to_block_index(i);
2463  for (unsigned int j = 0; j < fe.n_components(); ++j)
2464  {
2465  const unsigned int jb = fe.component_to_block_index(j);
2466  tables_by_block[0](ib, jb) |= table(i, j);
2467  }
2468  }
2469  }
2470  else
2471  {
2472  const hp::FECollection<dim> &fe_collection =
2473  dof_handler.get_fe_collection();
2474  tables_by_block.resize(fe_collection.size());
2475 
2476  for (unsigned int f = 0; f < fe_collection.size(); ++f)
2477  {
2478  const FiniteElement<dim, spacedim> &fe = fe_collection[f];
2479 
2480  const unsigned int nb = fe.n_blocks();
2481  tables_by_block[f].reinit(nb, nb);
2482  tables_by_block[f].fill(none);
2483  for (unsigned int i = 0; i < fe.n_components(); ++i)
2484  {
2485  const unsigned int ib = fe.component_to_block_index(i);
2486  for (unsigned int j = 0; j < fe.n_components(); ++j)
2487  {
2488  const unsigned int jb = fe.component_to_block_index(j);
2489  tables_by_block[f](ib, jb) |= table(i, j);
2490  }
2491  }
2492  }
2493  }
2494  }
2495 
2496 
2497 
2498  template <int dim, int spacedim>
2499  void
2501  const DoFHandler<dim, spacedim> &dof_handler,
2502  const unsigned int level,
2503  const std::vector<bool> &selected_dofs,
2504  const types::global_dof_index offset)
2505  {
2506  std::vector<types::global_dof_index> indices;
2507 
2508  unsigned int i = 0;
2509 
2510  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2511  if (cell->is_locally_owned_on_level())
2512  ++i;
2513  block_list.reinit(i,
2514  dof_handler.n_dofs(),
2515  dof_handler.get_fe().n_dofs_per_cell());
2516  i = 0;
2517  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2518  if (cell->is_locally_owned_on_level())
2519  {
2520  indices.resize(cell->get_fe().n_dofs_per_cell());
2521  cell->get_mg_dof_indices(indices);
2522 
2523  if (selected_dofs.size() != 0)
2524  AssertDimension(indices.size(), selected_dofs.size());
2525 
2526  for (types::global_dof_index j = 0; j < indices.size(); ++j)
2527  {
2528  if (selected_dofs.empty())
2529  block_list.add(i, indices[j] - offset);
2530  else
2531  {
2532  if (selected_dofs[j])
2533  block_list.add(i, indices[j] - offset);
2534  }
2535  }
2536  ++i;
2537  }
2538  }
2539 
2540 
2541  template <int dim, int spacedim>
2542  void
2544  const DoFHandler<dim, spacedim> &dof_handler,
2545  const unsigned int level,
2546  const bool interior_only)
2547  {
2548  const FiniteElement<dim> &fe = dof_handler.get_fe();
2549  block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level));
2550 
2551  std::vector<types::global_dof_index> indices;
2552  std::vector<bool> exclude;
2553 
2554  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2555  {
2556  indices.resize(cell->get_fe().n_dofs_per_cell());
2557  cell->get_mg_dof_indices(indices);
2558 
2559  if (interior_only)
2560  {
2561  // Exclude degrees of freedom on faces opposite to the vertex
2562  exclude.resize(fe.n_dofs_per_cell());
2563  std::fill(exclude.begin(), exclude.end(), false);
2564 
2565  for (const unsigned int face : cell->face_indices())
2566  if (cell->at_boundary(face) ||
2567  cell->neighbor(face)->level() != cell->level())
2568  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2569  exclude[fe.face_to_cell_index(i, face)] = true;
2570  for (types::global_dof_index j = 0; j < indices.size(); ++j)
2571  if (!exclude[j])
2572  block_list.add(0, indices[j]);
2573  }
2574  else
2575  {
2576  for (const auto index : indices)
2577  block_list.add(0, index);
2578  }
2579  }
2580  }
2581 
2582 
2583  template <int dim, int spacedim>
2584  void
2586  const DoFHandler<dim, spacedim> &dof_handler,
2587  const unsigned int level,
2588  const bool interior_dofs_only,
2589  const bool boundary_dofs)
2590  {
2591  Assert(level > 0 && level < dof_handler.get_triangulation().n_levels(),
2592  ExcIndexRange(level, 1, dof_handler.get_triangulation().n_levels()));
2593 
2594  std::vector<types::global_dof_index> indices;
2595  std::vector<bool> exclude;
2596 
2597  unsigned int block = 0;
2598  for (const auto &pcell : dof_handler.cell_iterators_on_level(level - 1))
2599  {
2600  if (pcell->is_active())
2601  continue;
2602 
2603  for (unsigned int child = 0; child < pcell->n_children(); ++child)
2604  {
2605  const auto cell = pcell->child(child);
2606 
2607  // For hp, only this line here would have to be replaced.
2608  const FiniteElement<dim> &fe = dof_handler.get_fe();
2609  const unsigned int n_dofs = fe.n_dofs_per_cell();
2610  indices.resize(n_dofs);
2611  exclude.resize(n_dofs);
2612  std::fill(exclude.begin(), exclude.end(), false);
2613  cell->get_mg_dof_indices(indices);
2614 
2615  if (interior_dofs_only)
2616  {
2617  // Eliminate dofs on faces of the child which are on faces
2618  // of the parent
2619  for (unsigned int d = 0; d < dim; ++d)
2620  {
2621  const unsigned int face =
2623  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2624  exclude[fe.face_to_cell_index(i, face)] = true;
2625  }
2626 
2627  // Now remove all degrees of freedom on the domain boundary
2628  // from the exclusion list
2629  if (boundary_dofs)
2630  for (const unsigned int face :
2632  if (cell->at_boundary(face))
2633  for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
2634  ++i)
2635  exclude[fe.face_to_cell_index(i, face)] = false;
2636  }
2637 
2638  for (unsigned int i = 0; i < n_dofs; ++i)
2639  if (!exclude[i])
2640  block_list.add(block, indices[i]);
2641  }
2642  ++block;
2643  }
2644  }
2645 
2646  template <int dim, int spacedim>
2647  std::vector<unsigned int>
2649  const DoFHandler<dim, spacedim> &dof_handler,
2650  const unsigned int level,
2651  const bool interior_only,
2652  const bool boundary_patches,
2653  const bool level_boundary_patches,
2654  const bool single_cell_patches,
2655  const bool invert_vertex_mapping)
2656  {
2657  const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
2658  BlockMask exclude_boundary_dofs = BlockMask(n_blocks, interior_only);
2659  return make_vertex_patches(block_list,
2660  dof_handler,
2661  level,
2662  exclude_boundary_dofs,
2663  boundary_patches,
2664  level_boundary_patches,
2665  single_cell_patches,
2666  invert_vertex_mapping);
2667  }
2668 
2669  template <int dim, int spacedim>
2670  std::vector<unsigned int>
2672  const DoFHandler<dim, spacedim> &dof_handler,
2673  const unsigned int level,
2674  const BlockMask &exclude_boundary_dofs,
2675  const bool boundary_patches,
2676  const bool level_boundary_patches,
2677  const bool single_cell_patches,
2678  const bool invert_vertex_mapping)
2679  {
2680  // Vector mapping from vertex index in the triangulation to consecutive
2681  // block indices on this level The number of cells at a vertex
2682  std::vector<unsigned int> vertex_cell_count(
2683  dof_handler.get_triangulation().n_vertices(), 0);
2684 
2685  // Is a vertex at the boundary?
2686  std::vector<bool> vertex_boundary(
2687  dof_handler.get_triangulation().n_vertices(), false);
2688 
2689  std::vector<unsigned int> vertex_mapping(
2690  dof_handler.get_triangulation().n_vertices(),
2692 
2693  // Estimate for the number of dofs at this point
2694  std::vector<unsigned int> vertex_dof_count(
2695  dof_handler.get_triangulation().n_vertices(), 0);
2696 
2697  // Identify all vertices active on this level and remember some data
2698  // about them
2699  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2700  for (const unsigned int v : cell->vertex_indices())
2701  {
2702  const unsigned int vg = cell->vertex_index(v);
2703  vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2704  ++vertex_cell_count[vg];
2705  for (unsigned int d = 0; d < dim; ++d)
2706  {
2707  const unsigned int face = GeometryInfo<dim>::vertex_to_face[v][d];
2708  if (cell->at_boundary(face))
2709  vertex_boundary[vg] = true;
2710  else if ((!level_boundary_patches) &&
2711  (cell->neighbor(face)->level() !=
2712  static_cast<int>(level)))
2713  vertex_boundary[vg] = true;
2714  }
2715  }
2716  // From now on, only vertices with positive dof count are "in".
2717 
2718  // Remove vertices at boundaries or in corners
2719  for (unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2720  if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2721  (!boundary_patches && vertex_boundary[vg]))
2722  vertex_dof_count[vg] = 0;
2723 
2724  // Create a mapping from all vertices to the ones used here
2725  unsigned int n_vertex_count = 0;
2726  for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2727  if (vertex_dof_count[vg] != 0)
2728  vertex_mapping[vg] = n_vertex_count++;
2729 
2730  // Compactify dof count
2731  for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2732  if (vertex_dof_count[vg] != 0)
2733  vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2734 
2735  // Now that we have all the data, we reduce it to the part we actually
2736  // want
2737  vertex_dof_count.resize(n_vertex_count);
2738 
2739  // At this point, the list of patches is ready. Now we enter the dofs
2740  // into the sparsity pattern.
2741  block_list.reinit(vertex_dof_count.size(),
2742  dof_handler.n_dofs(level),
2743  vertex_dof_count);
2744 
2745  std::vector<types::global_dof_index> indices;
2746  std::vector<bool> exclude;
2747 
2748  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2749  {
2750  const FiniteElement<dim> &fe = cell->get_fe();
2751  indices.resize(fe.n_dofs_per_cell());
2752  cell->get_mg_dof_indices(indices);
2753 
2754  for (const unsigned int v : cell->vertex_indices())
2755  {
2756  const unsigned int vg = cell->vertex_index(v);
2757  const unsigned int block = vertex_mapping[vg];
2758  if (block == numbers::invalid_unsigned_int)
2759  continue;
2760 
2761  // Collect excluded dofs for some block(s) if boundary dofs
2762  // for a block are decided to be excluded
2763  if (exclude_boundary_dofs.size() == 0 ||
2764  exclude_boundary_dofs.n_selected_blocks() != 0)
2765  {
2766  // Exclude degrees of freedom on faces opposite to the
2767  // vertex
2768  exclude.resize(fe.n_dofs_per_cell());
2769  std::fill(exclude.begin(), exclude.end(), false);
2770 
2771  for (unsigned int d = 0; d < dim; ++d)
2772  {
2773  const unsigned int a_face =
2775  const unsigned int face =
2777  for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
2778  {
2779  // For each dof, get the block it is in and decide to
2780  // exclude it or not
2781  if (exclude_boundary_dofs[fe.system_to_block_index(
2782  fe.face_to_cell_index(
2783  i, face))
2784  .first] == true)
2785  exclude[fe.face_to_cell_index(i, face)] = true;
2786  }
2787  }
2788  for (unsigned int j = 0; j < indices.size(); ++j)
2789  if (!exclude[j])
2790  block_list.add(block, indices[j]);
2791  }
2792  else
2793  {
2794  for (const auto index : indices)
2795  block_list.add(block, index);
2796  }
2797  }
2798  }
2799 
2800  if (invert_vertex_mapping)
2801  {
2802  // Compress vertex mapping
2803  unsigned int n_vertex_count = 0;
2804  for (unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2805  if (vertex_mapping[vg] != numbers::invalid_unsigned_int)
2806  vertex_mapping[n_vertex_count++] = vg;
2807 
2808  // Now we reduce it to the part we actually want
2809  vertex_mapping.resize(n_vertex_count);
2810  }
2811 
2812  return vertex_mapping;
2813  }
2814 
2815 
2816  template <int dim, int spacedim>
2817  unsigned int
2819  const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2820  &patch)
2821  {
2822  std::set<types::global_dof_index> dofs_on_patch;
2823  std::vector<types::global_dof_index> local_dof_indices;
2824 
2825  // loop over the cells in the patch and get the DoFs on each.
2826  // add all of them to a std::set which automatically makes sure
2827  // all duplicates are ignored
2828  for (unsigned int i = 0; i < patch.size(); ++i)
2829  {
2831  patch[i];
2832  Assert(cell->is_artificial() == false,
2833  ExcMessage("This function can not be called with cells that are "
2834  "not either locally owned or ghost cells."));
2835  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2836  cell->get_dof_indices(local_dof_indices);
2837  dofs_on_patch.insert(local_dof_indices.begin(),
2838  local_dof_indices.end());
2839  }
2840 
2841  // now return the number of DoFs (duplicates were ignored)
2842  return dofs_on_patch.size();
2843  }
2844 
2845 
2846 
2847  template <int dim, int spacedim>
2848  std::vector<types::global_dof_index>
2850  const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2851  &patch)
2852  {
2853  std::set<types::global_dof_index> dofs_on_patch;
2854  std::vector<types::global_dof_index> local_dof_indices;
2855 
2856  // loop over the cells in the patch and get the DoFs on each.
2857  // add all of them to a std::set which automatically makes sure
2858  // all duplicates are ignored
2859  for (unsigned int i = 0; i < patch.size(); ++i)
2860  {
2862  patch[i];
2863  Assert(cell->is_artificial() == false,
2864  ExcMessage("This function can not be called with cells that are "
2865  "not either locally owned or ghost cells."));
2866  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2867  cell->get_dof_indices(local_dof_indices);
2868  dofs_on_patch.insert(local_dof_indices.begin(),
2869  local_dof_indices.end());
2870  }
2871 
2872  Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2873  ExcInternalError());
2874 
2875  // return a vector with the content of the set above. copying
2876  // also ensures that we retain sortedness as promised in the
2877  // documentation and as necessary to retain the block structure
2878  // also on the local system
2879  return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2880  dofs_on_patch.end());
2881  }
2882 
2883 
2884 } // end of namespace DoFTools
2885 
2886 
2887 
2888 // explicit instantiations
2889 
2890 #include "dof_tools.inst"
2891 
2892 
2893 
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
std::vector< types::fe_index > get_active_fe_indices() const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
cell_iterator begin(const unsigned int level=0) const
const IndexSet & locally_owned_dofs() const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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, unsigned int > system_to_component_index(const unsigned int index) const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) 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
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
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
unsigned int component_to_block_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
unsigned int n_base_elements() const
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
size_type n_elements() const
Definition: index_set.h:1854
bool is_element(const size_type index) const
Definition: index_set.h:1814
void add_index(const size_type index)
Definition: index_set.h:1715
void fill_binary_vector(VectorType &vector) const
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1902
void compress() const
Definition: index_set.h:1704
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1751
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
Definition: vector.h:110
virtual size_type size() const override
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)
Definition: q_collection.h:216
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
unsigned int level
Definition: grid_out.cc:4617
unsigned int cell_index
Definition: grid_tools.cc:1192
IteratorRange< FilteredIterator< BaseIterator > > filter_iterators(IteratorRange< BaseIterator > i, const Predicate &p)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
@ 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)
Definition: dof_tools.cc:1803
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)
Definition: dof_tools.cc:1596
IndexSet dof_indices_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition: dof_tools.cc:1703
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={})
Definition: dof_tools.cc:796
void write_gnuplot_dof_support_point_info(std::ostream &out, const std::map< types::global_dof_index, Point< spacedim >> &support_points)
Definition: dof_tools.cc:2408
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1425
IndexSet extract_boundary_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask={}, const std::set< types::boundary_id > &boundary_ids={})
Definition: dof_tools.cc:586
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1136
std::map< typename DoFHandler< dim - 1, spacedim >::active_cell_iterator, std::pair< typename DoFHandler< dim, spacedim >::active_cell_iterator, unsigned int > > map_boundary_to_bulk_dof_iterators(const std::map< typename Triangulation< dim - 1, spacedim >::cell_iterator, typename Triangulation< dim, spacedim >::face_iterator > &c1_to_c0, const DoFHandler< dim, spacedim > &c0_dh, const DoFHandler< dim - 1, spacedim > &c1_dh)
Definition: dof_tools.cc:1360
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={})
Definition: dof_tools.cc:2293
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1412
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)
Definition: dof_tools.cc:2500
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)
Definition: dof_tools.cc:1044
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)
Definition: dof_tools.cc:2447
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
Definition: dof_tools.cc:1184
void extract_subdomain_dofs(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain_id, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:1013
IndexSet extract_locally_active_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
Definition: dof_tools.cc:1087
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)
Definition: dof_tools.cc:2648
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 >())
Definition: dof_tools.cc:2017
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={})
Definition: dof_tools.cc:1925
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)
Definition: dof_tools.cc:2585
void map_dof_to_boundary_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &mapping)
Definition: dof_tools.cc:2106
std::vector< IndexSet > locally_relevant_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1520
void make_single_patch(SparsityPattern &block_list, const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool interior_dofs_only=false)
Definition: dof_tools.cc:2543
std::vector< types::global_dof_index > get_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
Definition: dof_tools.cc:2849
std::vector< IndexSet > locally_owned_dofs_per_component(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &components={})
Definition: dof_tools.cc:439
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
Definition: dof_tools.cc:1686
unsigned int count_dofs_on_patch(const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &patch)
Definition: dof_tools.cc:2818
IndexSet extract_hanging_node_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1004
void extract_constant_modes(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool >> &constant_modes)
Definition: dof_tools.cc:1242
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
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)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition: mpi.cc:149
std::string compress(const std::string &input)
Definition: utilities.cc:390
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
unsigned int global_dof_index
Definition: types.h:82
unsigned int subdomain_id
Definition: types.h:44
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
const ::Triangulation< dim, spacedim > & tria
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:99