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