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