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