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