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