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