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