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