Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2846-g6fb608615f 2025-03-15 04:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
fe_poly_tensor.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
25
29
31#include <deal.II/grid/tria.h>
33
35
36namespace internal
37{
38 namespace FE_PolyTensor
39 {
40 namespace
41 {
42 template <int spacedim>
43 void
44 get_dof_sign_change_h_div(
45 const typename ::Triangulation<1, spacedim>::cell_iterator &,
47 const std::vector<MappingKind> &,
48 std::vector<double> &)
49 {
50 // Nothing to do in 1d.
51 }
52
53
54
55 // TODO: This function is not a consistent fix of the orientation issue
56 // like in 3d. It is rather kept not to break legacy behavior in 2d but
57 // should be replaced. See also the implementation of
58 // FE_RaviartThomas<dim>::initialize_quad_dof_index_permutation_and_sign_change()
59 // or other H(div) conforming elements such as FE_ABF<dim> and
60 // FE_BDM<dim>.
61 template <int spacedim>
62 void
63 get_dof_sign_change_h_div(
64 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
66 const std::vector<MappingKind> &mapping_kind,
67 std::vector<double> &face_sign)
68 {
69 const unsigned int dim = 2;
70 // const unsigned int spacedim = 2;
71
72 const CellId this_cell_id = cell->id();
73
74 for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
75 {
76 typename ::Triangulation<dim, spacedim>::face_iterator face =
77 cell->face(f);
78 if (!face->at_boundary())
79 {
80 const unsigned int nn = cell->neighbor_face_no(f);
81 const typename ::Triangulation<dim,
82 spacedim>::cell_iterator
83 neighbor_cell_at_face = cell->neighbor(f);
84 const CellId neighbor_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 < neighbor_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 template <int spacedim>
138 void
139 get_dof_sign_change_nedelec(
140 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
142 const std::vector<MappingKind> &mapping_kind,
143 std::vector<double> &line_dof_sign)
144 {
145 // The Nedelec finite elements in two spacial dimensions have two types
146 // of dofs: the line dofs and the quad dofs. The line dofs are
147 // associated with the edges. They are shared between the neighbouring
148 // cells and, consequently, need to be adjusted to compensate for a
149 // possible mismatch in the edge orientation. The quad dofs, they are
150 // associated with the interiors of the cells, are not shared between
151 // the cells in two spacial dimensions and need no adjustments.
152 //
153 // The Nedelec finite elements in two spacial dimensions have
154 // 2*(k+1)*(k+2) dofs per cell. All dofs are distributed between the
155 // line and quad dofs as the following:
156 //
157 // 4*(k+1) line dofs; (k+1) dofs per line.
158 // 2*k*(k+1) quad dofs.
159 //
160 // The dofs are indexed in the following order: first all line dofs,
161 // then all quad dofs.
162 //
163 // Here we adjust the line dofs. The sign of a line dof needs to be
164 // changed if the edge on which the dof resides points in the opposite
165 // direction.
166 const unsigned int k = fe.tensor_degree() - 1;
167
168 for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell; ++l)
169 if (cell->line_orientation(l) !=
171 mapping_kind[0] == mapping_nedelec)
172 {
173 if (k == 0)
174 {
175 // The lowest order element (k=0) is straightforward, because
176 // there is a single dof per edge, which needs to be flipped:
177 line_dof_sign[l] = -1.0;
178 }
179 else
180 {
181 // The case k > 0 is a bit more complicated. As we adjust
182 // only edge dofs in this function, we need to concern
183 // ourselves with the first 4*(k+1) entries in line_dof_sign
184 // vector ignoring the rest. There are (k+1) dofs per edge.
185 // Let us consider the local dof indices on one edge,
186 // local_line_dof = 0...k. The shape functions with even
187 // indices are asymmetric. The corresponding dofs need sign
188 // adjustment if the edge points in the opposite direction.
189 // The shape functions with odd indices are symmetric. The
190 // corresponding dofs need no sign adjustment even if the edge
191 // points in the opposite direction. In the current context
192 // the notion of symmetry of a shape function means that a
193 // shape function looks exactly the same if it is looked upon
194 // from the centers of the two neighbouring cells that share
195 // it.
196 for (unsigned int local_line_dof = 0;
197 local_line_dof < (k + 1);
198 local_line_dof++)
199 if (local_line_dof % 2 == 0)
200 line_dof_sign[local_line_dof + l * (k + 1)] = -1.0;
201 }
202 }
203 }
204
205 template <int spacedim>
206 void
207 get_dof_sign_change_nedelec(
208 const typename ::Triangulation<3, spacedim>::cell_iterator &cell,
210 const std::vector<MappingKind> &mapping_kind,
211 std::vector<double> &line_dof_sign)
212 {
213 // This function does half of the job - it adjusts the sign of the
214 // line (edge) dofs. In the three-dimensional space the quad (face) dofs
215 // need to be adjusted as well. The quad dofs are treated by
216 // FE_Nedelec<dim>::initialize_quad_dof_index_permutation_and_sign_change()
217 // in fe_nedelec.cc. The dofs associated with the interior of the cells,
218 // the hex dofs, need no adjustments as they are not shared between the
219 // neighboring cells.
220
221 const unsigned int k = fe.tensor_degree() - 1;
222 // The order of the Nedelec elements equals the tensor degree minus one,
223 // k = n - 1. In the three-dimensional space the Nedelec elements of the
224 // lowermost order, k = 0, have only 12 line (edge) dofs. The Nedelec
225 // elements of the higher orders, k > 0, have 3*(k+1)*(k+2)^2 dofs in
226 // total if dim=3. The dofs in a cell are distributed between lines
227 // (edges), quads (faces), and the hex (the interior of the cell) as the
228 // following:
229 //
230 // 12*(k+1) line dofs; (k+1) dofs per line.
231 // 2*6*k*(k+1) quad dofs; 2*k*(k+1) dofs per quad.
232 // 3*(k+2)^2*(k+1) hex dofs.
233 //
234 // The dofs are indexed in the following order: first all line dofs,
235 // then all quad dofs, and then all hex dofs.
236 //
237 // Here we adjust only the line (edge) dofs. The line dofs need only
238 // sign adjustment. That is, no permutation of the line dofs is needed.
239 for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
240 if (cell->line_orientation(l) !=
242 mapping_kind[0] == mapping_nedelec)
243 {
244 if (k == 0)
245 {
246 // The lowest order element (k=0) is straightforward, because
247 // there is a single dof per edge, which needs to be flipped:
248 line_dof_sign[l] = -1.0;
249 }
250 else
251 {
252 // The case k > 0 is a bit more complicated. As we adjust
253 // only edge dofs in this function, we need to concern
254 // ourselves with the first 12*(k+1) entries in line_dof_sign
255 // vector ignoring the rest. There are (k+1) dofs per edge.
256 // Let us consider the local dof indices on one edge,
257 // local_line_dof = 0...k. The shape functions with even
258 // indices are asymmetric. The corresponding dofs need sign
259 // adjustment if the edge points in the opposite direction.
260 // The shape functions with odd indices are symmetric. The
261 // corresponding dofs need no sign adjustment even if the edge
262 // points in the opposite direction. In the current context
263 // the notion of symmetry of a shape function means that a
264 // shape function looks exactly the same if it is looked upon
265 // from the centers of the two neighbouring cells that share
266 // it.
267 for (unsigned int local_line_dof = 0;
268 local_line_dof < (k + 1);
269 local_line_dof++)
270 if (local_line_dof % 2 == 0)
271 line_dof_sign[local_line_dof + l * (k + 1)] = -1.0;
272 }
273 }
274 }
275 } // namespace
276 } // namespace FE_PolyTensor
277} // namespace internal
278
279#ifndef DOXYGEN
280
281template <int dim, int spacedim>
283 const TensorPolynomialsBase<dim> &polynomials,
284 const FiniteElementData<dim> &fe_data,
285 const std::vector<bool> &restriction_is_additive_flags,
286 const std::vector<ComponentMask> &nonzero_components)
287 : FiniteElement<dim, spacedim>(fe_data,
288 restriction_is_additive_flags,
289 nonzero_components)
290 , mapping_kind({MappingKind::mapping_none})
291 , poly_space(polynomials.clone())
292{
293 cached_point[0] = -1;
294 // Set up the table converting
295 // components to base
296 // components. Since we have only
297 // one base element, everything
298 // remains zero except the
299 // component in the base, which is
300 // the component itself
301 for (unsigned int comp = 0; comp < this->n_components(); ++comp)
302 this->component_to_base_table[comp].first.second = comp;
303
304 if (dim == 3)
305 {
306 adjust_quad_dof_sign_for_face_orientation_table.resize(
307 this->n_unique_2d_subobjects());
308
309 for (unsigned int f = 0; f < this->n_unique_2d_subobjects(); ++f)
310 {
311 adjust_quad_dof_sign_for_face_orientation_table[f] =
312 Table<2, bool>(this->n_dofs_per_quad(f),
313 this->reference_cell().n_face_orientations(f));
314 adjust_quad_dof_sign_for_face_orientation_table[f].fill(false);
315 }
316 }
317}
318
319
320
321template <int dim, int spacedim>
323 : FiniteElement<dim, spacedim>(fe)
324 , mapping_kind(fe.mapping_kind)
325 , adjust_quad_dof_sign_for_face_orientation_table(
326 fe.adjust_quad_dof_sign_for_face_orientation_table)
327 , poly_space(fe.poly_space->clone())
328 , inverse_node_matrix(fe.inverse_node_matrix)
329{}
330
331
332
333template <int dim, int spacedim>
334bool
336{
337 return mapping_kind.size() == 1;
338}
339
340
341template <int dim, int spacedim>
342bool
344 const unsigned int index,
345 const unsigned int face,
346 const bool face_orientation,
347 const bool face_flip,
348 const bool face_rotation) const
349{
350 // do nothing in 1d and 2d
351 if (dim < 3)
352 return false;
353
354 // The exception are discontinuous
355 // elements for which there should be no
356 // face dofs anyway (i.e. dofs_per_quad==0
357 // in 3d), so we don't need the table, but
358 // the function should also not have been
359 // called
360 AssertIndexRange(index, this->n_dofs_per_quad(face));
361 Assert(adjust_quad_dof_sign_for_face_orientation_table
362 [this->n_unique_2d_subobjects() == 1 ? 0 : face]
363 .n_elements() ==
364 this->reference_cell().n_face_orientations(face) *
365 this->n_dofs_per_quad(face),
367
368 return adjust_quad_dof_sign_for_face_orientation_table
369 [this->n_unique_2d_subobjects() == 1 ? 0 : face](
370 index,
372 face_rotation,
373 face_flip));
374}
375
376
377template <int dim, int spacedim>
379FE_PolyTensor<dim, spacedim>::get_mapping_kind(const unsigned int i) const
380{
381 if (single_mapping_kind())
382 return mapping_kind[0];
383
384 AssertIndexRange(i, mapping_kind.size());
385 return mapping_kind[i];
386}
387
388
389
390template <int dim, int spacedim>
391double
393 const Point<dim> &) const
394
395{
397 return 0.;
398}
399
400
401
402template <int dim, int spacedim>
403double
405 const unsigned int i,
406 const Point<dim> &p,
407 const unsigned int component) const
408{
409 AssertIndexRange(i, this->n_dofs_per_cell());
410 AssertIndexRange(component, dim);
411
412 std::lock_guard<std::mutex> lock(cache_mutex);
413
414 if (cached_point != p || cached_values.empty())
415 {
416 cached_point = p;
417 cached_values.resize(poly_space->n());
418
419 std::vector<Tensor<4, dim>> dummy1;
420 std::vector<Tensor<5, dim>> dummy2;
421 poly_space->evaluate(
422 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
423 }
424
425 double s = 0;
426 if (inverse_node_matrix.n_cols() == 0)
427 return cached_values[i][component];
428 else
429 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
430 s += inverse_node_matrix(j, i) * cached_values[j][component];
431 return s;
432}
433
434
435
436template <int dim, int spacedim>
439 const Point<dim> &) const
440{
442 return Tensor<1, dim>();
443}
444
445
446
447template <int dim, int spacedim>
450 const unsigned int i,
451 const Point<dim> &p,
452 const unsigned int component) const
453{
454 AssertIndexRange(i, this->n_dofs_per_cell());
455 AssertIndexRange(component, dim);
456
457 std::lock_guard<std::mutex> lock(cache_mutex);
458
459 if (cached_point != p || cached_grads.empty())
460 {
461 cached_point = p;
462 cached_grads.resize(poly_space->n());
463
464 std::vector<Tensor<4, dim>> dummy1;
465 std::vector<Tensor<5, dim>> dummy2;
466 poly_space->evaluate(
467 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
468 }
469
471 if (inverse_node_matrix.n_cols() == 0)
472 return cached_grads[i][component];
473 else
474 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
475 s += inverse_node_matrix(j, i) * cached_grads[j][component];
476
477 return s;
478}
479
480
481
482template <int dim, int spacedim>
485 const Point<dim> &) const
486{
488 return Tensor<2, dim>();
489}
490
491
492
493template <int dim, int spacedim>
496 const unsigned int i,
497 const Point<dim> &p,
498 const unsigned int component) const
499{
500 AssertIndexRange(i, this->n_dofs_per_cell());
501 AssertIndexRange(component, dim);
502
503 std::lock_guard<std::mutex> lock(cache_mutex);
504
505 if (cached_point != p || cached_grad_grads.empty())
506 {
507 cached_point = p;
508 cached_grad_grads.resize(poly_space->n());
509
510 std::vector<Tensor<4, dim>> dummy1;
511 std::vector<Tensor<5, dim>> dummy2;
512 poly_space->evaluate(
513 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
514 }
515
517 if (inverse_node_matrix.n_cols() == 0)
518 return cached_grad_grads[i][component];
519 else
520 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
521 s += inverse_node_matrix(i, j) * cached_grad_grads[j][component];
522
523 return s;
524}
525
526
527//---------------------------------------------------------------------------
528// Fill data of FEValues
529//---------------------------------------------------------------------------
530
531template <int dim, int spacedim>
532void
535 const CellSimilarity::Similarity cell_similarity,
536 const Quadrature<dim> &quadrature,
537 const Mapping<dim, spacedim> &mapping,
538 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
540 &mapping_data,
541 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
543 spacedim>
544 &output_data) const
545{
546 // convert data object to internal
547 // data for this class. fails with
548 // an exception if that is not
549 // possible
550 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
552 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
553
554 const unsigned int n_q_points = quadrature.size();
555
556 Assert(!(fe_data.update_each & update_values) ||
557 fe_data.shape_values.size()[0] == this->n_dofs_per_cell(),
558 ExcDimensionMismatch(fe_data.shape_values.size()[0],
559 this->n_dofs_per_cell()));
560 Assert(!(fe_data.update_each & update_values) ||
561 fe_data.shape_values.size()[1] == n_q_points,
562 ExcDimensionMismatch(fe_data.shape_values.size()[1], n_q_points));
563
564 // TODO: The dof_sign_change only affects Nedelec elements and is not the
565 // correct thing on complicated meshes for higher order Nedelec elements.
566 // Something similar to FE_Q should be done to permute dofs and to change the
567 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
568 // class) is preferable.
569 std::fill(fe_data.dof_sign_change.begin(),
570 fe_data.dof_sign_change.end(),
571 1.0);
572 if (fe_data.update_each & update_values)
573 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
574 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
575
576 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
577 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
578 // Also nothing in 3d since we take care of it by using the
579 // adjust_quad_dof_sign_for_face_orientation_table.
580 if (fe_data.update_each & update_values)
581 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
582 *this,
583 this->mapping_kind,
584 fe_data.dof_sign_change);
585
586 // What is the first dof_index on a quad?
587 const unsigned int first_quad_index = this->get_first_quad_index();
588 // How many dofs per quad and how many quad dofs do we have at all?
589 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
590 const unsigned int n_quad_dofs =
591 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
592
593 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
594 ++dof_index)
595 {
596 /*
597 * This assumes that the dofs are ordered by first vertices, lines, quads
598 * and volume dofs. Note that in 2d this always gives false.
599 */
600 const bool is_quad_dof =
601 (dim == 2 ? false :
602 (first_quad_index <= dof_index) &&
603 (dof_index < first_quad_index + n_quad_dofs));
604
605 // TODO: This hack is not pretty and it is only here to handle the 2d
606 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
607 // handled by the
608 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
609 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
610 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
611 // it is handles with a table. This array is allocated in
612 // fe_poly_tensor.h.
613 double dof_sign = 1.0;
614 // under some circumstances fe_data.dof_sign_change is not allocated
615 if (fe_data.update_each & update_values)
616 dof_sign = fe_data.dof_sign_change[dof_index];
617
618 if (is_quad_dof)
619 {
620 /*
621 * Find the face belonging to this dof_index. This is integer
622 * division.
623 */
624 const unsigned int face_index_from_dof_index =
625 (dof_index - first_quad_index) / (n_dofs_per_quad);
626
627 const unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
628
629 // Correct the dof_sign if necessary
630 if (adjust_quad_dof_sign_for_face_orientation(
631 local_quad_dof_index,
632 face_index_from_dof_index,
633 cell->face_orientation(face_index_from_dof_index),
634 cell->face_flip(face_index_from_dof_index),
635 cell->face_rotation(face_index_from_dof_index)))
636 dof_sign = -1.0;
637 }
638
639 const MappingKind mapping_kind = get_mapping_kind(dof_index);
640
641 const unsigned int first =
642 output_data.shape_function_to_row_table
643 [dof_index * this->n_components() +
644 this->get_nonzero_components(dof_index).first_selected_component()];
645
646 // update the shape function values as necessary
647 //
648 // we only need to do this if the current cell is not a translation of
649 // the previous one; or, even if it is a translation, if we use
650 // mappings other than the standard mappings that require us to
651 // recompute values and derivatives because of possible sign changes
652 if (fe_data.update_each & update_values &&
653 ((cell_similarity != CellSimilarity::translation) ||
654 ((mapping_kind == mapping_piola) ||
655 (mapping_kind == mapping_raviart_thomas) ||
656 (mapping_kind == mapping_nedelec))))
657 {
658 switch (mapping_kind)
659 {
660 case mapping_none:
661 {
662 for (unsigned int k = 0; k < n_q_points; ++k)
663 for (unsigned int d = 0; d < dim; ++d)
664 output_data.shape_values(first + d, k) =
665 fe_data.shape_values[dof_index][k][d];
666 break;
667 }
668
671 {
672 mapping.transform(
673 make_array_view(fe_data.shape_values, dof_index),
674 mapping_kind,
675 mapping_internal,
676 make_array_view(fe_data.transformed_shape_values));
677
678 for (unsigned int k = 0; k < n_q_points; ++k)
679 for (unsigned int d = 0; d < dim; ++d)
680 output_data.shape_values(first + d, k) =
681 fe_data.transformed_shape_values[k][d];
682
683 break;
684 }
685
687 case mapping_piola:
688 {
689 mapping.transform(
690 make_array_view(fe_data.shape_values, dof_index),
692 mapping_internal,
693 make_array_view(fe_data.transformed_shape_values));
694 for (unsigned int k = 0; k < n_q_points; ++k)
695 for (unsigned int d = 0; d < dim; ++d)
696 output_data.shape_values(first + d, k) =
697 dof_sign * fe_data.transformed_shape_values[k][d];
698 break;
699 }
700
701 case mapping_nedelec:
702 {
703 mapping.transform(
704 make_array_view(fe_data.shape_values, dof_index),
706 mapping_internal,
707 make_array_view(fe_data.transformed_shape_values));
708
709 for (unsigned int k = 0; k < n_q_points; ++k)
710 for (unsigned int d = 0; d < dim; ++d)
711 output_data.shape_values(first + d, k) =
712 dof_sign * fe_data.transformed_shape_values[k][d];
713
714 break;
715 }
716
717 default:
719 }
720 }
721
722 // update gradients. apply the same logic as above
723 if (fe_data.update_each & update_gradients &&
724 ((cell_similarity != CellSimilarity::translation) ||
725 ((mapping_kind == mapping_piola) ||
726 (mapping_kind == mapping_raviart_thomas) ||
727 (mapping_kind == mapping_nedelec))))
728
729 {
730 switch (mapping_kind)
731 {
732 case mapping_none:
733 {
734 mapping.transform(
735 make_array_view(fe_data.shape_grads, dof_index),
737 mapping_internal,
738 make_array_view(fe_data.transformed_shape_grads));
739 for (unsigned int k = 0; k < n_q_points; ++k)
740 for (unsigned int d = 0; d < dim; ++d)
741 output_data.shape_gradients[first + d][k] =
742 fe_data.transformed_shape_grads[k][d];
743 break;
744 }
746 {
747 mapping.transform(
748 make_array_view(fe_data.shape_grads, dof_index),
750 mapping_internal,
751 make_array_view(fe_data.transformed_shape_grads));
752
753 for (unsigned int k = 0; k < n_q_points; ++k)
754 for (unsigned int d = 0; d < spacedim; ++d)
755 for (unsigned int n = 0; n < spacedim; ++n)
756 fe_data.transformed_shape_grads[k][d] -=
757 output_data.shape_values(first + n, k) *
758 mapping_data.jacobian_pushed_forward_grads[k][n][d];
759
760 for (unsigned int k = 0; k < n_q_points; ++k)
761 for (unsigned int d = 0; d < dim; ++d)
762 output_data.shape_gradients[first + d][k] =
763 fe_data.transformed_shape_grads[k][d];
764
765 break;
766 }
768 {
769 for (unsigned int k = 0; k < n_q_points; ++k)
770 fe_data.untransformed_shape_grads[k] =
771 fe_data.shape_grads[dof_index][k];
772 mapping.transform(
773 make_array_view(fe_data.untransformed_shape_grads),
775 mapping_internal,
776 make_array_view(fe_data.transformed_shape_grads));
777
778 for (unsigned int k = 0; k < n_q_points; ++k)
779 for (unsigned int d = 0; d < spacedim; ++d)
780 for (unsigned int n = 0; n < spacedim; ++n)
781 fe_data.transformed_shape_grads[k][d] +=
782 output_data.shape_values(first + n, k) *
783 mapping_data.jacobian_pushed_forward_grads[k][d][n];
784
785
786 for (unsigned int k = 0; k < n_q_points; ++k)
787 for (unsigned int d = 0; d < dim; ++d)
788 output_data.shape_gradients[first + d][k] =
789 fe_data.transformed_shape_grads[k][d];
790
791 break;
792 }
794 case mapping_piola:
795 {
796 for (unsigned int k = 0; k < n_q_points; ++k)
797 fe_data.untransformed_shape_grads[k] =
798 fe_data.shape_grads[dof_index][k];
799 mapping.transform(
800 make_array_view(fe_data.untransformed_shape_grads),
802 mapping_internal,
803 make_array_view(fe_data.transformed_shape_grads));
804
805 for (unsigned int k = 0; k < n_q_points; ++k)
806 for (unsigned int d = 0; d < spacedim; ++d)
807 for (unsigned int n = 0; n < spacedim; ++n)
808 fe_data.transformed_shape_grads[k][d] +=
809 (output_data.shape_values(first + n, k) *
810 mapping_data
812 (output_data.shape_values(first + d, k) *
813 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
814
815 for (unsigned int k = 0; k < n_q_points; ++k)
816 for (unsigned int d = 0; d < dim; ++d)
817 output_data.shape_gradients[first + d][k] =
818 dof_sign * fe_data.transformed_shape_grads[k][d];
819
820 break;
821 }
822
823 case mapping_nedelec:
824 {
825 // treat the gradients of
826 // this particular shape
827 // function at all
828 // q-points. if Dv is the
829 // gradient of the shape
830 // function on the unit
831 // cell, then
832 // (J^-T)Dv(J^-1) is the
833 // value we want to have on
834 // the real cell.
835 for (unsigned int k = 0; k < n_q_points; ++k)
836 fe_data.untransformed_shape_grads[k] =
837 fe_data.shape_grads[dof_index][k];
838
839 mapping.transform(
840 make_array_view(fe_data.untransformed_shape_grads),
842 mapping_internal,
843 make_array_view(fe_data.transformed_shape_grads));
844
845 for (unsigned int k = 0; k < n_q_points; ++k)
846 for (unsigned int d = 0; d < spacedim; ++d)
847 for (unsigned int n = 0; n < spacedim; ++n)
848 fe_data.transformed_shape_grads[k][d] -=
849 output_data.shape_values(first + n, k) *
850 mapping_data.jacobian_pushed_forward_grads[k][n][d];
851
852 for (unsigned int k = 0; k < n_q_points; ++k)
853 for (unsigned int d = 0; d < dim; ++d)
854 output_data.shape_gradients[first + d][k] =
855 dof_sign * fe_data.transformed_shape_grads[k][d];
856
857 break;
858 }
859
860 default:
862 }
863 }
864
865 // update hessians. apply the same logic as above
866 if (fe_data.update_each & update_hessians &&
867 ((cell_similarity != CellSimilarity::translation) ||
868 ((mapping_kind == mapping_piola) ||
869 (mapping_kind == mapping_raviart_thomas) ||
870 (mapping_kind == mapping_nedelec))))
871
872 {
873 switch (mapping_kind)
874 {
875 case mapping_none:
876 {
877 mapping.transform(
878 make_array_view(fe_data.shape_grad_grads, dof_index),
880 mapping_internal,
881 make_array_view(fe_data.transformed_shape_hessians));
882
883 for (unsigned int k = 0; k < n_q_points; ++k)
884 for (unsigned int d = 0; d < spacedim; ++d)
885 for (unsigned int n = 0; n < spacedim; ++n)
886 fe_data.transformed_shape_hessians[k][d] -=
887 output_data.shape_gradients[first + d][k][n] *
888 mapping_data.jacobian_pushed_forward_grads[k][n];
889
890 for (unsigned int k = 0; k < n_q_points; ++k)
891 for (unsigned int d = 0; d < dim; ++d)
892 output_data.shape_hessians[first + d][k] =
893 fe_data.transformed_shape_hessians[k][d];
894
895 break;
896 }
898 {
899 for (unsigned int k = 0; k < n_q_points; ++k)
900 fe_data.untransformed_shape_hessian_tensors[k] =
901 fe_data.shape_grad_grads[dof_index][k];
902
903 mapping.transform(
905 fe_data.untransformed_shape_hessian_tensors),
907 mapping_internal,
908 make_array_view(fe_data.transformed_shape_hessians));
909
910 for (unsigned int k = 0; k < n_q_points; ++k)
911 for (unsigned int d = 0; d < spacedim; ++d)
912 for (unsigned int n = 0; n < spacedim; ++n)
913 for (unsigned int i = 0; i < spacedim; ++i)
914 for (unsigned int j = 0; j < spacedim; ++j)
915 {
916 fe_data.transformed_shape_hessians[k][d][i][j] -=
917 (output_data.shape_values(first + n, k) *
918 mapping_data
920 [k][n][d][i][j]) +
921 (output_data.shape_gradients[first + d][k][n] *
922 mapping_data
923 .jacobian_pushed_forward_grads[k][n][i][j]) +
924 (output_data.shape_gradients[first + n][k][i] *
925 mapping_data
926 .jacobian_pushed_forward_grads[k][n][d][j]) +
927 (output_data.shape_gradients[first + n][k][j] *
928 mapping_data
929 .jacobian_pushed_forward_grads[k][n][i][d]);
930 }
931
932 for (unsigned int k = 0; k < n_q_points; ++k)
933 for (unsigned int d = 0; d < dim; ++d)
934 output_data.shape_hessians[first + d][k] =
935 fe_data.transformed_shape_hessians[k][d];
936
937 break;
938 }
940 {
941 for (unsigned int k = 0; k < n_q_points; ++k)
942 fe_data.untransformed_shape_hessian_tensors[k] =
943 fe_data.shape_grad_grads[dof_index][k];
944
945 mapping.transform(
947 fe_data.untransformed_shape_hessian_tensors),
949 mapping_internal,
950 make_array_view(fe_data.transformed_shape_hessians));
951
952 for (unsigned int k = 0; k < n_q_points; ++k)
953 for (unsigned int d = 0; d < spacedim; ++d)
954 for (unsigned int n = 0; n < spacedim; ++n)
955 for (unsigned int i = 0; i < spacedim; ++i)
956 for (unsigned int j = 0; j < spacedim; ++j)
957 {
958 fe_data.transformed_shape_hessians[k][d][i][j] +=
959 (output_data.shape_values(first + n, k) *
960 mapping_data
962 [k][d][n][i][j]) +
963 (output_data.shape_gradients[first + n][k][i] *
964 mapping_data
965 .jacobian_pushed_forward_grads[k][d][n][j]) +
966 (output_data.shape_gradients[first + n][k][j] *
967 mapping_data
968 .jacobian_pushed_forward_grads[k][d][i][n]) -
969 (output_data.shape_gradients[first + d][k][n] *
970 mapping_data
971 .jacobian_pushed_forward_grads[k][n][i][j]);
972 for (unsigned int m = 0; m < spacedim; ++m)
973 fe_data
974 .transformed_shape_hessians[k][d][i][j] -=
975 (mapping_data
976 .jacobian_pushed_forward_grads[k][d][i]
977 [m] *
978 mapping_data
979 .jacobian_pushed_forward_grads[k][m][n]
980 [j] *
981 output_data.shape_values(first + n, k)) +
982 (mapping_data
983 .jacobian_pushed_forward_grads[k][d][m]
984 [j] *
985 mapping_data
986 .jacobian_pushed_forward_grads[k][m][i]
987 [n] *
988 output_data.shape_values(first + n, k));
989 }
990
991 for (unsigned int k = 0; k < n_q_points; ++k)
992 for (unsigned int d = 0; d < dim; ++d)
993 output_data.shape_hessians[first + d][k] =
994 fe_data.transformed_shape_hessians[k][d];
995
996 break;
997 }
999 case mapping_piola:
1000 {
1001 for (unsigned int k = 0; k < n_q_points; ++k)
1002 fe_data.untransformed_shape_hessian_tensors[k] =
1003 fe_data.shape_grad_grads[dof_index][k];
1004
1005 mapping.transform(
1007 fe_data.untransformed_shape_hessian_tensors),
1009 mapping_internal,
1010 make_array_view(fe_data.transformed_shape_hessians));
1011
1012 for (unsigned int k = 0; k < n_q_points; ++k)
1013 for (unsigned int d = 0; d < spacedim; ++d)
1014 for (unsigned int n = 0; n < spacedim; ++n)
1015 for (unsigned int i = 0; i < spacedim; ++i)
1016 for (unsigned int j = 0; j < spacedim; ++j)
1017 {
1018 fe_data.transformed_shape_hessians[k][d][i][j] +=
1019 (output_data.shape_values(first + n, k) *
1020 mapping_data
1022 [k][d][n][i][j]) +
1023 (output_data.shape_gradients[first + n][k][i] *
1024 mapping_data
1025 .jacobian_pushed_forward_grads[k][d][n][j]) +
1026 (output_data.shape_gradients[first + n][k][j] *
1027 mapping_data
1028 .jacobian_pushed_forward_grads[k][d][i][n]) -
1029 (output_data.shape_gradients[first + d][k][n] *
1030 mapping_data
1031 .jacobian_pushed_forward_grads[k][n][i][j]);
1032
1033 fe_data.transformed_shape_hessians[k][d][i][j] -=
1034 (output_data.shape_values(first + d, k) *
1035 mapping_data
1037 [k][n][n][i][j]) +
1038 (output_data.shape_gradients[first + d][k][i] *
1039 mapping_data
1040 .jacobian_pushed_forward_grads[k][n][n][j]) +
1041 (output_data.shape_gradients[first + d][k][j] *
1042 mapping_data
1043 .jacobian_pushed_forward_grads[k][n][n][i]);
1044
1045 for (unsigned int m = 0; m < spacedim; ++m)
1046 {
1047 fe_data
1048 .transformed_shape_hessians[k][d][i][j] -=
1049 (mapping_data
1051 [m] *
1052 mapping_data
1054 [j] *
1055 output_data.shape_values(first + n, k)) +
1056 (mapping_data
1057 .jacobian_pushed_forward_grads[k][d][m]
1058 [j] *
1059 mapping_data
1060 .jacobian_pushed_forward_grads[k][m][i]
1061 [n] *
1062 output_data.shape_values(first + n, k));
1063
1064 fe_data
1065 .transformed_shape_hessians[k][d][i][j] +=
1066 (mapping_data
1068 [m] *
1069 mapping_data
1071 [j] *
1072 output_data.shape_values(first + d, k)) +
1073 (mapping_data
1074 .jacobian_pushed_forward_grads[k][n][m]
1075 [j] *
1076 mapping_data
1077 .jacobian_pushed_forward_grads[k][m][i]
1078 [n] *
1079 output_data.shape_values(first + d, k));
1080 }
1081 }
1082
1083 for (unsigned int k = 0; k < n_q_points; ++k)
1084 for (unsigned int d = 0; d < dim; ++d)
1085 output_data.shape_hessians[first + d][k] =
1086 dof_sign * fe_data.transformed_shape_hessians[k][d];
1087
1088 break;
1089 }
1090
1091 case mapping_nedelec:
1092 {
1093 for (unsigned int k = 0; k < n_q_points; ++k)
1094 fe_data.untransformed_shape_hessian_tensors[k] =
1095 fe_data.shape_grad_grads[dof_index][k];
1096
1097 mapping.transform(
1099 fe_data.untransformed_shape_hessian_tensors),
1101 mapping_internal,
1102 make_array_view(fe_data.transformed_shape_hessians));
1103
1104 for (unsigned int k = 0; k < n_q_points; ++k)
1105 for (unsigned int d = 0; d < spacedim; ++d)
1106 for (unsigned int n = 0; n < spacedim; ++n)
1107 for (unsigned int i = 0; i < spacedim; ++i)
1108 for (unsigned int j = 0; j < spacedim; ++j)
1109 {
1110 fe_data.transformed_shape_hessians[k][d][i][j] -=
1111 (output_data.shape_values(first + n, k) *
1112 mapping_data
1114 [k][n][d][i][j]) +
1115 (output_data.shape_gradients[first + d][k][n] *
1116 mapping_data
1117 .jacobian_pushed_forward_grads[k][n][i][j]) +
1118 (output_data.shape_gradients[first + n][k][i] *
1119 mapping_data
1120 .jacobian_pushed_forward_grads[k][n][d][j]) +
1121 (output_data.shape_gradients[first + n][k][j] *
1122 mapping_data
1123 .jacobian_pushed_forward_grads[k][n][i][d]);
1124 }
1125
1126 for (unsigned int k = 0; k < n_q_points; ++k)
1127 for (unsigned int d = 0; d < dim; ++d)
1128 output_data.shape_hessians[first + d][k] =
1129 dof_sign * fe_data.transformed_shape_hessians[k][d];
1130
1131 break;
1132 }
1133
1134 default:
1136 }
1137 }
1138
1139 // third derivatives are not implemented
1140 if (fe_data.update_each & update_3rd_derivatives &&
1141 ((cell_similarity != CellSimilarity::translation) ||
1142 ((mapping_kind == mapping_piola) ||
1143 (mapping_kind == mapping_raviart_thomas) ||
1144 (mapping_kind == mapping_nedelec))))
1145 {
1147 }
1148 }
1149}
1150
1151
1152
1153template <int dim, int spacedim>
1154void
1157 const unsigned int face_no,
1158 const hp::QCollection<dim - 1> &quadrature,
1159 const Mapping<dim, spacedim> &mapping,
1160 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1162 &mapping_data,
1163 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1165 spacedim>
1166 &output_data) const
1167{
1168 AssertDimension(quadrature.size(), 1);
1169
1170 // convert data object to internal
1171 // data for this class. fails with
1172 // an exception if that is not
1173 // possible
1174 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1176 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1177
1178 const unsigned int n_q_points = quadrature[0].size();
1179 // offset determines which data set
1180 // to take (all data sets for all
1181 // faces are stored contiguously)
1182
1183 const auto offset =
1185 face_no,
1186 cell->combined_face_orientation(
1187 face_no),
1188 n_q_points);
1189
1190 // TODO: Size assertions
1191
1192 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1193 // correct thing on complicated meshes for higher order Nedelec elements.
1194 // Something similar to FE_Q should be done to permute dofs and to change the
1195 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1196 // class) is preferable.
1197 std::fill(fe_data.dof_sign_change.begin(),
1198 fe_data.dof_sign_change.end(),
1199 1.0);
1200 if (fe_data.update_each & update_values)
1201 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
1202 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
1203
1204 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1205 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1206 // Also nothing in 3d since we take care of it by using the
1207 // adjust_quad_dof_sign_for_face_orientation_table.
1208 if (fe_data.update_each & update_values)
1209 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1210 *this,
1211 this->mapping_kind,
1212 fe_data.dof_sign_change);
1213
1214 // What is the first dof_index on a quad?
1215 const unsigned int first_quad_index = this->get_first_quad_index();
1216 // How many dofs per quad and how many quad dofs do we have at all?
1217 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1218 const unsigned int n_quad_dofs =
1219 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1220
1221 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1222 ++dof_index)
1223 {
1224 /*
1225 * This assumes that the dofs are ordered by first vertices, lines, quads
1226 * and volume dofs. Note that in 2d this always gives false.
1227 */
1228 const bool is_quad_dof =
1229 (dim == 2 ? false :
1230 (first_quad_index <= dof_index) &&
1231 (dof_index < first_quad_index + n_quad_dofs));
1232
1233 // TODO: This hack is not pretty and it is only here to handle the 2d
1234 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1235 // handled by the
1236 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1237 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1238 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1239 // it is handles with a table. This array is allocated in
1240 // fe_poly_tensor.h.
1241 double dof_sign = 1.0;
1242 // under some circumstances fe_data.dof_sign_change is not allocated
1243 if (fe_data.update_each & update_values)
1244 dof_sign = fe_data.dof_sign_change[dof_index];
1245
1246 if (is_quad_dof)
1247 {
1248 /*
1249 * Find the face belonging to this dof_index. This is integer
1250 * division.
1251 */
1252 unsigned int face_index_from_dof_index =
1253 (dof_index - first_quad_index) / (n_dofs_per_quad);
1254
1255 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1256
1257 // Correct the dof_sign if necessary
1258 if (adjust_quad_dof_sign_for_face_orientation(
1259 local_quad_dof_index,
1260 face_index_from_dof_index,
1261 cell->face_orientation(face_index_from_dof_index),
1262 cell->face_flip(face_index_from_dof_index),
1263 cell->face_rotation(face_index_from_dof_index)))
1264 dof_sign = -1.0;
1265 }
1266
1267 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1268
1269 const unsigned int first =
1270 output_data.shape_function_to_row_table
1271 [dof_index * this->n_components() +
1272 this->get_nonzero_components(dof_index).first_selected_component()];
1273
1274 if (fe_data.update_each & update_values)
1275 {
1276 switch (mapping_kind)
1277 {
1278 case mapping_none:
1279 {
1280 for (unsigned int k = 0; k < n_q_points; ++k)
1281 for (unsigned int d = 0; d < dim; ++d)
1282 output_data.shape_values(first + d, k) =
1283 fe_data.shape_values[dof_index][k + offset][d];
1284 break;
1285 }
1286
1287 case mapping_covariant:
1289 {
1291 transformed_shape_values =
1292 make_array_view(fe_data.transformed_shape_values,
1293 offset,
1294 n_q_points);
1295 mapping.transform(make_array_view(fe_data.shape_values,
1296 dof_index,
1297 offset,
1298 n_q_points),
1299 mapping_kind,
1300 mapping_internal,
1301 transformed_shape_values);
1302
1303 for (unsigned int k = 0; k < n_q_points; ++k)
1304 for (unsigned int d = 0; d < dim; ++d)
1305 output_data.shape_values(first + d, k) =
1306 transformed_shape_values[k][d];
1307
1308 break;
1309 }
1311 case mapping_piola:
1312 {
1314 transformed_shape_values =
1315 make_array_view(fe_data.transformed_shape_values,
1316 offset,
1317 n_q_points);
1318 mapping.transform(make_array_view(fe_data.shape_values,
1319 dof_index,
1320 offset,
1321 n_q_points),
1323 mapping_internal,
1324 transformed_shape_values);
1325 for (unsigned int k = 0; k < n_q_points; ++k)
1326 for (unsigned int d = 0; d < dim; ++d)
1327 output_data.shape_values(first + d, k) =
1328 dof_sign * transformed_shape_values[k][d];
1329 break;
1330 }
1331
1332 case mapping_nedelec:
1333 {
1335 transformed_shape_values =
1336 make_array_view(fe_data.transformed_shape_values,
1337 offset,
1338 n_q_points);
1339 mapping.transform(make_array_view(fe_data.shape_values,
1340 dof_index,
1341 offset,
1342 n_q_points),
1344 mapping_internal,
1345 transformed_shape_values);
1346
1347 for (unsigned int k = 0; k < n_q_points; ++k)
1348 for (unsigned int d = 0; d < dim; ++d)
1349 output_data.shape_values(first + d, k) =
1350 dof_sign * transformed_shape_values[k][d];
1351
1352 break;
1353 }
1354
1355 default:
1357 }
1358 }
1359
1360 if (fe_data.update_each & update_gradients)
1361 {
1362 switch (mapping_kind)
1363 {
1364 case mapping_none:
1365 {
1366 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1367 make_array_view(fe_data.transformed_shape_grads,
1368 offset,
1369 n_q_points);
1370 mapping.transform(make_array_view(fe_data.shape_grads,
1371 dof_index,
1372 offset,
1373 n_q_points),
1375 mapping_internal,
1376 transformed_shape_grads);
1377 for (unsigned int k = 0; k < n_q_points; ++k)
1378 for (unsigned int d = 0; d < dim; ++d)
1379 output_data.shape_gradients[first + d][k] =
1380 transformed_shape_grads[k][d];
1381 break;
1382 }
1383
1384 case mapping_covariant:
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(make_array_view(fe_data.shape_grads,
1391 dof_index,
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 transformed_shape_grads[k][d];
1409 break;
1410 }
1411
1413 {
1414 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1415 make_array_view(fe_data.transformed_shape_grads,
1416 offset,
1417 n_q_points);
1418 for (unsigned int k = 0; k < n_q_points; ++k)
1419 fe_data.untransformed_shape_grads[k + offset] =
1420 fe_data.shape_grads[dof_index][k + offset];
1421 mapping.transform(
1422 make_array_view(fe_data.untransformed_shape_grads,
1423 offset,
1424 n_q_points),
1426 mapping_internal,
1427 transformed_shape_grads);
1428
1429 for (unsigned int k = 0; k < n_q_points; ++k)
1430 for (unsigned int d = 0; d < spacedim; ++d)
1431 for (unsigned int n = 0; n < spacedim; ++n)
1432 transformed_shape_grads[k][d] +=
1433 output_data.shape_values(first + n, k) *
1434 mapping_data.jacobian_pushed_forward_grads[k][d][n];
1435
1436 for (unsigned int k = 0; k < n_q_points; ++k)
1437 for (unsigned int d = 0; d < dim; ++d)
1438 output_data.shape_gradients[first + d][k] =
1439 transformed_shape_grads[k][d];
1440
1441 break;
1442 }
1443
1445 case mapping_piola:
1446 {
1447 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1448 make_array_view(fe_data.transformed_shape_grads,
1449 offset,
1450 n_q_points);
1451 for (unsigned int k = 0; k < n_q_points; ++k)
1452 fe_data.untransformed_shape_grads[k + offset] =
1453 fe_data.shape_grads[dof_index][k + offset];
1454 mapping.transform(
1455 make_array_view(fe_data.untransformed_shape_grads,
1456 offset,
1457 n_q_points),
1459 mapping_internal,
1460 transformed_shape_grads);
1461
1462 for (unsigned int k = 0; k < n_q_points; ++k)
1463 for (unsigned int d = 0; d < spacedim; ++d)
1464 for (unsigned int n = 0; n < spacedim; ++n)
1465 transformed_shape_grads[k][d] +=
1466 (output_data.shape_values(first + n, k) *
1467 mapping_data
1469 (output_data.shape_values(first + d, k) *
1470 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
1471
1472 for (unsigned int k = 0; k < n_q_points; ++k)
1473 for (unsigned int d = 0; d < dim; ++d)
1474 output_data.shape_gradients[first + d][k] =
1475 dof_sign * transformed_shape_grads[k][d];
1476
1477 break;
1478 }
1479
1480 case mapping_nedelec:
1481 {
1482 // treat the gradients of
1483 // this particular shape
1484 // function at all
1485 // q-points. if Dv is the
1486 // gradient of the shape
1487 // function on the unit
1488 // cell, then
1489 // (J^-T)Dv(J^-1) is the
1490 // value we want to have on
1491 // the real cell.
1492 for (unsigned int k = 0; k < n_q_points; ++k)
1493 fe_data.untransformed_shape_grads[k + offset] =
1494 fe_data.shape_grads[dof_index][k + offset];
1495
1496 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1497 make_array_view(fe_data.transformed_shape_grads,
1498 offset,
1499 n_q_points);
1500 mapping.transform(
1501 make_array_view(fe_data.untransformed_shape_grads,
1502 offset,
1503 n_q_points),
1505 mapping_internal,
1506 transformed_shape_grads);
1507
1508 for (unsigned int k = 0; k < n_q_points; ++k)
1509 for (unsigned int d = 0; d < spacedim; ++d)
1510 for (unsigned int n = 0; n < spacedim; ++n)
1511 transformed_shape_grads[k][d] -=
1512 output_data.shape_values(first + n, k) *
1513 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1514
1515 for (unsigned int k = 0; k < n_q_points; ++k)
1516 for (unsigned int d = 0; d < dim; ++d)
1517 output_data.shape_gradients[first + d][k] =
1518 dof_sign * transformed_shape_grads[k][d];
1519
1520 break;
1521 }
1522
1523 default:
1525 }
1526 }
1527
1528 if (fe_data.update_each & update_hessians)
1529 {
1530 switch (mapping_kind)
1531 {
1532 case mapping_none:
1533 {
1535 transformed_shape_hessians =
1536 make_array_view(fe_data.transformed_shape_hessians,
1537 offset,
1538 n_q_points);
1539 mapping.transform(make_array_view(fe_data.shape_grad_grads,
1540 dof_index,
1541 offset,
1542 n_q_points),
1544 mapping_internal,
1545 transformed_shape_hessians);
1546
1547 for (unsigned int k = 0; k < n_q_points; ++k)
1548 for (unsigned int d = 0; d < spacedim; ++d)
1549 for (unsigned int n = 0; n < spacedim; ++n)
1550 transformed_shape_hessians[k][d] -=
1551 output_data.shape_gradients[first + d][k][n] *
1552 mapping_data.jacobian_pushed_forward_grads[k][n];
1553
1554 for (unsigned int k = 0; k < n_q_points; ++k)
1555 for (unsigned int d = 0; d < dim; ++d)
1556 output_data.shape_hessians[first + d][k] =
1557 transformed_shape_hessians[k][d];
1558
1559 break;
1560 }
1561 case mapping_covariant:
1562 {
1563 for (unsigned int k = 0; k < n_q_points; ++k)
1564 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1565 fe_data.shape_grad_grads[dof_index][k + offset];
1566
1568 transformed_shape_hessians =
1569 make_array_view(fe_data.transformed_shape_hessians,
1570 offset,
1571 n_q_points);
1572 mapping.transform(
1573 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1574 offset,
1575 n_q_points),
1577 mapping_internal,
1578 transformed_shape_hessians);
1579
1580 for (unsigned int k = 0; k < n_q_points; ++k)
1581 for (unsigned int d = 0; d < spacedim; ++d)
1582 for (unsigned int n = 0; n < spacedim; ++n)
1583 for (unsigned int i = 0; i < spacedim; ++i)
1584 for (unsigned int j = 0; j < spacedim; ++j)
1585 {
1586 transformed_shape_hessians[k][d][i][j] -=
1587 (output_data.shape_values(first + n, k) *
1588 mapping_data
1590 [k][n][d][i][j]) +
1591 (output_data.shape_gradients[first + d][k][n] *
1592 mapping_data
1593 .jacobian_pushed_forward_grads[k][n][i][j]) +
1594 (output_data.shape_gradients[first + n][k][i] *
1595 mapping_data
1596 .jacobian_pushed_forward_grads[k][n][d][j]) +
1597 (output_data.shape_gradients[first + n][k][j] *
1598 mapping_data
1599 .jacobian_pushed_forward_grads[k][n][i][d]);
1600 }
1601
1602 for (unsigned int k = 0; k < n_q_points; ++k)
1603 for (unsigned int d = 0; d < dim; ++d)
1604 output_data.shape_hessians[first + d][k] =
1605 transformed_shape_hessians[k][d];
1606
1607 break;
1608 }
1609
1611 {
1612 for (unsigned int k = 0; k < n_q_points; ++k)
1613 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1614 fe_data.shape_grad_grads[dof_index][k + offset];
1615
1617 transformed_shape_hessians =
1618 make_array_view(fe_data.transformed_shape_hessians,
1619 offset,
1620 n_q_points);
1621 mapping.transform(
1622 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1623 offset,
1624 n_q_points),
1626 mapping_internal,
1627 transformed_shape_hessians);
1628
1629 for (unsigned int k = 0; k < n_q_points; ++k)
1630 for (unsigned int d = 0; d < spacedim; ++d)
1631 for (unsigned int n = 0; n < spacedim; ++n)
1632 for (unsigned int i = 0; i < spacedim; ++i)
1633 for (unsigned int j = 0; j < spacedim; ++j)
1634 {
1635 transformed_shape_hessians[k][d][i][j] +=
1636 (output_data.shape_values(first + n, k) *
1637 mapping_data
1639 [k][d][n][i][j]) +
1640 (output_data.shape_gradients[first + n][k][i] *
1641 mapping_data
1642 .jacobian_pushed_forward_grads[k][d][n][j]) +
1643 (output_data.shape_gradients[first + n][k][j] *
1644 mapping_data
1645 .jacobian_pushed_forward_grads[k][d][i][n]) -
1646 (output_data.shape_gradients[first + d][k][n] *
1647 mapping_data
1648 .jacobian_pushed_forward_grads[k][n][i][j]);
1649 for (unsigned int m = 0; m < spacedim; ++m)
1650 transformed_shape_hessians[k][d][i][j] -=
1651 (mapping_data
1652 .jacobian_pushed_forward_grads[k][d][i]
1653 [m] *
1654 mapping_data
1655 .jacobian_pushed_forward_grads[k][m][n]
1656 [j] *
1657 output_data.shape_values(first + n, k)) +
1658 (mapping_data
1659 .jacobian_pushed_forward_grads[k][d][m]
1660 [j] *
1661 mapping_data
1662 .jacobian_pushed_forward_grads[k][m][i]
1663 [n] *
1664 output_data.shape_values(first + n, k));
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 transformed_shape_hessians[k][d];
1671
1672 break;
1673 }
1674
1676 case mapping_piola:
1677 {
1678 for (unsigned int k = 0; k < n_q_points; ++k)
1679 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1680 fe_data.shape_grad_grads[dof_index][k + offset];
1681
1683 transformed_shape_hessians =
1684 make_array_view(fe_data.transformed_shape_hessians,
1685 offset,
1686 n_q_points);
1687 mapping.transform(
1688 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1689 offset,
1690 n_q_points),
1692 mapping_internal,
1693 transformed_shape_hessians);
1694
1695 for (unsigned int k = 0; k < n_q_points; ++k)
1696 for (unsigned int d = 0; d < spacedim; ++d)
1697 for (unsigned int n = 0; n < spacedim; ++n)
1698 for (unsigned int i = 0; i < spacedim; ++i)
1699 for (unsigned int j = 0; j < spacedim; ++j)
1700 {
1701 transformed_shape_hessians[k][d][i][j] +=
1702 (output_data.shape_values(first + n, k) *
1703 mapping_data
1705 [k][d][n][i][j]) +
1706 (output_data.shape_gradients[first + n][k][i] *
1707 mapping_data
1708 .jacobian_pushed_forward_grads[k][d][n][j]) +
1709 (output_data.shape_gradients[first + n][k][j] *
1710 mapping_data
1711 .jacobian_pushed_forward_grads[k][d][i][n]) -
1712 (output_data.shape_gradients[first + d][k][n] *
1713 mapping_data
1714 .jacobian_pushed_forward_grads[k][n][i][j]);
1715
1716 transformed_shape_hessians[k][d][i][j] -=
1717 (output_data.shape_values(first + d, k) *
1718 mapping_data
1720 [k][n][n][i][j]) +
1721 (output_data.shape_gradients[first + d][k][i] *
1722 mapping_data
1723 .jacobian_pushed_forward_grads[k][n][n][j]) +
1724 (output_data.shape_gradients[first + d][k][j] *
1725 mapping_data
1726 .jacobian_pushed_forward_grads[k][n][n][i]);
1727
1728 for (unsigned int m = 0; m < spacedim; ++m)
1729 {
1730 transformed_shape_hessians[k][d][i][j] -=
1731 (mapping_data
1733 [m] *
1734 mapping_data
1736 [j] *
1737 output_data.shape_values(first + n, k)) +
1738 (mapping_data
1739 .jacobian_pushed_forward_grads[k][d][m]
1740 [j] *
1741 mapping_data
1742 .jacobian_pushed_forward_grads[k][m][i]
1743 [n] *
1744 output_data.shape_values(first + n, k));
1745
1746 transformed_shape_hessians[k][d][i][j] +=
1747 (mapping_data
1749 [m] *
1750 mapping_data
1752 [j] *
1753 output_data.shape_values(first + d, k)) +
1754 (mapping_data
1755 .jacobian_pushed_forward_grads[k][n][m]
1756 [j] *
1757 mapping_data
1758 .jacobian_pushed_forward_grads[k][m][i]
1759 [n] *
1760 output_data.shape_values(first + d, k));
1761 }
1762 }
1763
1764 for (unsigned int k = 0; k < n_q_points; ++k)
1765 for (unsigned int d = 0; d < dim; ++d)
1766 output_data.shape_hessians[first + d][k] =
1767 dof_sign * transformed_shape_hessians[k][d];
1768
1769 break;
1770 }
1771
1772 case mapping_nedelec:
1773 {
1774 for (unsigned int k = 0; k < n_q_points; ++k)
1775 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1776 fe_data.shape_grad_grads[dof_index][k + offset];
1777
1779 transformed_shape_hessians =
1780 make_array_view(fe_data.transformed_shape_hessians,
1781 offset,
1782 n_q_points);
1783 mapping.transform(
1784 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1785 offset,
1786 n_q_points),
1788 mapping_internal,
1789 transformed_shape_hessians);
1790
1791 for (unsigned int k = 0; k < n_q_points; ++k)
1792 for (unsigned int d = 0; d < spacedim; ++d)
1793 for (unsigned int n = 0; n < spacedim; ++n)
1794 for (unsigned int i = 0; i < spacedim; ++i)
1795 for (unsigned int j = 0; j < spacedim; ++j)
1796 {
1797 transformed_shape_hessians[k][d][i][j] -=
1798 (output_data.shape_values(first + n, k) *
1799 mapping_data
1801 [k][n][d][i][j]) +
1802 (output_data.shape_gradients[first + d][k][n] *
1803 mapping_data
1804 .jacobian_pushed_forward_grads[k][n][i][j]) +
1805 (output_data.shape_gradients[first + n][k][i] *
1806 mapping_data
1807 .jacobian_pushed_forward_grads[k][n][d][j]) +
1808 (output_data.shape_gradients[first + n][k][j] *
1809 mapping_data
1810 .jacobian_pushed_forward_grads[k][n][i][d]);
1811 }
1812
1813 for (unsigned int k = 0; k < n_q_points; ++k)
1814 for (unsigned int d = 0; d < dim; ++d)
1815 output_data.shape_hessians[first + d][k] =
1816 dof_sign * transformed_shape_hessians[k][d];
1817
1818 break;
1819 }
1820
1821 default:
1823 }
1824 }
1825
1826 // third derivatives are not implemented
1827 if (fe_data.update_each & update_3rd_derivatives)
1828 {
1830 }
1831 }
1832}
1833
1834
1835
1836template <int dim, int spacedim>
1837void
1840 const unsigned int face_no,
1841 const unsigned int sub_no,
1842 const Quadrature<dim - 1> &quadrature,
1843 const Mapping<dim, spacedim> &mapping,
1844 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1846 &mapping_data,
1847 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1849 spacedim>
1850 &output_data) const
1851{
1852 // convert data object to internal
1853 // data for this class. fails with
1854 // an exception if that is not
1855 // possible
1856 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1858 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1859
1860 const unsigned int n_q_points = quadrature.size();
1861
1862 // offset determines which data set
1863 // to take (all data sets for all
1864 // sub-faces are stored contiguously)
1865 const auto offset =
1867 face_no,
1868 sub_no,
1869 cell->combined_face_orientation(
1870 face_no),
1871 n_q_points,
1872 cell->subface_case(face_no));
1873
1874 // TODO: Size assertions
1875
1876 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1877 // correct thing on complicated meshes for higher order Nedelec elements.
1878 // Something similar to FE_Q should be done to permute dofs and to change the
1879 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1880 // class) is preferable.
1881 std::fill(fe_data.dof_sign_change.begin(),
1882 fe_data.dof_sign_change.end(),
1883 1.0);
1884 if (fe_data.update_each & update_values)
1885 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
1886 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
1887
1888 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1889 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1890 // Also nothing in 3d since we take care of it by using the
1891 // adjust_quad_dof_sign_for_face_orientation_table.
1892 if (fe_data.update_each & update_values)
1893 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1894 *this,
1895 this->mapping_kind,
1896 fe_data.dof_sign_change);
1897
1898 // What is the first dof_index on a quad?
1899 const unsigned int first_quad_index = this->get_first_quad_index();
1900 // How many dofs per quad and how many quad dofs do we have at all?
1901 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1902 const unsigned int n_quad_dofs =
1903 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1904
1905 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1906 ++dof_index)
1907 {
1908 /*
1909 * This assumes that the dofs are ordered by first vertices, lines, quads
1910 * and volume dofs. Note that in 2d this always gives false.
1911 */
1912 const bool is_quad_dof =
1913 (dim == 2 ? false :
1914 (first_quad_index <= dof_index) &&
1915 (dof_index < first_quad_index + n_quad_dofs));
1916
1917 // TODO: This hack is not pretty and it is only here to handle the 2d
1918 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1919 // handled by the
1920 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1921 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1922 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1923 // it is handles with a table. This array is allocated in
1924 // fe_poly_tensor.h.
1925 double dof_sign = 1.0;
1926 // under some circumstances fe_data.dof_sign_change is not allocated
1927 if (fe_data.update_each & update_values)
1928 dof_sign = fe_data.dof_sign_change[dof_index];
1929
1930 if (is_quad_dof)
1931 {
1932 /*
1933 * Find the face belonging to this dof_index. This is integer
1934 * division.
1935 */
1936 unsigned int face_index_from_dof_index =
1937 (dof_index - first_quad_index) / (n_dofs_per_quad);
1938
1939 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1940
1941 // Correct the dof_sign if necessary
1942 if (adjust_quad_dof_sign_for_face_orientation(
1943 local_quad_dof_index,
1944 face_index_from_dof_index,
1945 cell->face_orientation(face_index_from_dof_index),
1946 cell->face_flip(face_index_from_dof_index),
1947 cell->face_rotation(face_index_from_dof_index)))
1948 dof_sign = -1.0;
1949 }
1950
1951 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1952
1953 const unsigned int first =
1954 output_data.shape_function_to_row_table
1955 [dof_index * this->n_components() +
1956 this->get_nonzero_components(dof_index).first_selected_component()];
1957
1958 if (fe_data.update_each & update_values)
1959 {
1960 switch (mapping_kind)
1961 {
1962 case mapping_none:
1963 {
1964 for (unsigned int k = 0; k < n_q_points; ++k)
1965 for (unsigned int d = 0; d < dim; ++d)
1966 output_data.shape_values(first + d, k) =
1967 fe_data.shape_values[dof_index][k + offset][d];
1968 break;
1969 }
1970
1971 case mapping_covariant:
1973 {
1975 transformed_shape_values =
1976 make_array_view(fe_data.transformed_shape_values,
1977 offset,
1978 n_q_points);
1979 mapping.transform(make_array_view(fe_data.shape_values,
1980 dof_index,
1981 offset,
1982 n_q_points),
1983 mapping_kind,
1984 mapping_internal,
1985 transformed_shape_values);
1986
1987 for (unsigned int k = 0; k < n_q_points; ++k)
1988 for (unsigned int d = 0; d < dim; ++d)
1989 output_data.shape_values(first + d, k) =
1990 transformed_shape_values[k][d];
1991
1992 break;
1993 }
1994
1996 case mapping_piola:
1997 {
1999 transformed_shape_values =
2000 make_array_view(fe_data.transformed_shape_values,
2001 offset,
2002 n_q_points);
2003
2004 mapping.transform(make_array_view(fe_data.shape_values,
2005 dof_index,
2006 offset,
2007 n_q_points),
2009 mapping_internal,
2010 transformed_shape_values);
2011 for (unsigned int k = 0; k < n_q_points; ++k)
2012 for (unsigned int d = 0; d < dim; ++d)
2013 output_data.shape_values(first + d, k) =
2014 dof_sign * transformed_shape_values[k][d];
2015 break;
2016 }
2017
2018 case mapping_nedelec:
2019 {
2021 transformed_shape_values =
2022 make_array_view(fe_data.transformed_shape_values,
2023 offset,
2024 n_q_points);
2025
2026 mapping.transform(make_array_view(fe_data.shape_values,
2027 dof_index,
2028 offset,
2029 n_q_points),
2031 mapping_internal,
2032 transformed_shape_values);
2033
2034 for (unsigned int k = 0; k < n_q_points; ++k)
2035 for (unsigned int d = 0; d < dim; ++d)
2036 output_data.shape_values(first + d, k) =
2037 dof_sign * transformed_shape_values[k][d];
2038
2039 break;
2040 }
2041
2042 default:
2044 }
2045 }
2046
2047 if (fe_data.update_each & update_gradients)
2048 {
2049 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
2050 make_array_view(fe_data.transformed_shape_grads,
2051 offset,
2052 n_q_points);
2053 switch (mapping_kind)
2054 {
2055 case mapping_none:
2056 {
2057 mapping.transform(make_array_view(fe_data.shape_grads,
2058 dof_index,
2059 offset,
2060 n_q_points),
2062 mapping_internal,
2063 transformed_shape_grads);
2064 for (unsigned int k = 0; k < n_q_points; ++k)
2065 for (unsigned int d = 0; d < dim; ++d)
2066 output_data.shape_gradients[first + d][k] =
2067 transformed_shape_grads[k][d];
2068 break;
2069 }
2070
2071 case mapping_covariant:
2072 {
2073 mapping.transform(make_array_view(fe_data.shape_grads,
2074 dof_index,
2075 offset,
2076 n_q_points),
2078 mapping_internal,
2079 transformed_shape_grads);
2080
2081 for (unsigned int k = 0; k < n_q_points; ++k)
2082 for (unsigned int d = 0; d < spacedim; ++d)
2083 for (unsigned int n = 0; n < spacedim; ++n)
2084 transformed_shape_grads[k][d] -=
2085 output_data.shape_values(first + n, k) *
2086 mapping_data.jacobian_pushed_forward_grads[k][n][d];
2087
2088 for (unsigned int k = 0; k < n_q_points; ++k)
2089 for (unsigned int d = 0; d < dim; ++d)
2090 output_data.shape_gradients[first + d][k] =
2091 transformed_shape_grads[k][d];
2092
2093 break;
2094 }
2095
2097 {
2098 for (unsigned int k = 0; k < n_q_points; ++k)
2099 fe_data.untransformed_shape_grads[k + offset] =
2100 fe_data.shape_grads[dof_index][k + offset];
2101
2102 mapping.transform(
2103 make_array_view(fe_data.untransformed_shape_grads,
2104 offset,
2105 n_q_points),
2107 mapping_internal,
2108 transformed_shape_grads);
2109
2110 for (unsigned int k = 0; k < n_q_points; ++k)
2111 for (unsigned int d = 0; d < spacedim; ++d)
2112 for (unsigned int n = 0; n < spacedim; ++n)
2113 transformed_shape_grads[k][d] +=
2114 output_data.shape_values(first + n, k) *
2115 mapping_data.jacobian_pushed_forward_grads[k][d][n];
2116
2117 for (unsigned int k = 0; k < n_q_points; ++k)
2118 for (unsigned int d = 0; d < dim; ++d)
2119 output_data.shape_gradients[first + d][k] =
2120 transformed_shape_grads[k][d];
2121
2122 break;
2123 }
2124
2126 case mapping_piola:
2127 {
2128 for (unsigned int k = 0; k < n_q_points; ++k)
2129 fe_data.untransformed_shape_grads[k + offset] =
2130 fe_data.shape_grads[dof_index][k + offset];
2131
2132 mapping.transform(
2133 make_array_view(fe_data.untransformed_shape_grads,
2134 offset,
2135 n_q_points),
2137 mapping_internal,
2138 transformed_shape_grads);
2139
2140 for (unsigned int k = 0; k < n_q_points; ++k)
2141 for (unsigned int d = 0; d < spacedim; ++d)
2142 for (unsigned int n = 0; n < spacedim; ++n)
2143 transformed_shape_grads[k][d] +=
2144 (output_data.shape_values(first + n, k) *
2145 mapping_data
2147 (output_data.shape_values(first + d, k) *
2148 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
2149
2150 for (unsigned int k = 0; k < n_q_points; ++k)
2151 for (unsigned int d = 0; d < dim; ++d)
2152 output_data.shape_gradients[first + d][k] =
2153 dof_sign * transformed_shape_grads[k][d];
2154
2155 break;
2156 }
2157
2158 case mapping_nedelec:
2159 {
2160 // this particular shape
2161 // function at all
2162 // q-points. if Dv is the
2163 // gradient of the shape
2164 // function on the unit
2165 // cell, then
2166 // (J^-T)Dv(J^-1) is the
2167 // value we want to have on
2168 // the real cell.
2169 for (unsigned int k = 0; k < n_q_points; ++k)
2170 fe_data.untransformed_shape_grads[k + offset] =
2171 fe_data.shape_grads[dof_index][k + offset];
2172
2173 mapping.transform(
2174 make_array_view(fe_data.untransformed_shape_grads,
2175 offset,
2176 n_q_points),
2178 mapping_internal,
2179 transformed_shape_grads);
2180
2181 for (unsigned int k = 0; k < n_q_points; ++k)
2182 for (unsigned int d = 0; d < spacedim; ++d)
2183 for (unsigned int n = 0; n < spacedim; ++n)
2184 transformed_shape_grads[k][d] -=
2185 output_data.shape_values(first + n, k) *
2186 mapping_data.jacobian_pushed_forward_grads[k][n][d];
2187
2188 for (unsigned int k = 0; k < n_q_points; ++k)
2189 for (unsigned int d = 0; d < dim; ++d)
2190 output_data.shape_gradients[first + d][k] =
2191 dof_sign * transformed_shape_grads[k][d];
2192
2193 break;
2194 }
2195
2196 default:
2198 }
2199 }
2200
2201 if (fe_data.update_each & update_hessians)
2202 {
2203 switch (mapping_kind)
2204 {
2205 case mapping_none:
2206 {
2208 transformed_shape_hessians =
2209 make_array_view(fe_data.transformed_shape_hessians,
2210 offset,
2211 n_q_points);
2212 mapping.transform(make_array_view(fe_data.shape_grad_grads,
2213 dof_index,
2214 offset,
2215 n_q_points),
2217 mapping_internal,
2218 transformed_shape_hessians);
2219
2220 for (unsigned int k = 0; k < n_q_points; ++k)
2221 for (unsigned int d = 0; d < spacedim; ++d)
2222 for (unsigned int n = 0; n < spacedim; ++n)
2223 transformed_shape_hessians[k][d] -=
2224 output_data.shape_gradients[first + d][k][n] *
2225 mapping_data.jacobian_pushed_forward_grads[k][n];
2226
2227 for (unsigned int k = 0; k < n_q_points; ++k)
2228 for (unsigned int d = 0; d < dim; ++d)
2229 output_data.shape_hessians[first + d][k] =
2230 transformed_shape_hessians[k][d];
2231
2232 break;
2233 }
2234 case mapping_covariant:
2235 {
2236 for (unsigned int k = 0; k < n_q_points; ++k)
2237 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2238 fe_data.shape_grad_grads[dof_index][k + offset];
2239
2241 transformed_shape_hessians =
2242 make_array_view(fe_data.transformed_shape_hessians,
2243 offset,
2244 n_q_points);
2245 mapping.transform(
2246 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2247 offset,
2248 n_q_points),
2250 mapping_internal,
2251 transformed_shape_hessians);
2252
2253 for (unsigned int k = 0; k < n_q_points; ++k)
2254 for (unsigned int d = 0; d < spacedim; ++d)
2255 for (unsigned int n = 0; n < spacedim; ++n)
2256 for (unsigned int i = 0; i < spacedim; ++i)
2257 for (unsigned int j = 0; j < spacedim; ++j)
2258 {
2259 transformed_shape_hessians[k][d][i][j] -=
2260 (output_data.shape_values(first + n, k) *
2261 mapping_data
2263 [k][n][d][i][j]) +
2264 (output_data.shape_gradients[first + d][k][n] *
2265 mapping_data
2266 .jacobian_pushed_forward_grads[k][n][i][j]) +
2267 (output_data.shape_gradients[first + n][k][i] *
2268 mapping_data
2269 .jacobian_pushed_forward_grads[k][n][d][j]) +
2270 (output_data.shape_gradients[first + n][k][j] *
2271 mapping_data
2272 .jacobian_pushed_forward_grads[k][n][i][d]);
2273 }
2274
2275 for (unsigned int k = 0; k < n_q_points; ++k)
2276 for (unsigned int d = 0; d < dim; ++d)
2277 output_data.shape_hessians[first + d][k] =
2278 transformed_shape_hessians[k][d];
2279
2280 break;
2281 }
2282
2284 {
2285 for (unsigned int k = 0; k < n_q_points; ++k)
2286 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2287 fe_data.shape_grad_grads[dof_index][k + offset];
2288
2290 transformed_shape_hessians =
2291 make_array_view(fe_data.transformed_shape_hessians,
2292 offset,
2293 n_q_points);
2294 mapping.transform(
2295 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2296 offset,
2297 n_q_points),
2299 mapping_internal,
2300 transformed_shape_hessians);
2301
2302 for (unsigned int k = 0; k < n_q_points; ++k)
2303 for (unsigned int d = 0; d < spacedim; ++d)
2304 for (unsigned int n = 0; n < spacedim; ++n)
2305 for (unsigned int i = 0; i < spacedim; ++i)
2306 for (unsigned int j = 0; j < spacedim; ++j)
2307 {
2308 transformed_shape_hessians[k][d][i][j] +=
2309 (output_data.shape_values(first + n, k) *
2310 mapping_data
2312 [k][d][n][i][j]) +
2313 (output_data.shape_gradients[first + n][k][i] *
2314 mapping_data
2315 .jacobian_pushed_forward_grads[k][d][n][j]) +
2316 (output_data.shape_gradients[first + n][k][j] *
2317 mapping_data
2318 .jacobian_pushed_forward_grads[k][d][i][n]) -
2319 (output_data.shape_gradients[first + d][k][n] *
2320 mapping_data
2321 .jacobian_pushed_forward_grads[k][n][i][j]);
2322 for (unsigned int m = 0; m < spacedim; ++m)
2323 transformed_shape_hessians[k][d][i][j] -=
2324 (mapping_data
2325 .jacobian_pushed_forward_grads[k][d][i]
2326 [m] *
2327 mapping_data
2328 .jacobian_pushed_forward_grads[k][m][n]
2329 [j] *
2330 output_data.shape_values(first + n, k)) +
2331 (mapping_data
2332 .jacobian_pushed_forward_grads[k][d][m]
2333 [j] *
2334 mapping_data
2335 .jacobian_pushed_forward_grads[k][m][i]
2336 [n] *
2337 output_data.shape_values(first + n, k));
2338 }
2339
2340 for (unsigned int k = 0; k < n_q_points; ++k)
2341 for (unsigned int d = 0; d < dim; ++d)
2342 output_data.shape_hessians[first + d][k] =
2343 transformed_shape_hessians[k][d];
2344
2345 break;
2346 }
2347
2349 case mapping_piola:
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
2378 [k][d][n][i][j]) +
2379 (output_data.shape_gradients[first + n][k][i] *
2380 mapping_data
2381 .jacobian_pushed_forward_grads[k][d][n][j]) +
2382 (output_data.shape_gradients[first + n][k][j] *
2383 mapping_data
2384 .jacobian_pushed_forward_grads[k][d][i][n]) -
2385 (output_data.shape_gradients[first + d][k][n] *
2386 mapping_data
2387 .jacobian_pushed_forward_grads[k][n][i][j]);
2388
2389 transformed_shape_hessians[k][d][i][j] -=
2390 (output_data.shape_values(first + d, k) *
2391 mapping_data
2393 [k][n][n][i][j]) +
2394 (output_data.shape_gradients[first + d][k][i] *
2395 mapping_data
2396 .jacobian_pushed_forward_grads[k][n][n][j]) +
2397 (output_data.shape_gradients[first + d][k][j] *
2398 mapping_data
2399 .jacobian_pushed_forward_grads[k][n][n][i]);
2400 for (unsigned int m = 0; m < spacedim; ++m)
2401 {
2402 transformed_shape_hessians[k][d][i][j] -=
2403 (mapping_data
2405 [m] *
2406 mapping_data
2408 [j] *
2409 output_data.shape_values(first + n, k)) +
2410 (mapping_data
2411 .jacobian_pushed_forward_grads[k][d][m]
2412 [j] *
2413 mapping_data
2414 .jacobian_pushed_forward_grads[k][m][i]
2415 [n] *
2416 output_data.shape_values(first + n, k));
2417
2418 transformed_shape_hessians[k][d][i][j] +=
2419 (mapping_data
2421 [m] *
2422 mapping_data
2424 [j] *
2425 output_data.shape_values(first + d, k)) +
2426 (mapping_data
2427 .jacobian_pushed_forward_grads[k][n][m]
2428 [j] *
2429 mapping_data
2430 .jacobian_pushed_forward_grads[k][m][i]
2431 [n] *
2432 output_data.shape_values(first + d, k));
2433 }
2434 }
2435
2436 for (unsigned int k = 0; k < n_q_points; ++k)
2437 for (unsigned int d = 0; d < dim; ++d)
2438 output_data.shape_hessians[first + d][k] =
2439 dof_sign * transformed_shape_hessians[k][d];
2440
2441 break;
2442 }
2443
2444 case mapping_nedelec:
2445 {
2446 for (unsigned int k = 0; k < n_q_points; ++k)
2447 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2448 fe_data.shape_grad_grads[dof_index][k + offset];
2449
2451 transformed_shape_hessians =
2452 make_array_view(fe_data.transformed_shape_hessians,
2453 offset,
2454 n_q_points);
2455 mapping.transform(
2456 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2457 offset,
2458 n_q_points),
2460 mapping_internal,
2461 transformed_shape_hessians);
2462
2463 for (unsigned int k = 0; k < n_q_points; ++k)
2464 for (unsigned int d = 0; d < spacedim; ++d)
2465 for (unsigned int n = 0; n < spacedim; ++n)
2466 for (unsigned int i = 0; i < spacedim; ++i)
2467 for (unsigned int j = 0; j < spacedim; ++j)
2468 {
2469 transformed_shape_hessians[k][d][i][j] -=
2470 (output_data.shape_values(first + n, k) *
2471 mapping_data
2473 [k][n][d][i][j]) +
2474 (output_data.shape_gradients[first + d][k][n] *
2475 mapping_data
2476 .jacobian_pushed_forward_grads[k][n][i][j]) +
2477 (output_data.shape_gradients[first + n][k][i] *
2478 mapping_data
2479 .jacobian_pushed_forward_grads[k][n][d][j]) +
2480 (output_data.shape_gradients[first + n][k][j] *
2481 mapping_data
2482 .jacobian_pushed_forward_grads[k][n][i][d]);
2483 }
2484
2485 for (unsigned int k = 0; k < n_q_points; ++k)
2486 for (unsigned int d = 0; d < dim; ++d)
2487 output_data.shape_hessians[first + d][k] =
2488 dof_sign * transformed_shape_hessians[k][d];
2489
2490 break;
2491 }
2492
2493 default:
2495 }
2496 }
2497
2498 // third derivatives are not implemented
2499 if (fe_data.update_each & update_3rd_derivatives)
2500 {
2502 }
2503 }
2504}
2505
2506
2507
2508template <int dim, int spacedim>
2511 const UpdateFlags flags) const
2512{
2514
2515 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2516 {
2517 const MappingKind mapping_kind = get_mapping_kind(i);
2518
2519 switch (mapping_kind)
2520 {
2521 case mapping_none:
2522 {
2523 if (flags & update_values)
2524 out |= update_values;
2525
2526 if (flags & update_gradients)
2529
2530 if (flags & update_hessians)
2534 break;
2535 }
2537 case mapping_piola:
2538 {
2539 if (flags & update_values)
2540 out |= update_values | update_piola;
2541
2542 if (flags & update_gradients)
2547
2548 if (flags & update_hessians)
2553
2554 break;
2555 }
2556
2557
2559 {
2560 if (flags & update_values)
2561 out |= update_values | update_piola;
2562
2563 if (flags & update_gradients)
2568
2569 if (flags & update_hessians)
2574
2575 break;
2576 }
2577
2578 case mapping_nedelec:
2579 case mapping_covariant:
2580 {
2581 if (flags & update_values)
2583
2584 if (flags & update_gradients)
2588
2589 if (flags & update_hessians)
2594
2595 break;
2596 }
2597
2598 default:
2599 {
2601 }
2602 }
2603 }
2604
2605 return out;
2606}
2607
2608#endif
2609// explicit instantiations
2610#include "fe/fe_poly_tensor.inst"
2611
2612
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
virtual void fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
bool adjust_quad_dof_sign_for_face_orientation(const unsigned int index, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation) const
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
MappingKind get_mapping_kind(const unsigned int i) const
FE_PolyTensor(const TensorPolynomialsBase< dim > &polynomials, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
bool single_mapping_kind() const
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int tensor_degree() const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation) const
Abstract base class for mapping classes.
Definition mapping.h:320
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:113
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
unsigned int size() const
unsigned int size() const
Definition collection.h:316
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > first
Definition grid_out.cc:4632
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_covariant_transformation
Covariant transformation.
@ update_gradients
Shape function gradients.
@ update_default
No update.
@ update_piola
Values needed for Piola transform.
MappingKind
Definition mapping.h:81
@ mapping_piola
Definition mapping.h:116
@ mapping_covariant_gradient
Definition mapping.h:102
@ mapping_covariant
Definition mapping.h:91
@ mapping_nedelec
Definition mapping.h:131
@ mapping_contravariant
Definition mapping.h:96
@ mapping_raviart_thomas
Definition mapping.h:136
@ mapping_contravariant_hessian
Definition mapping.h:158
@ mapping_covariant_hessian
Definition mapping.h:152
@ mapping_contravariant_gradient
Definition mapping.h:108
@ mapping_none
Definition mapping.h:86
@ mapping_piola_gradient
Definition mapping.h:122
@ mapping_piola_hessian
Definition mapping.h:164
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)
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346