deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20: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\}}\)
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 Assert(fe_degree % 2 == 1,
68 ExcMessage("FE_Hermite only supports odd polynomial degrees."));
69 return (fe_degree == 0) ? 0 : (fe_degree - 1) / 2;
70 }
71
72
73
74 inline std::vector<unsigned int>
75 get_hermite_dpo_vector(const unsigned int dim, const unsigned int regularity)
76 {
77 std::vector<unsigned int> result(dim + 1, 0);
78 result[0] = Utilities::pow(regularity + 1, dim);
79
80 return result;
81 }
82
83
84
85 /*
86 * Renumbering function. Function needs different levels of for loop nesting
87 * for different values of dim, so different definitions are used for
88 * simplicity.
89 */
90 template <int dim>
91 void
92 hermite_hierarchic_to_lexicographic_numbering(const unsigned int regularity,
93 std::vector<unsigned int> &h2l);
94
95
96
97 template <>
98 void
100 const unsigned int regularity,
101 std::vector<unsigned int> &h2l)
102 {
103 const unsigned int node_dofs_1d = regularity + 1;
104
105 AssertDimension(h2l.size(), 2 * node_dofs_1d);
106
107 // Assign DOFs at vertices
108 for (unsigned int di = 0; di < 2; ++di)
109 for (unsigned int i = 0; i < node_dofs_1d; ++i)
110 h2l[i + di * node_dofs_1d] = i + di * node_dofs_1d;
111 }
112
113
114
115 template <>
116 void
118 const unsigned int regularity,
119 std::vector<unsigned int> &h2l)
120 {
121 const unsigned int node_dofs_1d = regularity + 1;
122 const unsigned int dim_dofs_1d = 2 * node_dofs_1d;
123 unsigned int offset = 0;
124
125 AssertDimension(h2l.size(), dim_dofs_1d * dim_dofs_1d);
126
127 // Assign DOFs at vertices
128 for (unsigned int di = 0; di < 2; ++di)
129 for (unsigned int dj = 0; dj < 2; ++dj)
130 {
131 for (unsigned int i = 0; i < node_dofs_1d; ++i)
132 for (unsigned int j = 0; j < node_dofs_1d; ++j)
133 h2l[j + i * node_dofs_1d + offset] =
134 j + i * dim_dofs_1d + (dj + di * dim_dofs_1d) * node_dofs_1d;
135
136 offset += node_dofs_1d * node_dofs_1d;
137 }
138 }
139
140
141
142 template <>
143 void
145 const unsigned int regularity,
146 std::vector<unsigned int> &h2l)
147 {
148 const unsigned int node_dofs_1d = regularity + 1;
149 const unsigned int node_dofs_2d = node_dofs_1d * node_dofs_1d;
150
151 const unsigned int dim_dofs_1d = 2 * node_dofs_1d;
152 const unsigned int dim_dofs_2d = dim_dofs_1d * dim_dofs_1d;
153
154 unsigned int offset = 0;
155
156 AssertDimension(h2l.size(), dim_dofs_2d * dim_dofs_1d);
157
158 // Assign DOFs at nodes
159 for (unsigned int di = 0; di < 2; ++di)
160 for (unsigned int dj = 0; dj < 2; ++dj)
161 for (unsigned int dk = 0; dk < 2; ++dk)
162 {
163 for (unsigned int i = 0; i < node_dofs_1d; ++i)
164 for (unsigned int j = 0; j < node_dofs_1d; ++j)
165 for (unsigned int k = 0; k < node_dofs_1d; ++k)
166 h2l[k + j * node_dofs_1d + i * node_dofs_2d + offset] =
167 k + j * dim_dofs_1d + i * dim_dofs_2d +
168 node_dofs_1d * (dk + dj * dim_dofs_1d + di * dim_dofs_2d);
169
170 offset += node_dofs_1d * node_dofs_2d;
171 }
172 }
173
174
175
176 template <int dim>
177 inline std::vector<unsigned int>
179 {
180 const std::vector<unsigned int> dpo =
181 get_hermite_dpo_vector(dim, regularity);
182 const ::FiniteElementData<dim> face_data(dpo, 1, 2 * regularity + 1);
183 std::vector<unsigned int> renumbering(face_data.dofs_per_cell);
184
185 hermite_hierarchic_to_lexicographic_numbering<dim>(regularity, renumbering);
186
187 return renumbering;
188 }
189
190
191
192 template <int dim>
193 std::vector<unsigned int>
195 {
197 hermite_hierarchic_to_lexicographic_numbering<dim>(regularity));
198 }
199
200
201
202 template <int dim>
203 inline std::vector<unsigned int>
205 const unsigned int regularity)
206 {
207 (void)regularity;
208 if constexpr (dim > 1)
209 return hermite_lexicographic_to_hierarchic_numbering<dim - 1>(regularity);
210 else
211 return std::vector<unsigned int>();
212 }
213
214
215
216 template <int dim>
218 get_hermite_polynomials(const unsigned int fe_degree)
219 {
220 const unsigned int regularity = get_regularity_from_degree(fe_degree);
221
222 TensorProductPolynomials<dim> polynomial_basis(
224
225 std::vector<unsigned int> renumber =
226 internal::hermite_hierarchic_to_lexicographic_numbering<dim>(regularity);
227 polynomial_basis.set_numbering(renumber);
228
229 return polynomial_basis;
230 }
231
232
233
240 template <int xdim, int xspacedim = xdim, typename xNumber = double>
242 {
243 public:
244 template <int spacedim, typename Number>
245 void
247 Rescaler<1, spacedim, Number> & /*rescaler*/,
248 const FE_Hermite<1, spacedim> &fe_herm,
249 const typename Mapping<1, spacedim>::InternalDataBase &mapping_data,
250 Table<2, Number> &value_list)
251 {
252 double cell_extent = 1.0;
253
254 // Check mapping_data is associated with a compatible mapping class
255 if (dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
256 &mapping_data) != nullptr)
257 {
258 const typename MappingCartesian<1>::InternalData *data =
259 dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
260 &mapping_data);
261 cell_extent = data->cell_extents[0];
262 }
263 else
265
266 const unsigned int regularity = fe_herm.get_regularity();
267 const unsigned int n_dofs_per_cell = fe_herm.n_dofs_per_cell();
268 const unsigned int n_q_points_out = value_list.size(1);
269 (void)n_dofs_per_cell;
270
271 AssertDimension(value_list.size(0), n_dofs_per_cell);
272 AssertDimension(n_dofs_per_cell, 2 * regularity + 2);
273
274 std::vector<unsigned int> l2h =
275 hermite_lexicographic_to_hierarchic_numbering<1>(regularity);
276
277 for (unsigned int q = 0; q < n_q_points_out; ++q)
278 {
279 double factor_1 = 1.0;
280
281 for (unsigned int d1 = 0, d2 = regularity + 1; d2 < n_dofs_per_cell;
282 ++d1, ++d2)
283 {
284 /*
285 * d1 is used to count over indices on the left and d2 counts
286 * over indices on the right. These variables are used
287 * to avoid the need to loop over vertices.
288 */
289 value_list(l2h[d1], q) *= factor_1;
290 value_list(l2h[d2], q) *= factor_1;
291
292 factor_1 *= cell_extent;
293 }
294 }
295 }
296
297
298
299 template <int spacedim, typename Number>
300 void
302 Rescaler<2, spacedim, Number> & /*rescaler*/,
303 const FE_Hermite<2, spacedim> &fe_herm,
304 const typename Mapping<2, spacedim>::InternalDataBase &mapping_data,
305 Table<2, Number> &value_list)
306 {
307 Tensor<1, 2> cell_extents;
308
309 // Check mapping_data is associated with a compatible mapping class
310 if (dynamic_cast<const typename MappingCartesian<2>::InternalData *>(
311 &mapping_data) != nullptr)
312 {
313 const typename MappingCartesian<2>::InternalData *data =
314 dynamic_cast<const typename MappingCartesian<2>::InternalData *>(
315 &mapping_data);
316 cell_extents = data->cell_extents;
317 }
318 else
320
321 const unsigned int regularity = fe_herm.get_regularity();
322 const unsigned int n_dofs_per_cell = fe_herm.n_dofs_per_cell();
323 const unsigned int n_dofs_per_dim = 2 * regularity + 2;
324 const unsigned int n_q_points_out = value_list.size(1);
325 (void)n_dofs_per_cell;
326
327 AssertDimension(value_list.size(0), n_dofs_per_cell);
328 AssertDimension(n_dofs_per_dim * n_dofs_per_dim, n_dofs_per_cell);
329
330 std::vector<unsigned int> l2h =
331 hermite_lexicographic_to_hierarchic_numbering<2>(regularity);
332
333 AssertDimension(l2h.size(), n_dofs_per_cell);
334
335 for (unsigned int q = 0; q < n_q_points_out; ++q)
336 {
337 double factor_2 = 1.0;
338
339 for (unsigned int d3 = 0, d4 = regularity + 1; d4 < n_dofs_per_dim;
340 ++d3, ++d4)
341 {
342 double factor_1 = factor_2;
343
344 for (unsigned int d1 = 0, d2 = regularity + 1;
345 d2 < n_dofs_per_dim;
346 ++d1, ++d2)
347 {
348 /*
349 * d1 and d2 represent "left" and "right" in the
350 * x-direction, d3 and d4 represent "bottom" and "top"
351 * in the y-direction. As before, this is to avoid looping
352 * over vertices.
353 */
354 value_list(l2h[d1 + d3 * n_dofs_per_dim], q) *= factor_1;
355 value_list(l2h[d2 + d3 * n_dofs_per_dim], q) *= factor_1;
356 value_list(l2h[d1 + d4 * n_dofs_per_dim], q) *= factor_1;
357 value_list(l2h[d2 + d4 * n_dofs_per_dim], q) *= factor_1;
358
359 factor_1 *= cell_extents[0];
360 }
361
362 factor_2 *= cell_extents[1];
363 }
364 }
365 }
366
367
368
369 template <int spacedim, typename Number>
370 void
372 Rescaler<3, spacedim, Number> & /*rescaler*/,
373 const FE_Hermite<3, spacedim> &fe_herm,
374 const typename Mapping<3, spacedim>::InternalDataBase &mapping_data,
375 Table<2, Number> &value_list)
376 {
377 Tensor<1, 3> cell_extents;
378
379 // Check mapping_data is associated with a compatible mapping class
380 if (dynamic_cast<const typename MappingCartesian<3>::InternalData *>(
381 &mapping_data) != nullptr)
382 {
383 const typename MappingCartesian<3>::InternalData *data =
384 dynamic_cast<const typename MappingCartesian<3>::InternalData *>(
385 &mapping_data);
386 cell_extents = data->cell_extents;
387 }
388 else
390
391 const unsigned int regularity = fe_herm.get_regularity();
392 const unsigned int n_dofs_per_cell = fe_herm.n_dofs_per_cell();
393 const unsigned int n_dofs_per_dim = 2 * regularity + 2;
394 const unsigned int n_dofs_per_quad = n_dofs_per_dim * n_dofs_per_dim;
395 const unsigned int n_q_points_out = value_list.size(1);
396 (void)n_dofs_per_cell;
397
398 AssertDimension(value_list.size(0), n_dofs_per_cell);
399 AssertDimension(Utilities::pow(n_dofs_per_dim, 3), n_dofs_per_cell);
400
401 std::vector<unsigned int> l2h =
402 hermite_lexicographic_to_hierarchic_numbering<3>(regularity);
403
404 for (unsigned int q = 0; q < n_q_points_out; ++q)
405 {
406 double factor_3 = 1.0;
407
408 for (unsigned int d5 = 0, d6 = regularity + 1; d6 < n_dofs_per_dim;
409 ++d5, ++d6)
410 {
411 double factor_2 = factor_3;
412
413 for (unsigned int d3 = 0, d4 = regularity + 1;
414 d4 < n_dofs_per_dim;
415 ++d3, ++d4)
416 {
417 double factor_1 = factor_2;
418
419 for (unsigned int d1 = 0, d2 = regularity + 1;
420 d2 < n_dofs_per_dim;
421 ++d1, ++d2)
422 {
423 /*
424 * d1, d2: "left" and "right" (x-direction)
425 * d3, d4: "bottom" and "top" (y-direction)
426 * d5, d6: "down" and "up" (z-direction)
427 * This avoids looping over vertices
428 */
429 value_list(
430 l2h[d1 + d3 * n_dofs_per_dim + d5 * n_dofs_per_quad],
431 q) *= factor_1;
432 value_list(
433 l2h[d2 + d3 * n_dofs_per_dim + d5 * n_dofs_per_quad],
434 q) *= factor_1;
435 value_list(
436 l2h[d1 + d4 * n_dofs_per_dim + d5 * n_dofs_per_quad],
437 q) *= factor_1;
438 value_list(
439 l2h[d2 + d4 * n_dofs_per_dim + d5 * n_dofs_per_quad],
440 q) *= factor_1;
441 value_list(
442 l2h[d1 + d3 * n_dofs_per_dim + d6 * n_dofs_per_quad],
443 q) *= factor_1;
444 value_list(
445 l2h[d2 + d3 * n_dofs_per_dim + d6 * n_dofs_per_quad],
446 q) *= factor_1;
447 value_list(
448 l2h[d1 + d4 * n_dofs_per_dim + d6 * n_dofs_per_quad],
449 q) *= factor_1;
450 value_list(
451 l2h[d2 + d4 * n_dofs_per_dim + d6 * n_dofs_per_quad],
452 q) *= factor_1;
453
454 factor_1 *= cell_extents[0];
455 }
456
457 factor_2 *= cell_extents[1];
458 }
459
460 factor_3 *= cell_extents[2];
461 }
462 }
463 }
464 }; // class Rescaler
465} // namespace internal
466
467
468
469// Constructors
470template <int dim, int spacedim>
471FE_Hermite<dim, spacedim>::FE_Hermite(const unsigned int fe_degree)
472 : FE_Poly<dim, spacedim>(
473 internal::get_hermite_polynomials<dim>(fe_degree),
474 FiniteElementData<dim>(internal::get_hermite_dpo_vector(
475 dim,
476 internal::get_regularity_from_degree(fe_degree)),
477 1,
478 std::max(1U, fe_degree),
479 ((fe_degree > 2) ? FiniteElementData<dim>::H2 :
480 FiniteElementData<dim>::H1)),
481 std::vector<bool>(Utilities::pow(std::max(2U, fe_degree + 1), dim),
482 false),
483 std::vector<ComponentMask>(Utilities::pow(std::max(2U, fe_degree + 1),
484 dim),
485 ComponentMask(1, true)))
486 , regularity(internal::get_regularity_from_degree(fe_degree))
487{
488 Assert((fe_degree % 2 == 1),
490 "ERROR: The current implementation of Hermite interpolation "
491 "polynomials is only defined for odd polynomial degrees. Running "
492 "in release mode will use a polynomial degree of max(1,fe_degree-1) "
493 "to protect against unexpected internal bugs."));
494}
495
496
497
498template <int dim, int spacedim>
499std::string
501{
502 std::ostringstream name_buffer;
503 name_buffer << "FE_Hermite<" << Utilities::dim_string(dim, spacedim) << ">("
504 << this->degree << ")";
505 return name_buffer.str();
506}
507
508
509
510template <int dim, int spacedim>
511std::unique_ptr<FiniteElement<dim, spacedim>>
513{
514 return std::make_unique<FE_Hermite<dim, spacedim>>(*this);
515}
516
517
518
519template <int dim, int spacedim>
522{
526 out |= update_rescale; // since we need to rescale values, gradients, ...
527 return out;
528}
529
530
531
538template <int dim, int spacedim>
539std::vector<std::pair<unsigned int, unsigned int>>
541 const FiniteElement<dim, spacedim> &fe_other) const
542{
543 if (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other) != nullptr)
544 {
545 // there should be exactly one single DoF of FE_Q_Base at a vertex, and it
546 // should have an identical value to the first Hermite DoF
547 return {{0U, 0U}};
548 }
549 else if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) !=
550 nullptr)
551 {
552 // there should be exactly one single DoF of FE_Q_Base at a vertex, and it
553 // should have an identical value to the first Hermite DoF
554 return {{0U, 0U}};
555 }
556 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
557 {
558 // the FE_Nothing has no degrees of freedom, so there are no
559 // equivalencies to be recorded
560 return {};
561 }
562 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
563 {
564 // if the other element has no elements on faces at all,
565 // then it would be impossible to enforce any kind of
566 // continuity even if we knew exactly what kind of element
567 // we have -- simply because the other element declares
568 // that it is discontinuous because it has no DoFs on
569 // its faces. in that case, just state that we have no
570 // constraints to declare
571 return {};
572 }
573 else if (const FE_Hermite<dim, spacedim> *fe_herm_other =
574 dynamic_cast<const FE_Hermite<dim, spacedim> *>(&fe_other))
575 {
577 return {};
578 }
579 else
580 {
582 return {};
583 }
584}
585
586
587
593template <int dim, int spacedim>
594std::vector<std::pair<unsigned int, unsigned int>>
596 const FiniteElement<dim, spacedim> &fe_other) const
597{
598 (void)fe_other;
599 return {};
600}
601
602
603
608template <int dim, int spacedim>
609std::vector<std::pair<unsigned int, unsigned int>>
611 const FiniteElement<dim, spacedim> &fe_other,
612 const unsigned int face_no) const
613{
614 (void)fe_other;
615 (void)face_no;
616 return {};
617}
618
619
620
621/*
622 * The layout of this function is largely copied directly from FE_Q,
623 * however FE_Hermite can behave significantly differently in terms
624 * of domination due to how the function space is defined */
625template <int dim, int spacedim>
628 const FiniteElement<dim, spacedim> &fe_other,
629 const unsigned int codim) const
630{
631 Assert(codim <= dim, ExcImpossibleInDim(dim));
632
633 if (codim > 0)
634 if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
635 // there are no requirements between continuous and discontinuous elements
637
638
639 // vertex/line/face domination
640 // (if fe_other is not derived from FE_DGQ)
641 // & cell domination
642 // ----------------------------------------
643 if (const FE_Hermite<dim, spacedim> *fe_hermite_other =
644 dynamic_cast<const FE_Hermite<dim, spacedim> *>(&fe_other))
645 {
646 if (this->degree < fe_hermite_other->degree)
648 else if (this->degree == fe_hermite_other->degree)
650 else
652 }
653 if (const FE_Q<dim, spacedim> *fe_q_other =
654 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
655 {
656 if (fe_q_other->degree == 1)
657 {
658 if (this->degree == 1)
660 else
662 }
663 else if (this->degree <= fe_q_other->degree)
665 else
667 }
668 else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
669 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
670 {
671 if (fe_p_other->degree == 1)
672 {
673 if (this->degree == 1)
675 else
677 }
678 else if (this->degree <= fe_p_other->degree)
680 else
682 }
683 else if (const FE_WedgeP<dim, spacedim> *fe_wp_other =
684 dynamic_cast<const FE_WedgeP<dim, spacedim> *>(&fe_other))
685 {
686 if (fe_wp_other->degree == 1)
687 {
688 if (this->degree == 1)
690 else
692 }
693 else if (this->degree <= fe_wp_other->degree)
695 else
697 }
698 else if (const FE_PyramidP<dim, spacedim> *fe_pp_other =
699 dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other))
700 {
701 if (fe_pp_other->degree == 1)
702 {
703 if (this->degree == 1)
705 else
707 }
708 else if (this->degree <= fe_pp_other->degree)
710 else
712 }
713 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
714 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
715 {
716 if (fe_nothing->is_dominating())
718 else
719 // the FE_Nothing has no degrees of freedom and it is typically used
720 // in a context where we don't require any continuity along the
721 // interface
723 }
724
727}
728
729
730
731template <int dim, int spacedim>
732std::vector<unsigned int>
734{
735 return internal::hermite_lexicographic_to_hierarchic_numbering<dim>(
736 this->regularity);
737}
738
739
740
741template <int dim, int spacedim>
742void
745 const CellSimilarity::Similarity cell_similarity,
746 const Quadrature<dim> & /*quadrature*/,
747 const Mapping<dim, spacedim> &mapping,
748 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
749 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
750 spacedim>
751 & /*mapping_data*/,
752 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
754 spacedim>
755 &output_data) const
756{
757 // Convert data object to internal data for this class.
758 // Fails with an exception if that is not possible.
759 Assert(
760 (dynamic_cast<const typename FE_Hermite<dim, spacedim>::InternalData *>(
761 &fe_internal) != nullptr),
763 const typename FE_Hermite<dim, spacedim>::InternalData &fe_data =
764 static_cast<const typename FE_Hermite<dim, spacedim>::InternalData &>(
765 fe_internal);
766
767 const UpdateFlags flags(fe_data.update_each);
768
769 // Transform values, gradients and higher derivatives. Values also need to
770 // be rescaled according the the nodal derivative they correspond to.
771 if ((flags & update_values) &&
772 (cell_similarity != CellSimilarity::translation))
773 {
775 for (unsigned int i = 0; i < output_data.shape_values.size(0); ++i)
776 for (unsigned int q = 0; q < output_data.shape_values.size(1); ++q)
777 output_data.shape_values(i, q) = fe_data.shape_values(i, q);
778 shape_fix.rescale_fe_hermite_values(shape_fix,
779 *this,
780 mapping_internal,
781 output_data.shape_values);
782 }
783
784 if ((flags & update_gradients) &&
785 (cell_similarity != CellSimilarity::translation))
786 {
787 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
788 mapping.transform(make_array_view(fe_data.shape_gradients, k),
790 mapping_internal,
791 make_array_view(output_data.shape_gradients, k));
792
794 grad_fix.rescale_fe_hermite_values(grad_fix,
795 *this,
796 mapping_internal,
797 output_data.shape_gradients);
798 }
799
800 if ((flags & update_hessians) &&
801 (cell_similarity != CellSimilarity::translation))
802 {
803 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
804 mapping.transform(make_array_view(fe_data.shape_hessians, k),
806 mapping_internal,
807 make_array_view(output_data.shape_hessians, k));
808
810 hessian_fix.rescale_fe_hermite_values(hessian_fix,
811 *this,
812 mapping_internal,
813 output_data.shape_hessians);
814 }
815
816 if ((flags & update_3rd_derivatives) &&
817 (cell_similarity != CellSimilarity::translation))
818 {
819 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
822 mapping_internal,
823 make_array_view(output_data.shape_3rd_derivatives,
824 k));
825
827 third_dev_fix.rescale_fe_hermite_values(
828 third_dev_fix,
829 *this,
830 mapping_internal,
831 output_data.shape_3rd_derivatives);
832 }
833}
834
835
836
837template <int dim, int spacedim>
838void
841 const unsigned int face_no,
842 const hp::QCollection<dim - 1> &quadrature,
843 const Mapping<dim, spacedim> &mapping,
844 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
845 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
846 spacedim>
847 &,
848 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
850 spacedim>
851 &output_data) const
852{
853 /*
854 * Convert data object to internal data for this class. Fails with
855 * an exception if that is not possible.
856 */
857 Assert(
858 (dynamic_cast<const typename FE_Hermite<dim, spacedim>::InternalData *>(
859 &fe_internal) != nullptr),
861 const typename FE_Hermite<dim, spacedim>::InternalData &fe_data =
862 static_cast<const typename FE_Hermite<dim, spacedim>::InternalData &>(
863 fe_internal);
864
865 Assert((dynamic_cast<
867 &mapping_internal) != nullptr),
869
870 AssertDimension(quadrature.size(), 1U);
871
872 /*
873 * offset determines which data set to take (all data sets for all
874 * faces are stored contiguously)
875 */
876 const typename QProjector<dim>::DataSetDescriptor offset =
878 ReferenceCells::get_hypercube<dim>(),
879 face_no,
880 cell->combined_face_orientation(face_no),
881 quadrature[0].size());
882
883 const UpdateFlags flags(fe_data.update_each);
884
885 // Transform values, gradients and higher derivatives.
886 if (flags & update_values)
887 {
888 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
889 for (unsigned int i = 0; i < quadrature[0].size(); ++i)
890 output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset];
891
893 shape_face_fix.rescale_fe_hermite_values(shape_face_fix,
894 *this,
895 mapping_internal,
896 output_data.shape_values);
897 }
898
899 if (flags & update_gradients)
900 {
901 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
903 k,
904 offset,
905 quadrature[0].size()),
907 mapping_internal,
908 make_array_view(output_data.shape_gradients, k));
909
911 grad_face_fix.rescale_fe_hermite_values(grad_face_fix,
912 *this,
913 mapping_internal,
914 output_data.shape_gradients);
915 }
916
917 if (flags & update_hessians)
918 {
919 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
921 k,
922 offset,
923 quadrature[0].size()),
925 mapping_internal,
926 make_array_view(output_data.shape_hessians, k));
927
929 hessian_face_fix.rescale_fe_hermite_values(hessian_face_fix,
930 *this,
931 mapping_internal,
932 output_data.shape_hessians);
933 }
934
935 if (flags & update_3rd_derivatives)
936 {
937 for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
939 k,
940 offset,
941 quadrature[0].size()),
943 mapping_internal,
944 make_array_view(output_data.shape_3rd_derivatives,
945 k));
946
948 shape_3rd_face_fix.rescale_fe_hermite_values(
949 shape_3rd_face_fix,
950 *this,
951 mapping_internal,
952 output_data.shape_3rd_derivatives);
953 }
954}
955
956
957
958// Explicit instantiations
959#include "fe_hermite.inst"
960
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
Table< 2, double > shape_values
Definition fe_poly.h:438
Table< 2, Tensor< 3, dim > > shape_3rd_derivatives
Definition fe_poly.h:471
Table< 2, Tensor< 2, dim > > shape_hessians
Definition fe_poly.h:460
Table< 2, Tensor< 1, dim > > shape_gradients
Definition fe_poly.h:449
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:315
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#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()
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:75
void hermite_hierarchic_to_lexicographic_numbering< 1 >(const unsigned int regularity, std::vector< unsigned int > &h2l)
Definition fe_hermite.cc:99
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.