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