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