Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20:02+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_poly_tensor.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 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
25
29
31#include <deal.II/grid/tria.h>
33
35
36namespace internal
37{
38 namespace FE_PolyTensor
39 {
40 namespace
41 {
42 template <int spacedim>
43 void
44 get_dof_sign_change_h_div(
45 const typename ::Triangulation<1, spacedim>::cell_iterator &,
47 const std::vector<MappingKind> &,
48 std::vector<double> &)
49 {
50 // Nothing to do in 1d.
51 }
52
53
54
55 // TODO: This function is not a consistent fix of the orientation issue
56 // like in 3d. It is rather kept not to break legacy behavior in 2d but
57 // should be replaced. See also the implementation of
58 // FE_RaviartThomas<dim>::initialize_quad_dof_index_permutation_and_sign_change()
59 // or other H(div) conforming elements such as FE_ABF<dim> and
60 // FE_BDM<dim>.
61 template <int spacedim>
62 void
63 get_dof_sign_change_h_div(
64 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
66 const std::vector<MappingKind> &mapping_kind,
67 std::vector<double> &face_sign)
68 {
69 const unsigned int dim = 2;
70 // const unsigned int spacedim = 2;
71
72 const CellId this_cell_id = cell->id();
73
74 for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
75 {
76 typename ::Triangulation<dim, spacedim>::face_iterator face =
77 cell->face(f);
78 if (!face->at_boundary())
79 {
80 const unsigned int nn = cell->neighbor_face_no(f);
81 const typename ::Triangulation<dim,
82 spacedim>::cell_iterator
83 neighbor_cell_at_face = cell->neighbor(f);
84 const CellId neigbor_cell_id = neighbor_cell_at_face->id();
85
86 // Only fix sign if the orientation is opposite and only do so
87 // on the face dofs on the cell with smaller cell_id.
88 if (((nn + f) % 2 == 0) && this_cell_id < neigbor_cell_id)
89 for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j)
90 {
91 const unsigned int cell_j = fe.face_to_cell_index(j, f);
92
93 Assert(f * fe.n_dofs_per_face(f) + j < face_sign.size(),
95 Assert(mapping_kind.size() == 1 ||
96 cell_j < mapping_kind.size(),
98
99 // TODO: This is probably only going to work for those
100 // elements for which all dofs are face dofs
101 if ((mapping_kind.size() > 1 ?
102 mapping_kind[cell_j] :
103 mapping_kind[0]) == mapping_raviart_thomas)
104 face_sign[f * fe.n_dofs_per_face(f) + j] = -1.0;
105 }
106 }
107 }
108 }
109
110
111
112 template <int spacedim>
113 void
114 get_dof_sign_change_h_div(
115 const typename ::Triangulation<3, spacedim>::cell_iterator
116 & /*cell*/,
117 const FiniteElement<3, spacedim> & /*fe*/,
118 const std::vector<MappingKind> & /*mapping_kind*/,
119 std::vector<double> & /*face_sign*/)
120 {
121 // Nothing to do. In 3d we take care of it through the
122 // adjust_quad_dof_sign_for_face_orientation_table
123 }
124
125 template <int spacedim>
126 void
127 get_dof_sign_change_nedelec(
128 const typename ::Triangulation<1, spacedim>::cell_iterator
129 & /*cell*/,
130 const FiniteElement<1, spacedim> & /*fe*/,
131 const std::vector<MappingKind> & /*mapping_kind*/,
132 std::vector<double> & /*line_dof_sign*/)
133 {
134 // nothing to do in 1d
135 }
136
137
138
139 template <int spacedim>
140 void
141 get_dof_sign_change_nedelec(
142 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
143 const FiniteElement<2, spacedim> & /*fe*/,
144 const std::vector<MappingKind> &mapping_kind,
145 std::vector<double> &line_dof_sign)
146 {
147 const unsigned int dim = 2;
148 // TODO: This fixes only lowest order
149 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
150 if (!(cell->line_orientation(l)) &&
151 mapping_kind[0] == mapping_nedelec)
152 line_dof_sign[l] = -1.0;
153 }
154
155
156 template <int spacedim>
157 void
158 get_dof_sign_change_nedelec(
159 const typename ::Triangulation<3, spacedim>::cell_iterator &cell,
160 const FiniteElement<3, spacedim> & /*fe*/,
161 const std::vector<MappingKind> &mapping_kind,
162 std::vector<double> &line_dof_sign)
163 {
164 const unsigned int dim = 3;
165 // TODO: This is probably only going to work for those elements for
166 // which all dofs are face dofs
167 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
168 if (!(cell->line_orientation(l)) &&
169 mapping_kind[0] == mapping_nedelec)
170 line_dof_sign[l] = -1.0;
171 }
172 } // namespace
173 } // namespace FE_PolyTensor
174} // namespace internal
175
176#ifndef DOXYGEN
177
178template <int dim, int spacedim>
180 const TensorPolynomialsBase<dim> &polynomials,
181 const FiniteElementData<dim> &fe_data,
182 const std::vector<bool> &restriction_is_additive_flags,
183 const std::vector<ComponentMask> &nonzero_components)
184 : FiniteElement<dim, spacedim>(fe_data,
185 restriction_is_additive_flags,
186 nonzero_components)
187 , mapping_kind({MappingKind::mapping_none})
188 , poly_space(polynomials.clone())
189{
190 cached_point[0] = -1;
191 // Set up the table converting
192 // components to base
193 // components. Since we have only
194 // one base element, everything
195 // remains zero except the
196 // component in the base, which is
197 // the component itself
198 for (unsigned int comp = 0; comp < this->n_components(); ++comp)
199 this->component_to_base_table[comp].first.second = comp;
200
201 if (dim == 3)
202 {
203 adjust_quad_dof_sign_for_face_orientation_table.resize(
204 this->n_unique_2d_subobjects());
205
206 for (unsigned int f = 0; f < this->n_unique_2d_subobjects(); ++f)
207 {
208 adjust_quad_dof_sign_for_face_orientation_table[f] =
209 Table<2, bool>(this->n_dofs_per_quad(f),
210 this->reference_cell().n_face_orientations(f));
211 adjust_quad_dof_sign_for_face_orientation_table[f].fill(false);
212 }
213 }
214}
215
216
217
218template <int dim, int spacedim>
220 : FiniteElement<dim, spacedim>(fe)
221 , mapping_kind(fe.mapping_kind)
222 , adjust_quad_dof_sign_for_face_orientation_table(
223 fe.adjust_quad_dof_sign_for_face_orientation_table)
224 , poly_space(fe.poly_space->clone())
225 , inverse_node_matrix(fe.inverse_node_matrix)
226{}
227
228
229
230template <int dim, int spacedim>
231bool
233{
234 return mapping_kind.size() == 1;
235}
236
237
238template <int dim, int spacedim>
239bool
241 const unsigned int index,
242 const unsigned int face,
243 const bool face_orientation,
244 const bool face_flip,
245 const bool face_rotation) const
246{
247 // do nothing in 1d and 2d
248 if (dim < 3)
249 return false;
250
251 // The exception are discontinuous
252 // elements for which there should be no
253 // face dofs anyway (i.e. dofs_per_quad==0
254 // in 3d), so we don't need the table, but
255 // the function should also not have been
256 // called
257 AssertIndexRange(index, this->n_dofs_per_quad(face));
258 Assert(adjust_quad_dof_sign_for_face_orientation_table
259 [this->n_unique_2d_subobjects() == 1 ? 0 : face]
260 .n_elements() ==
261 this->reference_cell().n_face_orientations(face) *
262 this->n_dofs_per_quad(face),
264
265 return adjust_quad_dof_sign_for_face_orientation_table
266 [this->n_unique_2d_subobjects() == 1 ? 0 : face](
267 index,
269 face_rotation,
270 face_flip));
271}
272
273
274template <int dim, int spacedim>
276FE_PolyTensor<dim, spacedim>::get_mapping_kind(const unsigned int i) const
277{
278 if (single_mapping_kind())
279 return mapping_kind[0];
280
281 AssertIndexRange(i, mapping_kind.size());
282 return mapping_kind[i];
283}
284
285
286
287template <int dim, int spacedim>
288double
290 const Point<dim> &) const
291
292{
294 return 0.;
295}
296
297
298
299template <int dim, int spacedim>
300double
302 const unsigned int i,
303 const Point<dim> &p,
304 const unsigned int component) const
305{
306 AssertIndexRange(i, this->n_dofs_per_cell());
307 AssertIndexRange(component, dim);
308
309 std::lock_guard<std::mutex> lock(cache_mutex);
310
311 if (cached_point != p || cached_values.empty())
312 {
313 cached_point = p;
314 cached_values.resize(poly_space->n());
315
316 std::vector<Tensor<4, dim>> dummy1;
317 std::vector<Tensor<5, dim>> dummy2;
318 poly_space->evaluate(
319 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
320 }
321
322 double s = 0;
323 if (inverse_node_matrix.n_cols() == 0)
324 return cached_values[i][component];
325 else
326 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
327 s += inverse_node_matrix(j, i) * cached_values[j][component];
328 return s;
329}
330
331
332
333template <int dim, int spacedim>
336 const Point<dim> &) const
337{
339 return Tensor<1, dim>();
340}
341
342
343
344template <int dim, int spacedim>
347 const unsigned int i,
348 const Point<dim> &p,
349 const unsigned int component) const
350{
351 AssertIndexRange(i, this->n_dofs_per_cell());
352 AssertIndexRange(component, dim);
353
354 std::lock_guard<std::mutex> lock(cache_mutex);
355
356 if (cached_point != p || cached_grads.empty())
357 {
358 cached_point = p;
359 cached_grads.resize(poly_space->n());
360
361 std::vector<Tensor<4, dim>> dummy1;
362 std::vector<Tensor<5, dim>> dummy2;
363 poly_space->evaluate(
364 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
365 }
366
368 if (inverse_node_matrix.n_cols() == 0)
369 return cached_grads[i][component];
370 else
371 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
372 s += inverse_node_matrix(j, i) * cached_grads[j][component];
373
374 return s;
375}
376
377
378
379template <int dim, int spacedim>
382 const Point<dim> &) const
383{
385 return Tensor<2, dim>();
386}
387
388
389
390template <int dim, int spacedim>
393 const unsigned int i,
394 const Point<dim> &p,
395 const unsigned int component) const
396{
397 AssertIndexRange(i, this->n_dofs_per_cell());
398 AssertIndexRange(component, dim);
399
400 std::lock_guard<std::mutex> lock(cache_mutex);
401
402 if (cached_point != p || cached_grad_grads.empty())
403 {
404 cached_point = p;
405 cached_grad_grads.resize(poly_space->n());
406
407 std::vector<Tensor<4, dim>> dummy1;
408 std::vector<Tensor<5, dim>> dummy2;
409 poly_space->evaluate(
410 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
411 }
412
414 if (inverse_node_matrix.n_cols() == 0)
415 return cached_grad_grads[i][component];
416 else
417 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
418 s += inverse_node_matrix(i, j) * cached_grad_grads[j][component];
419
420 return s;
421}
422
423
424//---------------------------------------------------------------------------
425// Fill data of FEValues
426//---------------------------------------------------------------------------
427
428template <int dim, int spacedim>
429void
432 const CellSimilarity::Similarity cell_similarity,
433 const Quadrature<dim> &quadrature,
434 const Mapping<dim, spacedim> &mapping,
435 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
437 &mapping_data,
438 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
440 spacedim>
441 &output_data) const
442{
443 // convert data object to internal
444 // data for this class. fails with
445 // an exception if that is not
446 // possible
447 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
449 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
450
451 const unsigned int n_q_points = quadrature.size();
452
453 Assert(!(fe_data.update_each & update_values) ||
454 fe_data.shape_values.size()[0] == this->n_dofs_per_cell(),
455 ExcDimensionMismatch(fe_data.shape_values.size()[0],
456 this->n_dofs_per_cell()));
457 Assert(!(fe_data.update_each & update_values) ||
458 fe_data.shape_values.size()[1] == n_q_points,
459 ExcDimensionMismatch(fe_data.shape_values.size()[1], n_q_points));
460
461 // TODO: The dof_sign_change only affects Nedelec elements and is not the
462 // correct thing on complicated meshes for higher order Nedelec elements.
463 // Something similar to FE_Q should be done to permute dofs and to change the
464 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
465 // class) is preferable.
466 std::fill(fe_data.dof_sign_change.begin(),
467 fe_data.dof_sign_change.end(),
468 1.0);
469 if (fe_data.update_each & update_values)
470 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
471 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
472
473 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
474 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
475 // Also nothing in 3d since we take care of it by using the
476 // adjust_quad_dof_sign_for_face_orientation_table.
477 if (fe_data.update_each & update_values)
478 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
479 *this,
480 this->mapping_kind,
481 fe_data.dof_sign_change);
482
483 // What is the first dof_index on a quad?
484 const unsigned int first_quad_index = this->get_first_quad_index();
485 // How many dofs per quad and how many quad dofs do we have at all?
486 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
487 const unsigned int n_quad_dofs =
488 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
489
490 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
491 ++dof_index)
492 {
493 /*
494 * This assumes that the dofs are ordered by first vertices, lines, quads
495 * and volume dofs. Note that in 2d this always gives false.
496 */
497 const bool is_quad_dof =
498 (dim == 2 ? false :
499 (first_quad_index <= dof_index) &&
500 (dof_index < first_quad_index + n_quad_dofs));
501
502 // TODO: This hack is not pretty and it is only here to handle the 2d
503 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
504 // handled by the
505 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
506 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
507 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
508 // it is handles with a table. This array is allocated in
509 // fe_poly_tensor.h.
510 double dof_sign = 1.0;
511 // under some circumstances fe_data.dof_sign_change is not allocated
512 if (fe_data.update_each & update_values)
513 dof_sign = fe_data.dof_sign_change[dof_index];
514
515 if (is_quad_dof)
516 {
517 /*
518 * Find the face belonging to this dof_index. This is integer
519 * division.
520 */
521 const unsigned int face_index_from_dof_index =
522 (dof_index - first_quad_index) / (n_dofs_per_quad);
523
524 const unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
525
526 // Correct the dof_sign if necessary
527 if (adjust_quad_dof_sign_for_face_orientation(
528 local_quad_dof_index,
529 face_index_from_dof_index,
530 cell->face_orientation(face_index_from_dof_index),
531 cell->face_flip(face_index_from_dof_index),
532 cell->face_rotation(face_index_from_dof_index)))
533 dof_sign = -1.0;
534 }
535
536 const MappingKind mapping_kind = get_mapping_kind(dof_index);
537
538 const unsigned int first =
539 output_data.shape_function_to_row_table
540 [dof_index * this->n_components() +
541 this->get_nonzero_components(dof_index).first_selected_component()];
542
543 // update the shape function values as necessary
544 //
545 // we only need to do this if the current cell is not a translation of
546 // the previous one; or, even if it is a translation, if we use
547 // mappings other than the standard mappings that require us to
548 // recompute values and derivatives because of possible sign changes
549 if (fe_data.update_each & update_values &&
550 ((cell_similarity != CellSimilarity::translation) ||
551 ((mapping_kind == mapping_piola) ||
552 (mapping_kind == mapping_raviart_thomas) ||
553 (mapping_kind == mapping_nedelec))))
554 {
555 switch (mapping_kind)
556 {
557 case mapping_none:
558 {
559 for (unsigned int k = 0; k < n_q_points; ++k)
560 for (unsigned int d = 0; d < dim; ++d)
561 output_data.shape_values(first + d, k) =
562 fe_data.shape_values[dof_index][k][d];
563 break;
564 }
565
568 {
569 mapping.transform(
570 make_array_view(fe_data.shape_values, dof_index),
571 mapping_kind,
572 mapping_internal,
573 make_array_view(fe_data.transformed_shape_values));
574
575 for (unsigned int k = 0; k < n_q_points; ++k)
576 for (unsigned int d = 0; d < dim; ++d)
577 output_data.shape_values(first + d, k) =
578 fe_data.transformed_shape_values[k][d];
579
580 break;
581 }
582
584 case mapping_piola:
585 {
586 mapping.transform(
587 make_array_view(fe_data.shape_values, dof_index),
589 mapping_internal,
590 make_array_view(fe_data.transformed_shape_values));
591 for (unsigned int k = 0; k < n_q_points; ++k)
592 for (unsigned int d = 0; d < dim; ++d)
593 output_data.shape_values(first + d, k) =
594 dof_sign * fe_data.transformed_shape_values[k][d];
595 break;
596 }
597
598 case mapping_nedelec:
599 {
600 mapping.transform(
601 make_array_view(fe_data.shape_values, dof_index),
603 mapping_internal,
604 make_array_view(fe_data.transformed_shape_values));
605
606 for (unsigned int k = 0; k < n_q_points; ++k)
607 for (unsigned int d = 0; d < dim; ++d)
608 output_data.shape_values(first + d, k) =
609 dof_sign * fe_data.transformed_shape_values[k][d];
610
611 break;
612 }
613
614 default:
616 }
617 }
618
619 // update gradients. apply the same logic as above
620 if (fe_data.update_each & update_gradients &&
621 ((cell_similarity != CellSimilarity::translation) ||
622 ((mapping_kind == mapping_piola) ||
623 (mapping_kind == mapping_raviart_thomas) ||
624 (mapping_kind == mapping_nedelec))))
625
626 {
627 switch (mapping_kind)
628 {
629 case mapping_none:
630 {
631 mapping.transform(
632 make_array_view(fe_data.shape_grads, dof_index),
634 mapping_internal,
635 make_array_view(fe_data.transformed_shape_grads));
636 for (unsigned int k = 0; k < n_q_points; ++k)
637 for (unsigned int d = 0; d < dim; ++d)
638 output_data.shape_gradients[first + d][k] =
639 fe_data.transformed_shape_grads[k][d];
640 break;
641 }
643 {
644 mapping.transform(
645 make_array_view(fe_data.shape_grads, dof_index),
647 mapping_internal,
648 make_array_view(fe_data.transformed_shape_grads));
649
650 for (unsigned int k = 0; k < n_q_points; ++k)
651 for (unsigned int d = 0; d < spacedim; ++d)
652 for (unsigned int n = 0; n < spacedim; ++n)
653 fe_data.transformed_shape_grads[k][d] -=
654 output_data.shape_values(first + n, k) *
655 mapping_data.jacobian_pushed_forward_grads[k][n][d];
656
657 for (unsigned int k = 0; k < n_q_points; ++k)
658 for (unsigned int d = 0; d < dim; ++d)
659 output_data.shape_gradients[first + d][k] =
660 fe_data.transformed_shape_grads[k][d];
661
662 break;
663 }
665 {
666 for (unsigned int k = 0; k < n_q_points; ++k)
667 fe_data.untransformed_shape_grads[k] =
668 fe_data.shape_grads[dof_index][k];
669 mapping.transform(
670 make_array_view(fe_data.untransformed_shape_grads),
672 mapping_internal,
673 make_array_view(fe_data.transformed_shape_grads));
674
675 for (unsigned int k = 0; k < n_q_points; ++k)
676 for (unsigned int d = 0; d < spacedim; ++d)
677 for (unsigned int n = 0; n < spacedim; ++n)
678 fe_data.transformed_shape_grads[k][d] +=
679 output_data.shape_values(first + n, k) *
680 mapping_data.jacobian_pushed_forward_grads[k][d][n];
681
682
683 for (unsigned int k = 0; k < n_q_points; ++k)
684 for (unsigned int d = 0; d < dim; ++d)
685 output_data.shape_gradients[first + d][k] =
686 fe_data.transformed_shape_grads[k][d];
687
688 break;
689 }
691 case mapping_piola:
692 {
693 for (unsigned int k = 0; k < n_q_points; ++k)
694 fe_data.untransformed_shape_grads[k] =
695 fe_data.shape_grads[dof_index][k];
696 mapping.transform(
697 make_array_view(fe_data.untransformed_shape_grads),
699 mapping_internal,
700 make_array_view(fe_data.transformed_shape_grads));
701
702 for (unsigned int k = 0; k < n_q_points; ++k)
703 for (unsigned int d = 0; d < spacedim; ++d)
704 for (unsigned int n = 0; n < spacedim; ++n)
705 fe_data.transformed_shape_grads[k][d] +=
706 (output_data.shape_values(first + n, k) *
707 mapping_data
709 (output_data.shape_values(first + d, k) *
710 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
711
712 for (unsigned int k = 0; k < n_q_points; ++k)
713 for (unsigned int d = 0; d < dim; ++d)
714 output_data.shape_gradients[first + d][k] =
715 dof_sign * fe_data.transformed_shape_grads[k][d];
716
717 break;
718 }
719
720 case mapping_nedelec:
721 {
722 // treat the gradients of
723 // this particular shape
724 // function at all
725 // q-points. if Dv is the
726 // gradient of the shape
727 // function on the unit
728 // cell, then
729 // (J^-T)Dv(J^-1) is the
730 // value we want to have on
731 // the real cell.
732 for (unsigned int k = 0; k < n_q_points; ++k)
733 fe_data.untransformed_shape_grads[k] =
734 fe_data.shape_grads[dof_index][k];
735
736 mapping.transform(
737 make_array_view(fe_data.untransformed_shape_grads),
739 mapping_internal,
740 make_array_view(fe_data.transformed_shape_grads));
741
742 for (unsigned int k = 0; k < n_q_points; ++k)
743 for (unsigned int d = 0; d < spacedim; ++d)
744 for (unsigned int n = 0; n < spacedim; ++n)
745 fe_data.transformed_shape_grads[k][d] -=
746 output_data.shape_values(first + n, k) *
747 mapping_data.jacobian_pushed_forward_grads[k][n][d];
748
749 for (unsigned int k = 0; k < n_q_points; ++k)
750 for (unsigned int d = 0; d < dim; ++d)
751 output_data.shape_gradients[first + d][k] =
752 dof_sign * fe_data.transformed_shape_grads[k][d];
753
754 break;
755 }
756
757 default:
759 }
760 }
761
762 // update hessians. apply the same logic as above
763 if (fe_data.update_each & update_hessians &&
764 ((cell_similarity != CellSimilarity::translation) ||
765 ((mapping_kind == mapping_piola) ||
766 (mapping_kind == mapping_raviart_thomas) ||
767 (mapping_kind == mapping_nedelec))))
768
769 {
770 switch (mapping_kind)
771 {
772 case mapping_none:
773 {
774 mapping.transform(
775 make_array_view(fe_data.shape_grad_grads, dof_index),
777 mapping_internal,
778 make_array_view(fe_data.transformed_shape_hessians));
779
780 for (unsigned int k = 0; k < n_q_points; ++k)
781 for (unsigned int d = 0; d < spacedim; ++d)
782 for (unsigned int n = 0; n < spacedim; ++n)
783 fe_data.transformed_shape_hessians[k][d] -=
784 output_data.shape_gradients[first + d][k][n] *
785 mapping_data.jacobian_pushed_forward_grads[k][n];
786
787 for (unsigned int k = 0; k < n_q_points; ++k)
788 for (unsigned int d = 0; d < dim; ++d)
789 output_data.shape_hessians[first + d][k] =
790 fe_data.transformed_shape_hessians[k][d];
791
792 break;
793 }
795 {
796 for (unsigned int k = 0; k < n_q_points; ++k)
797 fe_data.untransformed_shape_hessian_tensors[k] =
798 fe_data.shape_grad_grads[dof_index][k];
799
800 mapping.transform(
802 fe_data.untransformed_shape_hessian_tensors),
804 mapping_internal,
805 make_array_view(fe_data.transformed_shape_hessians));
806
807 for (unsigned int k = 0; k < n_q_points; ++k)
808 for (unsigned int d = 0; d < spacedim; ++d)
809 for (unsigned int n = 0; n < spacedim; ++n)
810 for (unsigned int i = 0; i < spacedim; ++i)
811 for (unsigned int j = 0; j < spacedim; ++j)
812 {
813 fe_data.transformed_shape_hessians[k][d][i][j] -=
814 (output_data.shape_values(first + n, k) *
815 mapping_data
817 [k][n][d][i][j]) +
818 (output_data.shape_gradients[first + d][k][n] *
819 mapping_data
820 .jacobian_pushed_forward_grads[k][n][i][j]) +
821 (output_data.shape_gradients[first + n][k][i] *
822 mapping_data
823 .jacobian_pushed_forward_grads[k][n][d][j]) +
824 (output_data.shape_gradients[first + n][k][j] *
825 mapping_data
826 .jacobian_pushed_forward_grads[k][n][i][d]);
827 }
828
829 for (unsigned int k = 0; k < n_q_points; ++k)
830 for (unsigned int d = 0; d < dim; ++d)
831 output_data.shape_hessians[first + d][k] =
832 fe_data.transformed_shape_hessians[k][d];
833
834 break;
835 }
837 {
838 for (unsigned int k = 0; k < n_q_points; ++k)
839 fe_data.untransformed_shape_hessian_tensors[k] =
840 fe_data.shape_grad_grads[dof_index][k];
841
842 mapping.transform(
844 fe_data.untransformed_shape_hessian_tensors),
846 mapping_internal,
847 make_array_view(fe_data.transformed_shape_hessians));
848
849 for (unsigned int k = 0; k < n_q_points; ++k)
850 for (unsigned int d = 0; d < spacedim; ++d)
851 for (unsigned int n = 0; n < spacedim; ++n)
852 for (unsigned int i = 0; i < spacedim; ++i)
853 for (unsigned int j = 0; j < spacedim; ++j)
854 {
855 fe_data.transformed_shape_hessians[k][d][i][j] +=
856 (output_data.shape_values(first + n, k) *
857 mapping_data
859 [k][d][n][i][j]) +
860 (output_data.shape_gradients[first + n][k][i] *
861 mapping_data
862 .jacobian_pushed_forward_grads[k][d][n][j]) +
863 (output_data.shape_gradients[first + n][k][j] *
864 mapping_data
865 .jacobian_pushed_forward_grads[k][d][i][n]) -
866 (output_data.shape_gradients[first + d][k][n] *
867 mapping_data
868 .jacobian_pushed_forward_grads[k][n][i][j]);
869 for (unsigned int m = 0; m < spacedim; ++m)
870 fe_data
871 .transformed_shape_hessians[k][d][i][j] -=
872 (mapping_data
873 .jacobian_pushed_forward_grads[k][d][i]
874 [m] *
875 mapping_data
876 .jacobian_pushed_forward_grads[k][m][n]
877 [j] *
878 output_data.shape_values(first + n, k)) +
879 (mapping_data
880 .jacobian_pushed_forward_grads[k][d][m]
881 [j] *
882 mapping_data
883 .jacobian_pushed_forward_grads[k][m][i]
884 [n] *
885 output_data.shape_values(first + n, k));
886 }
887
888 for (unsigned int k = 0; k < n_q_points; ++k)
889 for (unsigned int d = 0; d < dim; ++d)
890 output_data.shape_hessians[first + d][k] =
891 fe_data.transformed_shape_hessians[k][d];
892
893 break;
894 }
896 case mapping_piola:
897 {
898 for (unsigned int k = 0; k < n_q_points; ++k)
899 fe_data.untransformed_shape_hessian_tensors[k] =
900 fe_data.shape_grad_grads[dof_index][k];
901
902 mapping.transform(
904 fe_data.untransformed_shape_hessian_tensors),
906 mapping_internal,
907 make_array_view(fe_data.transformed_shape_hessians));
908
909 for (unsigned int k = 0; k < n_q_points; ++k)
910 for (unsigned int d = 0; d < spacedim; ++d)
911 for (unsigned int n = 0; n < spacedim; ++n)
912 for (unsigned int i = 0; i < spacedim; ++i)
913 for (unsigned int j = 0; j < spacedim; ++j)
914 {
915 fe_data.transformed_shape_hessians[k][d][i][j] +=
916 (output_data.shape_values(first + n, k) *
917 mapping_data
919 [k][d][n][i][j]) +
920 (output_data.shape_gradients[first + n][k][i] *
921 mapping_data
922 .jacobian_pushed_forward_grads[k][d][n][j]) +
923 (output_data.shape_gradients[first + n][k][j] *
924 mapping_data
925 .jacobian_pushed_forward_grads[k][d][i][n]) -
926 (output_data.shape_gradients[first + d][k][n] *
927 mapping_data
928 .jacobian_pushed_forward_grads[k][n][i][j]);
929
930 fe_data.transformed_shape_hessians[k][d][i][j] -=
931 (output_data.shape_values(first + d, k) *
932 mapping_data
934 [k][n][n][i][j]) +
935 (output_data.shape_gradients[first + d][k][i] *
936 mapping_data
937 .jacobian_pushed_forward_grads[k][n][n][j]) +
938 (output_data.shape_gradients[first + d][k][j] *
939 mapping_data
940 .jacobian_pushed_forward_grads[k][n][n][i]);
941
942 for (unsigned int m = 0; m < spacedim; ++m)
943 {
944 fe_data
945 .transformed_shape_hessians[k][d][i][j] -=
946 (mapping_data
948 [m] *
949 mapping_data
951 [j] *
952 output_data.shape_values(first + n, k)) +
953 (mapping_data
954 .jacobian_pushed_forward_grads[k][d][m]
955 [j] *
956 mapping_data
957 .jacobian_pushed_forward_grads[k][m][i]
958 [n] *
959 output_data.shape_values(first + n, k));
960
961 fe_data
962 .transformed_shape_hessians[k][d][i][j] +=
963 (mapping_data
965 [m] *
966 mapping_data
968 [j] *
969 output_data.shape_values(first + d, k)) +
970 (mapping_data
971 .jacobian_pushed_forward_grads[k][n][m]
972 [j] *
973 mapping_data
974 .jacobian_pushed_forward_grads[k][m][i]
975 [n] *
976 output_data.shape_values(first + d, k));
977 }
978 }
979
980 for (unsigned int k = 0; k < n_q_points; ++k)
981 for (unsigned int d = 0; d < dim; ++d)
982 output_data.shape_hessians[first + d][k] =
983 dof_sign * fe_data.transformed_shape_hessians[k][d];
984
985 break;
986 }
987
988 case mapping_nedelec:
989 {
990 for (unsigned int k = 0; k < n_q_points; ++k)
991 fe_data.untransformed_shape_hessian_tensors[k] =
992 fe_data.shape_grad_grads[dof_index][k];
993
994 mapping.transform(
996 fe_data.untransformed_shape_hessian_tensors),
998 mapping_internal,
999 make_array_view(fe_data.transformed_shape_hessians));
1000
1001 for (unsigned int k = 0; k < n_q_points; ++k)
1002 for (unsigned int d = 0; d < spacedim; ++d)
1003 for (unsigned int n = 0; n < spacedim; ++n)
1004 for (unsigned int i = 0; i < spacedim; ++i)
1005 for (unsigned int j = 0; j < spacedim; ++j)
1006 {
1007 fe_data.transformed_shape_hessians[k][d][i][j] -=
1008 (output_data.shape_values(first + n, k) *
1009 mapping_data
1011 [k][n][d][i][j]) +
1012 (output_data.shape_gradients[first + d][k][n] *
1013 mapping_data
1014 .jacobian_pushed_forward_grads[k][n][i][j]) +
1015 (output_data.shape_gradients[first + n][k][i] *
1016 mapping_data
1017 .jacobian_pushed_forward_grads[k][n][d][j]) +
1018 (output_data.shape_gradients[first + n][k][j] *
1019 mapping_data
1020 .jacobian_pushed_forward_grads[k][n][i][d]);
1021 }
1022
1023 for (unsigned int k = 0; k < n_q_points; ++k)
1024 for (unsigned int d = 0; d < dim; ++d)
1025 output_data.shape_hessians[first + d][k] =
1026 dof_sign * fe_data.transformed_shape_hessians[k][d];
1027
1028 break;
1029 }
1030
1031 default:
1033 }
1034 }
1035
1036 // third derivatives are not implemented
1037 if (fe_data.update_each & update_3rd_derivatives &&
1038 ((cell_similarity != CellSimilarity::translation) ||
1039 ((mapping_kind == mapping_piola) ||
1040 (mapping_kind == mapping_raviart_thomas) ||
1041 (mapping_kind == mapping_nedelec))))
1042 {
1044 }
1045 }
1046}
1047
1048
1049
1050template <int dim, int spacedim>
1051void
1054 const unsigned int face_no,
1055 const hp::QCollection<dim - 1> &quadrature,
1056 const Mapping<dim, spacedim> &mapping,
1057 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1059 &mapping_data,
1060 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1062 spacedim>
1063 &output_data) const
1064{
1065 AssertDimension(quadrature.size(), 1);
1066
1067 // convert data object to internal
1068 // data for this class. fails with
1069 // an exception if that is not
1070 // possible
1071 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1073 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1074
1075 const unsigned int n_q_points = quadrature[0].size();
1076 // offset determines which data set
1077 // to take (all data sets for all
1078 // faces are stored contiguously)
1079
1080 const auto offset =
1082 face_no,
1083 cell->face_orientation(face_no),
1084 cell->face_flip(face_no),
1085 cell->face_rotation(face_no),
1086 n_q_points);
1087
1088 // TODO: Size assertions
1089
1090 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1091 // correct thing on complicated meshes for higher order Nedelec elements.
1092 // Something similar to FE_Q should be done to permute dofs and to change the
1093 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1094 // class) is preferable.
1095 std::fill(fe_data.dof_sign_change.begin(),
1096 fe_data.dof_sign_change.end(),
1097 1.0);
1098 if (fe_data.update_each & update_values)
1099 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
1100 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
1101
1102 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1103 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1104 // Also nothing in 3d since we take care of it by using the
1105 // adjust_quad_dof_sign_for_face_orientation_table.
1106 if (fe_data.update_each & update_values)
1107 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1108 *this,
1109 this->mapping_kind,
1110 fe_data.dof_sign_change);
1111
1112 // What is the first dof_index on a quad?
1113 const unsigned int first_quad_index = this->get_first_quad_index();
1114 // How many dofs per quad and how many quad dofs do we have at all?
1115 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1116 const unsigned int n_quad_dofs =
1117 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1118
1119 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1120 ++dof_index)
1121 {
1122 /*
1123 * This assumes that the dofs are ordered by first vertices, lines, quads
1124 * and volume dofs. Note that in 2d this always gives false.
1125 */
1126 const bool is_quad_dof =
1127 (dim == 2 ? false :
1128 (first_quad_index <= dof_index) &&
1129 (dof_index < first_quad_index + n_quad_dofs));
1130
1131 // TODO: This hack is not pretty and it is only here to handle the 2d
1132 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1133 // handled by the
1134 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1135 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1136 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1137 // it is handles with a table. This array is allocated in
1138 // fe_poly_tensor.h.
1139 double dof_sign = 1.0;
1140 // under some circumstances fe_data.dof_sign_change is not allocated
1141 if (fe_data.update_each & update_values)
1142 dof_sign = fe_data.dof_sign_change[dof_index];
1143
1144 if (is_quad_dof)
1145 {
1146 /*
1147 * Find the face belonging to this dof_index. This is integer
1148 * division.
1149 */
1150 unsigned int face_index_from_dof_index =
1151 (dof_index - first_quad_index) / (n_dofs_per_quad);
1152
1153 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1154
1155 // Correct the dof_sign if necessary
1156 if (adjust_quad_dof_sign_for_face_orientation(
1157 local_quad_dof_index,
1158 face_index_from_dof_index,
1159 cell->face_orientation(face_index_from_dof_index),
1160 cell->face_flip(face_index_from_dof_index),
1161 cell->face_rotation(face_index_from_dof_index)))
1162 dof_sign = -1.0;
1163 }
1164
1165 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1166
1167 const unsigned int first =
1168 output_data.shape_function_to_row_table
1169 [dof_index * this->n_components() +
1170 this->get_nonzero_components(dof_index).first_selected_component()];
1171
1172 if (fe_data.update_each & update_values)
1173 {
1174 switch (mapping_kind)
1175 {
1176 case mapping_none:
1177 {
1178 for (unsigned int k = 0; k < n_q_points; ++k)
1179 for (unsigned int d = 0; d < dim; ++d)
1180 output_data.shape_values(first + d, k) =
1181 fe_data.shape_values[dof_index][k + offset][d];
1182 break;
1183 }
1184
1185 case mapping_covariant:
1187 {
1189 transformed_shape_values =
1190 make_array_view(fe_data.transformed_shape_values,
1191 offset,
1192 n_q_points);
1193 mapping.transform(make_array_view(fe_data.shape_values,
1194 dof_index,
1195 offset,
1196 n_q_points),
1197 mapping_kind,
1198 mapping_internal,
1199 transformed_shape_values);
1200
1201 for (unsigned int k = 0; k < n_q_points; ++k)
1202 for (unsigned int d = 0; d < dim; ++d)
1203 output_data.shape_values(first + d, k) =
1204 transformed_shape_values[k][d];
1205
1206 break;
1207 }
1209 case mapping_piola:
1210 {
1212 transformed_shape_values =
1213 make_array_view(fe_data.transformed_shape_values,
1214 offset,
1215 n_q_points);
1216 mapping.transform(make_array_view(fe_data.shape_values,
1217 dof_index,
1218 offset,
1219 n_q_points),
1221 mapping_internal,
1222 transformed_shape_values);
1223 for (unsigned int k = 0; k < n_q_points; ++k)
1224 for (unsigned int d = 0; d < dim; ++d)
1225 output_data.shape_values(first + d, k) =
1226 dof_sign * transformed_shape_values[k][d];
1227 break;
1228 }
1229
1230 case mapping_nedelec:
1231 {
1233 transformed_shape_values =
1234 make_array_view(fe_data.transformed_shape_values,
1235 offset,
1236 n_q_points);
1237 mapping.transform(make_array_view(fe_data.shape_values,
1238 dof_index,
1239 offset,
1240 n_q_points),
1242 mapping_internal,
1243 transformed_shape_values);
1244
1245 for (unsigned int k = 0; k < n_q_points; ++k)
1246 for (unsigned int d = 0; d < dim; ++d)
1247 output_data.shape_values(first + d, k) =
1248 dof_sign * transformed_shape_values[k][d];
1249
1250 break;
1251 }
1252
1253 default:
1255 }
1256 }
1257
1258 if (fe_data.update_each & update_gradients)
1259 {
1260 switch (mapping_kind)
1261 {
1262 case mapping_none:
1263 {
1264 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1265 make_array_view(fe_data.transformed_shape_grads,
1266 offset,
1267 n_q_points);
1268 mapping.transform(make_array_view(fe_data.shape_grads,
1269 dof_index,
1270 offset,
1271 n_q_points),
1273 mapping_internal,
1274 transformed_shape_grads);
1275 for (unsigned int k = 0; k < n_q_points; ++k)
1276 for (unsigned int d = 0; d < dim; ++d)
1277 output_data.shape_gradients[first + d][k] =
1278 transformed_shape_grads[k][d];
1279 break;
1280 }
1281
1282 case mapping_covariant:
1283 {
1284 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1285 make_array_view(fe_data.transformed_shape_grads,
1286 offset,
1287 n_q_points);
1288 mapping.transform(make_array_view(fe_data.shape_grads,
1289 dof_index,
1290 offset,
1291 n_q_points),
1293 mapping_internal,
1294 transformed_shape_grads);
1295
1296 for (unsigned int k = 0; k < n_q_points; ++k)
1297 for (unsigned int d = 0; d < spacedim; ++d)
1298 for (unsigned int n = 0; n < spacedim; ++n)
1299 transformed_shape_grads[k][d] -=
1300 output_data.shape_values(first + n, k) *
1301 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1302
1303 for (unsigned int k = 0; k < n_q_points; ++k)
1304 for (unsigned int d = 0; d < dim; ++d)
1305 output_data.shape_gradients[first + d][k] =
1306 transformed_shape_grads[k][d];
1307 break;
1308 }
1309
1311 {
1312 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1313 make_array_view(fe_data.transformed_shape_grads,
1314 offset,
1315 n_q_points);
1316 for (unsigned int k = 0; k < n_q_points; ++k)
1317 fe_data.untransformed_shape_grads[k + offset] =
1318 fe_data.shape_grads[dof_index][k + offset];
1319 mapping.transform(
1320 make_array_view(fe_data.untransformed_shape_grads,
1321 offset,
1322 n_q_points),
1324 mapping_internal,
1325 transformed_shape_grads);
1326
1327 for (unsigned int k = 0; k < n_q_points; ++k)
1328 for (unsigned int d = 0; d < spacedim; ++d)
1329 for (unsigned int n = 0; n < spacedim; ++n)
1330 transformed_shape_grads[k][d] +=
1331 output_data.shape_values(first + n, k) *
1332 mapping_data.jacobian_pushed_forward_grads[k][d][n];
1333
1334 for (unsigned int k = 0; k < n_q_points; ++k)
1335 for (unsigned int d = 0; d < dim; ++d)
1336 output_data.shape_gradients[first + d][k] =
1337 transformed_shape_grads[k][d];
1338
1339 break;
1340 }
1341
1343 case mapping_piola:
1344 {
1345 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1346 make_array_view(fe_data.transformed_shape_grads,
1347 offset,
1348 n_q_points);
1349 for (unsigned int k = 0; k < n_q_points; ++k)
1350 fe_data.untransformed_shape_grads[k + offset] =
1351 fe_data.shape_grads[dof_index][k + offset];
1352 mapping.transform(
1353 make_array_view(fe_data.untransformed_shape_grads,
1354 offset,
1355 n_q_points),
1357 mapping_internal,
1358 transformed_shape_grads);
1359
1360 for (unsigned int k = 0; k < n_q_points; ++k)
1361 for (unsigned int d = 0; d < spacedim; ++d)
1362 for (unsigned int n = 0; n < spacedim; ++n)
1363 transformed_shape_grads[k][d] +=
1364 (output_data.shape_values(first + n, k) *
1365 mapping_data
1367 (output_data.shape_values(first + d, k) *
1368 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
1369
1370 for (unsigned int k = 0; k < n_q_points; ++k)
1371 for (unsigned int d = 0; d < dim; ++d)
1372 output_data.shape_gradients[first + d][k] =
1373 dof_sign * transformed_shape_grads[k][d];
1374
1375 break;
1376 }
1377
1378 case mapping_nedelec:
1379 {
1380 // treat the gradients of
1381 // this particular shape
1382 // function at all
1383 // q-points. if Dv is the
1384 // gradient of the shape
1385 // function on the unit
1386 // cell, then
1387 // (J^-T)Dv(J^-1) is the
1388 // value we want to have on
1389 // the real cell.
1390 for (unsigned int k = 0; k < n_q_points; ++k)
1391 fe_data.untransformed_shape_grads[k + offset] =
1392 fe_data.shape_grads[dof_index][k + offset];
1393
1394 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1395 make_array_view(fe_data.transformed_shape_grads,
1396 offset,
1397 n_q_points);
1398 mapping.transform(
1399 make_array_view(fe_data.untransformed_shape_grads,
1400 offset,
1401 n_q_points),
1403 mapping_internal,
1404 transformed_shape_grads);
1405
1406 for (unsigned int k = 0; k < n_q_points; ++k)
1407 for (unsigned int d = 0; d < spacedim; ++d)
1408 for (unsigned int n = 0; n < spacedim; ++n)
1409 transformed_shape_grads[k][d] -=
1410 output_data.shape_values(first + n, k) *
1411 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1412
1413 for (unsigned int k = 0; k < n_q_points; ++k)
1414 for (unsigned int d = 0; d < dim; ++d)
1415 output_data.shape_gradients[first + d][k] =
1416 dof_sign * transformed_shape_grads[k][d];
1417
1418 break;
1419 }
1420
1421 default:
1423 }
1424 }
1425
1426 if (fe_data.update_each & update_hessians)
1427 {
1428 switch (mapping_kind)
1429 {
1430 case mapping_none:
1431 {
1433 transformed_shape_hessians =
1434 make_array_view(fe_data.transformed_shape_hessians,
1435 offset,
1436 n_q_points);
1437 mapping.transform(make_array_view(fe_data.shape_grad_grads,
1438 dof_index,
1439 offset,
1440 n_q_points),
1442 mapping_internal,
1443 transformed_shape_hessians);
1444
1445 for (unsigned int k = 0; k < n_q_points; ++k)
1446 for (unsigned int d = 0; d < spacedim; ++d)
1447 for (unsigned int n = 0; n < spacedim; ++n)
1448 transformed_shape_hessians[k][d] -=
1449 output_data.shape_gradients[first + d][k][n] *
1450 mapping_data.jacobian_pushed_forward_grads[k][n];
1451
1452 for (unsigned int k = 0; k < n_q_points; ++k)
1453 for (unsigned int d = 0; d < dim; ++d)
1454 output_data.shape_hessians[first + d][k] =
1455 transformed_shape_hessians[k][d];
1456
1457 break;
1458 }
1459 case mapping_covariant:
1460 {
1461 for (unsigned int k = 0; k < n_q_points; ++k)
1462 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1463 fe_data.shape_grad_grads[dof_index][k + offset];
1464
1466 transformed_shape_hessians =
1467 make_array_view(fe_data.transformed_shape_hessians,
1468 offset,
1469 n_q_points);
1470 mapping.transform(
1471 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1472 offset,
1473 n_q_points),
1475 mapping_internal,
1476 transformed_shape_hessians);
1477
1478 for (unsigned int k = 0; k < n_q_points; ++k)
1479 for (unsigned int d = 0; d < spacedim; ++d)
1480 for (unsigned int n = 0; n < spacedim; ++n)
1481 for (unsigned int i = 0; i < spacedim; ++i)
1482 for (unsigned int j = 0; j < spacedim; ++j)
1483 {
1484 transformed_shape_hessians[k][d][i][j] -=
1485 (output_data.shape_values(first + n, k) *
1486 mapping_data
1488 [k][n][d][i][j]) +
1489 (output_data.shape_gradients[first + d][k][n] *
1490 mapping_data
1491 .jacobian_pushed_forward_grads[k][n][i][j]) +
1492 (output_data.shape_gradients[first + n][k][i] *
1493 mapping_data
1494 .jacobian_pushed_forward_grads[k][n][d][j]) +
1495 (output_data.shape_gradients[first + n][k][j] *
1496 mapping_data
1497 .jacobian_pushed_forward_grads[k][n][i][d]);
1498 }
1499
1500 for (unsigned int k = 0; k < n_q_points; ++k)
1501 for (unsigned int d = 0; d < dim; ++d)
1502 output_data.shape_hessians[first + d][k] =
1503 transformed_shape_hessians[k][d];
1504
1505 break;
1506 }
1507
1509 {
1510 for (unsigned int k = 0; k < n_q_points; ++k)
1511 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1512 fe_data.shape_grad_grads[dof_index][k + offset];
1513
1515 transformed_shape_hessians =
1516 make_array_view(fe_data.transformed_shape_hessians,
1517 offset,
1518 n_q_points);
1519 mapping.transform(
1520 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1521 offset,
1522 n_q_points),
1524 mapping_internal,
1525 transformed_shape_hessians);
1526
1527 for (unsigned int k = 0; k < n_q_points; ++k)
1528 for (unsigned int d = 0; d < spacedim; ++d)
1529 for (unsigned int n = 0; n < spacedim; ++n)
1530 for (unsigned int i = 0; i < spacedim; ++i)
1531 for (unsigned int j = 0; j < spacedim; ++j)
1532 {
1533 transformed_shape_hessians[k][d][i][j] +=
1534 (output_data.shape_values(first + n, k) *
1535 mapping_data
1537 [k][d][n][i][j]) +
1538 (output_data.shape_gradients[first + n][k][i] *
1539 mapping_data
1540 .jacobian_pushed_forward_grads[k][d][n][j]) +
1541 (output_data.shape_gradients[first + n][k][j] *
1542 mapping_data
1543 .jacobian_pushed_forward_grads[k][d][i][n]) -
1544 (output_data.shape_gradients[first + d][k][n] *
1545 mapping_data
1546 .jacobian_pushed_forward_grads[k][n][i][j]);
1547 for (unsigned int m = 0; m < spacedim; ++m)
1548 transformed_shape_hessians[k][d][i][j] -=
1549 (mapping_data
1550 .jacobian_pushed_forward_grads[k][d][i]
1551 [m] *
1552 mapping_data
1553 .jacobian_pushed_forward_grads[k][m][n]
1554 [j] *
1555 output_data.shape_values(first + n, k)) +
1556 (mapping_data
1557 .jacobian_pushed_forward_grads[k][d][m]
1558 [j] *
1559 mapping_data
1560 .jacobian_pushed_forward_grads[k][m][i]
1561 [n] *
1562 output_data.shape_values(first + n, k));
1563 }
1564
1565 for (unsigned int k = 0; k < n_q_points; ++k)
1566 for (unsigned int d = 0; d < dim; ++d)
1567 output_data.shape_hessians[first + d][k] =
1568 transformed_shape_hessians[k][d];
1569
1570 break;
1571 }
1572
1574 case mapping_piola:
1575 {
1576 for (unsigned int k = 0; k < n_q_points; ++k)
1577 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1578 fe_data.shape_grad_grads[dof_index][k + offset];
1579
1581 transformed_shape_hessians =
1582 make_array_view(fe_data.transformed_shape_hessians,
1583 offset,
1584 n_q_points);
1585 mapping.transform(
1586 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1587 offset,
1588 n_q_points),
1590 mapping_internal,
1591 transformed_shape_hessians);
1592
1593 for (unsigned int k = 0; k < n_q_points; ++k)
1594 for (unsigned int d = 0; d < spacedim; ++d)
1595 for (unsigned int n = 0; n < spacedim; ++n)
1596 for (unsigned int i = 0; i < spacedim; ++i)
1597 for (unsigned int j = 0; j < spacedim; ++j)
1598 {
1599 transformed_shape_hessians[k][d][i][j] +=
1600 (output_data.shape_values(first + n, k) *
1601 mapping_data
1603 [k][d][n][i][j]) +
1604 (output_data.shape_gradients[first + n][k][i] *
1605 mapping_data
1606 .jacobian_pushed_forward_grads[k][d][n][j]) +
1607 (output_data.shape_gradients[first + n][k][j] *
1608 mapping_data
1609 .jacobian_pushed_forward_grads[k][d][i][n]) -
1610 (output_data.shape_gradients[first + d][k][n] *
1611 mapping_data
1612 .jacobian_pushed_forward_grads[k][n][i][j]);
1613
1614 transformed_shape_hessians[k][d][i][j] -=
1615 (output_data.shape_values(first + d, k) *
1616 mapping_data
1618 [k][n][n][i][j]) +
1619 (output_data.shape_gradients[first + d][k][i] *
1620 mapping_data
1621 .jacobian_pushed_forward_grads[k][n][n][j]) +
1622 (output_data.shape_gradients[first + d][k][j] *
1623 mapping_data
1624 .jacobian_pushed_forward_grads[k][n][n][i]);
1625
1626 for (unsigned int m = 0; m < spacedim; ++m)
1627 {
1628 transformed_shape_hessians[k][d][i][j] -=
1629 (mapping_data
1631 [m] *
1632 mapping_data
1634 [j] *
1635 output_data.shape_values(first + n, k)) +
1636 (mapping_data
1637 .jacobian_pushed_forward_grads[k][d][m]
1638 [j] *
1639 mapping_data
1640 .jacobian_pushed_forward_grads[k][m][i]
1641 [n] *
1642 output_data.shape_values(first + n, k));
1643
1644 transformed_shape_hessians[k][d][i][j] +=
1645 (mapping_data
1647 [m] *
1648 mapping_data
1650 [j] *
1651 output_data.shape_values(first + d, k)) +
1652 (mapping_data
1653 .jacobian_pushed_forward_grads[k][n][m]
1654 [j] *
1655 mapping_data
1656 .jacobian_pushed_forward_grads[k][m][i]
1657 [n] *
1658 output_data.shape_values(first + d, k));
1659 }
1660 }
1661
1662 for (unsigned int k = 0; k < n_q_points; ++k)
1663 for (unsigned int d = 0; d < dim; ++d)
1664 output_data.shape_hessians[first + d][k] =
1665 dof_sign * transformed_shape_hessians[k][d];
1666
1667 break;
1668 }
1669
1670 case mapping_nedelec:
1671 {
1672 for (unsigned int k = 0; k < n_q_points; ++k)
1673 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1674 fe_data.shape_grad_grads[dof_index][k + offset];
1675
1677 transformed_shape_hessians =
1678 make_array_view(fe_data.transformed_shape_hessians,
1679 offset,
1680 n_q_points);
1681 mapping.transform(
1682 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1683 offset,
1684 n_q_points),
1686 mapping_internal,
1687 transformed_shape_hessians);
1688
1689 for (unsigned int k = 0; k < n_q_points; ++k)
1690 for (unsigned int d = 0; d < spacedim; ++d)
1691 for (unsigned int n = 0; n < spacedim; ++n)
1692 for (unsigned int i = 0; i < spacedim; ++i)
1693 for (unsigned int j = 0; j < spacedim; ++j)
1694 {
1695 transformed_shape_hessians[k][d][i][j] -=
1696 (output_data.shape_values(first + n, k) *
1697 mapping_data
1699 [k][n][d][i][j]) +
1700 (output_data.shape_gradients[first + d][k][n] *
1701 mapping_data
1702 .jacobian_pushed_forward_grads[k][n][i][j]) +
1703 (output_data.shape_gradients[first + n][k][i] *
1704 mapping_data
1705 .jacobian_pushed_forward_grads[k][n][d][j]) +
1706 (output_data.shape_gradients[first + n][k][j] *
1707 mapping_data
1708 .jacobian_pushed_forward_grads[k][n][i][d]);
1709 }
1710
1711 for (unsigned int k = 0; k < n_q_points; ++k)
1712 for (unsigned int d = 0; d < dim; ++d)
1713 output_data.shape_hessians[first + d][k] =
1714 dof_sign * transformed_shape_hessians[k][d];
1715
1716 break;
1717 }
1718
1719 default:
1721 }
1722 }
1723
1724 // third derivatives are not implemented
1725 if (fe_data.update_each & update_3rd_derivatives)
1726 {
1728 }
1729 }
1730}
1731
1732
1733
1734template <int dim, int spacedim>
1735void
1738 const unsigned int face_no,
1739 const unsigned int sub_no,
1740 const Quadrature<dim - 1> &quadrature,
1741 const Mapping<dim, spacedim> &mapping,
1742 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1744 &mapping_data,
1745 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1747 spacedim>
1748 &output_data) const
1749{
1750 // convert data object to internal
1751 // data for this class. fails with
1752 // an exception if that is not
1753 // possible
1754 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1756 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1757
1758 const unsigned int n_q_points = quadrature.size();
1759
1760 // offset determines which data set
1761 // to take (all data sets for all
1762 // sub-faces are stored contiguously)
1763 const auto offset =
1765 face_no,
1766 sub_no,
1767 cell->face_orientation(face_no),
1768 cell->face_flip(face_no),
1769 cell->face_rotation(face_no),
1770 n_q_points,
1771 cell->subface_case(face_no));
1772
1773 // TODO: Size assertions
1774
1775 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1776 // correct thing on complicated meshes for higher order Nedelec elements.
1777 // Something similar to FE_Q should be done to permute dofs and to change the
1778 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1779 // class) is preferable.
1780 std::fill(fe_data.dof_sign_change.begin(),
1781 fe_data.dof_sign_change.end(),
1782 1.0);
1783 if (fe_data.update_each & update_values)
1784 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
1785 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
1786
1787 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1788 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1789 // Also nothing in 3d since we take care of it by using the
1790 // adjust_quad_dof_sign_for_face_orientation_table.
1791 if (fe_data.update_each & update_values)
1792 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1793 *this,
1794 this->mapping_kind,
1795 fe_data.dof_sign_change);
1796
1797 // What is the first dof_index on a quad?
1798 const unsigned int first_quad_index = this->get_first_quad_index();
1799 // How many dofs per quad and how many quad dofs do we have at all?
1800 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1801 const unsigned int n_quad_dofs =
1802 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1803
1804 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1805 ++dof_index)
1806 {
1807 /*
1808 * This assumes that the dofs are ordered by first vertices, lines, quads
1809 * and volume dofs. Note that in 2d this always gives false.
1810 */
1811 const bool is_quad_dof =
1812 (dim == 2 ? false :
1813 (first_quad_index <= dof_index) &&
1814 (dof_index < first_quad_index + n_quad_dofs));
1815
1816 // TODO: This hack is not pretty and it is only here to handle the 2d
1817 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1818 // handled by the
1819 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1820 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1821 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1822 // it is handles with a table. This array is allocated in
1823 // fe_poly_tensor.h.
1824 double dof_sign = 1.0;
1825 // under some circumstances fe_data.dof_sign_change is not allocated
1826 if (fe_data.update_each & update_values)
1827 dof_sign = fe_data.dof_sign_change[dof_index];
1828
1829 if (is_quad_dof)
1830 {
1831 /*
1832 * Find the face belonging to this dof_index. This is integer
1833 * division.
1834 */
1835 unsigned int face_index_from_dof_index =
1836 (dof_index - first_quad_index) / (n_dofs_per_quad);
1837
1838 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1839
1840 // Correct the dof_sign if necessary
1841 if (adjust_quad_dof_sign_for_face_orientation(
1842 local_quad_dof_index,
1843 face_index_from_dof_index,
1844 cell->face_orientation(face_index_from_dof_index),
1845 cell->face_flip(face_index_from_dof_index),
1846 cell->face_rotation(face_index_from_dof_index)))
1847 dof_sign = -1.0;
1848 }
1849
1850 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1851
1852 const unsigned int first =
1853 output_data.shape_function_to_row_table
1854 [dof_index * this->n_components() +
1855 this->get_nonzero_components(dof_index).first_selected_component()];
1856
1857 if (fe_data.update_each & update_values)
1858 {
1859 switch (mapping_kind)
1860 {
1861 case mapping_none:
1862 {
1863 for (unsigned int k = 0; k < n_q_points; ++k)
1864 for (unsigned int d = 0; d < dim; ++d)
1865 output_data.shape_values(first + d, k) =
1866 fe_data.shape_values[dof_index][k + offset][d];
1867 break;
1868 }
1869
1870 case mapping_covariant:
1872 {
1874 transformed_shape_values =
1875 make_array_view(fe_data.transformed_shape_values,
1876 offset,
1877 n_q_points);
1878 mapping.transform(make_array_view(fe_data.shape_values,
1879 dof_index,
1880 offset,
1881 n_q_points),
1882 mapping_kind,
1883 mapping_internal,
1884 transformed_shape_values);
1885
1886 for (unsigned int k = 0; k < n_q_points; ++k)
1887 for (unsigned int d = 0; d < dim; ++d)
1888 output_data.shape_values(first + d, k) =
1889 transformed_shape_values[k][d];
1890
1891 break;
1892 }
1893
1895 case mapping_piola:
1896 {
1898 transformed_shape_values =
1899 make_array_view(fe_data.transformed_shape_values,
1900 offset,
1901 n_q_points);
1902
1903 mapping.transform(make_array_view(fe_data.shape_values,
1904 dof_index,
1905 offset,
1906 n_q_points),
1908 mapping_internal,
1909 transformed_shape_values);
1910 for (unsigned int k = 0; k < n_q_points; ++k)
1911 for (unsigned int d = 0; d < dim; ++d)
1912 output_data.shape_values(first + d, k) =
1913 dof_sign * transformed_shape_values[k][d];
1914 break;
1915 }
1916
1917 case mapping_nedelec:
1918 {
1920 transformed_shape_values =
1921 make_array_view(fe_data.transformed_shape_values,
1922 offset,
1923 n_q_points);
1924
1925 mapping.transform(make_array_view(fe_data.shape_values,
1926 dof_index,
1927 offset,
1928 n_q_points),
1930 mapping_internal,
1931 transformed_shape_values);
1932
1933 for (unsigned int k = 0; k < n_q_points; ++k)
1934 for (unsigned int d = 0; d < dim; ++d)
1935 output_data.shape_values(first + d, k) =
1936 dof_sign * transformed_shape_values[k][d];
1937
1938 break;
1939 }
1940
1941 default:
1943 }
1944 }
1945
1946 if (fe_data.update_each & update_gradients)
1947 {
1948 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1949 make_array_view(fe_data.transformed_shape_grads,
1950 offset,
1951 n_q_points);
1952 switch (mapping_kind)
1953 {
1954 case mapping_none:
1955 {
1956 mapping.transform(make_array_view(fe_data.shape_grads,
1957 dof_index,
1958 offset,
1959 n_q_points),
1961 mapping_internal,
1962 transformed_shape_grads);
1963 for (unsigned int k = 0; k < n_q_points; ++k)
1964 for (unsigned int d = 0; d < dim; ++d)
1965 output_data.shape_gradients[first + d][k] =
1966 transformed_shape_grads[k][d];
1967 break;
1968 }
1969
1970 case mapping_covariant:
1971 {
1972 mapping.transform(make_array_view(fe_data.shape_grads,
1973 dof_index,
1974 offset,
1975 n_q_points),
1977 mapping_internal,
1978 transformed_shape_grads);
1979
1980 for (unsigned int k = 0; k < n_q_points; ++k)
1981 for (unsigned int d = 0; d < spacedim; ++d)
1982 for (unsigned int n = 0; n < spacedim; ++n)
1983 transformed_shape_grads[k][d] -=
1984 output_data.shape_values(first + n, k) *
1985 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1986
1987 for (unsigned int k = 0; k < n_q_points; ++k)
1988 for (unsigned int d = 0; d < dim; ++d)
1989 output_data.shape_gradients[first + d][k] =
1990 transformed_shape_grads[k][d];
1991
1992 break;
1993 }
1994
1996 {
1997 for (unsigned int k = 0; k < n_q_points; ++k)
1998 fe_data.untransformed_shape_grads[k + offset] =
1999 fe_data.shape_grads[dof_index][k + offset];
2000
2001 mapping.transform(
2002 make_array_view(fe_data.untransformed_shape_grads,
2003 offset,
2004 n_q_points),
2006 mapping_internal,
2007 transformed_shape_grads);
2008
2009 for (unsigned int k = 0; k < n_q_points; ++k)
2010 for (unsigned int d = 0; d < spacedim; ++d)
2011 for (unsigned int n = 0; n < spacedim; ++n)
2012 transformed_shape_grads[k][d] +=
2013 output_data.shape_values(first + n, k) *
2014 mapping_data.jacobian_pushed_forward_grads[k][d][n];
2015
2016 for (unsigned int k = 0; k < n_q_points; ++k)
2017 for (unsigned int d = 0; d < dim; ++d)
2018 output_data.shape_gradients[first + d][k] =
2019 transformed_shape_grads[k][d];
2020
2021 break;
2022 }
2023
2025 case mapping_piola:
2026 {
2027 for (unsigned int k = 0; k < n_q_points; ++k)
2028 fe_data.untransformed_shape_grads[k + offset] =
2029 fe_data.shape_grads[dof_index][k + offset];
2030
2031 mapping.transform(
2032 make_array_view(fe_data.untransformed_shape_grads,
2033 offset,
2034 n_q_points),
2036 mapping_internal,
2037 transformed_shape_grads);
2038
2039 for (unsigned int k = 0; k < n_q_points; ++k)
2040 for (unsigned int d = 0; d < spacedim; ++d)
2041 for (unsigned int n = 0; n < spacedim; ++n)
2042 transformed_shape_grads[k][d] +=
2043 (output_data.shape_values(first + n, k) *
2044 mapping_data
2046 (output_data.shape_values(first + d, k) *
2047 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
2048
2049 for (unsigned int k = 0; k < n_q_points; ++k)
2050 for (unsigned int d = 0; d < dim; ++d)
2051 output_data.shape_gradients[first + d][k] =
2052 dof_sign * transformed_shape_grads[k][d];
2053
2054 break;
2055 }
2056
2057 case mapping_nedelec:
2058 {
2059 // this particular shape
2060 // function at all
2061 // q-points. if Dv is the
2062 // gradient of the shape
2063 // function on the unit
2064 // cell, then
2065 // (J^-T)Dv(J^-1) is the
2066 // value we want to have on
2067 // the real cell.
2068 for (unsigned int k = 0; k < n_q_points; ++k)
2069 fe_data.untransformed_shape_grads[k + offset] =
2070 fe_data.shape_grads[dof_index][k + offset];
2071
2072 mapping.transform(
2073 make_array_view(fe_data.untransformed_shape_grads,
2074 offset,
2075 n_q_points),
2077 mapping_internal,
2078 transformed_shape_grads);
2079
2080 for (unsigned int k = 0; k < n_q_points; ++k)
2081 for (unsigned int d = 0; d < spacedim; ++d)
2082 for (unsigned int n = 0; n < spacedim; ++n)
2083 transformed_shape_grads[k][d] -=
2084 output_data.shape_values(first + n, k) *
2085 mapping_data.jacobian_pushed_forward_grads[k][n][d];
2086
2087 for (unsigned int k = 0; k < n_q_points; ++k)
2088 for (unsigned int d = 0; d < dim; ++d)
2089 output_data.shape_gradients[first + d][k] =
2090 dof_sign * transformed_shape_grads[k][d];
2091
2092 break;
2093 }
2094
2095 default:
2097 }
2098 }
2099
2100 if (fe_data.update_each & update_hessians)
2101 {
2102 switch (mapping_kind)
2103 {
2104 case mapping_none:
2105 {
2107 transformed_shape_hessians =
2108 make_array_view(fe_data.transformed_shape_hessians,
2109 offset,
2110 n_q_points);
2111 mapping.transform(make_array_view(fe_data.shape_grad_grads,
2112 dof_index,
2113 offset,
2114 n_q_points),
2116 mapping_internal,
2117 transformed_shape_hessians);
2118
2119 for (unsigned int k = 0; k < n_q_points; ++k)
2120 for (unsigned int d = 0; d < spacedim; ++d)
2121 for (unsigned int n = 0; n < spacedim; ++n)
2122 transformed_shape_hessians[k][d] -=
2123 output_data.shape_gradients[first + d][k][n] *
2124 mapping_data.jacobian_pushed_forward_grads[k][n];
2125
2126 for (unsigned int k = 0; k < n_q_points; ++k)
2127 for (unsigned int d = 0; d < dim; ++d)
2128 output_data.shape_hessians[first + d][k] =
2129 transformed_shape_hessians[k][d];
2130
2131 break;
2132 }
2133 case mapping_covariant:
2134 {
2135 for (unsigned int k = 0; k < n_q_points; ++k)
2136 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2137 fe_data.shape_grad_grads[dof_index][k + offset];
2138
2140 transformed_shape_hessians =
2141 make_array_view(fe_data.transformed_shape_hessians,
2142 offset,
2143 n_q_points);
2144 mapping.transform(
2145 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2146 offset,
2147 n_q_points),
2149 mapping_internal,
2150 transformed_shape_hessians);
2151
2152 for (unsigned int k = 0; k < n_q_points; ++k)
2153 for (unsigned int d = 0; d < spacedim; ++d)
2154 for (unsigned int n = 0; n < spacedim; ++n)
2155 for (unsigned int i = 0; i < spacedim; ++i)
2156 for (unsigned int j = 0; j < spacedim; ++j)
2157 {
2158 transformed_shape_hessians[k][d][i][j] -=
2159 (output_data.shape_values(first + n, k) *
2160 mapping_data
2162 [k][n][d][i][j]) +
2163 (output_data.shape_gradients[first + d][k][n] *
2164 mapping_data
2165 .jacobian_pushed_forward_grads[k][n][i][j]) +
2166 (output_data.shape_gradients[first + n][k][i] *
2167 mapping_data
2168 .jacobian_pushed_forward_grads[k][n][d][j]) +
2169 (output_data.shape_gradients[first + n][k][j] *
2170 mapping_data
2171 .jacobian_pushed_forward_grads[k][n][i][d]);
2172 }
2173
2174 for (unsigned int k = 0; k < n_q_points; ++k)
2175 for (unsigned int d = 0; d < dim; ++d)
2176 output_data.shape_hessians[first + d][k] =
2177 transformed_shape_hessians[k][d];
2178
2179 break;
2180 }
2181
2183 {
2184 for (unsigned int k = 0; k < n_q_points; ++k)
2185 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2186 fe_data.shape_grad_grads[dof_index][k + offset];
2187
2189 transformed_shape_hessians =
2190 make_array_view(fe_data.transformed_shape_hessians,
2191 offset,
2192 n_q_points);
2193 mapping.transform(
2194 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2195 offset,
2196 n_q_points),
2198 mapping_internal,
2199 transformed_shape_hessians);
2200
2201 for (unsigned int k = 0; k < n_q_points; ++k)
2202 for (unsigned int d = 0; d < spacedim; ++d)
2203 for (unsigned int n = 0; n < spacedim; ++n)
2204 for (unsigned int i = 0; i < spacedim; ++i)
2205 for (unsigned int j = 0; j < spacedim; ++j)
2206 {
2207 transformed_shape_hessians[k][d][i][j] +=
2208 (output_data.shape_values(first + n, k) *
2209 mapping_data
2211 [k][d][n][i][j]) +
2212 (output_data.shape_gradients[first + n][k][i] *
2213 mapping_data
2214 .jacobian_pushed_forward_grads[k][d][n][j]) +
2215 (output_data.shape_gradients[first + n][k][j] *
2216 mapping_data
2217 .jacobian_pushed_forward_grads[k][d][i][n]) -
2218 (output_data.shape_gradients[first + d][k][n] *
2219 mapping_data
2220 .jacobian_pushed_forward_grads[k][n][i][j]);
2221 for (unsigned int m = 0; m < spacedim; ++m)
2222 transformed_shape_hessians[k][d][i][j] -=
2223 (mapping_data
2224 .jacobian_pushed_forward_grads[k][d][i]
2225 [m] *
2226 mapping_data
2227 .jacobian_pushed_forward_grads[k][m][n]
2228 [j] *
2229 output_data.shape_values(first + n, k)) +
2230 (mapping_data
2231 .jacobian_pushed_forward_grads[k][d][m]
2232 [j] *
2233 mapping_data
2234 .jacobian_pushed_forward_grads[k][m][i]
2235 [n] *
2236 output_data.shape_values(first + n, k));
2237 }
2238
2239 for (unsigned int k = 0; k < n_q_points; ++k)
2240 for (unsigned int d = 0; d < dim; ++d)
2241 output_data.shape_hessians[first + d][k] =
2242 transformed_shape_hessians[k][d];
2243
2244 break;
2245 }
2246
2248 case mapping_piola:
2249 {
2250 for (unsigned int k = 0; k < n_q_points; ++k)
2251 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2252 fe_data.shape_grad_grads[dof_index][k + offset];
2253
2255 transformed_shape_hessians =
2256 make_array_view(fe_data.transformed_shape_hessians,
2257 offset,
2258 n_q_points);
2259 mapping.transform(
2260 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2261 offset,
2262 n_q_points),
2264 mapping_internal,
2265 transformed_shape_hessians);
2266
2267 for (unsigned int k = 0; k < n_q_points; ++k)
2268 for (unsigned int d = 0; d < spacedim; ++d)
2269 for (unsigned int n = 0; n < spacedim; ++n)
2270 for (unsigned int i = 0; i < spacedim; ++i)
2271 for (unsigned int j = 0; j < spacedim; ++j)
2272 {
2273 transformed_shape_hessians[k][d][i][j] +=
2274 (output_data.shape_values(first + n, k) *
2275 mapping_data
2277 [k][d][n][i][j]) +
2278 (output_data.shape_gradients[first + n][k][i] *
2279 mapping_data
2280 .jacobian_pushed_forward_grads[k][d][n][j]) +
2281 (output_data.shape_gradients[first + n][k][j] *
2282 mapping_data
2283 .jacobian_pushed_forward_grads[k][d][i][n]) -
2284 (output_data.shape_gradients[first + d][k][n] *
2285 mapping_data
2286 .jacobian_pushed_forward_grads[k][n][i][j]);
2287
2288 transformed_shape_hessians[k][d][i][j] -=
2289 (output_data.shape_values(first + d, k) *
2290 mapping_data
2292 [k][n][n][i][j]) +
2293 (output_data.shape_gradients[first + d][k][i] *
2294 mapping_data
2295 .jacobian_pushed_forward_grads[k][n][n][j]) +
2296 (output_data.shape_gradients[first + d][k][j] *
2297 mapping_data
2298 .jacobian_pushed_forward_grads[k][n][n][i]);
2299 for (unsigned int m = 0; m < spacedim; ++m)
2300 {
2301 transformed_shape_hessians[k][d][i][j] -=
2302 (mapping_data
2304 [m] *
2305 mapping_data
2307 [j] *
2308 output_data.shape_values(first + n, k)) +
2309 (mapping_data
2310 .jacobian_pushed_forward_grads[k][d][m]
2311 [j] *
2312 mapping_data
2313 .jacobian_pushed_forward_grads[k][m][i]
2314 [n] *
2315 output_data.shape_values(first + n, k));
2316
2317 transformed_shape_hessians[k][d][i][j] +=
2318 (mapping_data
2320 [m] *
2321 mapping_data
2323 [j] *
2324 output_data.shape_values(first + d, k)) +
2325 (mapping_data
2326 .jacobian_pushed_forward_grads[k][n][m]
2327 [j] *
2328 mapping_data
2329 .jacobian_pushed_forward_grads[k][m][i]
2330 [n] *
2331 output_data.shape_values(first + d, k));
2332 }
2333 }
2334
2335 for (unsigned int k = 0; k < n_q_points; ++k)
2336 for (unsigned int d = 0; d < dim; ++d)
2337 output_data.shape_hessians[first + d][k] =
2338 dof_sign * transformed_shape_hessians[k][d];
2339
2340 break;
2341 }
2342
2343 case mapping_nedelec:
2344 {
2345 for (unsigned int k = 0; k < n_q_points; ++k)
2346 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2347 fe_data.shape_grad_grads[dof_index][k + offset];
2348
2350 transformed_shape_hessians =
2351 make_array_view(fe_data.transformed_shape_hessians,
2352 offset,
2353 n_q_points);
2354 mapping.transform(
2355 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2356 offset,
2357 n_q_points),
2359 mapping_internal,
2360 transformed_shape_hessians);
2361
2362 for (unsigned int k = 0; k < n_q_points; ++k)
2363 for (unsigned int d = 0; d < spacedim; ++d)
2364 for (unsigned int n = 0; n < spacedim; ++n)
2365 for (unsigned int i = 0; i < spacedim; ++i)
2366 for (unsigned int j = 0; j < spacedim; ++j)
2367 {
2368 transformed_shape_hessians[k][d][i][j] -=
2369 (output_data.shape_values(first + n, k) *
2370 mapping_data
2372 [k][n][d][i][j]) +
2373 (output_data.shape_gradients[first + d][k][n] *
2374 mapping_data
2375 .jacobian_pushed_forward_grads[k][n][i][j]) +
2376 (output_data.shape_gradients[first + n][k][i] *
2377 mapping_data
2378 .jacobian_pushed_forward_grads[k][n][d][j]) +
2379 (output_data.shape_gradients[first + n][k][j] *
2380 mapping_data
2381 .jacobian_pushed_forward_grads[k][n][i][d]);
2382 }
2383
2384 for (unsigned int k = 0; k < n_q_points; ++k)
2385 for (unsigned int d = 0; d < dim; ++d)
2386 output_data.shape_hessians[first + d][k] =
2387 dof_sign * transformed_shape_hessians[k][d];
2388
2389 break;
2390 }
2391
2392 default:
2394 }
2395 }
2396
2397 // third derivatives are not implemented
2398 if (fe_data.update_each & update_3rd_derivatives)
2399 {
2401 }
2402 }
2403}
2404
2405
2406
2407template <int dim, int spacedim>
2410 const UpdateFlags flags) const
2411{
2413
2414 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2415 {
2416 const MappingKind mapping_kind = get_mapping_kind(i);
2417
2418 switch (mapping_kind)
2419 {
2420 case mapping_none:
2421 {
2422 if (flags & update_values)
2423 out |= update_values;
2424
2425 if (flags & update_gradients)
2428
2429 if (flags & update_hessians)
2433 break;
2434 }
2436 case mapping_piola:
2437 {
2438 if (flags & update_values)
2440
2441 if (flags & update_gradients)
2446
2447 if (flags & update_hessians)
2452
2453 break;
2454 }
2455
2456
2458 {
2459 if (flags & update_values)
2461
2462 if (flags & update_gradients)
2467
2468 if (flags & update_hessians)
2473
2474 break;
2475 }
2476
2477 case mapping_nedelec:
2478 case mapping_covariant:
2479 {
2480 if (flags & update_values)
2482
2483 if (flags & update_gradients)
2487
2488 if (flags & update_hessians)
2493
2494 break;
2495 }
2496
2497 default:
2498 {
2500 }
2501 }
2502 }
2503
2504 return out;
2505}
2506
2507#endif
2508// explicit instantiations
2509#include "fe_poly_tensor.inst"
2510
2511
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) 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 double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) 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 double shape_value(const unsigned int i, const Point< dim > &p) 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
bool adjust_quad_dof_sign_for_face_orientation(const unsigned int index, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation) 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
MappingKind get_mapping_kind(const unsigned int i) const
FE_PolyTensor(const TensorPolynomialsBase< dim > &polynomials, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
bool single_mapping_kind() const
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const unsigned char combined_orientation=ReferenceCell::default_combined_face_orientation()) const
Abstract base class for mapping classes.
Definition mapping.h:316
Definition point.h:111
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
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)
unsigned int size() const
unsigned int size() const
Definition collection.h:264
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > first
Definition grid_out.cc:4613
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_covariant_transformation
Covariant transformation.
@ update_gradients
Shape function gradients.
@ update_default
No update.
@ update_piola
Values needed for Piola transform.
MappingKind
Definition mapping.h:77
@ mapping_piola
Definition mapping.h:112
@ mapping_covariant_gradient
Definition mapping.h:98
@ mapping_covariant
Definition mapping.h:87
@ mapping_nedelec
Definition mapping.h:127
@ mapping_contravariant
Definition mapping.h:92
@ mapping_raviart_thomas
Definition mapping.h:132
@ mapping_contravariant_hessian
Definition mapping.h:154
@ mapping_covariant_hessian
Definition mapping.h:148
@ mapping_contravariant_gradient
Definition mapping.h:104
@ mapping_none
Definition mapping.h:82
@ mapping_piola_gradient
Definition mapping.h:118
@ mapping_piola_hessian
Definition mapping.h:160
#define DEAL_II_NOT_IMPLEMENTED()
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
spacedim & mapping
spacedim const Point< spacedim > & p
Definition grid_tools.h:980
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
const InputIterator OutputIterator out
Definition parallel.h:167