Reference documentation for deal.II version Git 54f37b890c 2020-02-28 19:52:28 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_system.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/memory_consumption.h>
17 #include <deal.II/base/qprojector.h>
18 #include <deal.II/base/quadrature.h>
19 #include <deal.II/base/std_cxx14/memory.h>
20 
21 #include <deal.II/dofs/dof_accessor.h>
22 
23 #include <deal.II/fe/fe_system.h>
24 #include <deal.II/fe/fe_tools.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
28 #include <deal.II/grid/tria.h>
29 #include <deal.II/grid/tria_iterator.h>
30 
31 #include <sstream>
32 
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 namespace
37 {
38  unsigned int
39  count_nonzeros(const std::vector<unsigned int> &vec)
40  {
41  return std::count_if(vec.begin(), vec.end(), [](const unsigned int i) {
42  return i > 0;
43  });
44  }
45 } // namespace
46 /* ----------------------- FESystem::InternalData ------------------- */
47 
48 
49 template <int dim, int spacedim>
51  const unsigned int n_base_elements)
52  : base_fe_datas(n_base_elements)
53  , base_fe_output_objects(n_base_elements)
54 {}
55 
56 
57 
58 template <int dim, int spacedim>
60 {
61  // delete pointers and set them to zero to avoid inadvertent use
62  for (unsigned int i = 0; i < base_fe_datas.size(); ++i)
63  base_fe_datas[i].reset();
64 }
65 
66 
67 template <int dim, int spacedim>
70  const unsigned int base_no) const
71 {
72  AssertIndexRange(base_no, base_fe_datas.size());
73  return *base_fe_datas[base_no];
74 }
75 
76 
77 
78 template <int dim, int spacedim>
79 void
81  const unsigned int base_no,
82  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> ptr)
83 {
84  AssertIndexRange(base_no, base_fe_datas.size());
85  base_fe_datas[base_no] = std::move(ptr);
86 }
87 
88 
89 
90 template <int dim, int spacedim>
93  const unsigned int base_no) const
94 {
96  return base_fe_output_objects[base_no];
97 }
98 
99 
100 
101 /* ---------------------------------- FESystem ------------------- */
102 
103 
104 template <int dim, int spacedim>
106 
107 
108 template <int dim, int spacedim>
110  const unsigned int n_elements)
111  : FiniteElement<dim, spacedim>(
112  FETools::Compositing::multiply_dof_numbers(&fe, n_elements),
113  FETools::Compositing::compute_restriction_is_additive_flags(&fe,
114  n_elements),
115  FETools::Compositing::compute_nonzero_components(&fe, n_elements))
116  , base_elements((n_elements > 0))
117 {
118  std::vector<const FiniteElement<dim, spacedim> *> fes;
119  fes.push_back(&fe);
120  std::vector<unsigned int> multiplicities;
121  multiplicities.push_back(n_elements);
122  initialize(fes, multiplicities);
123 }
124 
125 
126 
127 template <int dim, int spacedim>
129  const unsigned int n1,
130  const FiniteElement<dim, spacedim> &fe2,
131  const unsigned int n2)
132  : FiniteElement<dim, spacedim>(
133  FETools::Compositing::multiply_dof_numbers(&fe1, n1, &fe2, n2),
134  FETools::Compositing::compute_restriction_is_additive_flags(&fe1,
135  n1,
136  &fe2,
137  n2),
138  FETools::Compositing::compute_nonzero_components(&fe1, n1, &fe2, n2))
139  , base_elements((n1 > 0) + (n2 > 0))
140 {
141  std::vector<const FiniteElement<dim, spacedim> *> fes;
142  fes.push_back(&fe1);
143  fes.push_back(&fe2);
144  std::vector<unsigned int> multiplicities;
145  multiplicities.push_back(n1);
146  multiplicities.push_back(n2);
147  initialize(fes, multiplicities);
148 }
149 
150 
151 
152 template <int dim, int spacedim>
154  const unsigned int n1,
155  const FiniteElement<dim, spacedim> &fe2,
156  const unsigned int n2,
157  const FiniteElement<dim, spacedim> &fe3,
158  const unsigned int n3)
159  : FiniteElement<dim, spacedim>(
160  FETools::Compositing::multiply_dof_numbers(&fe1, n1, &fe2, n2, &fe3, n3),
161  FETools::Compositing::compute_restriction_is_additive_flags(&fe1,
162  n1,
163  &fe2,
164  n2,
165  &fe3,
166  n3),
167  FETools::Compositing::compute_nonzero_components(&fe1,
168  n1,
169  &fe2,
170  n2,
171  &fe3,
172  n3))
173  , base_elements((n1 > 0) + (n2 > 0) + (n3 > 0))
174 {
175  std::vector<const FiniteElement<dim, spacedim> *> fes;
176  fes.push_back(&fe1);
177  fes.push_back(&fe2);
178  fes.push_back(&fe3);
179  std::vector<unsigned int> multiplicities;
180  multiplicities.push_back(n1);
181  multiplicities.push_back(n2);
182  multiplicities.push_back(n3);
183  initialize(fes, multiplicities);
184 }
185 
186 
187 
188 template <int dim, int spacedim>
190  const unsigned int n1,
191  const FiniteElement<dim, spacedim> &fe2,
192  const unsigned int n2,
193  const FiniteElement<dim, spacedim> &fe3,
194  const unsigned int n3,
195  const FiniteElement<dim, spacedim> &fe4,
196  const unsigned int n4)
197  : FiniteElement<dim, spacedim>(
198  FETools::Compositing::multiply_dof_numbers(&fe1,
199  n1,
200  &fe2,
201  n2,
202  &fe3,
203  n3,
204  &fe4,
205  n4),
206  FETools::Compositing::compute_restriction_is_additive_flags(&fe1,
207  n1,
208  &fe2,
209  n2,
210  &fe3,
211  n3,
212  &fe4,
213  n4),
214  FETools::Compositing::compute_nonzero_components(&fe1,
215  n1,
216  &fe2,
217  n2,
218  &fe3,
219  n3,
220  &fe4,
221  n4))
222  , base_elements((n1 > 0) + (n2 > 0) + (n3 > 0) + (n4 > 0))
223 {
224  std::vector<const FiniteElement<dim, spacedim> *> fes;
225  fes.push_back(&fe1);
226  fes.push_back(&fe2);
227  fes.push_back(&fe3);
228  fes.push_back(&fe4);
229  std::vector<unsigned int> multiplicities;
230  multiplicities.push_back(n1);
231  multiplicities.push_back(n2);
232  multiplicities.push_back(n3);
233  multiplicities.push_back(n4);
234  initialize(fes, multiplicities);
235 }
236 
237 
238 
239 template <int dim, int spacedim>
241  const unsigned int n1,
242  const FiniteElement<dim, spacedim> &fe2,
243  const unsigned int n2,
244  const FiniteElement<dim, spacedim> &fe3,
245  const unsigned int n3,
246  const FiniteElement<dim, spacedim> &fe4,
247  const unsigned int n4,
248  const FiniteElement<dim, spacedim> &fe5,
249  const unsigned int n5)
250  : FiniteElement<dim, spacedim>(
251  FETools::Compositing::
252  multiply_dof_numbers(&fe1, n1, &fe2, n2, &fe3, n3, &fe4, n4, &fe5, n5),
253  FETools::Compositing::compute_restriction_is_additive_flags(&fe1,
254  n1,
255  &fe2,
256  n2,
257  &fe3,
258  n3,
259  &fe4,
260  n4,
261  &fe5,
262  n5),
263  FETools::Compositing::compute_nonzero_components(&fe1,
264  n1,
265  &fe2,
266  n2,
267  &fe3,
268  n3,
269  &fe4,
270  n4,
271  &fe5,
272  n5))
273  , base_elements((n1 > 0) + (n2 > 0) + (n3 > 0) + (n4 > 0) + (n5 > 0))
274 {
275  std::vector<const FiniteElement<dim, spacedim> *> fes;
276  fes.push_back(&fe1);
277  fes.push_back(&fe2);
278  fes.push_back(&fe3);
279  fes.push_back(&fe4);
280  fes.push_back(&fe5);
281  std::vector<unsigned int> multiplicities;
282  multiplicities.push_back(n1);
283  multiplicities.push_back(n2);
284  multiplicities.push_back(n3);
285  multiplicities.push_back(n4);
286  multiplicities.push_back(n5);
287  initialize(fes, multiplicities);
288 }
289 
290 
291 
292 template <int dim, int spacedim>
294  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
295  const std::vector<unsigned int> & multiplicities)
296  : FiniteElement<dim, spacedim>(
297  FETools::Compositing::multiply_dof_numbers(fes, multiplicities),
298  FETools::Compositing::compute_restriction_is_additive_flags(
299  fes,
300  multiplicities),
301  FETools::Compositing::compute_nonzero_components(fes, multiplicities))
302  , base_elements(count_nonzeros(multiplicities))
303 {
304  initialize(fes, multiplicities);
305 }
306 
307 
308 
309 template <int dim, int spacedim>
310 std::string
312 {
313  // note that the
314  // FETools::get_fe_by_name
315  // function depends on the
316  // particular format of the string
317  // this function returns, so they
318  // have to be kept in synch
319 
320  std::ostringstream namebuf;
321 
322  namebuf << "FESystem<" << Utilities::dim_string(dim, spacedim) << ">[";
323  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
324  {
325  namebuf << base_element(i).get_name();
326  if (this->element_multiplicity(i) != 1)
327  namebuf << '^' << this->element_multiplicity(i);
328  if (i != this->n_base_elements() - 1)
329  namebuf << '-';
330  }
331  namebuf << ']';
332 
333  return namebuf.str();
334 }
335 
336 
337 
338 template <int dim, int spacedim>
339 std::unique_ptr<FiniteElement<dim, spacedim>>
341 {
342  std::vector<const FiniteElement<dim, spacedim> *> fes;
343  std::vector<unsigned int> multiplicities;
344 
345  for (unsigned int i = 0; i < this->n_base_elements(); i++)
346  {
347  fes.push_back(&base_element(i));
348  multiplicities.push_back(this->element_multiplicity(i));
349  }
350  return std_cxx14::make_unique<FESystem<dim, spacedim>>(fes, multiplicities);
351 }
352 
353 
354 
355 template <int dim, int spacedim>
358  const unsigned int first_component,
359  const unsigned int n_selected_components) const
360 {
361  Assert(first_component + n_selected_components <= this->n_components(),
362  ExcMessage("Invalid arguments (not a part of this FiniteElement)."));
363 
364  const unsigned int base_index =
365  this->component_to_base_table[first_component].first.first;
366  const unsigned int component_in_base =
367  this->component_to_base_table[first_component].first.second;
368  const unsigned int base_components =
369  this->base_element(base_index).n_components();
370 
371  // Only select our child base_index if that is all the user wanted. Error
372  // handling will be done inside the recursion.
373  if (n_selected_components <= base_components)
374  return this->base_element(base_index)
375  .get_sub_fe(component_in_base, n_selected_components);
376 
377  Assert(n_selected_components == this->n_components(),
378  ExcMessage("You can not select a part of a FiniteElement."));
379  return *this;
380 }
381 
382 
383 
384 template <int dim, int spacedim>
385 double
387  const Point<dim> & p) const
388 {
390  Assert(this->is_primitive(i),
392  i)));
393 
394  return (base_element(this->system_to_base_table[i].first.first)
395  .shape_value(this->system_to_base_table[i].second, p));
396 }
397 
398 
399 
400 template <int dim, int spacedim>
401 double
403  const unsigned int i,
404  const Point<dim> & p,
405  const unsigned int component) const
406 {
408  AssertIndexRange(component, this->n_components());
409 
410  // if this value is supposed to be
411  // zero, then return right away...
412  if (this->nonzero_components[i][component] == false)
413  return 0;
414 
415  // ...otherwise: first find out to
416  // which of the base elements this
417  // desired component belongs, and
418  // which component within this base
419  // element it is
420  const unsigned int base = this->component_to_base_index(component).first;
421  const unsigned int component_in_base =
422  this->component_to_base_index(component).second;
423 
424  // then get value from base
425  // element. note that that will
426  // throw an error should the
427  // respective shape function not be
428  // primitive; thus, there is no
429  // need to check this here
430  return (base_element(base).shape_value_component(
431  this->system_to_base_table[i].second, p, component_in_base));
432 }
433 
434 
435 
436 template <int dim, int spacedim>
439  const Point<dim> & p) const
440 {
442  Assert(this->is_primitive(i),
444  i)));
445 
446  return (base_element(this->system_to_base_table[i].first.first)
447  .shape_grad(this->system_to_base_table[i].second, p));
448 }
449 
450 
451 
452 template <int dim, int spacedim>
455  const unsigned int i,
456  const Point<dim> & p,
457  const unsigned int component) const
458 {
460  AssertIndexRange(component, this->n_components());
461 
462  // if this value is supposed to be zero, then return right away...
463  if (this->nonzero_components[i][component] == false)
464  return Tensor<1, dim>();
465 
466  // ...otherwise: first find out to which of the base elements this desired
467  // component belongs, and which component within this base element it is
468  const unsigned int base = this->component_to_base_index(component).first;
469  const unsigned int component_in_base =
470  this->component_to_base_index(component).second;
471 
472  // then get value from base element. note that that will throw an error
473  // should the respective shape function not be primitive; thus, there is no
474  // need to check this here
475  return (base_element(base).shape_grad_component(
476  this->system_to_base_table[i].second, p, component_in_base));
477 }
478 
479 
480 
481 template <int dim, int spacedim>
484  const Point<dim> & p) const
485 {
487  Assert(this->is_primitive(i),
489  i)));
490 
491  return (base_element(this->system_to_base_table[i].first.first)
492  .shape_grad_grad(this->system_to_base_table[i].second, p));
493 }
494 
495 
496 
497 template <int dim, int spacedim>
500  const unsigned int i,
501  const Point<dim> & p,
502  const unsigned int component) const
503 {
505  AssertIndexRange(component, this->n_components());
506 
507  // if this value is supposed to be zero, then return right away...
508  if (this->nonzero_components[i][component] == false)
509  return Tensor<2, dim>();
510 
511  // ...otherwise: first find out to which of the base elements this desired
512  // component belongs, and which component within this base element it is
513  const unsigned int base = this->component_to_base_index(component).first;
514  const unsigned int component_in_base =
515  this->component_to_base_index(component).second;
516 
517  // then get value from base element. note that that will throw an error
518  // should the respective shape function not be primitive; thus, there is no
519  // need to check this here
521  this->system_to_base_table[i].second, p, component_in_base));
522 }
523 
524 
525 
526 template <int dim, int spacedim>
529  const Point<dim> & p) const
530 {
532  Assert(this->is_primitive(i),
534  i)));
535 
536  return (base_element(this->system_to_base_table[i].first.first)
537  .shape_3rd_derivative(this->system_to_base_table[i].second, p));
538 }
539 
540 
541 
542 template <int dim, int spacedim>
545  const unsigned int i,
546  const Point<dim> & p,
547  const unsigned int component) const
548 {
550  AssertIndexRange(component, this->n_components());
551 
552  // if this value is supposed to be zero, then return right away...
553  if (this->nonzero_components[i][component] == false)
554  return Tensor<3, dim>();
555 
556  // ...otherwise: first find out to which of the base elements this desired
557  // component belongs, and which component within this base element it is
558  const unsigned int base = this->component_to_base_index(component).first;
559  const unsigned int component_in_base =
560  this->component_to_base_index(component).second;
561 
562  // then get value from base element. note that that will throw an error
563  // should the respective shape function not be primitive; thus, there is no
564  // need to check this here
566  this->system_to_base_table[i].second, p, component_in_base));
567 }
568 
569 
570 
571 template <int dim, int spacedim>
574  const Point<dim> & p) const
575 {
577  Assert(this->is_primitive(i),
579  i)));
580 
581  return (base_element(this->system_to_base_table[i].first.first)
582  .shape_4th_derivative(this->system_to_base_table[i].second, p));
583 }
584 
585 
586 
587 template <int dim, int spacedim>
590  const unsigned int i,
591  const Point<dim> & p,
592  const unsigned int component) const
593 {
595  AssertIndexRange(component, this->n_components());
596 
597  // if this value is supposed to be zero, then return right away...
598  if (this->nonzero_components[i][component] == false)
599  return Tensor<4, dim>();
600 
601  // ...otherwise: first find out to which of the base elements this desired
602  // component belongs, and which component within this base element it is
603  const unsigned int base = this->component_to_base_index(component).first;
604  const unsigned int component_in_base =
605  this->component_to_base_index(component).second;
606 
607  // then get value from base element. note that that will throw an error
608  // should the respective shape function not be primitive; thus, there is no
609  // need to check this here
611  this->system_to_base_table[i].second, p, component_in_base));
612 }
613 
614 
615 
616 template <int dim, int spacedim>
617 void
619  const FiniteElement<dim, spacedim> &x_source_fe,
620  FullMatrix<double> & interpolation_matrix) const
621 {
622  // check that the size of the matrices is correct. for historical
623  // reasons, if you call matrix.reinit(8,0), it sets the sizes
624  // to m==n==0 internally. this may happen when we use a FE_Nothing,
625  // so write the test in a more lenient way
626  Assert((interpolation_matrix.m() == this->dofs_per_cell) ||
627  (x_source_fe.dofs_per_cell == 0),
628  ExcDimensionMismatch(interpolation_matrix.m(), this->dofs_per_cell));
629  Assert((interpolation_matrix.n() == x_source_fe.dofs_per_cell) ||
630  (this->dofs_per_cell == 0),
631  ExcDimensionMismatch(interpolation_matrix.m(),
632  x_source_fe.dofs_per_cell));
633 
634  // there are certain conditions that the two elements have to satisfy so
635  // that this can work.
636  //
637  // condition 1: the other element must also be a system element
638 
639  AssertThrow(
640  (x_source_fe.get_name().find("FESystem<") == 0) ||
641  (dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) != nullptr),
643 
644  // ok, source is a system element, so we may be able to do the work
645  const FESystem<dim, spacedim> &source_fe =
646  dynamic_cast<const FESystem<dim, spacedim> &>(x_source_fe);
647 
648  // condition 2: same number of basis elements
649  AssertThrow(
650  this->n_base_elements() == source_fe.n_base_elements(),
652 
653  // condition 3: same number of basis elements
654  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
655  AssertThrow(
656  this->element_multiplicity(i) == source_fe.element_multiplicity(i),
657  (typename FiniteElement<dim,
659 
660  // ok, so let's try whether it works:
661 
662  // first let's see whether all the basis elements actually generate their
663  // interpolation matrices. if we get past the following loop, then
664  // apparently none of the called base elements threw an exception, so we're
665  // fine continuing and assembling the one big matrix from the small ones of
666  // the base elements
667  std::vector<FullMatrix<double>> base_matrices(this->n_base_elements());
668  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
669  {
670  base_matrices[i].reinit(base_element(i).dofs_per_cell,
671  source_fe.base_element(i).dofs_per_cell);
672  base_element(i).get_interpolation_matrix(source_fe.base_element(i),
673  base_matrices[i]);
674  }
675 
676  // first clear big matrix, to make sure that entries that would couple
677  // different bases (or multiplicity indices) are really zero. then assign
678  // entries
679  interpolation_matrix = 0;
680  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
681  for (unsigned int j = 0; j < source_fe.dofs_per_cell; ++j)
682  if (this->system_to_base_table[i].first ==
683  source_fe.system_to_base_table[j].first)
684  interpolation_matrix(i, j) =
685  (base_matrices[this->system_to_base_table[i].first.first](
686  this->system_to_base_table[i].second,
687  source_fe.system_to_base_table[j].second));
688 }
689 
690 
691 
692 template <int dim, int spacedim>
693 const FullMatrix<double> &
695  const unsigned int child,
696  const RefinementCase<dim> &refinement_case) const
697 {
698  AssertIndexRange(refinement_case,
700  Assert(refinement_case != RefinementCase<dim>::no_refinement,
701  ExcMessage(
702  "Restriction matrices are only available for refined cells!"));
703  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
704 
705  // initialization upon first request
706  if (this->restriction[refinement_case - 1][child].n() == 0)
707  {
708  std::lock_guard<std::mutex> lock(this->mutex);
709 
710  // check if updated while waiting for lock
711  if (this->restriction[refinement_case - 1][child].n() ==
712  this->dofs_per_cell)
713  return this->restriction[refinement_case - 1][child];
714 
715  // Check if some of the matrices of the base elements are void.
716  bool do_restriction = true;
717 
718  // shortcut for accessing local restrictions further down
719  std::vector<const FullMatrix<double> *> base_matrices(
720  this->n_base_elements());
721 
722  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
723  {
724  base_matrices[i] =
725  &base_element(i).get_restriction_matrix(child, refinement_case);
726  if (base_matrices[i]->n() != base_element(i).dofs_per_cell)
727  do_restriction = false;
728  }
729  Assert(do_restriction,
731 
732  // if we did not encounter void matrices, initialize the matrix sizes
733  if (do_restriction)
734  {
736  this->dofs_per_cell);
737 
738  // distribute the matrices of the base finite elements to the
739  // matrices of this object. for this, loop over all degrees of
740  // freedom and take the respective entry of the underlying base
741  // element.
742  //
743  // note that we by definition of a base element, they are
744  // independent, i.e. do not couple. only DoFs that belong to the
745  // same instance of a base element may couple
746  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
747  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
748  {
749  // first find out to which base element indices i and j
750  // belong, and which instance thereof in case the base element
751  // has a multiplicity greater than one. if they should not
752  // happen to belong to the same instance of a base element,
753  // then they cannot couple, so go on with the next index
754  if (this->system_to_base_table[i].first !=
755  this->system_to_base_table[j].first)
756  continue;
757 
758  // so get the common base element and the indices therein:
759  const unsigned int base =
760  this->system_to_base_table[i].first.first;
761 
762  const unsigned int base_index_i =
763  this->system_to_base_table[i].second,
764  base_index_j =
765  this->system_to_base_table[j].second;
766 
767  // if we are sure that DoFs i and j may couple, then copy
768  // entries of the matrices:
769  restriction(i, j) =
770  (*base_matrices[base])(base_index_i, base_index_j);
771  }
772 
773  restriction.swap(const_cast<FullMatrix<double> &>(
774  this->restriction[refinement_case - 1][child]));
775  }
776  }
777 
778  return this->restriction[refinement_case - 1][child];
779 }
780 
781 
782 
783 template <int dim, int spacedim>
784 const FullMatrix<double> &
786  const unsigned int child,
787  const RefinementCase<dim> &refinement_case) const
788 {
789  AssertIndexRange(refinement_case,
791  Assert(refinement_case != RefinementCase<dim>::no_refinement,
792  ExcMessage(
793  "Restriction matrices are only available for refined cells!"));
794  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
795 
796  // initialization upon first request, construction completely analogous to
797  // restriction matrix
798  if (this->prolongation[refinement_case - 1][child].n() == 0)
799  {
800  std::lock_guard<std::mutex> lock(this->mutex);
801 
802  if (this->prolongation[refinement_case - 1][child].n() ==
803  this->dofs_per_cell)
804  return this->prolongation[refinement_case - 1][child];
805 
806  bool do_prolongation = true;
807  std::vector<const FullMatrix<double> *> base_matrices(
808  this->n_base_elements());
809  for (unsigned int i = 0; i < this->n_base_elements(); ++i)
810  {
811  base_matrices[i] =
812  &base_element(i).get_prolongation_matrix(child, refinement_case);
813  if (base_matrices[i]->n() != base_element(i).dofs_per_cell)
814  do_prolongation = false;
815  }
816  Assert(do_prolongation,
818 
819  if (do_prolongation)
820  {
821  FullMatrix<double> prolongate(this->dofs_per_cell,
822  this->dofs_per_cell);
823 
824  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
825  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
826  {
827  if (this->system_to_base_table[i].first !=
828  this->system_to_base_table[j].first)
829  continue;
830  const unsigned int base =
831  this->system_to_base_table[i].first.first;
832 
833  const unsigned int base_index_i =
834  this->system_to_base_table[i].second,
835  base_index_j =
836  this->system_to_base_table[j].second;
837  prolongate(i, j) =
838  (*base_matrices[base])(base_index_i, base_index_j);
839  }
840  prolongate.swap(const_cast<FullMatrix<double> &>(
841  this->prolongation[refinement_case - 1][child]));
842  }
843  }
844 
845  return this->prolongation[refinement_case - 1][child];
846 }
847 
848 
849 template <int dim, int spacedim>
850 unsigned int
851 FESystem<dim, spacedim>::face_to_cell_index(const unsigned int face_dof_index,
852  const unsigned int face,
853  const bool face_orientation,
854  const bool face_flip,
855  const bool face_rotation) const
856 {
857  // we need to ask the base elements how they want to translate
858  // the DoFs within their own numbering. thus, translate to
859  // the base element numbering and then back
860  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
861  face_base_index = this->face_system_to_base_index(face_dof_index);
862 
863  const unsigned int base_face_to_cell_index =
864  this->base_element(face_base_index.first.first)
865  .face_to_cell_index(face_base_index.second,
866  face,
867  face_orientation,
868  face_flip,
869  face_rotation);
870 
871  // it would be nice if we had a base_to_system_index function, but
872  // all that exists is a component_to_system_index function. we can't do
873  // this here because it won't work for non-primitive elements. consequently,
874  // simply do a loop over all dofs till we find whether it corresponds
875  // to the one we're interested in -- crude, maybe, but works for now
876  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int> target =
877  std::make_pair(face_base_index.first, base_face_to_cell_index);
878  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
879  if (this->system_to_base_index(i) == target)
880  return i;
881 
882  Assert(false, ExcInternalError());
884 }
885 
886 
887 
888 //---------------------------------------------------------------------------
889 // Data field initialization
890 //---------------------------------------------------------------------------
891 
892 
893 
894 template <int dim, int spacedim>
897 {
899  // generate maximal set of flags
900  // that are necessary
901  for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
902  out |= base_element(base_no).requires_update_flags(flags);
903  return out;
904 }
905 
906 
907 
908 template <int dim, int spacedim>
909 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
911  const UpdateFlags flags,
912  const Mapping<dim, spacedim> &mapping,
913  const Quadrature<dim> & quadrature,
915  spacedim>
916  & /*output_data*/) const
917 {
918  // create an internal data object and set the update flags we will need
919  // to deal with. the current object does not make use of these flags,
920  // but we need to nevertheless set them correctly since we look
921  // into the update_each flag of base elements in fill_fe_values,
922  // and so the current object's update_each flag needs to be
923  // correct in case the current FESystem is a base element for another,
924  // higher-level FESystem itself.
925  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
926  data_ptr = std_cxx14::make_unique<InternalData>(this->n_base_elements());
927  auto &data = dynamic_cast<InternalData &>(*data_ptr);
928  data.update_each = requires_update_flags(flags);
929 
930  // get data objects from each of the base elements and store
931  // them. one might think that doing this in parallel (over the
932  // base elements) would be a good idea, but this turns out to
933  // be wrong because we would then run these jobs on different
934  // threads/processors and this allocates memory in different
935  // NUMA domains; this has large detrimental effects when later
936  // writing into these objects in fill_fe_*_values. all of this
937  // is particularly true when using FEValues objects in
938  // WorkStream contexts where we explicitly make sure that
939  // every function only uses objects previously allocated
940  // in the same NUMA context and on the same thread as the
941  // function is called
942  for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
943  {
945  &base_fe_output_object = data.get_fe_output_object(base_no);
946  base_fe_output_object.initialize(
947  quadrature.size(),
948  base_element(base_no),
949  flags | base_element(base_no).requires_update_flags(flags));
950 
951  // let base objects produce their scratch objects. they may
952  // also at this time write into the output objects we provide
953  // for them; it would be nice if we could already copy something
954  // out of the base output object into the system output object,
955  // but we can't because we can't know what the elements already
956  // copied and/or will want to update on every cell
957  auto base_fe_data = base_element(base_no).get_data(flags,
958  mapping,
959  quadrature,
960  base_fe_output_object);
961 
962  data.set_fe_data(base_no, std::move(base_fe_data));
963  }
964 
965  return data_ptr;
966 }
967 
968 // The following function is a clone of get_data, with the exception
969 // that get_face_data of the base elements is called.
970 
971 template <int dim, int spacedim>
972 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
974  const UpdateFlags flags,
975  const Mapping<dim, spacedim> &mapping,
976  const Quadrature<dim - 1> & quadrature,
978  spacedim>
979  & /*output_data*/) const
980 {
981  // create an internal data object and set the update flags we will need
982  // to deal with. the current object does not make use of these flags,
983  // but we need to nevertheless set them correctly since we look
984  // into the update_each flag of base elements in fill_fe_values,
985  // and so the current object's update_each flag needs to be
986  // correct in case the current FESystem is a base element for another,
987  // higher-level FESystem itself.
988  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
989  data_ptr = std_cxx14::make_unique<InternalData>(this->n_base_elements());
990  auto &data = dynamic_cast<InternalData &>(*data_ptr);
991  data.update_each = requires_update_flags(flags);
992 
993  // get data objects from each of the base elements and store
994  // them. one might think that doing this in parallel (over the
995  // base elements) would be a good idea, but this turns out to
996  // be wrong because we would then run these jobs on different
997  // threads/processors and this allocates memory in different
998  // NUMA domains; this has large detrimental effects when later
999  // writing into these objects in fill_fe_*_values. all of this
1000  // is particularly true when using FEValues objects in
1001  // WorkStream contexts where we explicitly make sure that
1002  // every function only uses objects previously allocated
1003  // in the same NUMA context and on the same thread as the
1004  // function is called
1005  for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1006  {
1008  &base_fe_output_object = data.get_fe_output_object(base_no);
1009  base_fe_output_object.initialize(
1010  quadrature.size(),
1011  base_element(base_no),
1012  flags | base_element(base_no).requires_update_flags(flags));
1013 
1014  // let base objects produce their scratch objects. they may
1015  // also at this time write into the output objects we provide
1016  // for them; it would be nice if we could already copy something
1017  // out of the base output object into the system output object,
1018  // but we can't because we can't know what the elements already
1019  // copied and/or will want to update on every cell
1020  auto base_fe_data = base_element(base_no).get_face_data(
1021  flags, mapping, quadrature, base_fe_output_object);
1022 
1023  data.set_fe_data(base_no, std::move(base_fe_data));
1024  }
1025 
1026  return data_ptr;
1027 }
1028 
1029 
1030 
1031 // The following function is a clone of get_data, with the exception
1032 // that get_subface_data of the base elements is called.
1033 
1034 template <int dim, int spacedim>
1035 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1037  const UpdateFlags flags,
1038  const Mapping<dim, spacedim> &mapping,
1039  const Quadrature<dim - 1> & quadrature,
1041  spacedim>
1042  & /*output_data*/) const
1043 {
1044  // create an internal data object and set the update flags we will need
1045  // to deal with. the current object does not make use of these flags,
1046  // but we need to nevertheless set them correctly since we look
1047  // into the update_each flag of base elements in fill_fe_values,
1048  // and so the current object's update_each flag needs to be
1049  // correct in case the current FESystem is a base element for another,
1050  // higher-level FESystem itself.
1051  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1052  data_ptr = std_cxx14::make_unique<InternalData>(this->n_base_elements());
1053  auto &data = dynamic_cast<InternalData &>(*data_ptr);
1054 
1055  data.update_each = requires_update_flags(flags);
1056 
1057  // get data objects from each of the base elements and store
1058  // them. one might think that doing this in parallel (over the
1059  // base elements) would be a good idea, but this turns out to
1060  // be wrong because we would then run these jobs on different
1061  // threads/processors and this allocates memory in different
1062  // NUMA domains; this has large detrimental effects when later
1063  // writing into these objects in fill_fe_*_values. all of this
1064  // is particularly true when using FEValues objects in
1065  // WorkStream contexts where we explicitly make sure that
1066  // every function only uses objects previously allocated
1067  // in the same NUMA context and on the same thread as the
1068  // function is called
1069  for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1070  {
1072  &base_fe_output_object = data.get_fe_output_object(base_no);
1073  base_fe_output_object.initialize(
1074  quadrature.size(),
1075  base_element(base_no),
1076  flags | base_element(base_no).requires_update_flags(flags));
1077 
1078  // let base objects produce their scratch objects. they may
1079  // also at this time write into the output objects we provide
1080  // for them; it would be nice if we could already copy something
1081  // out of the base output object into the system output object,
1082  // but we can't because we can't know what the elements already
1083  // copied and/or will want to update on every cell
1084  auto base_fe_data = base_element(base_no).get_subface_data(
1085  flags, mapping, quadrature, base_fe_output_object);
1086 
1087  data.set_fe_data(base_no, std::move(base_fe_data));
1088  }
1089 
1090  return data_ptr;
1091 }
1092 
1093 
1094 
1095 template <int dim, int spacedim>
1096 void
1098  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1099  const CellSimilarity::Similarity cell_similarity,
1100  const Quadrature<dim> & quadrature,
1101  const Mapping<dim, spacedim> & mapping,
1102  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1103  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
1104  spacedim>
1105  & mapping_data,
1106  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1108  spacedim>
1109  &output_data) const
1110 {
1111  compute_fill(mapping,
1112  cell,
1115  quadrature,
1116  cell_similarity,
1117  mapping_internal,
1118  fe_internal,
1119  mapping_data,
1120  output_data);
1121 }
1122 
1123 
1124 
1125 template <int dim, int spacedim>
1126 void
1128  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1129  const unsigned int face_no,
1130  const Quadrature<dim - 1> & quadrature,
1131  const Mapping<dim, spacedim> & mapping,
1132  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1133  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
1134  spacedim>
1135  & mapping_data,
1136  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1138  spacedim>
1139  &output_data) const
1140 {
1141  compute_fill(mapping,
1142  cell,
1143  face_no,
1145  quadrature,
1147  mapping_internal,
1148  fe_internal,
1149  mapping_data,
1150  output_data);
1151 }
1152 
1153 
1154 
1155 template <int dim, int spacedim>
1156 void
1158  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1159  const unsigned int face_no,
1160  const unsigned int sub_no,
1161  const Quadrature<dim - 1> & quadrature,
1162  const Mapping<dim, spacedim> & mapping,
1163  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1164  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
1165  spacedim>
1166  & mapping_data,
1167  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1169  spacedim>
1170  &output_data) const
1171 {
1172  compute_fill(mapping,
1173  cell,
1174  face_no,
1175  sub_no,
1176  quadrature,
1178  mapping_internal,
1179  fe_internal,
1180  mapping_data,
1181  output_data);
1182 }
1183 
1184 
1185 
1186 template <int dim, int spacedim>
1187 template <int dim_1>
1188 void
1190  const Mapping<dim, spacedim> & mapping,
1191  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1192  const unsigned int face_no,
1193  const unsigned int sub_no,
1194  const Quadrature<dim_1> & quadrature,
1195  const CellSimilarity::Similarity cell_similarity,
1196  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1197  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1199  &mapping_data,
1201  &output_data) const
1202 {
1203  // convert data object to internal
1204  // data for this class. fails with
1205  // an exception if that is not
1206  // possible
1207  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1208  ExcInternalError());
1209  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1210 
1211  // Either dim_1==dim
1212  // (fill_fe_values) or dim_1==dim-1
1213  // (fill_fe_(sub)face_values)
1214  Assert(dim_1 == dim || dim_1 == dim - 1, ExcInternalError());
1215  const UpdateFlags flags = fe_data.update_each;
1216 
1217 
1218  // loop over the base elements, let them compute what they need to compute,
1219  // and then copy what is necessary.
1220  //
1221  // one may think that it would be a good idea to parallelize this over
1222  // base elements, but it turns out to be not worthwhile: doing so lets
1223  // multiple threads access data objects that were created by the current
1224  // thread, leading to many NUMA memory access inefficiencies. we specifically
1225  // want to avoid this if this class is called in a WorkStream context where
1226  // we very carefully allocate objects only on the thread where they
1227  // will actually be used; spawning new tasks here would be counterproductive
1230  for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1231  {
1232  const FiniteElement<dim, spacedim> &base_fe = base_element(base_no);
1233  typename FiniteElement<dim, spacedim>::InternalDataBase &base_fe_data =
1234  fe_data.get_fe_data(base_no);
1236  spacedim>
1237  &base_data = fe_data.get_fe_output_object(base_no);
1238 
1239  // fill_fe_face_values needs argument Quadrature<dim-1> for both cases
1240  // dim_1==dim-1 and dim_1=dim. Hence the following workaround
1241  const Quadrature<dim> * cell_quadrature = nullptr;
1242  const Quadrature<dim - 1> *face_quadrature = nullptr;
1243  const unsigned int n_q_points = quadrature.size();
1244 
1245  // static cast to the common base class of quadrature being either
1246  // Quadrature<dim> or Quadrature<dim-1>:
1247  const Subscriptor *quadrature_base_pointer = &quadrature;
1248 
1249  if (face_no == invalid_face_number)
1250  {
1251  Assert(dim_1 == dim, ExcDimensionMismatch(dim_1, dim));
1252  Assert(dynamic_cast<const Quadrature<dim> *>(
1253  quadrature_base_pointer) != nullptr,
1254  ExcInternalError());
1255 
1256  cell_quadrature =
1257  static_cast<const Quadrature<dim> *>(quadrature_base_pointer);
1258  }
1259  else
1260  {
1261  Assert(dim_1 == dim - 1, ExcDimensionMismatch(dim_1, dim - 1));
1262  Assert(dynamic_cast<const Quadrature<dim - 1> *>(
1263  quadrature_base_pointer) != nullptr,
1264  ExcInternalError());
1265 
1266  face_quadrature =
1267  static_cast<const Quadrature<dim - 1> *>(quadrature_base_pointer);
1268  }
1269 
1270 
1271  // Make sure that in the case of fill_fe_values the data is only
1272  // copied from base_data to data if base_data is changed. therefore
1273  // use fe_fe_data.current_update_flags()
1274  //
1275  // for the case of fill_fe_(sub)face_values the data needs to be
1276  // copied from base_data to data on each face, therefore use
1277  // base_fe_data.update_flags.
1278  if (face_no == invalid_face_number)
1279  base_fe.fill_fe_values(cell,
1280  cell_similarity,
1281  *cell_quadrature,
1282  mapping,
1283  mapping_internal,
1284  mapping_data,
1285  base_fe_data,
1286  base_data);
1287  else if (sub_no == invalid_face_number)
1288  base_fe.fill_fe_face_values(cell,
1289  face_no,
1290  *face_quadrature,
1291  mapping,
1292  mapping_internal,
1293  mapping_data,
1294  base_fe_data,
1295  base_data);
1296  else
1297  base_fe.fill_fe_subface_values(cell,
1298  face_no,
1299  sub_no,
1300  *face_quadrature,
1301  mapping,
1302  mapping_internal,
1303  mapping_data,
1304  base_fe_data,
1305  base_data);
1306 
1307  // now data has been generated, so copy it. we used to work by
1308  // looping over all base elements (i.e. this outer loop), then over
1309  // multiplicity, then over the shape functions from that base
1310  // element, but that requires that we can infer the global number of
1311  // a shape function from its number in the base element. for that we
1312  // had the component_to_system_table.
1313  //
1314  // however, this does of course no longer work since we have
1315  // non-primitive elements. so we go the other way round: loop over
1316  // all shape functions of the composed element, and here only treat
1317  // those shape functions that belong to a given base element
1318  // TODO: Introduce the needed table and loop only over base element
1319  // shape functions. This here is not efficient at all AND very bad style
1320  const UpdateFlags base_flags = base_fe_data.update_each;
1321 
1322  // some base element might involve values that depend on the shape
1323  // of the geometry, so we always need to copy the shape values around
1324  // also in case we detected a cell similarity (but no heavy work will
1325  // be done inside the individual elements in case we have a
1326  // translation and simple elements).
1327  for (unsigned int system_index = 0; system_index < this->dofs_per_cell;
1328  ++system_index)
1329  if (this->system_to_base_table[system_index].first.first == base_no)
1330  {
1331  const unsigned int base_index =
1332  this->system_to_base_table[system_index].second;
1333  Assert(base_index < base_fe.dofs_per_cell, ExcInternalError());
1334 
1335  // now copy. if the shape function is primitive, then there
1336  // is only one value to be copied, but for non-primitive
1337  // elements, there might be more values to be copied
1338  //
1339  // so, find out from which index to take this one value, and
1340  // to which index to put
1341  unsigned int out_index = 0;
1342  for (unsigned int i = 0; i < system_index; ++i)
1343  out_index += this->n_nonzero_components(i);
1344  unsigned int in_index = 0;
1345  for (unsigned int i = 0; i < base_index; ++i)
1346  in_index += base_fe.n_nonzero_components(i);
1347 
1348  // then loop over the number of components to be copied
1349  Assert(this->n_nonzero_components(system_index) ==
1350  base_fe.n_nonzero_components(base_index),
1351  ExcInternalError());
1352 
1353  if (base_flags & update_values)
1354  for (unsigned int s = 0;
1355  s < this->n_nonzero_components(system_index);
1356  ++s)
1357  for (unsigned int q = 0; q < n_q_points; ++q)
1358  output_data.shape_values[out_index + s][q] =
1359  base_data.shape_values(in_index + s, q);
1360 
1361  if (base_flags & update_gradients)
1362  for (unsigned int s = 0;
1363  s < this->n_nonzero_components(system_index);
1364  ++s)
1365  for (unsigned int q = 0; q < n_q_points; ++q)
1366  output_data.shape_gradients[out_index + s][q] =
1367  base_data.shape_gradients[in_index + s][q];
1368 
1369  if (base_flags & update_hessians)
1370  for (unsigned int s = 0;
1371  s < this->n_nonzero_components(system_index);
1372  ++s)
1373  for (unsigned int q = 0; q < n_q_points; ++q)
1374  output_data.shape_hessians[out_index + s][q] =
1375  base_data.shape_hessians[in_index + s][q];
1376 
1377  if (base_flags & update_3rd_derivatives)
1378  for (unsigned int s = 0;
1379  s < this->n_nonzero_components(system_index);
1380  ++s)
1381  for (unsigned int q = 0; q < n_q_points; ++q)
1382  output_data.shape_3rd_derivatives[out_index + s][q] =
1383  base_data.shape_3rd_derivatives[in_index + s][q];
1384  }
1385  }
1386 }
1387 
1388 
1389 template <int dim, int spacedim>
1390 void
1392 {
1393  // check whether all base elements implement their interface constraint
1394  // matrices. if this is not the case, then leave the interface costraints of
1395  // this composed element empty as well; however, the rest of the element is
1396  // usable
1397  for (unsigned int base = 0; base < this->n_base_elements(); ++base)
1398  if (base_element(base).constraints_are_implemented() == false)
1399  return;
1400 
1401  this->interface_constraints.TableBase<2, double>::reinit(
1402  this->interface_constraints_size());
1403 
1404  // the layout of the constraints matrix is described in the FiniteElement
1405  // class. you may want to look there first before trying to understand the
1406  // following, especially the mapping of the @p{m} index.
1407  //
1408  // in order to map it to the fe-system class, we have to know which base
1409  // element a degree of freedom within a vertex, line, etc belongs to. this
1410  // can be accomplished by the system_to_component_index function in
1411  // conjunction with the numbers first_{line,quad,...}_index
1412  for (unsigned int n = 0; n < this->interface_constraints.n(); ++n)
1413  for (unsigned int m = 0; m < this->interface_constraints.m(); ++m)
1414  {
1415  // for the pair (n,m) find out which base element they belong to and
1416  // the number therein
1417  //
1418  // first for the n index. this is simple since the n indices are in
1419  // the same order as they are usually on a face. note that for the
1420  // data type, first value in pair is (base element,instance of base
1421  // element), second is index within this instance
1422  const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
1423  n_index = this->face_system_to_base_table[n];
1424 
1425  // likewise for the m index. this is more complicated due to the
1426  // strange ordering we have for the dofs on the refined faces.
1427  std::pair<std::pair<unsigned int, unsigned int>, unsigned int> m_index;
1428  switch (dim)
1429  {
1430  case 1:
1431  {
1432  // we should never get here! (in 1d, the constraints matrix
1433  // should be of size zero)
1434  Assert(false, ExcInternalError());
1435  break;
1436  }
1437 
1438  case 2:
1439  {
1440  // the indices m=0..d_v-1 are from the center vertex. their
1441  // order is the same as for the first vertex of the whole cell,
1442  // so we can use the system_to_base_table variable (using the
1443  // face_s_t_base_t function would yield the same)
1444  if (m < this->dofs_per_vertex)
1445  m_index = this->system_to_base_table[m];
1446  else
1447  // then come the two sets of line indices
1448  {
1449  const unsigned int index_in_line =
1450  (m - this->dofs_per_vertex) % this->dofs_per_line;
1451  const unsigned int sub_line =
1452  (m - this->dofs_per_vertex) / this->dofs_per_line;
1453  Assert(sub_line < 2, ExcInternalError());
1454 
1455  // from this information, try to get base element and
1456  // instance of base element. we do so by constructing the
1457  // corresponding face index of m in the present element,
1458  // then use face_system_to_base_table
1459  const unsigned int tmp1 =
1460  2 * this->dofs_per_vertex + index_in_line;
1461  m_index.first = this->face_system_to_base_table[tmp1].first;
1462 
1463  // what we are still missing is the index of m within the
1464  // base elements interface_constraints table
1465  //
1466  // here, the second value of face_system_to_base_table can
1467  // help: it denotes the face index of that shape function
1468  // within the base element. since we know that it is a line
1469  // dof, we can construct the rest: tmp2 will denote the
1470  // index of this shape function among the line shape
1471  // functions:
1472  Assert(
1473  this->face_system_to_base_table[tmp1].second >=
1474  2 * base_element(m_index.first.first).dofs_per_vertex,
1475  ExcInternalError());
1476  const unsigned int tmp2 =
1477  this->face_system_to_base_table[tmp1].second -
1478  2 * base_element(m_index.first.first).dofs_per_vertex;
1479  Assert(tmp2 <
1480  base_element(m_index.first.first).dofs_per_line,
1481  ExcInternalError());
1482  m_index.second =
1483  base_element(m_index.first.first).dofs_per_vertex +
1484  base_element(m_index.first.first).dofs_per_line *
1485  sub_line +
1486  tmp2;
1487  }
1488  break;
1489  }
1490 
1491  case 3:
1492  {
1493  // same way as above, although a little more complicated...
1494 
1495  // the indices m=0..5*d_v-1 are from the center and the four
1496  // subline vertices. their order is the same as for the first
1497  // vertex of the whole cell, so we can use the simple arithmetic
1498  if (m < 5 * this->dofs_per_vertex)
1499  m_index = this->system_to_base_table[m];
1500  else
1501  // then come the 12 sets of line indices
1502  if (m < 5 * this->dofs_per_vertex + 12 * this->dofs_per_line)
1503  {
1504  // for the meaning of all this, see the 2d part
1505  const unsigned int index_in_line =
1506  (m - 5 * this->dofs_per_vertex) % this->dofs_per_line;
1507  const unsigned int sub_line =
1508  (m - 5 * this->dofs_per_vertex) / this->dofs_per_line;
1509  Assert(sub_line < 12, ExcInternalError());
1510 
1511  const unsigned int tmp1 =
1512  4 * this->dofs_per_vertex + index_in_line;
1513  m_index.first = this->face_system_to_base_table[tmp1].first;
1514 
1515  Assert(
1516  this->face_system_to_base_table[tmp1].second >=
1517  4 * base_element(m_index.first.first).dofs_per_vertex,
1518  ExcInternalError());
1519  const unsigned int tmp2 =
1520  this->face_system_to_base_table[tmp1].second -
1521  4 * base_element(m_index.first.first).dofs_per_vertex;
1522  Assert(tmp2 <
1523  base_element(m_index.first.first).dofs_per_line,
1524  ExcInternalError());
1525  m_index.second =
1526  5 * base_element(m_index.first.first).dofs_per_vertex +
1527  base_element(m_index.first.first).dofs_per_line *
1528  sub_line +
1529  tmp2;
1530  }
1531  else
1532  // on one of the four sub-quads
1533  {
1534  // for the meaning of all this, see the 2d part
1535  const unsigned int index_in_quad =
1536  (m - 5 * this->dofs_per_vertex -
1537  12 * this->dofs_per_line) %
1538  this->dofs_per_quad;
1539  Assert(index_in_quad < this->dofs_per_quad,
1540  ExcInternalError());
1541  const unsigned int sub_quad =
1542  ((m - 5 * this->dofs_per_vertex -
1543  12 * this->dofs_per_line) /
1544  this->dofs_per_quad);
1545  Assert(sub_quad < 4, ExcInternalError());
1546 
1547  const unsigned int tmp1 = 4 * this->dofs_per_vertex +
1548  4 * this->dofs_per_line +
1549  index_in_quad;
1550  Assert(tmp1 < this->face_system_to_base_table.size(),
1551  ExcInternalError());
1552  m_index.first = this->face_system_to_base_table[tmp1].first;
1553 
1554  Assert(
1555  this->face_system_to_base_table[tmp1].second >=
1556  4 * base_element(m_index.first.first).dofs_per_vertex +
1557  4 * base_element(m_index.first.first).dofs_per_line,
1558  ExcInternalError());
1559  const unsigned int tmp2 =
1560  this->face_system_to_base_table[tmp1].second -
1561  4 * base_element(m_index.first.first).dofs_per_vertex -
1562  4 * base_element(m_index.first.first).dofs_per_line;
1563  Assert(tmp2 <
1564  base_element(m_index.first.first).dofs_per_quad,
1565  ExcInternalError());
1566  m_index.second =
1567  5 * base_element(m_index.first.first).dofs_per_vertex +
1568  12 * base_element(m_index.first.first).dofs_per_line +
1569  base_element(m_index.first.first).dofs_per_quad *
1570  sub_quad +
1571  tmp2;
1572  }
1573 
1574  break;
1575  }
1576 
1577  default:
1578  Assert(false, ExcNotImplemented());
1579  }
1580 
1581  // now that we gathered all information: use it to build the
1582  // matrix. note that if n and m belong to different base elements or
1583  // instances, then there definitely will be no coupling
1584  if (n_index.first == m_index.first)
1585  this->interface_constraints(m, n) =
1586  (base_element(n_index.first.first)
1587  .constraints()(m_index.second, n_index.second));
1588  }
1589 }
1590 
1591 
1592 
1593 template <int dim, int spacedim>
1594 void
1596  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
1597  const std::vector<unsigned int> & multiplicities)
1598 {
1599  Assert(fes.size() == multiplicities.size(),
1600  ExcDimensionMismatch(fes.size(), multiplicities.size()));
1601  Assert(fes.size() > 0,
1602  ExcMessage("Need to pass at least one finite element."));
1603  Assert(count_nonzeros(multiplicities) > 0,
1604  ExcMessage("You only passed FiniteElements with multiplicity 0."));
1605 
1606  // Note that we need to skip every fe with multiplicity 0 in the following
1607  // block of code
1608 
1609  this->base_to_block_indices.reinit(0, 0);
1610 
1611  for (unsigned int i = 0; i < fes.size(); i++)
1612  if (multiplicities[i] > 0)
1613  this->base_to_block_indices.push_back(multiplicities[i]);
1614 
1615  {
1616  Threads::TaskGroup<> clone_base_elements;
1617 
1618  unsigned int ind = 0;
1619  for (unsigned int i = 0; i < fes.size(); i++)
1620  if (multiplicities[i] > 0)
1621  {
1622  clone_base_elements += Threads::new_task([&, i, ind]() {
1623  base_elements[ind] = {fes[i]->clone(), multiplicities[i]};
1624  });
1625  ++ind;
1626  }
1627  Assert(ind > 0, ExcInternalError());
1628 
1629  // wait for all of these clone operations to finish
1630  clone_base_elements.join_all();
1631  }
1632 
1633 
1634  {
1635  // If the system is not primitive, these have not been initialized by
1636  // FiniteElement
1637  this->system_to_component_table.resize(this->dofs_per_cell);
1638  this->face_system_to_component_table.resize(this->dofs_per_face);
1639 
1643  *this);
1644 
1648  *this);
1649  }
1650 
1651  // now initialize interface constraints, support points, and other tables.
1652  // (restriction and prolongation matrices are only built on demand.) do
1653  // this in parallel
1654 
1655  Threads::TaskGroup<> init_tasks;
1656 
1657  init_tasks +=
1658  Threads::new_task([&]() { this->build_interface_constraints(); });
1659 
1660  init_tasks += Threads::new_task([&]() {
1661  // if one of the base elements has no support points, then it makes no sense
1662  // to define support points for the composed element, so return an empty
1663  // array to demonstrate that fact. Note that we ignore FE_Nothing in this
1664  // logic.
1665  for (unsigned int base_el = 0; base_el < this->n_base_elements(); ++base_el)
1666  if (!base_element(base_el).has_support_points() &&
1667  base_element(base_el).dofs_per_cell != 0)
1668  {
1669  this->unit_support_points.resize(0);
1670  return;
1671  }
1672 
1673  // generate unit support points from unit support points of sub elements
1674  this->unit_support_points.resize(this->dofs_per_cell);
1675 
1676  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
1677  {
1678  const unsigned int base = this->system_to_base_table[i].first.first,
1679  base_index = this->system_to_base_table[i].second;
1680  Assert(base < this->n_base_elements(), ExcInternalError());
1681  Assert(base_index < base_element(base).unit_support_points.size(),
1682  ExcInternalError());
1683  this->unit_support_points[i] =
1684  base_element(base).unit_support_points[base_index];
1685  }
1686  });
1687 
1688  // initialize face support points (for dim==2,3). same procedure as above
1689  if (dim > 1)
1690  init_tasks += Threads::new_task([&]() {
1691  // if one of the base elements has no support points, then it makes no
1692  // sense to define support points for the composed element. In that case,
1693  // return an empty array to demonstrate that fact (note that we ask
1694  // whether the base element has no support points at all, not only none on
1695  // the face!)
1696  //
1697  // on the other hand, if there is an element that simply has no degrees of
1698  // freedom on the face at all, then we don't care whether it has support
1699  // points or not. this is, for example, the case for the stable Stokes
1700  // element Q(p)^dim \times DGP(p-1).
1701  for (unsigned int base_el = 0; base_el < this->n_base_elements();
1702  ++base_el)
1703  if (!base_element(base_el).has_support_points() &&
1704  (base_element(base_el).dofs_per_face > 0))
1705  {
1706  this->unit_face_support_points.resize(0);
1707  return;
1708  }
1709 
1710 
1711  // generate unit face support points from unit support points of sub
1712  // elements
1713  this->unit_face_support_points.resize(this->dofs_per_face);
1714 
1715  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
1716  {
1717  const unsigned int base_i =
1718  this->face_system_to_base_table[i].first.first;
1719  const unsigned int index_in_base =
1720  this->face_system_to_base_table[i].second;
1721 
1722  Assert(index_in_base <
1723  base_element(base_i).unit_face_support_points.size(),
1724  ExcInternalError());
1725 
1726  this->unit_face_support_points[i] =
1727  base_element(base_i).unit_face_support_points[index_in_base];
1728  }
1729  });
1730 
1731  // Initialize generalized support points and an (internal) index table
1732  init_tasks += Threads::new_task([&]() {
1733  // Iterate over all base elements, extract a representative set of
1734  // _unique_ generalized support points and store the information how
1735  // generalized support points of base elements are mapped to this list
1736  // of representatives. Complexity O(n^2), where n is the number of
1737  // generalized support points.
1738 
1740 
1741  for (unsigned int base = 0; base < this->n_base_elements(); ++base)
1742  {
1743  // If the current base element does not have generalized support
1744  // points, ignore it. Note that
1745  // * FESystem::convert_generalized_support_point_values_to_dof_values
1746  // will simply skip such non-interpolatory base elements by
1747  // assigning NaN to all dofs.
1748  // * If this routine does not pick up any generalized support
1749  // points the corresponding vector will be empty and
1750  // FiniteElement::has_generalized_support_points will return
1751  // false.
1752  if (!base_element(base).has_generalized_support_points())
1753  continue;
1754 
1755  for (const auto &point :
1756  base_element(base).get_generalized_support_points())
1757  {
1758  // Is point already an element of generalized_support_points?
1759  const auto p =
1760  std::find(std::begin(this->generalized_support_points),
1761  std::end(this->generalized_support_points),
1762  point);
1763 
1764  if (p == std::end(this->generalized_support_points))
1765  {
1766  // If no, update the table and add the point to the vector
1767  const auto n = this->generalized_support_points.size();
1768  generalized_support_points_index_table[base].push_back(n);
1769  this->generalized_support_points.push_back(point);
1770  }
1771  else
1772  {
1773  // If yes, just add the correct index to the table.
1774  const auto n = p - std::begin(this->generalized_support_points);
1775  generalized_support_points_index_table[base].push_back(n);
1776  }
1777  }
1778  }
1779 
1780 #ifdef DEBUG
1781  // check generalized_support_points_index_table for consistency
1782  for (unsigned int i = 0; i < base_elements.size(); ++i)
1783  {
1784  if (!base_element(i).has_generalized_support_points())
1785  continue;
1786 
1787  const auto &points =
1788  base_elements[i].first->get_generalized_support_points();
1789  for (unsigned int j = 0; j < points.size(); ++j)
1790  {
1791  const auto n = generalized_support_points_index_table[i][j];
1792  Assert(this->generalized_support_points[n] == points[j],
1793  ExcInternalError());
1794  }
1795  }
1796 #endif /* DEBUG */
1797  });
1798 
1799  // initialize quad dof index permutation in 3d and higher
1800  if (dim >= 3)
1801  init_tasks += Threads::new_task([&]() {
1802  // the array into which we want to write should have the correct size
1803  // already.
1805  .n_elements() == 8 * this->dofs_per_quad,
1806  ExcInternalError());
1807 
1808  // to obtain the shifts for this composed element, copy the shift
1809  // information of the base elements
1810  unsigned int index = 0;
1811  for (unsigned int b = 0; b < this->n_base_elements(); ++b)
1812  {
1813  const Table<2, int> &temp =
1814  this->base_element(b)
1815  .adjust_quad_dof_index_for_face_orientation_table;
1816  for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
1817  {
1818  for (unsigned int i = 0; i < temp.size(0); ++i)
1819  for (unsigned int j = 0; j < 8; ++j)
1821  index + i, j) = temp(i, j);
1822  index += temp.size(0);
1823  }
1824  }
1825  Assert(index == this->dofs_per_quad, ExcInternalError());
1826 
1827  // additionally compose the permutation information for lines
1829  this->dofs_per_line,
1830  ExcInternalError());
1831  index = 0;
1832  for (unsigned int b = 0; b < this->n_base_elements(); ++b)
1833  {
1834  const std::vector<int> &temp2 =
1835  this->base_element(b)
1836  .adjust_line_dof_index_for_line_orientation_table;
1837  for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
1838  {
1839  std::copy(
1840  temp2.begin(),
1841  temp2.end(),
1843  index);
1844  index += temp2.size();
1845  }
1846  }
1847  Assert(index == this->dofs_per_line, ExcInternalError());
1848  });
1849 
1850  // wait for all of this to finish
1851  init_tasks.join_all();
1852 }
1853 
1854 
1855 
1856 template <int dim, int spacedim>
1857 bool
1859 {
1860  for (unsigned int b = 0; b < this->n_base_elements(); ++b)
1861  if (base_element(b).hp_constraints_are_implemented() == false)
1862  return false;
1863 
1864  return true;
1865 }
1866 
1867 
1868 
1869 template <int dim, int spacedim>
1870 void
1872  const FiniteElement<dim, spacedim> &x_source_fe,
1873  FullMatrix<double> & interpolation_matrix) const
1874 {
1875  Assert(interpolation_matrix.n() == this->dofs_per_face,
1876  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
1877  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
1878  ExcDimensionMismatch(interpolation_matrix.m(),
1879  x_source_fe.dofs_per_face));
1880 
1881  // since dofs for each base are independent, we only have to stack things up
1882  // from base element to base element
1883  //
1884  // the problem is that we have to work with two FEs (this and
1885  // fe_other). only deal with the case that both are FESystems and that they
1886  // both have the same number of bases (counting multiplicity) each of which
1887  // match in their number of components. this covers
1888  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
1889  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
1890  // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
1891  if (const auto *fe_other_system =
1892  dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe))
1893  {
1894  // clear matrix, since we will not get to set all elements
1895  interpolation_matrix = 0;
1896 
1897  // loop over all the base elements of this and the other element, counting
1898  // their multiplicities
1899  unsigned int base_index = 0, base_index_other = 0;
1900  unsigned int multiplicity = 0, multiplicity_other = 0;
1901 
1902  FullMatrix<double> base_to_base_interpolation;
1903 
1904  while (true)
1905  {
1906  const FiniteElement<dim, spacedim> &base = base_element(base_index),
1907  &base_other =
1908  fe_other_system->base_element(
1909  base_index_other);
1910 
1911  Assert(base.n_components() == base_other.n_components(),
1912  ExcNotImplemented());
1913 
1914  // get the interpolation from the bases
1915  base_to_base_interpolation.reinit(base_other.dofs_per_face,
1916  base.dofs_per_face);
1917  base.get_face_interpolation_matrix(base_other,
1918  base_to_base_interpolation);
1919 
1920  // now translate entries. we'd like to have something like
1921  // face_base_to_system_index, but that doesn't exist. rather, all we
1922  // have is the reverse. well, use that then
1923  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
1924  if (this->face_system_to_base_index(i).first ==
1925  std::make_pair(base_index, multiplicity))
1926  for (unsigned int j = 0; j < fe_other_system->dofs_per_face; ++j)
1927  if (fe_other_system->face_system_to_base_index(j).first ==
1928  std::make_pair(base_index_other, multiplicity_other))
1929  interpolation_matrix(j, i) = base_to_base_interpolation(
1930  fe_other_system->face_system_to_base_index(j).second,
1931  this->face_system_to_base_index(i).second);
1932 
1933  // advance to the next base element for this and the other fe_system;
1934  // see if we can simply advance the multiplicity by one, or if have to
1935  // move on to the next base element
1936  ++multiplicity;
1937  if (multiplicity == this->element_multiplicity(base_index))
1938  {
1939  multiplicity = 0;
1940  ++base_index;
1941  }
1942  ++multiplicity_other;
1943  if (multiplicity_other ==
1944  fe_other_system->element_multiplicity(base_index_other))
1945  {
1946  multiplicity_other = 0;
1947  ++base_index_other;
1948  }
1949 
1950  // see if we have reached the end of the present element. if so, we
1951  // should have reached the end of the other one as well
1952  if (base_index == this->n_base_elements())
1953  {
1954  Assert(base_index_other == fe_other_system->n_base_elements(),
1955  ExcInternalError());
1956  break;
1957  }
1958 
1959  // if we haven't reached the end of this element, we shouldn't have
1960  // reached the end of the other one either
1961  Assert(base_index_other != fe_other_system->n_base_elements(),
1962  ExcInternalError());
1963  }
1964  }
1965  else
1966  {
1967  // repeat the cast to make the exception message more useful
1968  AssertThrow(
1969  (dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) !=
1970  nullptr),
1971  (typename FiniteElement<dim,
1972  spacedim>::ExcInterpolationNotImplemented()));
1973  }
1974 }
1975 
1976 
1977 
1978 template <int dim, int spacedim>
1979 void
1981  const FiniteElement<dim, spacedim> &x_source_fe,
1982  const unsigned int subface,
1983  FullMatrix<double> & interpolation_matrix) const
1984 {
1985  AssertThrow(
1986  (x_source_fe.get_name().find("FESystem<") == 0) ||
1987  (dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) != nullptr),
1989 
1990  Assert(interpolation_matrix.n() == this->dofs_per_face,
1991  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
1992  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
1993  ExcDimensionMismatch(interpolation_matrix.m(),
1994  x_source_fe.dofs_per_face));
1995 
1996  // since dofs for each base are independent, we only have to stack things up
1997  // from base element to base element
1998  //
1999  // the problem is that we have to work with two FEs (this and
2000  // fe_other). only deal with the case that both are FESystems and that they
2001  // both have the same number of bases (counting multiplicity) each of which
2002  // match in their number of components. this covers
2003  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2004  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2005  // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2006  const FESystem<dim, spacedim> *fe_other_system =
2007  dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe);
2008  if (fe_other_system != nullptr)
2009  {
2010  // clear matrix, since we will not get to set all elements
2011  interpolation_matrix = 0;
2012 
2013  // loop over all the base elements of this and the other element, counting
2014  // their multiplicities
2015  unsigned int base_index = 0, base_index_other = 0;
2016  unsigned int multiplicity = 0, multiplicity_other = 0;
2017 
2018  FullMatrix<double> base_to_base_interpolation;
2019 
2020  while (true)
2021  {
2022  const FiniteElement<dim, spacedim> &base = base_element(base_index),
2023  &base_other =
2024  fe_other_system->base_element(
2025  base_index_other);
2026 
2027  Assert(base.n_components() == base_other.n_components(),
2028  ExcNotImplemented());
2029 
2030  // get the interpolation from the bases
2031  base_to_base_interpolation.reinit(base_other.dofs_per_face,
2032  base.dofs_per_face);
2033  base.get_subface_interpolation_matrix(base_other,
2034  subface,
2035  base_to_base_interpolation);
2036 
2037  // now translate entries. we'd like to have something like
2038  // face_base_to_system_index, but that doesn't exist. rather, all we
2039  // have is the reverse. well, use that then
2040  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
2041  if (this->face_system_to_base_index(i).first ==
2042  std::make_pair(base_index, multiplicity))
2043  for (unsigned int j = 0; j < fe_other_system->dofs_per_face; ++j)
2044  if (fe_other_system->face_system_to_base_index(j).first ==
2045  std::make_pair(base_index_other, multiplicity_other))
2046  interpolation_matrix(j, i) = base_to_base_interpolation(
2047  fe_other_system->face_system_to_base_index(j).second,
2048  this->face_system_to_base_index(i).second);
2049 
2050  // advance to the next base element for this and the other fe_system;
2051  // see if we can simply advance the multiplicity by one, or if have to
2052  // move on to the next base element
2053  ++multiplicity;
2054  if (multiplicity == this->element_multiplicity(base_index))
2055  {
2056  multiplicity = 0;
2057  ++base_index;
2058  }
2059  ++multiplicity_other;
2060  if (multiplicity_other ==
2061  fe_other_system->element_multiplicity(base_index_other))
2062  {
2063  multiplicity_other = 0;
2064  ++base_index_other;
2065  }
2066 
2067  // see if we have reached the end of the present element. if so, we
2068  // should have reached the end of the other one as well
2069  if (base_index == this->n_base_elements())
2070  {
2071  Assert(base_index_other == fe_other_system->n_base_elements(),
2072  ExcInternalError());
2073  break;
2074  }
2075 
2076  // if we haven't reached the end of this element, we shouldn't have
2077  // reached the end of the other one either
2078  Assert(base_index_other != fe_other_system->n_base_elements(),
2079  ExcInternalError());
2080  }
2081  }
2082  else
2083  {
2084  // we should have caught this at the start, but check again anyway
2085  Assert(
2086  fe_other_system != nullptr,
2087  (typename FiniteElement<dim,
2088  spacedim>::ExcInterpolationNotImplemented()));
2089  }
2090 }
2091 
2092 
2093 
2094 template <int dim, int spacedim>
2095 template <int structdim>
2096 std::vector<std::pair<unsigned int, unsigned int>>
2098  const FiniteElement<dim, spacedim> &fe_other) const
2099 {
2100  // since dofs on each subobject (vertex, line, ...) are ordered such that
2101  // first come all from the first base element all multiplicities, then
2102  // second base element all multiplicities, etc., we simply have to stack all
2103  // the identities after each other
2104  //
2105  // the problem is that we have to work with two FEs (this and
2106  // fe_other). only deal with the case that both are FESystems and that they
2107  // both have the same number of bases (counting multiplicity) each of which
2108  // match in their number of components. this covers
2109  // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2110  // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2111  // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2112  if (const FESystem<dim, spacedim> *fe_other_system =
2113  dynamic_cast<const FESystem<dim, spacedim> *>(&fe_other))
2114  {
2115  // loop over all the base elements of this and the other element,
2116  // counting their multiplicities
2117  unsigned int base_index = 0, base_index_other = 0;
2118  unsigned int multiplicity = 0, multiplicity_other = 0;
2119 
2120  // we also need to keep track of the number of dofs already treated for
2121  // each of the elements
2122  unsigned int dof_offset = 0, dof_offset_other = 0;
2123 
2124  std::vector<std::pair<unsigned int, unsigned int>> identities;
2125 
2126  while (true)
2127  {
2128  const FiniteElement<dim, spacedim> &base = base_element(base_index),
2129  &base_other =
2130  fe_other_system->base_element(
2131  base_index_other);
2132 
2133  Assert(base.n_components() == base_other.n_components(),
2134  ExcNotImplemented());
2135 
2136  // now translate the identities returned by the base elements to the
2137  // indices of this system element
2138  std::vector<std::pair<unsigned int, unsigned int>> base_identities;
2139  switch (structdim)
2140  {
2141  case 0:
2142  base_identities = base.hp_vertex_dof_identities(base_other);
2143  break;
2144  case 1:
2145  base_identities = base.hp_line_dof_identities(base_other);
2146  break;
2147  case 2:
2148  base_identities = base.hp_quad_dof_identities(base_other);
2149  break;
2150  default:
2151  Assert(false, ExcNotImplemented());
2152  }
2153 
2154  for (const auto &base_identity : base_identities)
2155  identities.emplace_back(base_identity.first + dof_offset,
2156  base_identity.second + dof_offset_other);
2157 
2158  // record the dofs treated above as already taken care of
2159  dof_offset += base.template n_dofs_per_object<structdim>();
2160  dof_offset_other +=
2161  base_other.template n_dofs_per_object<structdim>();
2162 
2163  // advance to the next base element for this and the other
2164  // fe_system; see if we can simply advance the multiplicity by one,
2165  // or if have to move on to the next base element
2166  ++multiplicity;
2167  if (multiplicity == this->element_multiplicity(base_index))
2168  {
2169  multiplicity = 0;
2170  ++base_index;
2171  }
2172  ++multiplicity_other;
2173  if (multiplicity_other ==
2174  fe_other_system->element_multiplicity(base_index_other))
2175  {
2176  multiplicity_other = 0;
2177  ++base_index_other;
2178  }
2179 
2180  // see if we have reached the end of the present element. if so, we
2181  // should have reached the end of the other one as well
2182  if (base_index == this->n_base_elements())
2183  {
2184  Assert(base_index_other == fe_other_system->n_base_elements(),
2185  ExcInternalError());
2186  break;
2187  }
2188 
2189  // if we haven't reached the end of this element, we shouldn't have
2190  // reached the end of the other one either
2191  Assert(base_index_other != fe_other_system->n_base_elements(),
2192  ExcInternalError());
2193  }
2194 
2195  return identities;
2196  }
2197  else
2198  {
2199  Assert(false, ExcNotImplemented());
2200  return std::vector<std::pair<unsigned int, unsigned int>>();
2201  }
2202 }
2203 
2204 
2205 
2206 template <int dim, int spacedim>
2207 std::vector<std::pair<unsigned int, unsigned int>>
2209  const FiniteElement<dim, spacedim> &fe_other) const
2210 {
2211  return hp_object_dof_identities<0>(fe_other);
2212 }
2213 
2214 template <int dim, int spacedim>
2215 std::vector<std::pair<unsigned int, unsigned int>>
2217  const FiniteElement<dim, spacedim> &fe_other) const
2218 {
2219  return hp_object_dof_identities<1>(fe_other);
2220 }
2221 
2222 
2223 
2224 template <int dim, int spacedim>
2225 std::vector<std::pair<unsigned int, unsigned int>>
2227  const FiniteElement<dim, spacedim> &fe_other) const
2228 {
2229  return hp_object_dof_identities<2>(fe_other);
2230 }
2231 
2232 
2233 
2234 template <int dim, int spacedim>
2237  const FiniteElement<dim, spacedim> &fe_other,
2238  const unsigned int codim) const
2239 {
2240  Assert(codim <= dim, ExcImpossibleInDim(dim));
2241 
2242  // vertex/line/face/cell domination
2243  // --------------------------------
2244  // at present all we can do is to compare with other FESystems that have the
2245  // same number of components and bases
2246  if (const FESystem<dim, spacedim> *fe_sys_other =
2247  dynamic_cast<const FESystem<dim, spacedim> *>(&fe_other))
2248  {
2249  Assert(this->n_components() == fe_sys_other->n_components(),
2250  ExcNotImplemented());
2251  Assert(this->n_base_elements() == fe_sys_other->n_base_elements(),
2252  ExcNotImplemented());
2253 
2256 
2257  // loop over all base elements and do some sanity checks
2258  for (unsigned int b = 0; b < this->n_base_elements(); ++b)
2259  {
2260  Assert(this->base_element(b).n_components() ==
2261  fe_sys_other->base_element(b).n_components(),
2262  ExcNotImplemented());
2263  Assert(this->element_multiplicity(b) ==
2264  fe_sys_other->element_multiplicity(b),
2265  ExcNotImplemented());
2266 
2267  // for this pair of base elements, check who dominates and combine
2268  // with previous result
2269  const FiniteElementDomination::Domination base_domination =
2270  (this->base_element(b).compare_for_domination(
2271  fe_sys_other->base_element(b), codim));
2272  domination = domination & base_domination;
2273  }
2274 
2275  return domination;
2276  }
2277 
2278  Assert(false, ExcNotImplemented());
2280 }
2281 
2282 
2283 
2284 template <int dim, int spacedim>
2286 FESystem<dim, spacedim>::base_element(const unsigned int index) const
2287 {
2288  AssertIndexRange(index, base_elements.size());
2289  return *base_elements[index].first;
2290 }
2291 
2292 
2293 
2294 template <int dim, int spacedim>
2295 bool
2297  const unsigned int shape_index,
2298  const unsigned int face_index) const
2299 {
2300  return (base_element(this->system_to_base_index(shape_index).first.first)
2301  .has_support_on_face(this->system_to_base_index(shape_index).second,
2302  face_index));
2303 }
2304 
2305 
2306 
2307 template <int dim, int spacedim>
2308 Point<dim>
2309 FESystem<dim, spacedim>::unit_support_point(const unsigned int index) const
2310 {
2311  AssertIndexRange(index, this->dofs_per_cell);
2312  Assert((this->unit_support_points.size() == this->dofs_per_cell) ||
2313  (this->unit_support_points.size() == 0),
2315 
2316  // let's see whether we have the information pre-computed
2317  if (this->unit_support_points.size() != 0)
2318  return this->unit_support_points[index];
2319  else
2320  // no. ask the base element whether it would like to provide this
2321  // information
2322  return (base_element(this->system_to_base_index(index).first.first)
2323  .unit_support_point(this->system_to_base_index(index).second));
2324 }
2325 
2326 
2327 
2328 template <int dim, int spacedim>
2329 Point<dim - 1>
2331 {
2332  AssertIndexRange(index, this->dofs_per_face);
2333  Assert((this->unit_face_support_points.size() == this->dofs_per_face) ||
2334  (this->unit_face_support_points.size() == 0),
2336 
2337  // let's see whether we have the information pre-computed
2338  if (this->unit_face_support_points.size() != 0)
2339  return this->unit_face_support_points[index];
2340  else
2341  // no. ask the base element whether it would like to provide this
2342  // information
2343  return (base_element(this->face_system_to_base_index(index).first.first)
2344  .unit_face_support_point(
2345  this->face_system_to_base_index(index).second));
2346 }
2347 
2348 
2349 
2350 template <int dim, int spacedim>
2351 std::pair<Table<2, bool>, std::vector<unsigned int>>
2353 {
2354  // Note that this->n_components() is actually only an estimate of how many
2355  // constant modes we will need. There might be more than one such mode
2356  // (e.g. FE_Q_DG0).
2357  Table<2, bool> constant_modes(this->n_components(), this->dofs_per_cell);
2358  std::vector<unsigned int> components;
2359  for (unsigned int i = 0; i < base_elements.size(); ++i)
2360  {
2361  const std::pair<Table<2, bool>, std::vector<unsigned int>> base_table =
2362  base_elements[i].first->get_constant_modes();
2363  AssertDimension(base_table.first.n_rows(), base_table.second.size());
2364  const unsigned int element_multiplicity = this->element_multiplicity(i);
2365 
2366  // there might be more than one constant mode for some scalar elements,
2367  // so make sure the table actually fits: Create a new table with more
2368  // rows
2369  const unsigned int comp = components.size();
2370  if (constant_modes.n_rows() <
2371  comp + base_table.first.n_rows() * element_multiplicity)
2372  {
2373  Table<2, bool> new_constant_modes(comp + base_table.first.n_rows() *
2375  constant_modes.n_cols());
2376  for (unsigned int r = 0; r < comp; ++r)
2377  for (unsigned int c = 0; c < this->dofs_per_cell; ++c)
2378  new_constant_modes(r, c) = constant_modes(r, c);
2379  constant_modes.swap(new_constant_modes);
2380  }
2381 
2382  // next, fill the constant modes from the individual components as well
2383  // as the component numbers corresponding to the constant mode rows
2384  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
2385  {
2386  std::pair<std::pair<unsigned int, unsigned int>, unsigned int> ind =
2387  this->system_to_base_index(k);
2388  if (ind.first.first == i)
2389  for (unsigned int c = 0; c < base_table.first.n_rows(); ++c)
2390  constant_modes(comp +
2391  ind.first.second * base_table.first.n_rows() + c,
2392  k) = base_table.first(c, ind.second);
2393  }
2394  for (unsigned int r = 0; r < element_multiplicity; ++r)
2395  for (const unsigned int c : base_table.second)
2396  components.push_back(
2397  comp + r * this->base_elements[i].first->n_components() + c);
2398  }
2399  AssertDimension(components.size(), constant_modes.n_rows());
2400  return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
2401  components);
2402 }
2403 
2404 
2405 
2406 template <int dim, int spacedim>
2407 void
2409  const std::vector<Vector<double>> &point_values,
2410  std::vector<double> & dof_values) const
2411 {
2413  ExcMessage("The FESystem does not have generalized support points"));
2414 
2415  AssertDimension(point_values.size(),
2416  this->get_generalized_support_points().size());
2417  AssertDimension(dof_values.size(), this->dofs_per_cell);
2418 
2419  std::vector<double> base_dof_values;
2420  std::vector<Vector<double>> base_point_values;
2421 
2422  // loop over all base elements (respecting multiplicity) and let them do
2423  // the work on their share of the input argument
2424 
2425  unsigned int current_vector_component = 0;
2426  for (unsigned int base = 0; base < base_elements.size(); ++base)
2427  {
2428  // We need access to the base_element, its multiplicity, the
2429  // number of generalized support points (n_base_points) and the
2430  // number of components we're dealing with.
2431  const auto & base_element = this->base_element(base);
2432  const unsigned int multiplicity = this->element_multiplicity(base);
2433  const unsigned int n_base_dofs = base_element.dofs_per_cell;
2434  const unsigned int n_base_components = base_element.n_components();
2435 
2436  // If the number of base degrees of freedom is zero, there is nothing
2437  // to do, skip the rest of the body in this case and continue with
2438  // the next element
2439  if (n_base_dofs == 0)
2440  {
2441  current_vector_component += multiplicity * n_base_components;
2442  continue;
2443  }
2444 
2445  if (base_element.has_generalized_support_points())
2446  {
2447  const unsigned int n_base_points =
2448  base_element.get_generalized_support_points().size();
2449 
2450  base_dof_values.resize(n_base_dofs);
2451  base_point_values.resize(n_base_points);
2452 
2453  for (unsigned int m = 0; m < multiplicity;
2454  ++m, current_vector_component += n_base_components)
2455  {
2456  // populate base_point_values for a recursive call to
2457  // convert_generalized_support_point_values_to_dof_values
2458  for (unsigned int j = 0; j < base_point_values.size(); ++j)
2459  {
2460  base_point_values[j].reinit(n_base_components, false);
2461 
2462  const auto n =
2464 
2465  // we have to extract the correct slice out of the global
2466  // vector of values:
2467  const auto begin =
2468  std::begin(point_values[n]) + current_vector_component;
2469  const auto end = begin + n_base_components;
2470  std::copy(begin, end, std::begin(base_point_values[j]));
2471  }
2472 
2473  base_element
2474  .convert_generalized_support_point_values_to_dof_values(
2475  base_point_values, base_dof_values);
2476 
2477  // Finally put these dof values back into global dof values
2478  // vector.
2479 
2480  // To do this, we could really use a base_to_system_index()
2481  // function, but that doesn't exist -- just do it by using the
2482  // reverse table -- the amount of work done here is not worth
2483  // trying to optimizing this.
2484  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
2485  if (this->system_to_base_index(i).first ==
2486  std::make_pair(base, m))
2487  dof_values[i] =
2488  base_dof_values[this->system_to_base_index(i).second];
2489  } /*for*/
2490  }
2491  else
2492  {
2493  // If the base element is non-interpolatory, assign NaN to all
2494  // DoFs associated to it.
2495 
2496  // To do this, we could really use a base_to_system_index()
2497  // function, but that doesn't exist -- just do it by using the
2498  // reverse table -- the amount of work done here is not worth
2499  // trying to optimizing this.
2500  for (unsigned int m = 0; m < multiplicity; ++m)
2501  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
2502  if (this->system_to_base_index(i).first ==
2503  std::make_pair(base, m))
2504  dof_values[i] = std::numeric_limits<double>::signaling_NaN();
2505 
2506  current_vector_component += multiplicity * n_base_components;
2507  }
2508  } /*for*/
2509 }
2510 
2511 
2512 
2513 template <int dim, int spacedim>
2514 std::size_t
2516 {
2517  // neglect size of data stored in @p{base_elements} due to some problems
2518  // with the compiler. should be neglectable after all, considering the size
2519  // of the data of the subelements
2521  sizeof(base_elements));
2522  for (unsigned int i = 0; i < base_elements.size(); ++i)
2524  return mem;
2525 }
2526 
2527 
2528 
2529 // explicit instantiations
2530 #include "fe_system.inst"
2531 
2532 DEAL_II_NAMESPACE_CLOSE
~InternalData() override
Definition: fe_system.cc:59
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
virtual std::size_t memory_consumption() const override
Definition: fe_system.cc:2515
Shape function values.
virtual Point< dim > unit_support_point(const unsigned int index) const override
Definition: fe_system.cc:2309
size_type m() const
std::vector< internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > > base_fe_output_objects
Definition: fe_system.h:1234
static const unsigned int invalid_unsigned_int
Definition: types.h:187
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_system.cc:402
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2442
void build_interface_constraints()
Definition: fe_system.cc:1391
const unsigned int components
Definition: fe_base.h:292
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe.cc:878
void initialize(const std::vector< const FiniteElement< dim, spacedim > *> &fes, const std::vector< unsigned int > &multiplicities)
Definition: fe_system.cc:1595
virtual std::string get_name() const override
Definition: fe_system.cc:311
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2493
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_system.cc:2216
void swap(TableBase< N, T > &v)
std::vector< std::pair< unsigned int, unsigned int > > hp_object_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe_system.cc:2097
FullMatrix< double > interface_constraints
Definition: fe.h:2468
unsigned int n_nonzero_components(const unsigned int i) const
Definition: fe.h:3291
virtual Tensor< 4, dim > shape_4th_derivative(const unsigned int i, const Point< dim > &p) const override
Definition: fe_system.cc:573
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_system.cc:2352
Task< RT > new_task(const std::function< RT()> &function)
const std::vector< Point< dim > > & get_generalized_support_points() const
Definition: fe.cc:1016
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double >> &support_point_values, std::vector< double > &dof_values) const override
Definition: fe_system.cc:2408
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_system.cc:2296
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
const unsigned int dofs_per_quad
Definition: fe_base.h:237
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:933
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_system.cc:454
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const override
Definition: fe_system.cc:2286
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition: fe.h:2602
bool has_generalized_support_points() const
Definition: fe.cc:1029
std::vector< Point< dim - 1 > > unit_face_support_points
Definition: fe.h:2487
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
Definition: fe_system.cc:386
bool is_primitive() const
Definition: fe.h:3301
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 override
Definition: fe_system.cc:851
static ::ExceptionBase & ExcInterpolationNotImplemented()
virtual const FiniteElement< dim, spacedim > & get_sub_fe(const unsigned int first_component, const unsigned int n_selected_components) const override
Definition: fe_system.cc:357
virtual Tensor< 4, dim > shape_4th_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_system.cc:589
const unsigned int dofs_per_line
Definition: fe_base.h:231
InternalData(const unsigned int n_base_elements)
Definition: fe_system.cc:50
virtual Tensor< 3, dim > shape_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
Definition: fe_system.cc:528
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2480
std::vector< std::vector< std::size_t > > generalized_support_points_index_table
Definition: fe_system.h:1138
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_system.cc:340
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: fe_system.cc:896
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
virtual std::size_t memory_consumption() const
Definition: fe.cc:1204
size_type n() const
No update.
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2456
void build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int >> &system_to_base_table, std::vector< std::pair< unsigned int, unsigned int >> &system_to_component_table, std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int >> &component_to_base_table, const FiniteElement< dim, spacedim > &finite_element, const bool do_tensor_product=true)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
Definition: fe.h:3231
Third derivatives of shape functions.
std::mutex mutex
Definition: fe_system.h:1240
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1411
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3133
UpdateFlags
static const unsigned int invalid_face_number
Definition: fe_system.h:1113
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_system.cc:785
Abstract base class for mapping classes.
Definition: mapping.h:302
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
Definition: fe.h:44
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
Definition: fe.cc:894
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:3001
void set_fe_data(const unsigned int base_no, std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase >)
Definition: fe_system.cc:80
virtual std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim - 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
Definition: fe_system.cc:1036
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
Definition: fe_system.cc:483
virtual std::string get_name() const =0
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:911
Second derivatives of shape functions.
virtual bool hp_constraints_are_implemented() const override
Definition: fe_system.cc:1858
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:514
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
std::vector< std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > > base_fe_datas
Definition: fe_system.h:1222
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_system.cc:2208
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_system.cc:1980
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_system.cc:2236
virtual Tensor< 3, dim > shape_3rd_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_system.cc:544
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_system.cc:618
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:922
size_type size(const unsigned int i) const
Table< 2, int > adjust_quad_dof_index_for_face_orientation_table
Definition: fe.h:2516
std::pair< std::pair< unsigned int, unsigned int >, unsigned int > face_system_to_base_index(const unsigned int index) const
Definition: fe.h:3210
unsigned int n_components() const
std::pair< std::pair< unsigned int, unsigned int >, unsigned int > system_to_base_index(const unsigned int index) const
Definition: fe.h:3199
FiniteElement< dim, spacedim >::InternalDataBase & get_fe_data(const unsigned int base_no) const
Definition: fe_system.cc:69
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_system.cc:694
Shape function gradients.
std::vector< std::pair< unsigned int, unsigned int > > system_to_component_table
Definition: fe.h:2536
FESystem(const FiniteElement< dim, spacedim > &fe, const unsigned int n_elements)
Definition: fe_system.cc:109
void push_back(const size_type size)
void build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int >> &face_system_to_base_table, std::vector< std::pair< unsigned int, unsigned int >> &face_system_to_component_table, const FiniteElement< dim, spacedim > &finite_element, const bool do_tensor_product=true)
const unsigned int dofs_per_face
Definition: fe_base.h:275
const std::vector< ComponentMask > nonzero_components
Definition: fe.h:2618
internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > & get_fe_output_object(const unsigned int base_no) const
Definition: fe_system.cc:92
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > system_to_base_table
Definition: fe.h:2567
static ::ExceptionBase & ExcNotImplemented()
void compute_fill(const Mapping< dim, spacedim > &mapping, const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim_1 > &quadrature, const CellSimilarity::Similarity cell_similarity, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_data, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe_system.cc:1189
const unsigned int dofs_per_vertex
Definition: fe_base.h:225
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
Definition: fe_system.cc:438
std::vector< std::pair< unsigned int, unsigned int > > face_system_to_component_table
Definition: fe.h:2548
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_system.cc:1871
BlockIndices base_to_block_indices
Definition: fe.h:2579
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_system.cc:499
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > face_system_to_base_table
Definition: fe.h:2573
unsigned int n_base_elements() const
Definition: fe.h:3124
virtual void fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
std::vector< int > adjust_line_dof_index_for_line_orientation_table
Definition: fe.h:2531
TableIndices< 2 > interface_constraints_size() const
Definition: fe.cc:839
virtual std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
Definition: fe_system.cc:910
virtual std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim - 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
Definition: fe_system.cc:973
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_system.cc:2226
std::vector< std::pair< std::unique_ptr< const FiniteElement< dim, spacedim > >, unsigned int > > base_elements
Definition: fe_system.h:1124
UpdateFlags update_each
Definition: fe.h:715
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
virtual Point< dim - 1 > unit_face_support_point(const unsigned int index) const override
Definition: fe_system.cc:2330
static ::ExceptionBase & ExcInternalError()