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