Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_hermite.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 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
15#include <deal.II/base/config.h>
16
19#include <deal.II/base/table.h>
21
23
24#include <deal.II/fe/fe_dgq.h>
25#include <deal.II/fe/fe_face.h>
29#include <deal.II/fe/fe_q.h>
32#include <deal.II/fe/fe_tools.h>
36
38
47#include <deal.II/lac/vector.h>
48
50
51#include <algorithm>
52#include <cmath>
53#include <iostream>
54#include <iterator>
55#include <memory>
56#include <sstream>
57
59
60
61
62namespace internal
63{
64 inline unsigned int
65 get_regularity_from_degree(const unsigned int fe_degree)
66 {
67 return (fe_degree == 0) ? 0 : (fe_degree - 1) / 2;
68 }
69
70
71
72 inline std::vector<unsigned int>
73 get_hermite_dpo_vector(const unsigned int dim, const unsigned int regularity)
74 {
75 std::vector<unsigned int> result(dim + 1, 0);
76 result[0] = Utilities::pow(regularity + 1, dim);
77
78 return result;
79 }
80
81
82
83 /*
84 * Renumbering function. Function needs different levels of for loop nesting
85 * for different values of dim, so different definitions are used for
86 * simplicity.
87 */
88 template <int dim>
89 void
90 hermite_hierarchic_to_lexicographic_numbering(const unsigned int regularity,
91 std::vector<unsigned int> &h2l);
92
93
94
95 template <>
96 void
98 const unsigned int regularity,
99 std::vector<unsigned int> &h2l)
100 {
101 const unsigned int node_dofs_1d = regularity + 1;
102
103 AssertDimension(h2l.size(), 2 * node_dofs_1d);
104
105 // Assign DOFs at vertices
106 for (unsigned int di = 0; di < 2; ++di)
107 for (unsigned int i = 0; i < node_dofs_1d; ++i)
108 h2l[i + di * node_dofs_1d] = i + di * node_dofs_1d;
109 }
110
111
112
113 template <>
114 void
116 const unsigned int regularity,
117 std::vector<unsigned int> &h2l)
118 {
119 const unsigned int node_dofs_1d = regularity + 1;
120 const unsigned int dim_dofs_1d = 2 * node_dofs_1d;
121 unsigned int offset = 0;
122
123 AssertDimension(h2l.size(), dim_dofs_1d * dim_dofs_1d);
124
125 // Assign DOFs at vertices
126 for (unsigned int di = 0; di < 2; ++di)
127 for (unsigned int dj = 0; dj < 2; ++dj)
128 {
129 for (unsigned int i = 0; i < node_dofs_1d; ++i)
130 for (unsigned int j = 0; j < node_dofs_1d; ++j)
131 h2l[j + i * node_dofs_1d + offset] =
132 j + i * dim_dofs_1d + (dj + di * dim_dofs_1d) * node_dofs_1d;
133
134 offset += node_dofs_1d * node_dofs_1d;
135 }
136 }
137
138
139
140 template <>
141 void
143 const unsigned int regularity,
144 std::vector<unsigned int> &h2l)
145 {
146 const unsigned int node_dofs_1d = regularity + 1;
147 const unsigned int node_dofs_2d = node_dofs_1d * node_dofs_1d;
148
149 const unsigned int dim_dofs_1d = 2 * node_dofs_1d;
150 const unsigned int dim_dofs_2d = dim_dofs_1d * dim_dofs_1d;
151
152 unsigned int offset = 0;
153
154 AssertDimension(h2l.size(), dim_dofs_2d * dim_dofs_1d);
155
156 // Assign DOFs at nodes
157 for (unsigned int di = 0; di < 2; ++di)
158 for (unsigned int dj = 0; dj < 2; ++dj)
159 for (unsigned int dk = 0; dk < 2; ++dk)
160 {
161 for (unsigned int i = 0; i < node_dofs_1d; ++i)
162 for (unsigned int j = 0; j < node_dofs_1d; ++j)
163 for (unsigned int k = 0; k < node_dofs_1d; ++k)
164 h2l[k + j * node_dofs_1d + i * node_dofs_2d + offset] =
165 k + j * dim_dofs_1d + i * dim_dofs_2d +
166 node_dofs_1d * (dk + dj * dim_dofs_1d + di * dim_dofs_2d);
167
168 offset += node_dofs_1d * node_dofs_2d;
169 }
170 }
171
172
173
174 template <int dim>
175 inline std::vector<unsigned int>
177 {
178 const std::vector<unsigned int> dpo =
179 get_hermite_dpo_vector(dim, regularity);
180 const ::FiniteElementData<dim> face_data(dpo, 1, 2 * regularity + 1);
181 std::vector<unsigned int> renumbering(face_data.dofs_per_cell);
182
184
185 return renumbering;
186 }
187
188
189
190 template <int dim>
191 std::vector<unsigned int>
197
198
199
200 template <int dim>
201 inline std::vector<unsigned int>
203 const unsigned int regularity)
204 {
205 (void)regularity;
206 if constexpr (dim > 1)
208 else
209 return std::vector<unsigned int>();
210 }
211
212
213
214 template <int dim>
216 get_hermite_polynomials(const unsigned int fe_degree)
217 {
218 const unsigned int regularity = get_regularity_from_degree(fe_degree);
219
220 TensorProductPolynomials<dim> polynomial_basis(
222
223 std::vector<unsigned int> renumber =
225 polynomial_basis.set_numbering(renumber);
226
227 return polynomial_basis;
228 }
229
230
231
238 template <int xdim, int xspacedim = xdim, typename xNumber = double>
240 {
241 public:
242 template <int spacedim, typename Number>
243 void
245 Rescaler<1, spacedim, Number> & /*rescaler*/,
246 const FE_Hermite<1, spacedim> &fe_herm,
247 const typename Mapping<1, spacedim>::InternalDataBase &mapping_data,
248 Table<2, Number> &value_list)
249 {
250 double cell_extent = 1.0;
251
252 // Check mapping_data is associated with a compatible mapping class
253 if (dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
254 &mapping_data) != nullptr)
255 {
256 const typename MappingCartesian<1>::InternalData *data =
257 dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
258 &mapping_data);
259 cell_extent = data->cell_extents[0];
260 }
261 else
263
264 const unsigned int regularity = fe_herm.get_regularity();
265 const unsigned int n_dofs_per_cell = fe_herm.n_dofs_per_cell();
266 const unsigned int n_q_points_out = value_list.size(1);
267 (void)n_dofs_per_cell;
268
269 AssertDimension(value_list.size(0), n_dofs_per_cell);
270 AssertDimension(n_dofs_per_cell, 2 * regularity + 2);
271
272 std::vector<unsigned int> l2h =
274
275 for (unsigned int q = 0; q < n_q_points_out; ++q)
276 {
277 double factor_1 = 1.0;
278
279 for (unsigned int d1 = 0, d2 = regularity + 1; d2 < n_dofs_per_cell;
280 ++d1, ++d2)
281 {
282 /*
283 * d1 is used to count over indices on the left and d2 counts
284 * over indices on the right. These variables are used
285 * to avoid the need to loop over vertices.
286 */
287 value_list(l2h[d1], q) *= factor_1;
288 value_list(l2h[d2], q) *= factor_1;
289
290 factor_1 *= cell_extent;
291 }
292 }
293 }
294
295
296
297 template <int spacedim, typename Number>
298 void
300 Rescaler<2, spacedim, Number> & /*rescaler*/,
301 const FE_Hermite<2, spacedim> &fe_herm,
302 const typename Mapping<2, spacedim>::InternalDataBase &mapping_data,
303 Table<2, Number> &value_list)
304 {
305 Tensor<1, 2> cell_extents;
306
307 // Check mapping_data is associated with a compatible mapping class
308 if (dynamic_cast<const typename MappingCartesian<2>::InternalData *>(
309 &mapping_data) != nullptr)
310 {
311 const typename MappingCartesian<2>::InternalData *data =
312 dynamic_cast<const typename MappingCartesian<2>::InternalData *>(
313 &mapping_data);
314 cell_extents = data->cell_extents;
315 }
316 else
318
319 const unsigned int regularity = fe_herm.get_regularity();
320 const unsigned int n_dofs_per_cell = fe_herm.n_dofs_per_cell();
321 const unsigned int n_dofs_per_dim = 2 * regularity + 2;
322 const unsigned int n_q_points_out = value_list.size(1);
323 (void)n_dofs_per_cell;
324
325 AssertDimension(value_list.size(0), n_dofs_per_cell);
326 AssertDimension(n_dofs_per_dim * n_dofs_per_dim, n_dofs_per_cell);
327
328 std::vector<unsigned int> l2h =
330
331 AssertDimension(l2h.size(), n_dofs_per_cell);
332
333 for (unsigned int q = 0; q < n_q_points_out; ++q)
334 {
335 double factor_2 = 1.0;
336
337 for (unsigned int d3 = 0, d4 = regularity + 1; d4 < n_dofs_per_dim;
338 ++d3, ++d4)
339 {
340 double factor_1 = factor_2;
341
342 for (unsigned int d1 = 0, d2 = regularity + 1;
343 d2 < n_dofs_per_dim;
344 ++d1, ++d2)
345 {
346 /*
347 * d1 and d2 represent "left" and "right" in the
348 * x-direction, d3 and d4 represent "bottom" and "top"
349 * in the y-direction. As before, this is to avoid looping
350 * over vertices.
351 */
352 value_list(l2h[d1 + d3 * n_dofs_per_dim], q) *= factor_1;
353 value_list(l2h[d2 + d3 * n_dofs_per_dim], q) *= factor_1;
354 value_list(l2h[d1 + d4 * n_dofs_per_dim], q) *= factor_1;
355 value_list(l2h[d2 + d4 * n_dofs_per_dim], q) *= factor_1;
356
357 factor_1 *= cell_extents[0];
358 }
359
360 factor_2 *= cell_extents[1];
361 }
362 }
363 }
364
365
366
367 template <int spacedim, typename Number>
368 void
370 Rescaler<3, spacedim, Number> & /*rescaler*/,
371 const FE_Hermite<3, spacedim> &fe_herm,
372 const typename Mapping<3, spacedim>::InternalDataBase &mapping_data,
373 Table<2, Number> &value_list)
374 {
375 Tensor<1, 3> cell_extents;
376
377 // Check mapping_data is associated with a compatible mapping class
378 if (dynamic_cast<const typename MappingCartesian<3>::InternalData *>(
379 &mapping_data) != nullptr)
380 {
381 const typename MappingCartesian<3>::InternalData *data =
382 dynamic_cast<const typename MappingCartesian<3>::InternalData *>(
383 &mapping_data);
384 cell_extents = data->cell_extents;
385 }
386 else
388
389 const unsigned int regularity = fe_herm.get_regularity();
390 const unsigned int n_dofs_per_cell = fe_herm.n_dofs_per_cell();
391 const unsigned int n_dofs_per_dim = 2 * regularity + 2;
392 const unsigned int n_dofs_per_quad = n_dofs_per_dim * n_dofs_per_dim;
393 const unsigned int n_q_points_out = value_list.size(1);
394 (void)n_dofs_per_cell;
395
396 AssertDimension(value_list.size(0), n_dofs_per_cell);
397 AssertDimension(Utilities::pow(n_dofs_per_dim, 3), n_dofs_per_cell);
398
399 std::vector<unsigned int> l2h =
401
402 for (unsigned int q = 0; q < n_q_points_out; ++q)
403 {
404 double factor_3 = 1.0;
405
406 for (unsigned int d5 = 0, d6 = regularity + 1; d6 < n_dofs_per_dim;
407 ++d5, ++d6)
408 {
409 double factor_2 = factor_3;
410
411 for (unsigned int d3 = 0, d4 = regularity + 1;
412 d4 < n_dofs_per_dim;
413 ++d3, ++d4)
414 {
415 double factor_1 = factor_2;
416
417 for (unsigned int d1 = 0, d2 = regularity + 1;
418 d2 < n_dofs_per_dim;
419 ++d1, ++d2)
420 {
421 /*
422 * d1, d2: "left" and "right" (x-direction)
423 * d3, d4: "bottom" and "top" (y-direction)
424 * d5, d6: "down" and "up" (z-direction)
425 * This avoids looping over vertices
426 */
427 value_list(
428 l2h[d1 + d3 * n_dofs_per_dim + d5 * n_dofs_per_quad],
429 q) *= factor_1;
430 value_list(
431 l2h[d2 + d3 * n_dofs_per_dim + d5 * n_dofs_per_quad],
432 q) *= factor_1;
433 value_list(
434 l2h[d1 + d4 * n_dofs_per_dim + d5 * n_dofs_per_quad],
435 q) *= factor_1;
436 value_list(
437 l2h[d2 + d4 * n_dofs_per_dim + d5 * n_dofs_per_quad],
438 q) *= factor_1;
439 value_list(
440 l2h[d1 + d3 * n_dofs_per_dim + d6 * n_dofs_per_quad],
441 q) *= factor_1;
442 value_list(
443 l2h[d2 + d3 * n_dofs_per_dim + d6 * n_dofs_per_quad],
444 q) *= factor_1;
445 value_list(
446 l2h[d1 + d4 * n_dofs_per_dim + d6 * n_dofs_per_quad],
447 q) *= factor_1;
448 value_list(
449 l2h[d2 + d4 * n_dofs_per_dim + d6 * n_dofs_per_quad],
450 q) *= factor_1;
451
452 factor_1 *= cell_extents[0];
453 }
454
455 factor_2 *= cell_extents[1];
456 }
457
458 factor_3 *= cell_extents[2];
459 }
460 }
461 }
462 }; // class Rescaler
463} // namespace internal
464
465
466
467// Constructors
468template <int dim, int spacedim>
469FE_Hermite<dim, spacedim>::FE_Hermite(const unsigned int fe_degree)
470 : FE_Poly<dim, spacedim>(
471 internal::get_hermite_polynomials<dim>(fe_degree),
472 FiniteElementData<dim>(internal::get_hermite_dpo_vector(
473 dim,
474 internal::get_regularity_from_degree(fe_degree)),
475 1,
476 std::max(1U, fe_degree),
477 ((fe_degree > 2) ? FiniteElementData<dim>::H2 :
478 FiniteElementData<dim>::H1)),
479 std::vector<bool>(Utilities::pow(std::max(2U, fe_degree + 1), dim),
480 false),
481 std::vector<ComponentMask>(Utilities::pow(std::max(2U, fe_degree + 1),
482 dim),
483 ComponentMask(1, true)))
484 , regularity(internal::get_regularity_from_degree(fe_degree))
485{
486 Assert((fe_degree % 2 == 1),
488 "ERROR: The current implementation of Hermite interpolation "
489 "polynomials is only defined for odd polynomial degrees. Running "
490 "in release mode will use a polynomial degree of max(1,fe_degree-1) "
491 "to protect against unexpected internal bugs."));
492}
493
494
495
496template <int dim, int spacedim>
497std::string
499{
500 std::ostringstream name_buffer;
501 name_buffer << "FE_Hermite<" << Utilities::dim_string(dim, spacedim) << ">("
502 << this->degree << ")";
503 return name_buffer.str();
504}
505
506
507
508template <int dim, int spacedim>
509std::unique_ptr<FiniteElement<dim, spacedim>>
511{
512 return std::make_unique<FE_Hermite<dim, spacedim>>(*this);
513}
514
515
516
517template <int dim, int spacedim>
520{
524 out |= update_rescale; // since we need to rescale values, gradients, ...
525 return out;
526}
527
528
529
536template <int dim, int spacedim>
537std::vector<std::pair<unsigned int, unsigned int>>
539 const FiniteElement<dim, spacedim> &fe_other) const
540{
541 if (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other) != nullptr)
542 {
543 // there should be exactly one single DoF of FE_Q_Base at a vertex, and it
544 // should have an identical value to the first Hermite DoF
545 return {{0U, 0U}};
546 }
547 else if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) !=
548 nullptr)
549 {
550 // there should be exactly one single DoF of FE_Q_Base at a vertex, and it
551 // should have an identical value to the first Hermite DoF
552 return {{0U, 0U}};
553 }
554 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
555 {
556 // the FE_Nothing has no degrees of freedom, so there are no
557 // equivalencies to be recorded
558 return {};
559 }
560 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
561 {
562 // if the other element has no elements on faces at all,
563 // then it would be impossible to enforce any kind of
564 // continuity even if we knew exactly what kind of element
565 // we have -- simply because the other element declares
566 // that it is discontinuous because it has no DoFs on
567 // its faces. in that case, just state that we have no
568 // constraints to declare
569 return {};
570 }
571 else if (const FE_Hermite<dim, spacedim> *fe_herm_other =
572 dynamic_cast<const FE_Hermite<dim, spacedim> *>(&fe_other))
573 {
575 return {};
576 }
577 else
578 {
580 return {};
581 }
582}
583
584
585
591template <int dim, int spacedim>
592std::vector<std::pair<unsigned int, unsigned int>>
594 const FiniteElement<dim, spacedim> &fe_other) const
595{
596 (void)fe_other;
597 return {};
598}
599
600
601
606template <int dim, int spacedim>
607std::vector<std::pair<unsigned int, unsigned int>>
609 const FiniteElement<dim, spacedim> &fe_other,
610 const unsigned int face_no) const
611{
612 (void)fe_other;
613 (void)face_no;
614 return {};
615}
616
617
618
619/*
620 * The layout of this function is largely copied directly from FE_Q,
621 * however FE_Hermite can behave significantly differently in terms
622 * of domination due to how the function space is defined */
623template <int dim, int spacedim>
626 const FiniteElement<dim, spacedim> &fe_other,
627 const unsigned int codim) const
628{
629 Assert(codim <= dim, ExcImpossibleInDim(dim));
630
631 if (codim > 0)
632 if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
633 // there are no requirements between continuous and discontinuous elements
635
636
637 // vertex/line/face domination
638 // (if fe_other is not derived from FE_DGQ)
639 // & cell domination
640 // ----------------------------------------
641 if (const FE_Hermite<dim, spacedim> *fe_hermite_other =
642 dynamic_cast<const FE_Hermite<dim, spacedim> *>(&fe_other))
643 {
644 if (this->degree < fe_hermite_other->degree)
646 else if (this->degree == fe_hermite_other->degree)
648 else
650 }
651 if (const FE_Q<dim, spacedim> *fe_q_other =
652 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
653 {
654 if (fe_q_other->degree == 1)
655 {
656 if (this->degree == 1)
658 else
660 }
661 else if (this->degree <= fe_q_other->degree)
663 else
665 }
666 else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
667 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
668 {
669 if (fe_p_other->degree == 1)
670 {
671 if (this->degree == 1)
673 else
675 }
676 else if (this->degree <= fe_p_other->degree)
678 else
680 }
681 else if (const FE_WedgeP<dim, spacedim> *fe_wp_other =
682 dynamic_cast<const FE_WedgeP<dim, spacedim> *>(&fe_other))
683 {
684 if (fe_wp_other->degree == 1)
685 {
686 if (this->degree == 1)
688 else
690 }
691 else if (this->degree <= fe_wp_other->degree)
693 else
695 }
696 else if (const FE_PyramidP<dim, spacedim> *fe_pp_other =
697 dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other))
698 {
699 if (fe_pp_other->degree == 1)
700 {
701 if (this->degree == 1)
703 else
705 }
706 else if (this->degree <= fe_pp_other->degree)
708 else
710 }
711 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
712 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
713 {
714 if (fe_nothing->is_dominating())
716 else
717 // the FE_Nothing has no degrees of freedom and it is typically used
718 // in a context where we don't require any continuity along the
719 // interface
721 }
722
725}
726
727
728
729template <int dim, int spacedim>
730std::vector<unsigned int>
736
737
738
739template <int dim, int spacedim>
740void
743 const CellSimilarity::Similarity cell_similarity,
744 const Quadrature<dim> & /*quadrature*/,
745 const Mapping<dim, spacedim> &mapping,
746 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
747 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
748 spacedim>
749 & /*mapping_data*/,
750 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
752 spacedim>
753 &output_data) const
754{
755 // Convert data object to internal data for this class.
756 // Fails with an exception if that is not possible.
757 Assert(
758 (dynamic_cast<const typename FE_Hermite<dim, spacedim>::InternalData *>(
759 &fe_internal) != nullptr),
761 const typename FE_Hermite<dim, spacedim>::InternalData &fe_data =
762 static_cast<const typename FE_Hermite<dim, spacedim>::InternalData &>(
763 fe_internal);
764
765 const UpdateFlags flags(fe_data.update_each);
766
767 // Transform values, gradients and higher derivatives. Values also need to
768 // be rescaled according the the nodal derivative they correspond to.
769 if ((flags & update_values) &&
770 (cell_similarity != CellSimilarity::translation))
771 {
773 for (unsigned int i = 0; i < output_data.shape_values.size(0); ++i)
774 for (unsigned int q = 0; q < output_data.shape_values.size(1); ++q)
775 output_data.shape_values(i, q) = fe_data.shape_values(i, q);
776 shape_fix.rescale_fe_hermite_values(shape_fix,
777 *this,
778 mapping_internal,
779 output_data.shape_values);
780 }
781
782 if ((flags & update_gradients) &&
783 (cell_similarity != CellSimilarity::translation))
784 {
785 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
786 mapping.transform(make_array_view(fe_data.shape_gradients, k),
788 mapping_internal,
789 make_array_view(output_data.shape_gradients, k));
790
792 grad_fix.rescale_fe_hermite_values(grad_fix,
793 *this,
794 mapping_internal,
795 output_data.shape_gradients);
796 }
797
798 if ((flags & update_hessians) &&
799 (cell_similarity != CellSimilarity::translation))
800 {
801 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
802 mapping.transform(make_array_view(fe_data.shape_hessians, k),
804 mapping_internal,
805 make_array_view(output_data.shape_hessians, k));
806
808 hessian_fix.rescale_fe_hermite_values(hessian_fix,
809 *this,
810 mapping_internal,
811 output_data.shape_hessians);
812 }
813
814 if ((flags & update_3rd_derivatives) &&
815 (cell_similarity != CellSimilarity::translation))
816 {
817 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
818 mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k),
820 mapping_internal,
821 make_array_view(output_data.shape_3rd_derivatives,
822 k));
823
825 third_dev_fix.rescale_fe_hermite_values(
826 third_dev_fix,
827 *this,
828 mapping_internal,
829 output_data.shape_3rd_derivatives);
830 }
831}
832
833
834
835template <int dim, int spacedim>
836void
839 const unsigned int face_no,
840 const hp::QCollection<dim - 1> &quadrature,
841 const Mapping<dim, spacedim> &mapping,
842 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
843 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
844 spacedim>
845 &,
846 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
848 spacedim>
849 &output_data) const
850{
851 /*
852 * Convert data object to internal data for this class. Fails with
853 * an exception if that is not possible.
854 */
855 Assert(
856 (dynamic_cast<const typename FE_Hermite<dim, spacedim>::InternalData *>(
857 &fe_internal) != nullptr),
859 const typename FE_Hermite<dim, spacedim>::InternalData &fe_data =
860 static_cast<const typename FE_Hermite<dim, spacedim>::InternalData &>(
861 fe_internal);
862
863 Assert((dynamic_cast<
865 &mapping_internal) != nullptr),
867
868 AssertDimension(quadrature.size(), 1U);
869
870 /*
871 * offset determines which data set to take (all data sets for all
872 * faces are stored contiguously)
873 */
874 const typename QProjector<dim>::DataSetDescriptor offset =
877 face_no,
878 cell->combined_face_orientation(face_no),
879 quadrature[0].size());
880
881 const UpdateFlags flags(fe_data.update_each);
882
883 // Transform values, gradients and higher derivatives.
884 if (flags & update_values)
885 {
886 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
887 for (unsigned int i = 0; i < quadrature[0].size(); ++i)
888 output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset];
889
891 shape_face_fix.rescale_fe_hermite_values(shape_face_fix,
892 *this,
893 mapping_internal,
894 output_data.shape_values);
895 }
896
897 if (flags & update_gradients)
898 {
899 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
900 mapping.transform(make_array_view(fe_data.shape_gradients,
901 k,
902 offset,
903 quadrature[0].size()),
905 mapping_internal,
906 make_array_view(output_data.shape_gradients, k));
907
909 grad_face_fix.rescale_fe_hermite_values(grad_face_fix,
910 *this,
911 mapping_internal,
912 output_data.shape_gradients);
913 }
914
915 if (flags & update_hessians)
916 {
917 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
918 mapping.transform(make_array_view(fe_data.shape_hessians,
919 k,
920 offset,
921 quadrature[0].size()),
923 mapping_internal,
924 make_array_view(output_data.shape_hessians, k));
925
927 hessian_face_fix.rescale_fe_hermite_values(hessian_face_fix,
928 *this,
929 mapping_internal,
930 output_data.shape_hessians);
931 }
932
933 if (flags & update_3rd_derivatives)
934 {
935 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
936 mapping.transform(make_array_view(fe_data.shape_3rd_derivatives,
937 k,
938 offset,
939 quadrature[0].size()),
941 mapping_internal,
942 make_array_view(output_data.shape_3rd_derivatives,
943 k));
944
946 shape_3rd_face_fix.rescale_fe_hermite_values(
947 shape_3rd_face_fix,
948 *this,
949 mapping_internal,
950 output_data.shape_3rd_derivatives);
951 }
952}
953
954
955
956// Explicit instantiations
957#include "fe_hermite.inst"
958
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
unsigned int get_regularity() const
Definition fe_hermite.h:271
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &other_fe, const unsigned int codim) 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 UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual std::string get_name() const override
std::vector< unsigned int > get_lexicographic_to_hierarchic_numbering() const
FE_Hermite(const unsigned int fe_degree)
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
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_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition fe_q.h:554
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_unique_faces() const
Abstract base class for mapping classes.
Definition mapping.h:318
virtual void transform(const ArrayView< const Tensor< 1, dim > > &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim > > &output) const =0
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int regularity)
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
void set_numbering(const std::vector< unsigned int > &renumber)
unsigned int size() const
Definition collection.h:308
void rescale_fe_hermite_values(Rescaler< 1, spacedim, Number > &, const FE_Hermite< 1, spacedim > &fe_herm, const typename Mapping< 1, spacedim >::InternalDataBase &mapping_data, Table< 2, Number > &value_list)
void rescale_fe_hermite_values(Rescaler< 2, spacedim, Number > &, const FE_Hermite< 2, spacedim > &fe_herm, const typename Mapping< 2, spacedim >::InternalDataBase &mapping_data, Table< 2, Number > &value_list)
void rescale_fe_hermite_values(Rescaler< 3, spacedim, Number > &, const FE_Hermite< 3, spacedim > &fe_herm, const typename Mapping< 3, spacedim >::InternalDataBase &mapping_data, Table< 2, Number > &value_list)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
UpdateFlags
@ update_rescale
@ 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.
@ mapping_covariant_gradient
Definition mapping.h:100
@ mapping_covariant
Definition mapping.h:89
@ mapping_covariant_hessian
Definition mapping.h:150
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
constexpr const ReferenceCell & get_hypercube()
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:555
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1699
std::vector< unsigned int > get_hermite_dpo_vector(const unsigned int dim, const unsigned int regularity)
Definition fe_hermite.cc:73
void hermite_hierarchic_to_lexicographic_numbering< 1 >(const unsigned int regularity, std::vector< unsigned int > &h2l)
Definition fe_hermite.cc:97
std::vector< unsigned int > hermite_face_lexicographic_to_hierarchic_numbering(const unsigned int regularity)
std::vector< unsigned int > hermite_lexicographic_to_hierarchic_numbering(const unsigned int regularity)
void hermite_hierarchic_to_lexicographic_numbering(const unsigned int regularity, std::vector< unsigned int > &h2l)
void hermite_hierarchic_to_lexicographic_numbering< 2 >(const unsigned int regularity, std::vector< unsigned int > &h2l)
void hermite_hierarchic_to_lexicographic_numbering< 3 >(const unsigned int regularity, std::vector< unsigned int > &h2l)
TensorProductPolynomials< dim > get_hermite_polynomials(const unsigned int fe_degree)
unsigned int get_regularity_from_degree(const unsigned int fe_degree)
Definition fe_hermite.cc:65
STL namespace.