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