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