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