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