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