Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+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, GeometryInfo<dim>::n_children(refinement_case));
896
897 // initialization upon first request
898 if (this->restriction[refinement_case - 1][child].n() == 0)
899 {
900 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
901
902 // check if updated while waiting for lock
903 if (this->restriction[refinement_case - 1][child].n() ==
904 this->n_dofs_per_cell())
905 return this->restriction[refinement_case - 1][child];
906
907 // shortcut for accessing local restrictions further down
908 std::vector<const FullMatrix<double> *> base_matrices(
909 this->n_base_elements());
910
911 for (unsigned int i = 0; i < this->n_base_elements(); ++i)
912 {
913 base_matrices[i] =
914 &base_element(i).get_restriction_matrix(child, refinement_case);
915
916 Assert(base_matrices[i]->n() == base_element(i).n_dofs_per_cell(),
918 }
919
920 FullMatrix<double> restriction(this->n_dofs_per_cell(),
921 this->n_dofs_per_cell());
922
923 // distribute the matrices of the base finite elements to the
924 // matrices of this object. for this, loop over all degrees of
925 // freedom and take the respective entry of the underlying base
926 // element.
927 //
928 // note that we by definition of a base element, they are
929 // independent, i.e. do not couple. only DoFs that belong to the
930 // same instance of a base element may couple
931 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
932 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
933 {
934 // first find out to which base element indices i and j
935 // belong, and which instance thereof in case the base element
936 // has a multiplicity greater than one. if they should not
937 // happen to belong to the same instance of a base element,
938 // then they cannot couple, so go on with the next index
939 if (this->system_to_base_table[i].first !=
940 this->system_to_base_table[j].first)
941 continue;
942
943 // so get the common base element and the indices therein:
944 const unsigned int base = this->system_to_base_table[i].first.first;
945
946 const unsigned int base_index_i =
947 this->system_to_base_table[i].second,
948 base_index_j =
949 this->system_to_base_table[j].second;
950
951 // if we are sure that DoFs i and j may couple, then copy
952 // entries of the matrices:
953 restriction(i, j) =
954 (*base_matrices[base])(base_index_i, base_index_j);
955 }
956
957 const_cast<FullMatrix<double> &>(
958 this->restriction[refinement_case - 1][child]) = std::move(restriction);
959 }
960
961 return this->restriction[refinement_case - 1][child];
962}
963
964
965
966template <int dim, int spacedim>
967const FullMatrix<double> &
969 const unsigned int child,
970 const RefinementCase<dim> &refinement_case) const
971{
972 AssertIndexRange(refinement_case,
976 "Restriction matrices are only available for refined cells!"));
977 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
978
979 // initialization upon first request, construction completely analogous to
980 // restriction matrix
981 if (this->prolongation[refinement_case - 1][child].n() == 0)
982 {
983 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
984
985 if (this->prolongation[refinement_case - 1][child].n() ==
986 this->n_dofs_per_cell())
987 return this->prolongation[refinement_case - 1][child];
988
989 std::vector<const FullMatrix<double> *> base_matrices(
990 this->n_base_elements());
991 for (unsigned int i = 0; i < this->n_base_elements(); ++i)
992 {
993 base_matrices[i] =
994 &base_element(i).get_prolongation_matrix(child, refinement_case);
995
996 Assert(base_matrices[i]->n() == base_element(i).n_dofs_per_cell(),
998 }
999
1000 FullMatrix<double> prolongate(this->n_dofs_per_cell(),
1001 this->n_dofs_per_cell());
1002
1003 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1004 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
1005 {
1006 if (this->system_to_base_table[i].first !=
1007 this->system_to_base_table[j].first)
1008 continue;
1009 const unsigned int base = this->system_to_base_table[i].first.first;
1010
1011 const unsigned int base_index_i =
1012 this->system_to_base_table[i].second,
1013 base_index_j =
1014 this->system_to_base_table[j].second;
1015 prolongate(i, j) =
1016 (*base_matrices[base])(base_index_i, base_index_j);
1017 }
1018
1019 const_cast<FullMatrix<double> &>(
1020 this->prolongation[refinement_case - 1][child]) = std::move(prolongate);
1021 }
1022
1023 return this->prolongation[refinement_case - 1][child];
1024}
1025
1026
1027template <int dim, int spacedim>
1028unsigned int
1030 const unsigned int face_dof_index,
1031 const unsigned int face,
1032 const unsigned char combined_orientation) const
1033{
1034 // we need to ask the base elements how they want to translate
1035 // the DoFs within their own numbering. thus, translate to
1036 // the base element numbering and then back
1037 const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
1038 face_base_index = this->face_system_to_base_index(face_dof_index, face);
1039
1040 const unsigned int base_face_to_cell_index =
1041 this->base_element(face_base_index.first.first)
1042 .face_to_cell_index(face_base_index.second, face, combined_orientation);
1043
1044 // it would be nice if we had a base_to_system_index function, but
1045 // all that exists is a component_to_system_index function. we can't do
1046 // this here because it won't work for non-primitive elements. consequently,
1047 // simply do a loop over all dofs till we find whether it corresponds
1048 // to the one we're interested in -- crude, maybe, but works for now
1049 const std::pair<std::pair<unsigned int, unsigned int>, unsigned int> target =
1050 std::make_pair(face_base_index.first, base_face_to_cell_index);
1051 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1052 if (this->system_to_base_index(i) == target)
1053 return i;
1054
1057}
1058
1059
1060
1061//---------------------------------------------------------------------------
1062// Data field initialization
1063//---------------------------------------------------------------------------
1064
1065
1066
1067template <int dim, int spacedim>
1070{
1072 // generate maximal set of flags
1073 // that are necessary
1074 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1075 out |= base_element(base_no).requires_update_flags(flags);
1076 return out;
1077}
1078
1079
1080
1081template <int dim, int spacedim>
1082std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1084 const UpdateFlags flags,
1085 const Mapping<dim, spacedim> &mapping,
1086 const Quadrature<dim> &quadrature,
1088 spacedim>
1089 & /*output_data*/) const
1090{
1091 // create an internal data object and set the update flags we will need
1092 // to deal with. the current object does not make use of these flags,
1093 // but we need to nevertheless set them correctly since we look
1094 // into the update_each flag of base elements in fill_fe_values,
1095 // and so the current object's update_each flag needs to be
1096 // correct in case the current FESystem is a base element for another,
1097 // higher-level FESystem itself.
1098 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1099 data_ptr = std::make_unique<InternalData>(this->n_base_elements());
1100 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1101 data.update_each = requires_update_flags(flags);
1102
1103 // get data objects from each of the base elements and store
1104 // them. one might think that doing this in parallel (over the
1105 // base elements) would be a good idea, but this turns out to
1106 // be wrong because we would then run these jobs on different
1107 // threads/processors and this allocates memory in different
1108 // NUMA domains; this has large detrimental effects when later
1109 // writing into these objects in fill_fe_*_values. all of this
1110 // is particularly true when using FEValues objects in
1111 // WorkStream contexts where we explicitly make sure that
1112 // every function only uses objects previously allocated
1113 // in the same NUMA context and on the same thread as the
1114 // function is called
1115 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1116 {
1118 &base_fe_output_object = data.get_fe_output_object(base_no);
1119 base_fe_output_object.initialize(
1120 quadrature.size(),
1121 base_element(base_no),
1122 flags | base_element(base_no).requires_update_flags(flags));
1123
1124 // let base objects produce their scratch objects. they may
1125 // also at this time write into the output objects we provide
1126 // for them; it would be nice if we could already copy something
1127 // out of the base output object into the system output object,
1128 // but we can't because we can't know what the elements already
1129 // copied and/or will want to update on every cell
1130 auto base_fe_data = base_element(base_no).get_data(flags,
1131 mapping,
1132 quadrature,
1133 base_fe_output_object);
1134
1135 data.set_fe_data(base_no, std::move(base_fe_data));
1136 }
1137
1138 return data_ptr;
1139}
1140
1141// The following function is a clone of get_data, with the exception
1142// that get_face_data of the base elements is called.
1143
1144template <int dim, int spacedim>
1145std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1147 const UpdateFlags flags,
1148 const Mapping<dim, spacedim> &mapping,
1149 const hp::QCollection<dim - 1> &quadrature,
1151 spacedim>
1152 & /*output_data*/) const
1153{
1154 // create an internal data object and set the update flags we will need
1155 // to deal with. the current object does not make use of these flags,
1156 // but we need to nevertheless set them correctly since we look
1157 // into the update_each flag of base elements in fill_fe_values,
1158 // and so the current object's update_each flag needs to be
1159 // correct in case the current FESystem is a base element for another,
1160 // higher-level FESystem itself.
1161 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1162 data_ptr = std::make_unique<InternalData>(this->n_base_elements());
1163 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1164 data.update_each = requires_update_flags(flags);
1165
1166 // get data objects from each of the base elements and store
1167 // them. one might think that doing this in parallel (over the
1168 // base elements) would be a good idea, but this turns out to
1169 // be wrong because we would then run these jobs on different
1170 // threads/processors and this allocates memory in different
1171 // NUMA domains; this has large detrimental effects when later
1172 // writing into these objects in fill_fe_*_values. all of this
1173 // is particularly true when using FEValues objects in
1174 // WorkStream contexts where we explicitly make sure that
1175 // every function only uses objects previously allocated
1176 // in the same NUMA context and on the same thread as the
1177 // function is called
1178 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1179 {
1181 &base_fe_output_object = data.get_fe_output_object(base_no);
1182 base_fe_output_object.initialize(
1183 quadrature.max_n_quadrature_points(),
1184 base_element(base_no),
1185 flags | base_element(base_no).requires_update_flags(flags));
1186
1187 // let base objects produce their scratch objects. they may
1188 // also at this time write into the output objects we provide
1189 // for them; it would be nice if we could already copy something
1190 // out of the base output object into the system output object,
1191 // but we can't because we can't know what the elements already
1192 // copied and/or will want to update on every cell
1193 auto base_fe_data = base_element(base_no).get_face_data(
1194 flags, mapping, quadrature, base_fe_output_object);
1195
1196 data.set_fe_data(base_no, std::move(base_fe_data));
1197 }
1198
1199 return data_ptr;
1200}
1201
1202
1203
1204// The following function is a clone of get_data, with the exception
1205// that get_subface_data of the base elements is called.
1206
1207template <int dim, int spacedim>
1208std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1210 const UpdateFlags flags,
1211 const Mapping<dim, spacedim> &mapping,
1212 const Quadrature<dim - 1> &quadrature,
1214 spacedim>
1215 & /*output_data*/) const
1216{
1217 // create an internal data object and set the update flags we will need
1218 // to deal with. the current object does not make use of these flags,
1219 // but we need to nevertheless set them correctly since we look
1220 // into the update_each flag of base elements in fill_fe_values,
1221 // and so the current object's update_each flag needs to be
1222 // correct in case the current FESystem is a base element for another,
1223 // higher-level FESystem itself.
1224 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1225 data_ptr = std::make_unique<InternalData>(this->n_base_elements());
1226 auto &data = dynamic_cast<InternalData &>(*data_ptr);
1227
1228 data.update_each = requires_update_flags(flags);
1229
1230 // get data objects from each of the base elements and store
1231 // them. one might think that doing this in parallel (over the
1232 // base elements) would be a good idea, but this turns out to
1233 // be wrong because we would then run these jobs on different
1234 // threads/processors and this allocates memory in different
1235 // NUMA domains; this has large detrimental effects when later
1236 // writing into these objects in fill_fe_*_values. all of this
1237 // is particularly true when using FEValues objects in
1238 // WorkStream contexts where we explicitly make sure that
1239 // every function only uses objects previously allocated
1240 // in the same NUMA context and on the same thread as the
1241 // function is called
1242 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1243 {
1245 &base_fe_output_object = data.get_fe_output_object(base_no);
1246 base_fe_output_object.initialize(
1247 quadrature.size(),
1248 base_element(base_no),
1249 flags | base_element(base_no).requires_update_flags(flags));
1250
1251 // let base objects produce their scratch objects. they may
1252 // also at this time write into the output objects we provide
1253 // for them; it would be nice if we could already copy something
1254 // out of the base output object into the system output object,
1255 // but we can't because we can't know what the elements already
1256 // copied and/or will want to update on every cell
1257 auto base_fe_data = base_element(base_no).get_subface_data(
1258 flags, mapping, quadrature, base_fe_output_object);
1259
1260 data.set_fe_data(base_no, std::move(base_fe_data));
1261 }
1262
1263 return data_ptr;
1264}
1265
1266
1267
1268template <int dim, int spacedim>
1269void
1272 const CellSimilarity::Similarity cell_similarity,
1273 const Quadrature<dim> &quadrature,
1274 const Mapping<dim, spacedim> &mapping,
1275 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1277 &mapping_data,
1278 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1280 spacedim>
1281 &output_data) const
1282{
1283 compute_fill(mapping,
1284 cell,
1285 invalid_face_number,
1286 invalid_face_number,
1287 quadrature,
1288 cell_similarity,
1289 mapping_internal,
1290 fe_internal,
1291 mapping_data,
1292 output_data);
1293}
1294
1295
1296
1297template <int dim, int spacedim>
1298void
1301 const unsigned int face_no,
1302 const hp::QCollection<dim - 1> &quadrature,
1303 const Mapping<dim, spacedim> &mapping,
1304 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1306 &mapping_data,
1307 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1309 spacedim>
1310 &output_data) const
1311{
1312 compute_fill(mapping,
1313 cell,
1314 face_no,
1315 invalid_face_number,
1316 quadrature,
1318 mapping_internal,
1319 fe_internal,
1320 mapping_data,
1321 output_data);
1322}
1323
1324
1325
1326template <int dim, int spacedim>
1327void
1330 const unsigned int face_no,
1331 const unsigned int sub_no,
1332 const Quadrature<dim - 1> &quadrature,
1333 const Mapping<dim, spacedim> &mapping,
1334 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1336 &mapping_data,
1337 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1339 spacedim>
1340 &output_data) const
1341{
1342 compute_fill(mapping,
1343 cell,
1344 face_no,
1345 sub_no,
1346 quadrature,
1348 mapping_internal,
1349 fe_internal,
1350 mapping_data,
1351 output_data);
1352}
1353
1354
1355
1356template <int dim, int spacedim>
1357template <class Q_or_QC>
1358void
1360 const Mapping<dim, spacedim> &mapping,
1362 const unsigned int face_no,
1363 const unsigned int sub_no,
1364 const Q_or_QC &quadrature,
1365 const CellSimilarity::Similarity cell_similarity,
1366 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1367 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1369 &mapping_data,
1371 &output_data) const
1372{
1373 // convert data object to internal
1374 // data for this class. fails with
1375 // an exception if that is not
1376 // possible
1377 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1379 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1380
1381 const UpdateFlags flags = fe_data.update_each;
1382
1383
1384 // loop over the base elements, let them compute what they need to compute,
1385 // and then copy what is necessary.
1386 //
1387 // one may think that it would be a good idea to parallelize this over
1388 // base elements, but it turns out to be not worthwhile: doing so lets
1389 // multiple threads access data objects that were created by the current
1390 // thread, leading to many NUMA memory access inefficiencies. we specifically
1391 // want to avoid this if this class is called in a WorkStream context where
1392 // we very carefully allocate objects only on the thread where they
1393 // will actually be used; spawning new tasks here would be counterproductive
1396 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1397 {
1398 const FiniteElement<dim, spacedim> &base_fe = base_element(base_no);
1399 typename FiniteElement<dim, spacedim>::InternalDataBase &base_fe_data =
1400 fe_data.get_fe_data(base_no);
1402 spacedim>
1403 &base_data = fe_data.get_fe_output_object(base_no);
1404
1405 // If we have mixed meshes we need to support a QCollection here, hence
1406 // this pointer casting workaround:
1407 const Quadrature<dim> *cell_quadrature = nullptr;
1408 const hp::QCollection<dim - 1> *face_quadrature = nullptr;
1409 const Quadrature<dim - 1> *sub_face_quadrature = nullptr;
1410 unsigned int n_q_points = numbers::invalid_unsigned_int;
1411
1412 // static cast through the common base class:
1413 if (face_no == invalid_face_number)
1414 {
1415 cell_quadrature =
1416 dynamic_cast<const Quadrature<dim> *>(&quadrature);
1417 Assert(cell_quadrature != nullptr, ExcInternalError());
1418 n_q_points = cell_quadrature->size();
1419 }
1420 else if (sub_no == invalid_face_number)
1421 {
1422 // If we don't have wedges or pyramids then there should only be one
1423 // quadrature rule here
1424 face_quadrature =
1425 dynamic_cast<const hp::QCollection<dim - 1> *>(&quadrature);
1426 Assert(face_quadrature != nullptr, ExcInternalError());
1427
1428 n_q_points =
1429 (*face_quadrature)[face_quadrature->size() == 1 ? 0 : face_no]
1430 .size();
1431 }
1432 else
1433 {
1434 sub_face_quadrature =
1435 dynamic_cast<const Quadrature<dim - 1> *>(&quadrature);
1436 Assert(sub_face_quadrature != nullptr, ExcInternalError());
1437
1438 n_q_points = sub_face_quadrature->size();
1439 }
1441
1442
1443 // Make sure that in the case of fill_fe_values the data is only
1444 // copied from base_data to data if base_data is changed. therefore
1445 // use fe_fe_data.current_update_flags()
1446 //
1447 // for the case of fill_fe_(sub)face_values the data needs to be
1448 // copied from base_data to data on each face, therefore use
1449 // base_fe_data.update_flags.
1450 if (face_no == invalid_face_number)
1451 base_fe.fill_fe_values(cell,
1452 cell_similarity,
1453 *cell_quadrature,
1454 mapping,
1455 mapping_internal,
1456 mapping_data,
1457 base_fe_data,
1458 base_data);
1459 else if (sub_no == invalid_face_number)
1460 base_fe.fill_fe_face_values(cell,
1461 face_no,
1462 *face_quadrature,
1463 mapping,
1464 mapping_internal,
1465 mapping_data,
1466 base_fe_data,
1467 base_data);
1468 else
1469 base_fe.fill_fe_subface_values(cell,
1470 face_no,
1471 sub_no,
1472 *sub_face_quadrature,
1473 mapping,
1474 mapping_internal,
1475 mapping_data,
1476 base_fe_data,
1477 base_data);
1478
1479 // now data has been generated, so copy it. This procedure is different
1480 // for primitive and non-primitive base elements, so at this point we
1481 // dispatch to helper functions.
1482 const UpdateFlags base_flags = base_fe_data.update_each;
1483
1484 if (base_fe.is_primitive())
1485 {
1487 *this,
1488 base_no,
1489 base_flags,
1490 primitive_offset_tables[base_no],
1491 base_data,
1492 output_data);
1493 }
1494 else
1495 {
1497 *this,
1498 base_no,
1499 n_q_points,
1500 base_flags,
1501 nonprimitive_offset_tables[base_no],
1502 base_data,
1503 output_data);
1504 }
1505 }
1506}
1507
1508
1509template <int dim, int spacedim>
1510void
1512{
1513 // check whether all base elements implement their interface constraint
1514 // matrices. if this is not the case, then leave the interface constraints of
1515 // this composed element empty as well; however, the rest of the element is
1516 // usable
1517 for (unsigned int base = 0; base < this->n_base_elements(); ++base)
1518 if (base_element(base).constraints_are_implemented() == false)
1519 return;
1520
1521 // TODO: the implementation makes the assumption that all faces have the
1522 // same number of dofs
1523 AssertDimension(this->n_unique_faces(), 1);
1524 const unsigned int face_no = 0;
1525
1526 this->interface_constraints.TableBase<2, double>::reinit(
1527 this->interface_constraints_size());
1528
1529 // the layout of the constraints matrix is described in the FiniteElement
1530 // class. you may want to look there first before trying to understand the
1531 // following, especially the mapping of the @p{m} index.
1532 //
1533 // in order to map it to the fe-system class, we have to know which base
1534 // element a degree of freedom within a vertex, line, etc belongs to. this
1535 // can be accomplished by the system_to_component_index function in
1536 // conjunction with the numbers first_{line,quad,...}_index
1537 for (unsigned int n = 0; n < this->interface_constraints.n(); ++n)
1538 for (unsigned int m = 0; m < this->interface_constraints.m(); ++m)
1539 {
1540 // for the pair (n,m) find out which base element they belong to and
1541 // the number therein
1542 //
1543 // first for the n index. this is simple since the n indices are in
1544 // the same order as they are usually on a face. note that for the
1545 // data type, first value in pair is (base element,instance of base
1546 // element), second is index within this instance
1547 const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
1548 n_index = this->face_system_to_base_table[face_no][n];
1549
1550 // likewise for the m index. this is more complicated due to the
1551 // strange ordering we have for the dofs on the refined faces.
1552 std::pair<std::pair<unsigned int, unsigned int>, unsigned int> m_index;
1553 switch (dim)
1554 {
1555 case 1:
1556 {
1557 // we should never get here! (in 1d, the constraints matrix
1558 // should be of size zero)
1560 break;
1561 }
1562
1563 case 2:
1564 {
1565 // the indices m=0..d_v-1 are from the center vertex. their
1566 // order is the same as for the first vertex of the whole cell,
1567 // so we can use the system_to_base_table variable (using the
1568 // face_s_t_base_t function would yield the same)
1569 if (m < this->n_dofs_per_vertex())
1570 m_index = this->system_to_base_table[m];
1571 else
1572 // then come the two sets of line indices
1573 {
1574 const unsigned int index_in_line =
1575 (m - this->n_dofs_per_vertex()) % this->n_dofs_per_line();
1576 const unsigned int sub_line =
1577 (m - this->n_dofs_per_vertex()) / this->n_dofs_per_line();
1578 Assert(sub_line < 2, ExcInternalError());
1579
1580 // from this information, try to get base element and
1581 // instance of base element. we do so by constructing the
1582 // corresponding face index of m in the present element,
1583 // then use face_system_to_base_table
1584 const unsigned int tmp1 =
1585 2 * this->n_dofs_per_vertex() + index_in_line;
1586 m_index.first =
1587 this->face_system_to_base_table[face_no][tmp1].first;
1588
1589 // what we are still missing is the index of m within the
1590 // base elements interface_constraints table
1591 //
1592 // here, the second value of face_system_to_base_table can
1593 // help: it denotes the face index of that shape function
1594 // within the base element. since we know that it is a line
1595 // dof, we can construct the rest: tmp2 will denote the
1596 // index of this shape function among the line shape
1597 // functions:
1598 Assert(
1599 this->face_system_to_base_table[face_no][tmp1].second >=
1600 2 *
1601 base_element(m_index.first.first).n_dofs_per_vertex(),
1603 const unsigned int tmp2 =
1604 this->face_system_to_base_table[face_no][tmp1].second -
1605 2 * base_element(m_index.first.first).n_dofs_per_vertex();
1606 Assert(tmp2 < base_element(m_index.first.first)
1607 .n_dofs_per_line(),
1609 m_index.second =
1610 base_element(m_index.first.first).n_dofs_per_vertex() +
1611 base_element(m_index.first.first).n_dofs_per_line() *
1612 sub_line +
1613 tmp2;
1614 }
1615 break;
1616 }
1617
1618 case 3:
1619 {
1620 // same way as above, although a little more complicated...
1621
1622 // the indices m=0..5*d_v-1 are from the center and the four
1623 // subline vertices. their order is the same as for the first
1624 // vertex of the whole cell, so we can use the simple arithmetic
1625 if (m < 5 * this->n_dofs_per_vertex())
1626 m_index = this->system_to_base_table[m];
1627 else
1628 // then come the 12 sets of line indices
1629 if (m < 5 * this->n_dofs_per_vertex() +
1630 12 * this->n_dofs_per_line())
1631 {
1632 // for the meaning of all this, see the 2d part
1633 const unsigned int index_in_line =
1634 (m - 5 * this->n_dofs_per_vertex()) %
1635 this->n_dofs_per_line();
1636 const unsigned int sub_line =
1637 (m - 5 * this->n_dofs_per_vertex()) /
1638 this->n_dofs_per_line();
1639 Assert(sub_line < 12, ExcInternalError());
1640
1641 const unsigned int tmp1 =
1642 4 * this->n_dofs_per_vertex() + index_in_line;
1643 m_index.first =
1644 this->face_system_to_base_table[face_no][tmp1].first;
1645
1646 Assert(
1647 this->face_system_to_base_table[face_no][tmp1].second >=
1648 4 * base_element(m_index.first.first)
1649 .n_dofs_per_vertex(),
1651 const unsigned int tmp2 =
1652 this->face_system_to_base_table[face_no][tmp1].second -
1653 4 *
1654 base_element(m_index.first.first).n_dofs_per_vertex();
1655 Assert(tmp2 < base_element(m_index.first.first)
1656 .n_dofs_per_line(),
1658 m_index.second =
1659 5 * base_element(m_index.first.first)
1660 .n_dofs_per_vertex() +
1661 base_element(m_index.first.first).n_dofs_per_line() *
1662 sub_line +
1663 tmp2;
1664 }
1665 else
1666 // on one of the four sub-quads
1667 {
1668 // for the meaning of all this, see the 2d part
1669 const unsigned int index_in_quad =
1670 (m - 5 * this->n_dofs_per_vertex() -
1671 12 * this->n_dofs_per_line()) %
1672 this->n_dofs_per_quad(face_no);
1673 Assert(index_in_quad < this->n_dofs_per_quad(face_no),
1675 const unsigned int sub_quad =
1676 ((m - 5 * this->n_dofs_per_vertex() -
1677 12 * this->n_dofs_per_line()) /
1678 this->n_dofs_per_quad(face_no));
1679 Assert(sub_quad < 4, ExcInternalError());
1680
1681 const unsigned int tmp1 = 4 * this->n_dofs_per_vertex() +
1682 4 * this->n_dofs_per_line() +
1683 index_in_quad;
1684 Assert(tmp1 <
1685 this->face_system_to_base_table[face_no].size(),
1687 m_index.first =
1688 this->face_system_to_base_table[face_no][tmp1].first;
1689
1690 Assert(
1691 this->face_system_to_base_table[face_no][tmp1].second >=
1692 4 * base_element(m_index.first.first)
1693 .n_dofs_per_vertex() +
1694 4 * base_element(m_index.first.first)
1695 .n_dofs_per_line(),
1697 const unsigned int tmp2 =
1698 this->face_system_to_base_table[face_no][tmp1].second -
1699 4 * base_element(m_index.first.first)
1700 .n_dofs_per_vertex() -
1701 4 * base_element(m_index.first.first).n_dofs_per_line();
1702 Assert(tmp2 < base_element(m_index.first.first)
1703 .n_dofs_per_quad(face_no),
1705 m_index.second =
1706 5 * base_element(m_index.first.first)
1707 .n_dofs_per_vertex() +
1708 12 *
1709 base_element(m_index.first.first).n_dofs_per_line() +
1710 base_element(m_index.first.first)
1711 .n_dofs_per_quad(face_no) *
1712 sub_quad +
1713 tmp2;
1714 }
1715
1716 break;
1717 }
1718
1719 default:
1721 }
1722
1723 // now that we gathered all information: use it to build the
1724 // matrix. note that if n and m belong to different base elements or
1725 // instances, then there definitely will be no coupling
1726 if (n_index.first == m_index.first)
1727 this->interface_constraints(m, n) =
1728 (base_element(n_index.first.first)
1729 .constraints()(m_index.second, n_index.second));
1730 }
1731}
1732
1733
1734
1735template <int dim, int spacedim>
1736void
1738 const std::vector<const FiniteElement<dim, spacedim> *> &fes,
1739 const std::vector<unsigned int> &multiplicities)
1740{
1741 Assert(fes.size() == multiplicities.size(),
1742 ExcDimensionMismatch(fes.size(), multiplicities.size()));
1743 Assert(fes.size() > 0,
1744 ExcMessage("Need to pass at least one finite element."));
1745 Assert(count_nonzeros(multiplicities) > 0,
1746 ExcMessage("You only passed FiniteElements with multiplicity 0."));
1747
1748 // Note that we need to skip every FE with multiplicity 0 in the following
1749 // block of code
1750
1751 this->base_to_block_indices.reinit(0, 0);
1752
1753 for (unsigned int i = 0; i < fes.size(); ++i)
1754 if (multiplicities[i] > 0)
1755 this->base_to_block_indices.push_back(multiplicities[i]);
1756
1757 {
1758 Threads::TaskGroup<> clone_base_elements;
1759
1760 unsigned int ind = 0;
1761 for (unsigned int i = 0; i < fes.size(); ++i)
1762 if (multiplicities[i] > 0)
1763 {
1764 clone_base_elements += Threads::new_task([&, i, ind]() {
1765 base_elements[ind] = {fes[i]->clone(), multiplicities[i]};
1766 });
1767 ++ind;
1768 }
1769 Assert(ind > 0, ExcInternalError());
1770
1771 // wait for all of these clone operations to finish
1772 clone_base_elements.join_all();
1773 }
1774
1775
1776 {
1777 // If the system is not primitive, these have not been initialized by
1778 // FiniteElement
1779 this->system_to_component_table.resize(this->n_dofs_per_cell());
1780
1781 FETools::Compositing::build_cell_tables(this->system_to_base_table,
1782 this->system_to_component_table,
1783 this->component_to_base_table,
1784 *this);
1785
1786 this->face_system_to_component_table.resize(this->n_unique_faces());
1787
1788 for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
1789 {
1790 this->face_system_to_component_table[face_no].resize(
1791 this->n_dofs_per_face(face_no));
1792
1794 this->face_system_to_base_table[face_no],
1795 this->face_system_to_component_table[face_no],
1796 *this,
1797 true,
1798 face_no);
1799 }
1800 }
1801
1802 // now initialize interface constraints, support points, and other tables.
1803 // (restriction and prolongation matrices are only built on demand.) do
1804 // this in parallel
1805
1806 Threads::TaskGroup<> init_tasks;
1807
1808 init_tasks +=
1809 Threads::new_task([&]() { this->build_interface_constraints(); });
1810
1811 init_tasks += Threads::new_task([&]() {
1812 // if one of the base elements has no support points, then it makes no sense
1813 // to define support points for the composed element, so return an empty
1814 // array to demonstrate that fact. Note that we ignore FE_Nothing in this
1815 // logic.
1816 for (unsigned int base_el = 0; base_el < this->n_base_elements(); ++base_el)
1817 if (!base_element(base_el).has_support_points() &&
1818 base_element(base_el).n_dofs_per_cell() != 0)
1819 {
1820 this->unit_support_points.resize(0);
1821 return;
1822 }
1823
1824 // generate unit support points from unit support points of sub elements
1825 this->unit_support_points.resize(this->n_dofs_per_cell());
1826
1827 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1828 {
1829 const unsigned int base = this->system_to_base_table[i].first.first,
1830 base_index = this->system_to_base_table[i].second;
1831 Assert(base < this->n_base_elements(), ExcInternalError());
1832 Assert(base_index < base_element(base).unit_support_points.size(),
1834 this->unit_support_points[i] =
1835 base_element(base).unit_support_points[base_index];
1836 }
1837 });
1838
1839 init_tasks += Threads::new_task([&]() {
1840 primitive_offset_tables.resize(this->n_base_elements());
1841
1842 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1843 if (base_element(base_no).is_primitive())
1844 primitive_offset_tables[base_no] =
1846 });
1847
1848 init_tasks += Threads::new_task([&]() {
1849 nonprimitive_offset_tables.resize(this->n_base_elements());
1850
1851 for (unsigned int base_no = 0; base_no < this->n_base_elements(); ++base_no)
1852 if (!base_element(base_no).is_primitive())
1853 nonprimitive_offset_tables[base_no] =
1855 });
1856
1857 // initialize face support points (for dim==2,3). same procedure as above
1858 if (dim > 1)
1859 init_tasks += Threads::new_task([&]() {
1860 for (unsigned int face_no = 0; face_no < this->n_unique_faces();
1861 ++face_no)
1862 {
1863 // if one of the base elements has no support points, then it makes
1864 // no sense to define support points for the composed element. In
1865 // that case, return an empty array to demonstrate that fact (note
1866 // that we ask whether the base element has no support points at
1867 // all, not only none on the face!)
1868 //
1869 // on the other hand, if there is an element that simply has no
1870 // degrees of freedom on the face at all, then we don't care whether
1871 // it has support points or not. this is, for example, the case for
1872 // the stable Stokes element Q(p)^dim \times DGP(p-1).
1873 bool flag_has_no_support_points = false;
1874
1875 for (unsigned int base_el = 0; base_el < this->n_base_elements();
1876 ++base_el)
1877 if (!base_element(base_el).has_support_points() &&
1878 (base_element(base_el).n_dofs_per_face(face_no) > 0))
1879 {
1880 this->unit_face_support_points[face_no].resize(0);
1881 flag_has_no_support_points = true;
1882 break;
1883 }
1884
1885
1886 if (flag_has_no_support_points)
1887 continue;
1888
1889 // generate unit face support points from unit support points of sub
1890 // elements
1891 this->unit_face_support_points[face_no].resize(
1892 this->n_dofs_per_face(face_no));
1893
1894 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
1895 {
1896 const unsigned int base_i =
1897 this->face_system_to_base_table[face_no][i].first.first;
1898 const unsigned int index_in_base =
1899 this->face_system_to_base_table[face_no][i].second;
1900
1901 Assert(
1902 index_in_base <
1903 base_element(base_i).unit_face_support_points[face_no].size(),
1905
1906 this->unit_face_support_points[face_no][i] =
1907 base_element(base_i)
1908 .unit_face_support_points[face_no][index_in_base];
1909 }
1910 }
1911 });
1912
1913 // Initialize generalized support points and an (internal) index table
1914 init_tasks += Threads::new_task([&]() {
1915 // Iterate over all base elements, extract a representative set of
1916 // _unique_ generalized support points and store the information how
1917 // generalized support points of base elements are mapped to this list
1918 // of representatives. Complexity O(n^2), where n is the number of
1919 // generalized support points.
1920
1921 generalized_support_points_index_table.resize(this->n_base_elements());
1922
1923 for (unsigned int base = 0; base < this->n_base_elements(); ++base)
1924 {
1925 // If the current base element does not have generalized support
1926 // points, ignore it. Note that
1927 // * FESystem::convert_generalized_support_point_values_to_dof_values
1928 // will simply skip such non-interpolatory base elements by
1929 // assigning NaN to all dofs.
1930 // * If this routine does not pick up any generalized support
1931 // points the corresponding vector will be empty and
1932 // FiniteElement::has_generalized_support_points will return
1933 // false.
1934 if (!base_element(base).has_generalized_support_points())
1935 continue;
1936
1937 for (const auto &point :
1938 base_element(base).get_generalized_support_points())
1939 {
1940 // Is point already an element of generalized_support_points?
1941 const auto p =
1942 std::find(std::begin(this->generalized_support_points),
1943 std::end(this->generalized_support_points),
1944 point);
1945
1946 if (p == std::end(this->generalized_support_points))
1947 {
1948 // If no, update the table and add the point to the vector
1949 const auto n = this->generalized_support_points.size();
1950 generalized_support_points_index_table[base].push_back(n);
1951 this->generalized_support_points.push_back(point);
1952 }
1953 else
1954 {
1955 // If yes, just add the correct index to the table.
1956 const auto n = p - std::begin(this->generalized_support_points);
1957 generalized_support_points_index_table[base].push_back(n);
1958 }
1959 }
1960 }
1961
1962# ifdef DEBUG
1963 // check generalized_support_points_index_table for consistency
1964 for (unsigned int i = 0; i < base_elements.size(); ++i)
1965 {
1966 if (!base_element(i).has_generalized_support_points())
1967 continue;
1968
1969 const auto &points =
1970 base_elements[i].first->get_generalized_support_points();
1971 for (unsigned int j = 0; j < points.size(); ++j)
1972 {
1973 const auto n = generalized_support_points_index_table[i][j];
1974 Assert(this->generalized_support_points[n] == points[j],
1976 }
1977 }
1978# endif /* DEBUG */
1979 });
1980
1981 // initialize quad dof index permutation in 3d and higher
1982 if (dim >= 3)
1983 init_tasks += Threads::new_task([&]() {
1984 for (unsigned int face_no = 0; face_no < this->n_unique_faces();
1985 ++face_no)
1986 {
1987 // the array into which we want to write should have the correct size
1988 // already.
1989 Assert(this->adjust_quad_dof_index_for_face_orientation_table[face_no]
1990 .n_elements() ==
1991 this->reference_cell().n_face_orientations(face_no) *
1992 this->n_dofs_per_quad(face_no),
1994
1995 // to obtain the shifts for this composed element, copy the shift
1996 // information of the base elements
1997 unsigned int index = 0;
1998 for (unsigned int b = 0; b < this->n_base_elements(); ++b)
1999 {
2000 const Table<2, int> &temp =
2001 this->base_element(b)
2002 .adjust_quad_dof_index_for_face_orientation_table[face_no];
2003 for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
2004 {
2005 for (unsigned int i = 0; i < temp.size(0); ++i)
2006 for (unsigned int j = 0;
2007 j <
2008 this->reference_cell().n_face_orientations(face_no);
2009 ++j)
2010 this->adjust_quad_dof_index_for_face_orientation_table
2011 [face_no](index + i, j) = temp(i, j);
2012 index += temp.size(0);
2013 }
2014 }
2015 Assert(index == this->n_dofs_per_quad(face_no), ExcInternalError());
2016 }
2017
2018 // additionally compose the permutation information for lines
2019 Assert(this->adjust_line_dof_index_for_line_orientation_table.size() ==
2020 this->n_dofs_per_line(),
2022 unsigned int index = 0;
2023 for (unsigned int b = 0; b < this->n_base_elements(); ++b)
2024 {
2025 const std::vector<int> &temp2 =
2026 this->base_element(b)
2027 .adjust_line_dof_index_for_line_orientation_table;
2028 for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
2029 {
2030 std::copy(
2031 temp2.begin(),
2032 temp2.end(),
2033 this->adjust_line_dof_index_for_line_orientation_table.begin() +
2034 index);
2035 index += temp2.size();
2036 }
2037 }
2038 Assert(index == this->n_dofs_per_line(), ExcInternalError());
2039 });
2040
2041 // wait for all of this to finish
2042 init_tasks.join_all();
2043}
2044
2045
2046
2047template <int dim, int spacedim>
2048bool
2050{
2051 for (unsigned int b = 0; b < this->n_base_elements(); ++b)
2052 if (base_element(b).hp_constraints_are_implemented() == false)
2053 return false;
2054
2055 return true;
2056}
2057
2058
2059
2060template <int dim, int spacedim>
2061void
2063 const FiniteElement<dim, spacedim> &x_source_fe,
2064 FullMatrix<double> &interpolation_matrix,
2065 const unsigned int face_no) const
2066{
2067 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
2068 ExcDimensionMismatch(interpolation_matrix.n(),
2069 this->n_dofs_per_face(face_no)));
2070 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
2071 ExcDimensionMismatch(interpolation_matrix.m(),
2072 x_source_fe.n_dofs_per_face(face_no)));
2073
2074 // since dofs for each base are independent, we only have to stack things up
2075 // from base element to base element
2076 //
2077 // the problem is that we have to work with two FEs (this and
2078 // fe_other). only deal with the case that both are FESystems and that they
2079 // both have the same number of bases (counting multiplicity) each of which
2080 // match in their number of components. this covers
2081 // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2082 // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2083 // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2084 if (const auto *fe_other_system =
2085 dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe))
2086 {
2087 // clear matrix, since we will not get to set all elements
2088 interpolation_matrix = 0;
2089
2090 // loop over all the base elements of this and the other element, counting
2091 // their multiplicities
2092 unsigned int base_index = 0, base_index_other = 0;
2093 unsigned int multiplicity = 0, multiplicity_other = 0;
2094
2095 FullMatrix<double> base_to_base_interpolation;
2096
2097 while (true)
2098 {
2099 const FiniteElement<dim, spacedim> &base = base_element(base_index),
2100 &base_other =
2101 fe_other_system->base_element(
2102 base_index_other);
2103
2104 Assert(base.n_components() == base_other.n_components(),
2106
2107 // get the interpolation from the bases
2108 base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no),
2109 base.n_dofs_per_face(face_no));
2110 base.get_face_interpolation_matrix(base_other,
2111 base_to_base_interpolation,
2112 face_no);
2113
2114 // now translate entries. we'd like to have something like
2115 // face_base_to_system_index, but that doesn't exist. rather, all we
2116 // have is the reverse. well, use that then
2117 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
2118 if (this->face_system_to_base_index(i, face_no).first ==
2119 std::make_pair(base_index, multiplicity))
2120 for (unsigned int j = 0;
2121 j < fe_other_system->n_dofs_per_face(face_no);
2122 ++j)
2123 if (fe_other_system->face_system_to_base_index(j, face_no)
2124 .first ==
2125 std::make_pair(base_index_other, multiplicity_other))
2126 interpolation_matrix(j, i) = base_to_base_interpolation(
2127 fe_other_system->face_system_to_base_index(j, face_no)
2128 .second,
2129 this->face_system_to_base_index(i, face_no).second);
2130
2131 // advance to the next base element for this and the other fe_system;
2132 // see if we can simply advance the multiplicity by one, or if have to
2133 // move on to the next base element
2134 ++multiplicity;
2135 if (multiplicity == this->element_multiplicity(base_index))
2136 {
2137 multiplicity = 0;
2138 ++base_index;
2139 }
2140 ++multiplicity_other;
2141 if (multiplicity_other ==
2142 fe_other_system->element_multiplicity(base_index_other))
2143 {
2144 multiplicity_other = 0;
2145 ++base_index_other;
2146 }
2147
2148 // see if we have reached the end of the present element. if so, we
2149 // should have reached the end of the other one as well
2150 if (base_index == this->n_base_elements())
2151 {
2152 Assert(base_index_other == fe_other_system->n_base_elements(),
2154 break;
2155 }
2156
2157 // if we haven't reached the end of this element, we shouldn't have
2158 // reached the end of the other one either
2159 Assert(base_index_other != fe_other_system->n_base_elements(),
2161 }
2162 }
2163 else
2164 {
2165 // repeat the cast to make the exception message more useful
2167 (dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) !=
2168 nullptr),
2169 (typename FiniteElement<dim,
2170 spacedim>::ExcInterpolationNotImplemented()));
2171 }
2172}
2173
2174
2175
2176template <int dim, int spacedim>
2177void
2179 const FiniteElement<dim, spacedim> &x_source_fe,
2180 const unsigned int subface,
2181 FullMatrix<double> &interpolation_matrix,
2182 const unsigned int face_no) const
2183{
2185 (x_source_fe.get_name().find("FESystem<") == 0) ||
2186 (dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) != nullptr),
2188
2189 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
2190 ExcDimensionMismatch(interpolation_matrix.n(),
2191 this->n_dofs_per_face(face_no)));
2192 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
2193 ExcDimensionMismatch(interpolation_matrix.m(),
2194 x_source_fe.n_dofs_per_face(face_no)));
2195
2196 // since dofs for each base are independent, we only have to stack things up
2197 // from base element to base element
2198 //
2199 // the problem is that we have to work with two FEs (this and
2200 // fe_other). only deal with the case that both are FESystems and that they
2201 // both have the same number of bases (counting multiplicity) each of which
2202 // match in their number of components. this covers
2203 // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2204 // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2205 // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2206 const FESystem<dim, spacedim> *fe_other_system =
2207 dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe);
2208 if (fe_other_system != nullptr)
2209 {
2210 // clear matrix, since we will not get to set all elements
2211 interpolation_matrix = 0;
2212
2213 // loop over all the base elements of this and the other element, counting
2214 // their multiplicities
2215 unsigned int base_index = 0, base_index_other = 0;
2216 unsigned int multiplicity = 0, multiplicity_other = 0;
2217
2218 FullMatrix<double> base_to_base_interpolation;
2219
2220 while (true)
2221 {
2222 const FiniteElement<dim, spacedim> &base = base_element(base_index),
2223 &base_other =
2224 fe_other_system->base_element(
2225 base_index_other);
2226
2227 Assert(base.n_components() == base_other.n_components(),
2229
2230 // get the interpolation from the bases
2231 base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no),
2232 base.n_dofs_per_face(face_no));
2233 base.get_subface_interpolation_matrix(base_other,
2234 subface,
2235 base_to_base_interpolation,
2236 face_no);
2237
2238 // now translate entries. we'd like to have something like
2239 // face_base_to_system_index, but that doesn't exist. rather, all we
2240 // have is the reverse. well, use that then
2241 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
2242 if (this->face_system_to_base_index(i, face_no).first ==
2243 std::make_pair(base_index, multiplicity))
2244 for (unsigned int j = 0;
2245 j < fe_other_system->n_dofs_per_face(face_no);
2246 ++j)
2247 if (fe_other_system->face_system_to_base_index(j, face_no)
2248 .first ==
2249 std::make_pair(base_index_other, multiplicity_other))
2250 interpolation_matrix(j, i) = base_to_base_interpolation(
2251 fe_other_system->face_system_to_base_index(j, face_no)
2252 .second,
2253 this->face_system_to_base_index(i, face_no).second);
2254
2255 // advance to the next base element for this and the other fe_system;
2256 // see if we can simply advance the multiplicity by one, or if have to
2257 // move on to the next base element
2258 ++multiplicity;
2259 if (multiplicity == this->element_multiplicity(base_index))
2260 {
2261 multiplicity = 0;
2262 ++base_index;
2263 }
2264 ++multiplicity_other;
2265 if (multiplicity_other ==
2266 fe_other_system->element_multiplicity(base_index_other))
2267 {
2268 multiplicity_other = 0;
2269 ++base_index_other;
2270 }
2271
2272 // see if we have reached the end of the present element. if so, we
2273 // should have reached the end of the other one as well
2274 if (base_index == this->n_base_elements())
2275 {
2276 Assert(base_index_other == fe_other_system->n_base_elements(),
2278 break;
2279 }
2280
2281 // if we haven't reached the end of this element, we shouldn't have
2282 // reached the end of the other one either
2283 Assert(base_index_other != fe_other_system->n_base_elements(),
2285 }
2286 }
2287 else
2288 {
2289 // we should have caught this at the start, but check again anyway
2290 Assert(
2291 fe_other_system != nullptr,
2292 (typename FiniteElement<dim,
2293 spacedim>::ExcInterpolationNotImplemented()));
2294 }
2295}
2296
2297
2298
2299template <int dim, int spacedim>
2300template <int structdim>
2301std::vector<std::pair<unsigned int, unsigned int>>
2303 const FiniteElement<dim, spacedim> &fe_other,
2304 const unsigned int face_no) const
2305{
2306 // since dofs on each subobject (vertex, line, ...) are ordered such that
2307 // first come all from the first base element all multiplicities, then
2308 // second base element all multiplicities, etc., we simply have to stack all
2309 // the identities after each other
2310 //
2311 // the problem is that we have to work with two FEs (this and
2312 // fe_other). only deal with the case that both are FESystems and that they
2313 // both have the same number of bases (counting multiplicity) each of which
2314 // match in their number of components. this covers
2315 // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not
2316 // FESystem(FE_Q(p),1,FE_Q(q),2) vs
2317 // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1)
2318 if (const FESystem<dim, spacedim> *fe_other_system =
2319 dynamic_cast<const FESystem<dim, spacedim> *>(&fe_other))
2320 {
2321 // loop over all the base elements of this and the other element,
2322 // counting their multiplicities
2323 unsigned int base_index = 0, base_index_other = 0;
2324 unsigned int multiplicity = 0, multiplicity_other = 0;
2325
2326 // we also need to keep track of the number of dofs already treated for
2327 // each of the elements
2328 unsigned int dof_offset = 0, dof_offset_other = 0;
2329
2330 std::vector<std::pair<unsigned int, unsigned int>> identities;
2331
2332 while (true)
2333 {
2334 const FiniteElement<dim, spacedim> &base = base_element(base_index),
2335 &base_other =
2336 fe_other_system->base_element(
2337 base_index_other);
2338
2339 Assert(base.n_components() == base_other.n_components(),
2341
2342 // now translate the identities returned by the base elements to the
2343 // indices of this system element
2344 std::vector<std::pair<unsigned int, unsigned int>> base_identities;
2345 switch (structdim)
2346 {
2347 case 0:
2348 base_identities = base.hp_vertex_dof_identities(base_other);
2349 break;
2350 case 1:
2351 base_identities = base.hp_line_dof_identities(base_other);
2352 break;
2353 case 2:
2354 base_identities =
2355 base.hp_quad_dof_identities(base_other, face_no);
2356 break;
2357 default:
2359 }
2360
2361 for (const auto &base_identity : base_identities)
2362 identities.emplace_back(base_identity.first + dof_offset,
2363 base_identity.second + dof_offset_other);
2364
2365 // record the dofs treated above as already taken care of
2366 dof_offset += base.template n_dofs_per_object<structdim>();
2367 dof_offset_other +=
2368 base_other.template n_dofs_per_object<structdim>();
2369
2370 // advance to the next base element for this and the other
2371 // fe_system; see if we can simply advance the multiplicity by one,
2372 // or if have to move on to the next base element
2373 ++multiplicity;
2374 if (multiplicity == this->element_multiplicity(base_index))
2375 {
2376 multiplicity = 0;
2377 ++base_index;
2378 }
2379 ++multiplicity_other;
2380 if (multiplicity_other ==
2381 fe_other_system->element_multiplicity(base_index_other))
2382 {
2383 multiplicity_other = 0;
2384 ++base_index_other;
2385 }
2386
2387 // see if we have reached the end of the present element. if so, we
2388 // should have reached the end of the other one as well
2389 if (base_index == this->n_base_elements())
2390 {
2391 Assert(base_index_other == fe_other_system->n_base_elements(),
2393 break;
2394 }
2395
2396 // if we haven't reached the end of this element, we shouldn't have
2397 // reached the end of the other one either
2398 Assert(base_index_other != fe_other_system->n_base_elements(),
2400 }
2401
2402 return identities;
2403 }
2404 else
2405 {
2407 return std::vector<std::pair<unsigned int, unsigned int>>();
2408 }
2409}
2410
2411
2412
2413template <int dim, int spacedim>
2414std::vector<std::pair<unsigned int, unsigned int>>
2416 const FiniteElement<dim, spacedim> &fe_other) const
2417{
2418 return hp_object_dof_identities<0>(fe_other);
2419}
2420
2421template <int dim, int spacedim>
2422std::vector<std::pair<unsigned int, unsigned int>>
2424 const FiniteElement<dim, spacedim> &fe_other) const
2425{
2426 return hp_object_dof_identities<1>(fe_other);
2427}
2428
2429
2430
2431template <int dim, int spacedim>
2432std::vector<std::pair<unsigned int, unsigned int>>
2434 const FiniteElement<dim, spacedim> &fe_other,
2435 const unsigned int face_no) const
2436{
2437 return hp_object_dof_identities<2>(fe_other, face_no);
2438}
2439
2440
2441
2442template <int dim, int spacedim>
2445 const FiniteElement<dim, spacedim> &fe_other,
2446 const unsigned int codim) const
2447{
2448 Assert(codim <= dim, ExcImpossibleInDim(dim));
2449
2450 // vertex/line/face/cell domination
2451 // --------------------------------
2452 // at present all we can do is to compare with other FESystems that have the
2453 // same number of components and bases
2454 if (const FESystem<dim, spacedim> *fe_sys_other =
2455 dynamic_cast<const FESystem<dim, spacedim> *>(&fe_other))
2456 {
2457 Assert(this->n_components() == fe_sys_other->n_components(),
2459 Assert(this->n_base_elements() == fe_sys_other->n_base_elements(),
2461
2464
2465 // loop over all base elements and do some sanity checks
2466 for (unsigned int b = 0; b < this->n_base_elements(); ++b)
2467 {
2468 Assert(this->base_element(b).n_components() ==
2469 fe_sys_other->base_element(b).n_components(),
2471 Assert(this->element_multiplicity(b) ==
2472 fe_sys_other->element_multiplicity(b),
2474
2475 // for this pair of base elements, check who dominates and combine
2476 // with previous result
2477 const FiniteElementDomination::Domination base_domination =
2478 (this->base_element(b).compare_for_domination(
2479 fe_sys_other->base_element(b), codim));
2480 domination = domination & base_domination;
2481 }
2482
2483 return domination;
2484 }
2485
2488}
2489
2490
2491
2492template <int dim, int spacedim>
2494FESystem<dim, spacedim>::base_element(const unsigned int index) const
2495{
2496 AssertIndexRange(index, base_elements.size());
2497 return *base_elements[index].first;
2498}
2499
2500
2501
2502template <int dim, int spacedim>
2503bool
2505 const unsigned int shape_index,
2506 const unsigned int face_index) const
2507{
2508 return (base_element(this->system_to_base_index(shape_index).first.first)
2509 .has_support_on_face(this->system_to_base_index(shape_index).second,
2510 face_index));
2511}
2512
2513
2514
2515template <int dim, int spacedim>
2517FESystem<dim, spacedim>::unit_support_point(const unsigned int index) const
2518{
2519 AssertIndexRange(index, this->n_dofs_per_cell());
2520 Assert((this->unit_support_points.size() == this->n_dofs_per_cell()) ||
2521 (this->unit_support_points.empty()),
2523
2524 // let's see whether we have the information pre-computed
2525 if (this->unit_support_points.size() != 0)
2526 return this->unit_support_points[index];
2527 else
2528 // no. ask the base element whether it would like to provide this
2529 // information
2530 return (base_element(this->system_to_base_index(index).first.first)
2531 .unit_support_point(this->system_to_base_index(index).second));
2532}
2533
2534
2535
2536template <int dim, int spacedim>
2537Point<dim - 1>
2539 const unsigned int index,
2540 const unsigned int face_no) const
2541{
2542 AssertIndexRange(index, this->n_dofs_per_face(face_no));
2543 Assert(
2544 (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
2545 .size() == this->n_dofs_per_face(face_no)) ||
2546 (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
2547 .empty()),
2548 (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints()));
2549
2550 // let's see whether we have the information pre-computed
2551 if (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
2552 .size() != 0)
2553 return this
2554 ->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
2555 [index];
2556 else
2557 // no. ask the base element whether it would like to provide this
2558 // information
2559 return (
2560 base_element(this->face_system_to_base_index(index, face_no).first.first)
2561 .unit_face_support_point(
2562 this->face_system_to_base_index(index, face_no).second, face_no));
2563}
2564
2565
2566
2567template <int dim, int spacedim>
2568std::pair<Table<2, bool>, std::vector<unsigned int>>
2570{
2571 // Note that this->n_components() is actually only an estimate of how many
2572 // constant modes we will need. There might be more than one such mode
2573 // (e.g. FE_Q_DG0).
2574 Table<2, bool> constant_modes(this->n_components(), this->n_dofs_per_cell());
2575 std::vector<unsigned int> components;
2576 for (unsigned int i = 0; i < base_elements.size(); ++i)
2577 {
2578 const std::pair<Table<2, bool>, std::vector<unsigned int>> base_table =
2579 base_elements[i].first->get_constant_modes();
2580 AssertDimension(base_table.first.n_rows(), base_table.second.size());
2581 const unsigned int element_multiplicity = this->element_multiplicity(i);
2582
2583 // there might be more than one constant mode for some scalar elements,
2584 // so make sure the table actually fits: Create a new table with more
2585 // rows
2586 const unsigned int comp = components.size();
2587 if (constant_modes.n_rows() <
2588 comp + base_table.first.n_rows() * element_multiplicity)
2589 {
2590 Table<2, bool> new_constant_modes(comp + base_table.first.n_rows() *
2591 element_multiplicity,
2592 constant_modes.n_cols());
2593 for (unsigned int r = 0; r < comp; ++r)
2594 for (unsigned int c = 0; c < this->n_dofs_per_cell(); ++c)
2595 new_constant_modes(r, c) = constant_modes(r, c);
2596
2597 constant_modes = std::move(new_constant_modes);
2598 }
2599
2600 // next, fill the constant modes from the individual components as well
2601 // as the component numbers corresponding to the constant mode rows
2602 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
2603 {
2604 std::pair<std::pair<unsigned int, unsigned int>, unsigned int> ind =
2605 this->system_to_base_index(k);
2606 if (ind.first.first == i)
2607 for (unsigned int c = 0; c < base_table.first.n_rows(); ++c)
2608 constant_modes(comp +
2609 ind.first.second * base_table.first.n_rows() + c,
2610 k) = base_table.first(c, ind.second);
2611 }
2612 for (unsigned int r = 0; r < element_multiplicity; ++r)
2613 for (const unsigned int c : base_table.second)
2614 components.push_back(
2615 comp + r * this->base_elements[i].first->n_components() + c);
2616 }
2617 AssertDimension(components.size(), constant_modes.n_rows());
2618 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
2619 components);
2620}
2621
2622
2623
2624template <int dim, int spacedim>
2625void
2627 const std::vector<Vector<double>> &point_values,
2628 std::vector<double> &dof_values) const
2629{
2630 Assert(this->has_generalized_support_points(),
2631 ExcMessage("The FESystem does not have generalized support points"));
2632
2633 AssertDimension(point_values.size(),
2634 this->get_generalized_support_points().size());
2635 AssertDimension(dof_values.size(), this->n_dofs_per_cell());
2636
2637 std::vector<double> base_dof_values;
2638 std::vector<Vector<double>> base_point_values;
2639
2640 // loop over all base elements (respecting multiplicity) and let them do
2641 // the work on their share of the input argument
2642
2643 unsigned int current_vector_component = 0;
2644 for (unsigned int base = 0; base < base_elements.size(); ++base)
2645 {
2646 // We need access to the base_element, its multiplicity, the
2647 // number of generalized support points (n_base_points) and the
2648 // number of components we're dealing with.
2649 const auto &base_element = this->base_element(base);
2650 const unsigned int multiplicity = this->element_multiplicity(base);
2651 const unsigned int n_base_dofs = base_element.n_dofs_per_cell();
2652 const unsigned int n_base_components = base_element.n_components();
2653
2654 // If the number of base degrees of freedom is zero, there is nothing
2655 // to do, skip the rest of the body in this case and continue with
2656 // the next element
2657 if (n_base_dofs == 0)
2658 {
2659 current_vector_component += multiplicity * n_base_components;
2660 continue;
2661 }
2662
2663 if (base_element.has_generalized_support_points())
2664 {
2665 const unsigned int n_base_points =
2666 base_element.get_generalized_support_points().size();
2667
2668 base_dof_values.resize(n_base_dofs);
2669 base_point_values.resize(n_base_points);
2670
2671 for (unsigned int m = 0; m < multiplicity;
2672 ++m, current_vector_component += n_base_components)
2673 {
2674 // populate base_point_values for a recursive call to
2675 // convert_generalized_support_point_values_to_dof_values
2676 for (unsigned int j = 0; j < base_point_values.size(); ++j)
2677 {
2678 base_point_values[j].reinit(n_base_components, false);
2679
2680 const auto n =
2681 generalized_support_points_index_table[base][j];
2682
2683 // we have to extract the correct slice out of the global
2684 // vector of values:
2685 const auto *const begin =
2686 std::begin(point_values[n]) + current_vector_component;
2687 const auto *const end = begin + n_base_components;
2688 std::copy(begin, end, std::begin(base_point_values[j]));
2689 }
2690
2691 base_element
2692 .convert_generalized_support_point_values_to_dof_values(
2693 base_point_values, base_dof_values);
2694
2695 // Finally put these dof values back into global dof values
2696 // vector.
2697
2698 // To do this, we could really use a base_to_system_index()
2699 // function, but that doesn't exist -- just do it by using the
2700 // reverse table -- the amount of work done here is not worth
2701 // trying to optimizing this.
2702 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2703 if (this->system_to_base_index(i).first ==
2704 std::make_pair(base, m))
2705 dof_values[i] =
2706 base_dof_values[this->system_to_base_index(i).second];
2707 } /*for*/
2708 }
2709 else
2710 {
2711 // If the base element is non-interpolatory, assign NaN to all
2712 // DoFs associated to it.
2713
2714 // To do this, we could really use a base_to_system_index()
2715 // function, but that doesn't exist -- just do it by using the
2716 // reverse table -- the amount of work done here is not worth
2717 // trying to optimizing this.
2718 for (unsigned int m = 0; m < multiplicity; ++m)
2719 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2720 if (this->system_to_base_index(i).first ==
2721 std::make_pair(base, m))
2722 dof_values[i] = std::numeric_limits<double>::signaling_NaN();
2723
2724 current_vector_component += multiplicity * n_base_components;
2725 }
2726 } /*for*/
2727}
2728
2729
2730
2731template <int dim, int spacedim>
2732std::size_t
2734{
2735 // neglect size of data stored in @p{base_elements} due to some problems
2736 // with the compiler. should be neglectable after all, considering the size
2737 // of the data of the subelements
2739 sizeof(base_elements));
2740 for (unsigned int i = 0; i < base_elements.size(); ++i)
2741 mem += MemoryConsumption::memory_consumption(*base_elements[i].first);
2742 return mem;
2743}
2744
2745#endif
2746
2747// explicit instantiations
2748#include "fe_system.inst"
2749
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:2521
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