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