Reference documentation for deal.II version GIT 9fcc068bd8 2023-06-06 14:45:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_collection.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
18 
21 
22 #include <limits>
23 #include <set>
24 
25 
26 
28 
29 namespace hp
30 {
31  template <int dim, int spacedim>
33  {
34  set_default_hierarchy();
35  }
36 
37 
38 
39  template <int dim, int spacedim>
42  : FECollection()
43  {
44  push_back(fe);
45  }
46 
47 
48 
49  template <int dim, int spacedim>
51  const std::vector<const FiniteElement<dim, spacedim> *> &fes)
52  : FECollection()
53  {
54  Assert(fes.size() > 0,
55  ExcMessage("Need to pass at least one finite element."));
56 
57  for (unsigned int i = 0; i < fes.size(); ++i)
58  push_back(*fes[i]);
59  }
60 
61 
62 
63  template <int dim, int spacedim>
64  void
66  const FiniteElement<dim, spacedim> &new_fe)
67  {
68  // check that the new element has the right number of components. only check
69  // with the first element, since all the other elements have already passed
70  // the test against the first element
71  Assert(this->size() == 0 ||
72  new_fe.n_components() == this->operator[](0).n_components(),
73  ExcMessage("All elements inside a collection need to have the "
74  "same number of vector components!"));
75 
76  Collection<FiniteElement<dim, spacedim>>::push_back(new_fe.clone());
77  }
78 
79 
80 
81  template <int dim, int spacedim>
84  {
85  Assert(this->size() > 0, ExcNoFiniteElements());
86 
87  // Since we can only add elements to an FECollection, we are safe comparing
88  // the sizes of this object and the MappingCollection. One circumstance that
89  // might lead to their sizes diverging is this:
90  // - An FECollection is created and then this function is called. The shared
91  // map is now initialized.
92  // - A second FECollection is made as a copy of this one. The shared map is
93  // not recreated.
94  // - The second FECollection is then resized by adding a new FE. The shared
95  // map is thus invalid for the second instance.
96  if (!reference_cell_default_linear_mapping ||
97  reference_cell_default_linear_mapping->size() != this->size())
98  {
99  auto &this_nc = const_cast<FECollection<dim, spacedim> &>(*this);
100 
102  std::make_shared<MappingCollection<dim, spacedim>>();
103 
104  for (const auto &fe : *this)
105  this_nc.reference_cell_default_linear_mapping->push_back(
106  fe.reference_cell()
107  .template get_default_linear_mapping<dim, spacedim>());
108  }
109 
110  return *reference_cell_default_linear_mapping;
111  }
112 
113 
114 
115  template <int dim, int spacedim>
116  std::set<unsigned int>
118  const std::set<unsigned int> &fes,
119  const unsigned int codim) const
120  {
121 #ifdef DEBUG
122  // Validate user inputs.
123  Assert(codim <= dim, ExcImpossibleInDim(dim));
124  Assert(this->size() > 0, ExcEmptyObject());
125  for (const auto &fe : fes)
126  AssertIndexRange(fe, this->size());
127 #endif
128 
129  // Check if any element of this FECollection is able to dominate all
130  // elements of @p fes. If one was found, we add it to the set of
131  // dominating elements.
132  std::set<unsigned int> dominating_fes;
133  for (unsigned int current_fe = 0; current_fe < this->size(); ++current_fe)
134  {
135  // Check if current_fe can dominate all elements in @p fes.
138  for (const auto &other_fe : fes)
139  domination =
140  domination &
141  this->operator[](current_fe)
142  .compare_for_domination(this->operator[](other_fe), codim);
143 
144  // If current_fe dominates, add it to the set.
147  /*covers cases like {Q2,Q3,Q1,Q1} with fes={2,3}*/))
148  dominating_fes.insert(current_fe);
149  }
150  return dominating_fes;
151  }
152 
153 
154 
155  template <int dim, int spacedim>
156  std::set<unsigned int>
158  const std::set<unsigned int> &fes,
159  const unsigned int codim) const
160  {
161 #ifdef DEBUG
162  // Validate user inputs.
163  Assert(codim <= dim, ExcImpossibleInDim(dim));
164  Assert(this->size() > 0, ExcEmptyObject());
165  for (const auto &fe : fes)
166  AssertIndexRange(fe, this->size());
167 #endif
168 
169  // Check if any element of this FECollection is dominated by all
170  // elements of @p fes. If one was found, we add it to the set of
171  // dominated elements.
172  std::set<unsigned int> dominated_fes;
173  for (unsigned int current_fe = 0; current_fe < this->size(); ++current_fe)
174  {
175  // Check if current_fe is dominated by all other elements in @p fes.
178  for (const auto &other_fe : fes)
179  domination =
180  domination &
181  this->operator[](current_fe)
182  .compare_for_domination(this->operator[](other_fe), codim);
183 
184  // If current_fe is dominated, add it to the set.
187  /*covers cases like {Q2,Q3,Q1,Q1} with fes={2,3}*/))
188  dominated_fes.insert(current_fe);
189  }
190  return dominated_fes;
191  }
192 
193 
194 
195  template <int dim, int spacedim>
196  unsigned int
198  const std::set<unsigned int> &fes,
199  const unsigned int codim) const
200  {
201  // If the set of elements contains only a single element,
202  // then this very element is considered to be the dominating one.
203  if (fes.size() == 1)
204  return *fes.begin();
205 
206 #ifdef DEBUG
207  // Validate user inputs.
208  Assert(codim <= dim, ExcImpossibleInDim(dim));
209  Assert(this->size() > 0, ExcEmptyObject());
210  for (const auto &fe : fes)
211  AssertIndexRange(fe, this->size());
212 #endif
213 
214  // There may also be others, in which case we'll check if any of these
215  // elements is able to dominate all others. If one was found, we stop
216  // looking further and return the dominating element.
217  for (const auto &current_fe : fes)
218  {
219  // Check if current_fe can dominate all elements in @p fes.
222  for (const auto &other_fe : fes)
223  if (current_fe != other_fe)
224  domination =
225  domination &
226  this->operator[](current_fe)
227  .compare_for_domination(this->operator[](other_fe), codim);
228 
229  // If current_fe dominates, return its index.
232  /*covers cases like {Q2,Q3,Q1,Q1} with fes={2,3}*/))
233  return current_fe;
234  }
235 
236  // If we couldn't find the dominating object, return an invalid one.
238  }
239 
240 
241 
242  template <int dim, int spacedim>
243  unsigned int
245  const std::set<unsigned int> &fes,
246  const unsigned int codim) const
247  {
248  // If the set of elements contains only a single element,
249  // then this very element is considered to be the dominated one.
250  if (fes.size() == 1)
251  return *fes.begin();
252 
253 #ifdef DEBUG
254  // Validate user inputs.
255  Assert(codim <= dim, ExcImpossibleInDim(dim));
256  Assert(this->size() > 0, ExcEmptyObject());
257  for (const auto &fe : fes)
258  AssertIndexRange(fe, this->size());
259 #endif
260 
261  // There may also be others, in which case we'll check if any of these
262  // elements is dominated by all others. If one was found, we stop
263  // looking further and return the dominated element.
264  for (const auto &current_fe : fes)
265  {
266  // Check if current_fe is dominated by all other elements in @p fes.
269  for (const auto &other_fe : fes)
270  if (current_fe != other_fe)
271  domination =
272  domination &
273  this->operator[](current_fe)
274  .compare_for_domination(this->operator[](other_fe), codim);
275 
276  // If current_fe is dominated, return its index.
279  /*covers cases like {Q2,Q3,Q1,Q1} with fes={2,3}*/))
280  return current_fe;
281  }
282 
283  // If we couldn't find the dominated object, return an invalid one.
285  }
286 
287 
288 
289  template <int dim, int spacedim>
290  unsigned int
292  const std::set<unsigned int> &fes,
293  const unsigned int codim) const
294  {
295  unsigned int fe_index = find_dominating_fe(fes, codim);
296 
298  {
299  const std::set<unsigned int> dominating_fes =
300  find_common_fes(fes, codim);
301  fe_index = find_dominated_fe(dominating_fes, codim);
302  }
303 
304  return fe_index;
305  }
306 
307 
308 
309  template <int dim, int spacedim>
310  unsigned int
312  const std::set<unsigned int> &fes,
313  const unsigned int codim) const
314  {
315  unsigned int fe_index = find_dominated_fe(fes, codim);
316 
318  {
319  const std::set<unsigned int> dominated_fes =
320  find_enclosing_fes(fes, codim);
321  fe_index = find_dominating_fe(dominated_fes, codim);
322  }
323 
324  return fe_index;
325  }
326 
327 
328 
329  namespace
330  {
335  std::vector<std::map<unsigned int, unsigned int>>
336  compute_hp_dof_identities(
337  const std::set<unsigned int> &fes,
338  const std::function<std::vector<std::pair<unsigned int, unsigned int>>(
339  const unsigned int,
340  const unsigned int)> & query_identities)
341  {
342  // Let's deal with the easy cases first. If the set of fe indices is empty
343  // or has only one entry, then there are no identities:
344  if (fes.size() <= 1)
345  return {};
346 
347  // If the set has two entries, then the
348  // FiniteElement::hp_*_dof_identities() function directly returns what we
349  // need. We just need to prefix its output with the respective fe indices:
350  if (fes.size() == 2)
351  {
352  const unsigned int fe_index_1 = *fes.begin();
353  const unsigned int fe_index_2 = *(++fes.begin());
354  const auto reduced_identities =
355  query_identities(fe_index_1, fe_index_2);
356 
357  std::vector<std::map<unsigned int, unsigned int>> complete_identities;
358 
359  for (const auto &reduced_identity : reduced_identities)
360  {
361  // Each identity returned by query_identities() is a pair of
362  // dof indices. Prefix each with its fe index and put the result
363  // into a vector
364  std::map<unsigned int, unsigned int> complete_identity = {
365  {fe_index_1, reduced_identity.first},
366  {fe_index_2, reduced_identity.second}};
367  complete_identities.emplace_back(std::move(complete_identity));
368  }
369 
370  return complete_identities;
371  }
372 
373  // Now for the general case of three or more elements:
374  //
375  // Consider all degrees of freedom of the identified elements (represented
376  // via (fe_index,dof_index) pairs) as the nodes in a graph. Then draw
377  // edges for all DoFs that are identified based on what the elements
378  // selected in the argument say. Let us first build this graph, where we
379  // only store the edges of the graph, and as a consequence ignore nodes
380  // (DoFs) that simply don't show up at all in any of the identities:
381  using Node = std::pair<unsigned int, unsigned int>;
382  using Edge = std::pair<Node, Node>;
383  using Graph = std::set<Edge>;
384 
385  Graph identities_graph;
386  for (const unsigned int fe_index_1 : fes)
387  for (const unsigned int fe_index_2 : fes)
388  if (fe_index_1 != fe_index_2)
389  for (const auto &identity :
390  query_identities(fe_index_1, fe_index_2))
391  identities_graph.emplace(Node(fe_index_1, identity.first),
392  Node(fe_index_2, identity.second));
393 
394 #ifdef DEBUG
395  // Now verify that indeed the graph is symmetric: If one element
396  // declares that certain ones of its DoFs are to be unified with those
397  // of the other, then the other one should agree with this. As a
398  // consequence of this test succeeding, we know that the graph is actually
399  // undirected.
400  for (const auto &edge : identities_graph)
401  Assert(identities_graph.find({edge.second, edge.first}) !=
402  identities_graph.end(),
403  ExcInternalError());
404 #endif
405 
406  // The next step is that we ought to verify that if there is an identity
407  // between (fe1,dof1) and (fe2,dof2), as well as with (fe2,dof2) and
408  // (fe3,dof3), then there should also be an identity between (fe1,dof1)
409  // and (fe3,dof3). The same logic should apply to chains of four
410  // identities.
411  //
412  // This means that the graph we have built above must be composed of a
413  // collection of complete sub-graphs (complete = each possible edge in the
414  // sub-graph exists) -- or, using a different term, that the graph
415  // consists of a number of "cliques". Each of these cliques is then one
416  // extended identity between two or more DoFs, and these are the ones that
417  // we will want to return.
418  //
419  // To ascertain that this is true, and to figure out what we want to
420  // return, we need to divide the graph into its sub-graphs and then ensure
421  // that each sub-graph is indeed a clique. This is slightly cumbersome,
422  // but can be done as follows:
423  // - pick one edge 'e' of G
424  // - add e=(n1,n2) to the sub-graph SG
425  // - set N={n1,n2}
426  // - loop over the remainder of the edges 'e' of the graph:
427  // - if 'e' has one or both nodes in N:
428  // - add 'e' to SG and
429  // - add its two nodes to N (they may already be in there)
430  // - remove 'e' from G
431  //
432  // In general, this may not find the entire sub-graph if the edges are
433  // stored in random order. For example, if the graph consisted of the
434  // following edges in this order:
435  // (a,b)
436  // (c,d)
437  // (a,c)
438  // (a,d)
439  // (b,c)
440  // (b,d)
441  // then the graph itself is a clique, but the algorithm outlined above
442  // would skip the edge (c,d) because neither of the nodes are already
443  // in the set N which at that point is still (a,b).
444  //
445  // But, we store the graph with a std::set, which stored edges in sorted
446  // order where the order is the lexicographic order of nodes. This ensures
447  // that we really capture all edges that correspond to a sub-graph (but
448  // we will assert this as well).
449  //
450  // (For programmatic reasons, we skip the removal of 'e' from G in a first
451  // run through because it modifies the graph and thus invalidates
452  // iterators. But because SG stores all of these edges, we can remove them
453  // all from G after collecting the edges in SG.)
454  std::vector<std::map<unsigned int, unsigned int>> identities;
455  while (identities_graph.size() > 0)
456  {
457  Graph sub_graph; // SG
458  std::set<Node> sub_graph_nodes; // N
459 
460  sub_graph.emplace(*identities_graph.begin());
461  sub_graph_nodes.emplace(identities_graph.begin()->first);
462  sub_graph_nodes.emplace(identities_graph.begin()->second);
463 
464  for (const Edge &e : identities_graph)
465  if ((sub_graph_nodes.find(e.first) != sub_graph_nodes.end()) ||
466  (sub_graph_nodes.find(e.second) != sub_graph_nodes.end()))
467  {
468  sub_graph.insert(e);
469  sub_graph_nodes.insert(e.first);
470  sub_graph_nodes.insert(e.second);
471  }
472 
473  // We have now obtained a sub-graph from the overall graph.
474  // Now remove it from the bigger graph
475  for (const Edge &e : sub_graph)
476  identities_graph.erase(e);
477 
478 #ifdef DEBUG
479  // There are three checks we ought to perform:
480  // - That the sub-graph is undirected, i.e. that every edge appears
481  // in both directions
482  for (const auto &edge : sub_graph)
483  Assert(sub_graph.find({edge.second, edge.first}) != sub_graph.end(),
484  ExcInternalError());
485 
486  // - None of the nodes in the sub-graph should have appeared in
487  // any of the other sub-graphs. If they did, then we have a bug
488  // in extracting sub-graphs. This is actually more easily checked
489  // the other way around: none of the nodes of the sub-graph we
490  // just extracted should be in any of the edges of the *remaining*
491  // graph
492  for (const Node &n : sub_graph_nodes)
493  for (const Edge &e : identities_graph)
494  Assert((n != e.first) && (n != e.second), ExcInternalError());
495  // - Second, the sub-graph we just extracted needs to be complete,
496  // i.e.,
497  // be a "clique". We check this by counting how many edges it has.
498  // for 'n' nodes in 'N', we need to have n*(n-1) edges (we store
499  // both directed edges).
500  Assert(sub_graph.size() ==
501  sub_graph_nodes.size() * (sub_graph_nodes.size() - 1),
502  ExcInternalError());
503 #endif
504 
505  // At this point we're sure that we have extracted a complete
506  // sub-graph ("clique"). The DoFs involved are all identical then, and
507  // we will store this identity so we can return it later.
508  //
509  // The sub-graph is given as a set of Node objects, which is just
510  // a collection of (fe_index,dof_index) pairs. Because each
511  // fe_index can only appear once there, a better data structure
512  // is a std::map from fe_index to dof_index, which can conveniently
513  // be initialized from a range of iterators to pairs:
514  identities.emplace_back(sub_graph_nodes.begin(),
515  sub_graph_nodes.end());
516  Assert(identities.back().size() == sub_graph_nodes.size(),
517  ExcInternalError());
518  }
519 
520  return identities;
521  }
522  } // namespace
523 
524 
525 
526  template <int dim, int spacedim>
527  std::vector<std::map<unsigned int, unsigned int>>
529  const std::set<unsigned int> &fes) const
530  {
531  auto query_vertex_dof_identities = [this](const unsigned int fe_index_1,
532  const unsigned int fe_index_2) {
533  return (*this)[fe_index_1].hp_vertex_dof_identities((*this)[fe_index_2]);
534  };
535  return compute_hp_dof_identities(fes, query_vertex_dof_identities);
536  }
537 
538 
539 
540  template <int dim, int spacedim>
541  std::vector<std::map<unsigned int, unsigned int>>
543  const std::set<unsigned int> &fes) const
544  {
545  auto query_line_dof_identities = [this](const unsigned int fe_index_1,
546  const unsigned int fe_index_2) {
547  return (*this)[fe_index_1].hp_line_dof_identities((*this)[fe_index_2]);
548  };
549  return compute_hp_dof_identities(fes, query_line_dof_identities);
550  }
551 
552 
553 
554  template <int dim, int spacedim>
555  std::vector<std::map<unsigned int, unsigned int>>
557  const std::set<unsigned int> &fes,
558  const unsigned int face_no) const
559  {
560  auto query_quad_dof_identities = [this,
561  face_no](const unsigned int fe_index_1,
562  const unsigned int fe_index_2) {
563  return (*this)[fe_index_1].hp_quad_dof_identities((*this)[fe_index_2],
564  face_no);
565  };
566  return compute_hp_dof_identities(fes, query_quad_dof_identities);
567  }
568 
569 
570 
571  template <int dim, int spacedim>
572  void
574  const std::function<
575  unsigned int(const typename hp::FECollection<dim, spacedim> &,
576  const unsigned int)> &next,
577  const std::function<
578  unsigned int(const typename hp::FECollection<dim, spacedim> &,
579  const unsigned int)> &prev)
580  {
581  // copy hierarchy functions
582  hierarchy_next = next;
583  hierarchy_prev = prev;
584  }
585 
586 
587 
588  template <int dim, int spacedim>
589  void
591  {
592  // establish hierarchy corresponding to order of indices
593  set_hierarchy(&DefaultHierarchy::next_index,
594  &DefaultHierarchy::previous_index);
595  }
596 
597 
598 
599  template <int dim, int spacedim>
600  std::vector<unsigned int>
602  const unsigned int fe_index) const
603  {
604  AssertIndexRange(fe_index, this->size());
605 
606  std::deque<unsigned int> sequence = {fe_index};
607 
608  // get predecessors
609  {
610  unsigned int front = sequence.front();
611  unsigned int previous;
612  while ((previous = previous_in_hierarchy(front)) != front)
613  {
614  sequence.push_front(previous);
615  front = previous;
616 
617  Assert(sequence.size() <= this->size(),
618  ExcMessage(
619  "The registered hierarchy is not terminated: "
620  "previous_in_hierarchy() does not stop at a final index."));
621  }
622  }
623 
624  // get successors
625  {
626  unsigned int back = sequence.back();
627  unsigned int next;
628  while ((next = next_in_hierarchy(back)) != back)
629  {
630  sequence.push_back(next);
631  back = next;
632 
633  Assert(sequence.size() <= this->size(),
634  ExcMessage(
635  "The registered hierarchy is not terminated: "
636  "next_in_hierarchy() does not stop at a final index."));
637  }
638  }
639 
640  return {sequence.begin(), sequence.end()};
641  }
642 
643 
644 
645  template <int dim, int spacedim>
646  unsigned int
648  const unsigned int fe_index) const
649  {
650  AssertIndexRange(fe_index, this->size());
651 
652  const unsigned int new_fe_index = hierarchy_next(*this, fe_index);
653  AssertIndexRange(new_fe_index, this->size());
654 
655  return new_fe_index;
656  }
657 
658 
659 
660  template <int dim, int spacedim>
661  unsigned int
663  const unsigned int fe_index) const
664  {
665  AssertIndexRange(fe_index, this->size());
666 
667  const unsigned int new_fe_index = hierarchy_prev(*this, fe_index);
668  AssertIndexRange(new_fe_index, this->size());
669 
670  return new_fe_index;
671  }
672 
673 
674 
675  template <int dim, int spacedim>
678  const FEValuesExtractors::Scalar &scalar) const
679  {
680  Assert(this->size() > 0,
681  ExcMessage("This collection contains no finite element."));
682 
683  // get the mask from the first element of the collection
684  const ComponentMask mask = (*this)[0].component_mask(scalar);
685 
686  // but then also verify that the other elements of the collection
687  // would return the same mask
688  for (unsigned int c = 1; c < this->size(); ++c)
689  Assert(mask == (*this)[c].component_mask(scalar), ExcInternalError());
690 
691  return mask;
692  }
693 
694 
695  template <int dim, int spacedim>
698  const FEValuesExtractors::Vector &vector) const
699  {
700  Assert(this->size() > 0,
701  ExcMessage("This collection contains no finite element."));
702 
703  // get the mask from the first element of the collection
704  const ComponentMask mask = (*this)[0].component_mask(vector);
705 
706  // but then also verify that the other elements of the collection
707  // would return the same mask
708  for (unsigned int c = 1; c < this->size(); ++c)
709  Assert(mask == (*this)[c].component_mask(vector), ExcInternalError());
710 
711  return mask;
712  }
713 
714 
715  template <int dim, int spacedim>
718  const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
719  {
720  Assert(this->size() > 0,
721  ExcMessage("This collection contains no finite element."));
722 
723  // get the mask from the first element of the collection
724  const ComponentMask mask = (*this)[0].component_mask(sym_tensor);
725 
726  // but then also verify that the other elements of the collection
727  // would return the same mask
728  for (unsigned int c = 1; c < this->size(); ++c)
729  Assert(mask == (*this)[c].component_mask(sym_tensor), ExcInternalError());
730 
731  return mask;
732  }
733 
734 
735  template <int dim, int spacedim>
738  {
739  Assert(this->size() > 0,
740  ExcMessage("This collection contains no finite element."));
741 
742  // get the mask from the first element of the collection
743  const ComponentMask mask = (*this)[0].component_mask(block_mask);
744 
745  // but then also verify that the other elements of the collection
746  // would return the same mask
747  for (unsigned int c = 1; c < this->size(); ++c)
748  Assert(mask == (*this)[c].component_mask(block_mask),
749  ExcMessage("Not all elements of this collection agree on what "
750  "the appropriate mask should be."));
751 
752  return mask;
753  }
754 
755 
756  template <int dim, int spacedim>
757  BlockMask
759  const FEValuesExtractors::Scalar &scalar) const
760  {
761  Assert(this->size() > 0,
762  ExcMessage("This collection contains no finite element."));
763 
764  // get the mask from the first element of the collection
765  const BlockMask mask = (*this)[0].block_mask(scalar);
766 
767  // but then also verify that the other elements of the collection
768  // would return the same mask
769  for (unsigned int c = 1; c < this->size(); ++c)
770  Assert(mask == (*this)[c].block_mask(scalar),
771  ExcMessage("Not all elements of this collection agree on what "
772  "the appropriate mask should be."));
773 
774  return mask;
775  }
776 
777 
778  template <int dim, int spacedim>
779  BlockMask
781  const FEValuesExtractors::Vector &vector) const
782  {
783  Assert(this->size() > 0,
784  ExcMessage("This collection contains no finite element."));
785 
786  // get the mask from the first element of the collection
787  const BlockMask mask = (*this)[0].block_mask(vector);
788 
789  // but then also verify that the other elements of the collection
790  // would return the same mask
791  for (unsigned int c = 1; c < this->size(); ++c)
792  Assert(mask == (*this)[c].block_mask(vector),
793  ExcMessage("Not all elements of this collection agree on what "
794  "the appropriate mask should be."));
795 
796  return mask;
797  }
798 
799 
800  template <int dim, int spacedim>
801  BlockMask
803  const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
804  {
805  Assert(this->size() > 0,
806  ExcMessage("This collection contains no finite element."));
807 
808  // get the mask from the first element of the collection
809  const BlockMask mask = (*this)[0].block_mask(sym_tensor);
810 
811  // but then also verify that the other elements of the collection
812  // would return the same mask
813  for (unsigned int c = 1; c < this->size(); ++c)
814  Assert(mask == (*this)[c].block_mask(sym_tensor),
815  ExcMessage("Not all elements of this collection agree on what "
816  "the appropriate mask should be."));
817 
818  return mask;
819  }
820 
821 
822 
823  template <int dim, int spacedim>
824  BlockMask
826  const ComponentMask &component_mask) const
827  {
828  Assert(this->size() > 0,
829  ExcMessage("This collection contains no finite element."));
830 
831  // get the mask from the first element of the collection
832  const BlockMask mask = (*this)[0].block_mask(component_mask);
833 
834  // but then also verify that the other elements of the collection
835  // would return the same mask
836  for (unsigned int c = 1; c < this->size(); ++c)
837  Assert(mask == (*this)[c].block_mask(component_mask),
838  ExcMessage("Not all elements of this collection agree on what "
839  "the appropriate mask should be."));
840 
841  return mask;
842  }
843 
844 
845 
846  template <int dim, int spacedim>
847  unsigned int
849  {
850  Assert(this->size() > 0, ExcNoFiniteElements());
851 
852  const unsigned int nb = this->operator[](0).n_blocks();
853  for (unsigned int i = 1; i < this->size(); ++i)
854  Assert(this->operator[](i).n_blocks() == nb,
855  ExcMessage("Not all finite elements in this collection have "
856  "the same number of components."));
857 
858  return nb;
859  }
860 } // namespace hp
861 
862 
863 
864 // explicit instantiations
865 #include "fe_collection.inst"
866 
867 
unsigned int n_components() const
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const =0
CollectionIterator< T > begin() const
Definition: collection.h:283
std::vector< std::map< unsigned int, unsigned int > > hp_vertex_dof_identities(const std::set< unsigned int > &fes) const
void push_back(const FiniteElement< dim, spacedim > &new_fe)
std::shared_ptr< MappingCollection< dim, spacedim > > reference_cell_default_linear_mapping
std::vector< std::map< unsigned int, unsigned int > > hp_line_dof_identities(const std::set< unsigned int > &fes) const
std::vector< std::map< unsigned int, unsigned int > > hp_quad_dof_identities(const std::set< unsigned int > &fes, const unsigned int face_no=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcEmptyObject()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
static ::ExceptionBase & ExcMessage(std::string arg1)
adjacency_list< vecS, vecS, undirectedS, property< vertex_color_t, default_color_type, property< vertex_degree_t, int > >> Graph
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Definition: hp.h:118
const types::fe_index invalid_fe_index
Definition: types.h:228
unsigned short int fe_index
Definition: types.h:60