Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08: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 // TODO: The same 'legacy' comments for 2d apply here as well: these classes
1184 // do not handle non-standard orientations in 2d in a way consistent with the
1185 // rest of the library, but are consistent with themselves (see, e.g., the
1186 // fe_conformity_dim_2 tests).
1187 //
1188 // In this case: all of this code was written assuming that QProjector assumed
1189 // that all faces were in the default orientation in 2d, but contains special
1190 // workarounds in case that isn't the case. Hence, to keep those workarounds
1191 // working, we still assume that all faces are in the default orientation.
1192 const auto offset = QProjector<dim>::DataSetDescriptor::face(
1193 this->reference_cell(),
1194 face_no,
1196 cell->combined_face_orientation(face_no),
1197 n_q_points);
1198
1199 // TODO: Size assertions
1200
1201 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1202 // correct thing on complicated meshes for higher order Nedelec elements.
1203 // Something similar to FE_Q should be done to permute dofs and to change the
1204 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1205 // class) is preferable.
1206 std::fill(fe_data.dof_sign_change.begin(),
1207 fe_data.dof_sign_change.end(),
1208 1.0);
1209 if (fe_data.update_each & update_values)
1210 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
1211 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
1212
1213 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1214 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1215 // Also nothing in 3d since we take care of it by using the
1216 // adjust_quad_dof_sign_for_face_orientation_table.
1217 if (fe_data.update_each & update_values)
1218 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1219 *this,
1220 this->mapping_kind,
1221 fe_data.dof_sign_change);
1222
1223 // What is the first dof_index on a quad?
1224 const unsigned int first_quad_index = this->get_first_quad_index();
1225 // How many dofs per quad and how many quad dofs do we have at all?
1226 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1227 const unsigned int n_quad_dofs =
1228 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1229
1230 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1231 ++dof_index)
1232 {
1233 /*
1234 * This assumes that the dofs are ordered by first vertices, lines, quads
1235 * and volume dofs. Note that in 2d this always gives false.
1236 */
1237 const bool is_quad_dof =
1238 (dim == 2 ? false :
1239 (first_quad_index <= dof_index) &&
1240 (dof_index < first_quad_index + n_quad_dofs));
1241
1242 // TODO: This hack is not pretty and it is only here to handle the 2d
1243 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1244 // handled by the
1245 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1246 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1247 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1248 // it is handles with a table. This array is allocated in
1249 // fe_poly_tensor.h.
1250 double dof_sign = 1.0;
1251 // under some circumstances fe_data.dof_sign_change is not allocated
1252 if (fe_data.update_each & update_values)
1253 dof_sign = fe_data.dof_sign_change[dof_index];
1254
1255 if (is_quad_dof)
1256 {
1257 /*
1258 * Find the face belonging to this dof_index. This is integer
1259 * division.
1260 */
1261 unsigned int face_index_from_dof_index =
1262 (dof_index - first_quad_index) / (n_dofs_per_quad);
1263
1264 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1265
1266 // Correct the dof_sign if necessary
1267 if (adjust_quad_dof_sign_for_face_orientation(
1268 local_quad_dof_index,
1269 face_index_from_dof_index,
1270 cell->face_orientation(face_index_from_dof_index),
1271 cell->face_flip(face_index_from_dof_index),
1272 cell->face_rotation(face_index_from_dof_index)))
1273 dof_sign = -1.0;
1274 }
1275
1276 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1277
1278 const unsigned int first =
1279 output_data.shape_function_to_row_table
1280 [dof_index * this->n_components() +
1281 this->get_nonzero_components(dof_index).first_selected_component()];
1282
1283 if (fe_data.update_each & update_values)
1284 {
1285 switch (mapping_kind)
1286 {
1287 case mapping_none:
1288 {
1289 for (unsigned int k = 0; k < n_q_points; ++k)
1290 for (unsigned int d = 0; d < dim; ++d)
1291 output_data.shape_values(first + d, k) =
1292 fe_data.shape_values[dof_index][k + offset][d];
1293 break;
1294 }
1295
1296 case mapping_covariant:
1298 {
1300 transformed_shape_values =
1301 make_array_view(fe_data.transformed_shape_values,
1302 offset,
1303 n_q_points);
1304 mapping.transform(make_array_view(fe_data.shape_values,
1305 dof_index,
1306 offset,
1307 n_q_points),
1308 mapping_kind,
1309 mapping_internal,
1310 transformed_shape_values);
1311
1312 for (unsigned int k = 0; k < n_q_points; ++k)
1313 for (unsigned int d = 0; d < dim; ++d)
1314 output_data.shape_values(first + d, k) =
1315 transformed_shape_values[k][d];
1316
1317 break;
1318 }
1320 case mapping_piola:
1321 {
1323 transformed_shape_values =
1324 make_array_view(fe_data.transformed_shape_values,
1325 offset,
1326 n_q_points);
1327 mapping.transform(make_array_view(fe_data.shape_values,
1328 dof_index,
1329 offset,
1330 n_q_points),
1332 mapping_internal,
1333 transformed_shape_values);
1334 for (unsigned int k = 0; k < n_q_points; ++k)
1335 for (unsigned int d = 0; d < dim; ++d)
1336 output_data.shape_values(first + d, k) =
1337 dof_sign * transformed_shape_values[k][d];
1338 break;
1339 }
1340
1341 case mapping_nedelec:
1342 {
1344 transformed_shape_values =
1345 make_array_view(fe_data.transformed_shape_values,
1346 offset,
1347 n_q_points);
1348 mapping.transform(make_array_view(fe_data.shape_values,
1349 dof_index,
1350 offset,
1351 n_q_points),
1353 mapping_internal,
1354 transformed_shape_values);
1355
1356 for (unsigned int k = 0; k < n_q_points; ++k)
1357 for (unsigned int d = 0; d < dim; ++d)
1358 output_data.shape_values(first + d, k) =
1359 dof_sign * transformed_shape_values[k][d];
1360
1361 break;
1362 }
1363
1364 default:
1366 }
1367 }
1368
1369 if (fe_data.update_each & update_gradients)
1370 {
1371 switch (mapping_kind)
1372 {
1373 case mapping_none:
1374 {
1375 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1376 make_array_view(fe_data.transformed_shape_grads,
1377 offset,
1378 n_q_points);
1379 mapping.transform(make_array_view(fe_data.shape_grads,
1380 dof_index,
1381 offset,
1382 n_q_points),
1384 mapping_internal,
1385 transformed_shape_grads);
1386 for (unsigned int k = 0; k < n_q_points; ++k)
1387 for (unsigned int d = 0; d < dim; ++d)
1388 output_data.shape_gradients[first + d][k] =
1389 transformed_shape_grads[k][d];
1390 break;
1391 }
1392
1393 case mapping_covariant:
1394 {
1395 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1396 make_array_view(fe_data.transformed_shape_grads,
1397 offset,
1398 n_q_points);
1399 mapping.transform(make_array_view(fe_data.shape_grads,
1400 dof_index,
1401 offset,
1402 n_q_points),
1404 mapping_internal,
1405 transformed_shape_grads);
1406
1407 for (unsigned int k = 0; k < n_q_points; ++k)
1408 for (unsigned int d = 0; d < spacedim; ++d)
1409 for (unsigned int n = 0; n < spacedim; ++n)
1410 transformed_shape_grads[k][d] -=
1411 output_data.shape_values(first + n, k) *
1412 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1413
1414 for (unsigned int k = 0; k < n_q_points; ++k)
1415 for (unsigned int d = 0; d < dim; ++d)
1416 output_data.shape_gradients[first + d][k] =
1417 transformed_shape_grads[k][d];
1418 break;
1419 }
1420
1422 {
1423 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1424 make_array_view(fe_data.transformed_shape_grads,
1425 offset,
1426 n_q_points);
1427 for (unsigned int k = 0; k < n_q_points; ++k)
1428 fe_data.untransformed_shape_grads[k + offset] =
1429 fe_data.shape_grads[dof_index][k + offset];
1430 mapping.transform(
1431 make_array_view(fe_data.untransformed_shape_grads,
1432 offset,
1433 n_q_points),
1435 mapping_internal,
1436 transformed_shape_grads);
1437
1438 for (unsigned int k = 0; k < n_q_points; ++k)
1439 for (unsigned int d = 0; d < spacedim; ++d)
1440 for (unsigned int n = 0; n < spacedim; ++n)
1441 transformed_shape_grads[k][d] +=
1442 output_data.shape_values(first + n, k) *
1443 mapping_data.jacobian_pushed_forward_grads[k][d][n];
1444
1445 for (unsigned int k = 0; k < n_q_points; ++k)
1446 for (unsigned int d = 0; d < dim; ++d)
1447 output_data.shape_gradients[first + d][k] =
1448 transformed_shape_grads[k][d];
1449
1450 break;
1451 }
1452
1454 case mapping_piola:
1455 {
1456 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1457 make_array_view(fe_data.transformed_shape_grads,
1458 offset,
1459 n_q_points);
1460 for (unsigned int k = 0; k < n_q_points; ++k)
1461 fe_data.untransformed_shape_grads[k + offset] =
1462 fe_data.shape_grads[dof_index][k + offset];
1463 mapping.transform(
1464 make_array_view(fe_data.untransformed_shape_grads,
1465 offset,
1466 n_q_points),
1468 mapping_internal,
1469 transformed_shape_grads);
1470
1471 for (unsigned int k = 0; k < n_q_points; ++k)
1472 for (unsigned int d = 0; d < spacedim; ++d)
1473 for (unsigned int n = 0; n < spacedim; ++n)
1474 transformed_shape_grads[k][d] +=
1475 (output_data.shape_values(first + n, k) *
1476 mapping_data
1478 (output_data.shape_values(first + d, k) *
1479 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
1480
1481 for (unsigned int k = 0; k < n_q_points; ++k)
1482 for (unsigned int d = 0; d < dim; ++d)
1483 output_data.shape_gradients[first + d][k] =
1484 dof_sign * transformed_shape_grads[k][d];
1485
1486 break;
1487 }
1488
1489 case mapping_nedelec:
1490 {
1491 // treat the gradients of
1492 // this particular shape
1493 // function at all
1494 // q-points. if Dv is the
1495 // gradient of the shape
1496 // function on the unit
1497 // cell, then
1498 // (J^-T)Dv(J^-1) is the
1499 // value we want to have on
1500 // the real cell.
1501 for (unsigned int k = 0; k < n_q_points; ++k)
1502 fe_data.untransformed_shape_grads[k + offset] =
1503 fe_data.shape_grads[dof_index][k + offset];
1504
1505 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1506 make_array_view(fe_data.transformed_shape_grads,
1507 offset,
1508 n_q_points);
1509 mapping.transform(
1510 make_array_view(fe_data.untransformed_shape_grads,
1511 offset,
1512 n_q_points),
1514 mapping_internal,
1515 transformed_shape_grads);
1516
1517 for (unsigned int k = 0; k < n_q_points; ++k)
1518 for (unsigned int d = 0; d < spacedim; ++d)
1519 for (unsigned int n = 0; n < spacedim; ++n)
1520 transformed_shape_grads[k][d] -=
1521 output_data.shape_values(first + n, k) *
1522 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1523
1524 for (unsigned int k = 0; k < n_q_points; ++k)
1525 for (unsigned int d = 0; d < dim; ++d)
1526 output_data.shape_gradients[first + d][k] =
1527 dof_sign * transformed_shape_grads[k][d];
1528
1529 break;
1530 }
1531
1532 default:
1534 }
1535 }
1536
1537 if (fe_data.update_each & update_hessians)
1538 {
1539 switch (mapping_kind)
1540 {
1541 case mapping_none:
1542 {
1544 transformed_shape_hessians =
1545 make_array_view(fe_data.transformed_shape_hessians,
1546 offset,
1547 n_q_points);
1548 mapping.transform(make_array_view(fe_data.shape_grad_grads,
1549 dof_index,
1550 offset,
1551 n_q_points),
1553 mapping_internal,
1554 transformed_shape_hessians);
1555
1556 for (unsigned int k = 0; k < n_q_points; ++k)
1557 for (unsigned int d = 0; d < spacedim; ++d)
1558 for (unsigned int n = 0; n < spacedim; ++n)
1559 transformed_shape_hessians[k][d] -=
1560 output_data.shape_gradients[first + d][k][n] *
1561 mapping_data.jacobian_pushed_forward_grads[k][n];
1562
1563 for (unsigned int k = 0; k < n_q_points; ++k)
1564 for (unsigned int d = 0; d < dim; ++d)
1565 output_data.shape_hessians[first + d][k] =
1566 transformed_shape_hessians[k][d];
1567
1568 break;
1569 }
1570 case mapping_covariant:
1571 {
1572 for (unsigned int k = 0; k < n_q_points; ++k)
1573 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1574 fe_data.shape_grad_grads[dof_index][k + offset];
1575
1577 transformed_shape_hessians =
1578 make_array_view(fe_data.transformed_shape_hessians,
1579 offset,
1580 n_q_points);
1581 mapping.transform(
1582 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1583 offset,
1584 n_q_points),
1586 mapping_internal,
1587 transformed_shape_hessians);
1588
1589 for (unsigned int k = 0; k < n_q_points; ++k)
1590 for (unsigned int d = 0; d < spacedim; ++d)
1591 for (unsigned int n = 0; n < spacedim; ++n)
1592 for (unsigned int i = 0; i < spacedim; ++i)
1593 for (unsigned int j = 0; j < spacedim; ++j)
1594 {
1595 transformed_shape_hessians[k][d][i][j] -=
1596 (output_data.shape_values(first + n, k) *
1597 mapping_data
1599 [k][n][d][i][j]) +
1600 (output_data.shape_gradients[first + d][k][n] *
1601 mapping_data
1602 .jacobian_pushed_forward_grads[k][n][i][j]) +
1603 (output_data.shape_gradients[first + n][k][i] *
1604 mapping_data
1605 .jacobian_pushed_forward_grads[k][n][d][j]) +
1606 (output_data.shape_gradients[first + n][k][j] *
1607 mapping_data
1608 .jacobian_pushed_forward_grads[k][n][i][d]);
1609 }
1610
1611 for (unsigned int k = 0; k < n_q_points; ++k)
1612 for (unsigned int d = 0; d < dim; ++d)
1613 output_data.shape_hessians[first + d][k] =
1614 transformed_shape_hessians[k][d];
1615
1616 break;
1617 }
1618
1620 {
1621 for (unsigned int k = 0; k < n_q_points; ++k)
1622 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1623 fe_data.shape_grad_grads[dof_index][k + offset];
1624
1626 transformed_shape_hessians =
1627 make_array_view(fe_data.transformed_shape_hessians,
1628 offset,
1629 n_q_points);
1630 mapping.transform(
1631 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1632 offset,
1633 n_q_points),
1635 mapping_internal,
1636 transformed_shape_hessians);
1637
1638 for (unsigned int k = 0; k < n_q_points; ++k)
1639 for (unsigned int d = 0; d < spacedim; ++d)
1640 for (unsigned int n = 0; n < spacedim; ++n)
1641 for (unsigned int i = 0; i < spacedim; ++i)
1642 for (unsigned int j = 0; j < spacedim; ++j)
1643 {
1644 transformed_shape_hessians[k][d][i][j] +=
1645 (output_data.shape_values(first + n, k) *
1646 mapping_data
1648 [k][d][n][i][j]) +
1649 (output_data.shape_gradients[first + n][k][i] *
1650 mapping_data
1651 .jacobian_pushed_forward_grads[k][d][n][j]) +
1652 (output_data.shape_gradients[first + n][k][j] *
1653 mapping_data
1654 .jacobian_pushed_forward_grads[k][d][i][n]) -
1655 (output_data.shape_gradients[first + d][k][n] *
1656 mapping_data
1657 .jacobian_pushed_forward_grads[k][n][i][j]);
1658 for (unsigned int m = 0; m < spacedim; ++m)
1659 transformed_shape_hessians[k][d][i][j] -=
1660 (mapping_data
1661 .jacobian_pushed_forward_grads[k][d][i]
1662 [m] *
1663 mapping_data
1664 .jacobian_pushed_forward_grads[k][m][n]
1665 [j] *
1666 output_data.shape_values(first + n, k)) +
1667 (mapping_data
1668 .jacobian_pushed_forward_grads[k][d][m]
1669 [j] *
1670 mapping_data
1671 .jacobian_pushed_forward_grads[k][m][i]
1672 [n] *
1673 output_data.shape_values(first + n, k));
1674 }
1675
1676 for (unsigned int k = 0; k < n_q_points; ++k)
1677 for (unsigned int d = 0; d < dim; ++d)
1678 output_data.shape_hessians[first + d][k] =
1679 transformed_shape_hessians[k][d];
1680
1681 break;
1682 }
1683
1685 case mapping_piola:
1686 {
1687 for (unsigned int k = 0; k < n_q_points; ++k)
1688 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1689 fe_data.shape_grad_grads[dof_index][k + offset];
1690
1692 transformed_shape_hessians =
1693 make_array_view(fe_data.transformed_shape_hessians,
1694 offset,
1695 n_q_points);
1696 mapping.transform(
1697 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1698 offset,
1699 n_q_points),
1701 mapping_internal,
1702 transformed_shape_hessians);
1703
1704 for (unsigned int k = 0; k < n_q_points; ++k)
1705 for (unsigned int d = 0; d < spacedim; ++d)
1706 for (unsigned int n = 0; n < spacedim; ++n)
1707 for (unsigned int i = 0; i < spacedim; ++i)
1708 for (unsigned int j = 0; j < spacedim; ++j)
1709 {
1710 transformed_shape_hessians[k][d][i][j] +=
1711 (output_data.shape_values(first + n, k) *
1712 mapping_data
1714 [k][d][n][i][j]) +
1715 (output_data.shape_gradients[first + n][k][i] *
1716 mapping_data
1717 .jacobian_pushed_forward_grads[k][d][n][j]) +
1718 (output_data.shape_gradients[first + n][k][j] *
1719 mapping_data
1720 .jacobian_pushed_forward_grads[k][d][i][n]) -
1721 (output_data.shape_gradients[first + d][k][n] *
1722 mapping_data
1723 .jacobian_pushed_forward_grads[k][n][i][j]);
1724
1725 transformed_shape_hessians[k][d][i][j] -=
1726 (output_data.shape_values(first + d, k) *
1727 mapping_data
1729 [k][n][n][i][j]) +
1730 (output_data.shape_gradients[first + d][k][i] *
1731 mapping_data
1732 .jacobian_pushed_forward_grads[k][n][n][j]) +
1733 (output_data.shape_gradients[first + d][k][j] *
1734 mapping_data
1735 .jacobian_pushed_forward_grads[k][n][n][i]);
1736
1737 for (unsigned int m = 0; m < spacedim; ++m)
1738 {
1739 transformed_shape_hessians[k][d][i][j] -=
1740 (mapping_data
1742 [m] *
1743 mapping_data
1745 [j] *
1746 output_data.shape_values(first + n, k)) +
1747 (mapping_data
1748 .jacobian_pushed_forward_grads[k][d][m]
1749 [j] *
1750 mapping_data
1751 .jacobian_pushed_forward_grads[k][m][i]
1752 [n] *
1753 output_data.shape_values(first + n, k));
1754
1755 transformed_shape_hessians[k][d][i][j] +=
1756 (mapping_data
1758 [m] *
1759 mapping_data
1761 [j] *
1762 output_data.shape_values(first + d, k)) +
1763 (mapping_data
1764 .jacobian_pushed_forward_grads[k][n][m]
1765 [j] *
1766 mapping_data
1767 .jacobian_pushed_forward_grads[k][m][i]
1768 [n] *
1769 output_data.shape_values(first + d, k));
1770 }
1771 }
1772
1773 for (unsigned int k = 0; k < n_q_points; ++k)
1774 for (unsigned int d = 0; d < dim; ++d)
1775 output_data.shape_hessians[first + d][k] =
1776 dof_sign * transformed_shape_hessians[k][d];
1777
1778 break;
1779 }
1780
1781 case mapping_nedelec:
1782 {
1783 for (unsigned int k = 0; k < n_q_points; ++k)
1784 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1785 fe_data.shape_grad_grads[dof_index][k + offset];
1786
1788 transformed_shape_hessians =
1789 make_array_view(fe_data.transformed_shape_hessians,
1790 offset,
1791 n_q_points);
1792 mapping.transform(
1793 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1794 offset,
1795 n_q_points),
1797 mapping_internal,
1798 transformed_shape_hessians);
1799
1800 for (unsigned int k = 0; k < n_q_points; ++k)
1801 for (unsigned int d = 0; d < spacedim; ++d)
1802 for (unsigned int n = 0; n < spacedim; ++n)
1803 for (unsigned int i = 0; i < spacedim; ++i)
1804 for (unsigned int j = 0; j < spacedim; ++j)
1805 {
1806 transformed_shape_hessians[k][d][i][j] -=
1807 (output_data.shape_values(first + n, k) *
1808 mapping_data
1810 [k][n][d][i][j]) +
1811 (output_data.shape_gradients[first + d][k][n] *
1812 mapping_data
1813 .jacobian_pushed_forward_grads[k][n][i][j]) +
1814 (output_data.shape_gradients[first + n][k][i] *
1815 mapping_data
1816 .jacobian_pushed_forward_grads[k][n][d][j]) +
1817 (output_data.shape_gradients[first + n][k][j] *
1818 mapping_data
1819 .jacobian_pushed_forward_grads[k][n][i][d]);
1820 }
1821
1822 for (unsigned int k = 0; k < n_q_points; ++k)
1823 for (unsigned int d = 0; d < dim; ++d)
1824 output_data.shape_hessians[first + d][k] =
1825 dof_sign * transformed_shape_hessians[k][d];
1826
1827 break;
1828 }
1829
1830 default:
1832 }
1833 }
1834
1835 // third derivatives are not implemented
1836 if (fe_data.update_each & update_3rd_derivatives)
1837 {
1839 }
1840 }
1841}
1842
1843
1844
1845template <int dim, int spacedim>
1846void
1849 const unsigned int face_no,
1850 const unsigned int sub_no,
1851 const Quadrature<dim - 1> &quadrature,
1852 const Mapping<dim, spacedim> &mapping,
1853 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1855 &mapping_data,
1856 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1858 spacedim>
1859 &output_data) const
1860{
1861 // convert data object to internal
1862 // data for this class. fails with
1863 // an exception if that is not
1864 // possible
1865 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1867 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1868
1869 const unsigned int n_q_points = quadrature.size();
1870
1871 // offset determines which data set
1872 // to take (all data sets for all
1873 // sub-faces are stored contiguously)
1874 const auto offset =
1876 face_no,
1877 sub_no,
1878 cell->combined_face_orientation(
1879 face_no),
1880 n_q_points,
1881 cell->subface_case(face_no));
1882
1883 // TODO: Size assertions
1884
1885 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1886 // correct thing on complicated meshes for higher order Nedelec elements.
1887 // Something similar to FE_Q should be done to permute dofs and to change the
1888 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1889 // class) is preferable.
1890 std::fill(fe_data.dof_sign_change.begin(),
1891 fe_data.dof_sign_change.end(),
1892 1.0);
1893 if (fe_data.update_each & update_values)
1894 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
1895 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
1896
1897 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1898 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1899 // Also nothing in 3d since we take care of it by using the
1900 // adjust_quad_dof_sign_for_face_orientation_table.
1901 if (fe_data.update_each & update_values)
1902 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1903 *this,
1904 this->mapping_kind,
1905 fe_data.dof_sign_change);
1906
1907 // What is the first dof_index on a quad?
1908 const unsigned int first_quad_index = this->get_first_quad_index();
1909 // How many dofs per quad and how many quad dofs do we have at all?
1910 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1911 const unsigned int n_quad_dofs =
1912 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1913
1914 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1915 ++dof_index)
1916 {
1917 /*
1918 * This assumes that the dofs are ordered by first vertices, lines, quads
1919 * and volume dofs. Note that in 2d this always gives false.
1920 */
1921 const bool is_quad_dof =
1922 (dim == 2 ? false :
1923 (first_quad_index <= dof_index) &&
1924 (dof_index < first_quad_index + n_quad_dofs));
1925
1926 // TODO: This hack is not pretty and it is only here to handle the 2d
1927 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1928 // handled by the
1929 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1930 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1931 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1932 // it is handles with a table. This array is allocated in
1933 // fe_poly_tensor.h.
1934 double dof_sign = 1.0;
1935 // under some circumstances fe_data.dof_sign_change is not allocated
1936 if (fe_data.update_each & update_values)
1937 dof_sign = fe_data.dof_sign_change[dof_index];
1938
1939 if (is_quad_dof)
1940 {
1941 /*
1942 * Find the face belonging to this dof_index. This is integer
1943 * division.
1944 */
1945 unsigned int face_index_from_dof_index =
1946 (dof_index - first_quad_index) / (n_dofs_per_quad);
1947
1948 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1949
1950 // Correct the dof_sign if necessary
1951 if (adjust_quad_dof_sign_for_face_orientation(
1952 local_quad_dof_index,
1953 face_index_from_dof_index,
1954 cell->face_orientation(face_index_from_dof_index),
1955 cell->face_flip(face_index_from_dof_index),
1956 cell->face_rotation(face_index_from_dof_index)))
1957 dof_sign = -1.0;
1958 }
1959
1960 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1961
1962 const unsigned int first =
1963 output_data.shape_function_to_row_table
1964 [dof_index * this->n_components() +
1965 this->get_nonzero_components(dof_index).first_selected_component()];
1966
1967 if (fe_data.update_each & update_values)
1968 {
1969 switch (mapping_kind)
1970 {
1971 case mapping_none:
1972 {
1973 for (unsigned int k = 0; k < n_q_points; ++k)
1974 for (unsigned int d = 0; d < dim; ++d)
1975 output_data.shape_values(first + d, k) =
1976 fe_data.shape_values[dof_index][k + offset][d];
1977 break;
1978 }
1979
1980 case mapping_covariant:
1982 {
1984 transformed_shape_values =
1985 make_array_view(fe_data.transformed_shape_values,
1986 offset,
1987 n_q_points);
1988 mapping.transform(make_array_view(fe_data.shape_values,
1989 dof_index,
1990 offset,
1991 n_q_points),
1992 mapping_kind,
1993 mapping_internal,
1994 transformed_shape_values);
1995
1996 for (unsigned int k = 0; k < n_q_points; ++k)
1997 for (unsigned int d = 0; d < dim; ++d)
1998 output_data.shape_values(first + d, k) =
1999 transformed_shape_values[k][d];
2000
2001 break;
2002 }
2003
2005 case mapping_piola:
2006 {
2008 transformed_shape_values =
2009 make_array_view(fe_data.transformed_shape_values,
2010 offset,
2011 n_q_points);
2012
2013 mapping.transform(make_array_view(fe_data.shape_values,
2014 dof_index,
2015 offset,
2016 n_q_points),
2018 mapping_internal,
2019 transformed_shape_values);
2020 for (unsigned int k = 0; k < n_q_points; ++k)
2021 for (unsigned int d = 0; d < dim; ++d)
2022 output_data.shape_values(first + d, k) =
2023 dof_sign * transformed_shape_values[k][d];
2024 break;
2025 }
2026
2027 case mapping_nedelec:
2028 {
2030 transformed_shape_values =
2031 make_array_view(fe_data.transformed_shape_values,
2032 offset,
2033 n_q_points);
2034
2035 mapping.transform(make_array_view(fe_data.shape_values,
2036 dof_index,
2037 offset,
2038 n_q_points),
2040 mapping_internal,
2041 transformed_shape_values);
2042
2043 for (unsigned int k = 0; k < n_q_points; ++k)
2044 for (unsigned int d = 0; d < dim; ++d)
2045 output_data.shape_values(first + d, k) =
2046 dof_sign * transformed_shape_values[k][d];
2047
2048 break;
2049 }
2050
2051 default:
2053 }
2054 }
2055
2056 if (fe_data.update_each & update_gradients)
2057 {
2058 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
2059 make_array_view(fe_data.transformed_shape_grads,
2060 offset,
2061 n_q_points);
2062 switch (mapping_kind)
2063 {
2064 case mapping_none:
2065 {
2066 mapping.transform(make_array_view(fe_data.shape_grads,
2067 dof_index,
2068 offset,
2069 n_q_points),
2071 mapping_internal,
2072 transformed_shape_grads);
2073 for (unsigned int k = 0; k < n_q_points; ++k)
2074 for (unsigned int d = 0; d < dim; ++d)
2075 output_data.shape_gradients[first + d][k] =
2076 transformed_shape_grads[k][d];
2077 break;
2078 }
2079
2080 case mapping_covariant:
2081 {
2082 mapping.transform(make_array_view(fe_data.shape_grads,
2083 dof_index,
2084 offset,
2085 n_q_points),
2087 mapping_internal,
2088 transformed_shape_grads);
2089
2090 for (unsigned int k = 0; k < n_q_points; ++k)
2091 for (unsigned int d = 0; d < spacedim; ++d)
2092 for (unsigned int n = 0; n < spacedim; ++n)
2093 transformed_shape_grads[k][d] -=
2094 output_data.shape_values(first + n, k) *
2095 mapping_data.jacobian_pushed_forward_grads[k][n][d];
2096
2097 for (unsigned int k = 0; k < n_q_points; ++k)
2098 for (unsigned int d = 0; d < dim; ++d)
2099 output_data.shape_gradients[first + d][k] =
2100 transformed_shape_grads[k][d];
2101
2102 break;
2103 }
2104
2106 {
2107 for (unsigned int k = 0; k < n_q_points; ++k)
2108 fe_data.untransformed_shape_grads[k + offset] =
2109 fe_data.shape_grads[dof_index][k + offset];
2110
2111 mapping.transform(
2112 make_array_view(fe_data.untransformed_shape_grads,
2113 offset,
2114 n_q_points),
2116 mapping_internal,
2117 transformed_shape_grads);
2118
2119 for (unsigned int k = 0; k < n_q_points; ++k)
2120 for (unsigned int d = 0; d < spacedim; ++d)
2121 for (unsigned int n = 0; n < spacedim; ++n)
2122 transformed_shape_grads[k][d] +=
2123 output_data.shape_values(first + n, k) *
2124 mapping_data.jacobian_pushed_forward_grads[k][d][n];
2125
2126 for (unsigned int k = 0; k < n_q_points; ++k)
2127 for (unsigned int d = 0; d < dim; ++d)
2128 output_data.shape_gradients[first + d][k] =
2129 transformed_shape_grads[k][d];
2130
2131 break;
2132 }
2133
2135 case mapping_piola:
2136 {
2137 for (unsigned int k = 0; k < n_q_points; ++k)
2138 fe_data.untransformed_shape_grads[k + offset] =
2139 fe_data.shape_grads[dof_index][k + offset];
2140
2141 mapping.transform(
2142 make_array_view(fe_data.untransformed_shape_grads,
2143 offset,
2144 n_q_points),
2146 mapping_internal,
2147 transformed_shape_grads);
2148
2149 for (unsigned int k = 0; k < n_q_points; ++k)
2150 for (unsigned int d = 0; d < spacedim; ++d)
2151 for (unsigned int n = 0; n < spacedim; ++n)
2152 transformed_shape_grads[k][d] +=
2153 (output_data.shape_values(first + n, k) *
2154 mapping_data
2156 (output_data.shape_values(first + d, k) *
2157 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
2158
2159 for (unsigned int k = 0; k < n_q_points; ++k)
2160 for (unsigned int d = 0; d < dim; ++d)
2161 output_data.shape_gradients[first + d][k] =
2162 dof_sign * transformed_shape_grads[k][d];
2163
2164 break;
2165 }
2166
2167 case mapping_nedelec:
2168 {
2169 // this particular shape
2170 // function at all
2171 // q-points. if Dv is the
2172 // gradient of the shape
2173 // function on the unit
2174 // cell, then
2175 // (J^-T)Dv(J^-1) is the
2176 // value we want to have on
2177 // the real cell.
2178 for (unsigned int k = 0; k < n_q_points; ++k)
2179 fe_data.untransformed_shape_grads[k + offset] =
2180 fe_data.shape_grads[dof_index][k + offset];
2181
2182 mapping.transform(
2183 make_array_view(fe_data.untransformed_shape_grads,
2184 offset,
2185 n_q_points),
2187 mapping_internal,
2188 transformed_shape_grads);
2189
2190 for (unsigned int k = 0; k < n_q_points; ++k)
2191 for (unsigned int d = 0; d < spacedim; ++d)
2192 for (unsigned int n = 0; n < spacedim; ++n)
2193 transformed_shape_grads[k][d] -=
2194 output_data.shape_values(first + n, k) *
2195 mapping_data.jacobian_pushed_forward_grads[k][n][d];
2196
2197 for (unsigned int k = 0; k < n_q_points; ++k)
2198 for (unsigned int d = 0; d < dim; ++d)
2199 output_data.shape_gradients[first + d][k] =
2200 dof_sign * transformed_shape_grads[k][d];
2201
2202 break;
2203 }
2204
2205 default:
2207 }
2208 }
2209
2210 if (fe_data.update_each & update_hessians)
2211 {
2212 switch (mapping_kind)
2213 {
2214 case mapping_none:
2215 {
2217 transformed_shape_hessians =
2218 make_array_view(fe_data.transformed_shape_hessians,
2219 offset,
2220 n_q_points);
2221 mapping.transform(make_array_view(fe_data.shape_grad_grads,
2222 dof_index,
2223 offset,
2224 n_q_points),
2226 mapping_internal,
2227 transformed_shape_hessians);
2228
2229 for (unsigned int k = 0; k < n_q_points; ++k)
2230 for (unsigned int d = 0; d < spacedim; ++d)
2231 for (unsigned int n = 0; n < spacedim; ++n)
2232 transformed_shape_hessians[k][d] -=
2233 output_data.shape_gradients[first + d][k][n] *
2234 mapping_data.jacobian_pushed_forward_grads[k][n];
2235
2236 for (unsigned int k = 0; k < n_q_points; ++k)
2237 for (unsigned int d = 0; d < dim; ++d)
2238 output_data.shape_hessians[first + d][k] =
2239 transformed_shape_hessians[k][d];
2240
2241 break;
2242 }
2243 case mapping_covariant:
2244 {
2245 for (unsigned int k = 0; k < n_q_points; ++k)
2246 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2247 fe_data.shape_grad_grads[dof_index][k + offset];
2248
2250 transformed_shape_hessians =
2251 make_array_view(fe_data.transformed_shape_hessians,
2252 offset,
2253 n_q_points);
2254 mapping.transform(
2255 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2256 offset,
2257 n_q_points),
2259 mapping_internal,
2260 transformed_shape_hessians);
2261
2262 for (unsigned int k = 0; k < n_q_points; ++k)
2263 for (unsigned int d = 0; d < spacedim; ++d)
2264 for (unsigned int n = 0; n < spacedim; ++n)
2265 for (unsigned int i = 0; i < spacedim; ++i)
2266 for (unsigned int j = 0; j < spacedim; ++j)
2267 {
2268 transformed_shape_hessians[k][d][i][j] -=
2269 (output_data.shape_values(first + n, k) *
2270 mapping_data
2272 [k][n][d][i][j]) +
2273 (output_data.shape_gradients[first + d][k][n] *
2274 mapping_data
2275 .jacobian_pushed_forward_grads[k][n][i][j]) +
2276 (output_data.shape_gradients[first + n][k][i] *
2277 mapping_data
2278 .jacobian_pushed_forward_grads[k][n][d][j]) +
2279 (output_data.shape_gradients[first + n][k][j] *
2280 mapping_data
2281 .jacobian_pushed_forward_grads[k][n][i][d]);
2282 }
2283
2284 for (unsigned int k = 0; k < n_q_points; ++k)
2285 for (unsigned int d = 0; d < dim; ++d)
2286 output_data.shape_hessians[first + d][k] =
2287 transformed_shape_hessians[k][d];
2288
2289 break;
2290 }
2291
2293 {
2294 for (unsigned int k = 0; k < n_q_points; ++k)
2295 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2296 fe_data.shape_grad_grads[dof_index][k + offset];
2297
2299 transformed_shape_hessians =
2300 make_array_view(fe_data.transformed_shape_hessians,
2301 offset,
2302 n_q_points);
2303 mapping.transform(
2304 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2305 offset,
2306 n_q_points),
2308 mapping_internal,
2309 transformed_shape_hessians);
2310
2311 for (unsigned int k = 0; k < n_q_points; ++k)
2312 for (unsigned int d = 0; d < spacedim; ++d)
2313 for (unsigned int n = 0; n < spacedim; ++n)
2314 for (unsigned int i = 0; i < spacedim; ++i)
2315 for (unsigned int j = 0; j < spacedim; ++j)
2316 {
2317 transformed_shape_hessians[k][d][i][j] +=
2318 (output_data.shape_values(first + n, k) *
2319 mapping_data
2321 [k][d][n][i][j]) +
2322 (output_data.shape_gradients[first + n][k][i] *
2323 mapping_data
2324 .jacobian_pushed_forward_grads[k][d][n][j]) +
2325 (output_data.shape_gradients[first + n][k][j] *
2326 mapping_data
2327 .jacobian_pushed_forward_grads[k][d][i][n]) -
2328 (output_data.shape_gradients[first + d][k][n] *
2329 mapping_data
2330 .jacobian_pushed_forward_grads[k][n][i][j]);
2331 for (unsigned int m = 0; m < spacedim; ++m)
2332 transformed_shape_hessians[k][d][i][j] -=
2333 (mapping_data
2334 .jacobian_pushed_forward_grads[k][d][i]
2335 [m] *
2336 mapping_data
2337 .jacobian_pushed_forward_grads[k][m][n]
2338 [j] *
2339 output_data.shape_values(first + n, k)) +
2340 (mapping_data
2341 .jacobian_pushed_forward_grads[k][d][m]
2342 [j] *
2343 mapping_data
2344 .jacobian_pushed_forward_grads[k][m][i]
2345 [n] *
2346 output_data.shape_values(first + n, k));
2347 }
2348
2349 for (unsigned int k = 0; k < n_q_points; ++k)
2350 for (unsigned int d = 0; d < dim; ++d)
2351 output_data.shape_hessians[first + d][k] =
2352 transformed_shape_hessians[k][d];
2353
2354 break;
2355 }
2356
2358 case mapping_piola:
2359 {
2360 for (unsigned int k = 0; k < n_q_points; ++k)
2361 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2362 fe_data.shape_grad_grads[dof_index][k + offset];
2363
2365 transformed_shape_hessians =
2366 make_array_view(fe_data.transformed_shape_hessians,
2367 offset,
2368 n_q_points);
2369 mapping.transform(
2370 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2371 offset,
2372 n_q_points),
2374 mapping_internal,
2375 transformed_shape_hessians);
2376
2377 for (unsigned int k = 0; k < n_q_points; ++k)
2378 for (unsigned int d = 0; d < spacedim; ++d)
2379 for (unsigned int n = 0; n < spacedim; ++n)
2380 for (unsigned int i = 0; i < spacedim; ++i)
2381 for (unsigned int j = 0; j < spacedim; ++j)
2382 {
2383 transformed_shape_hessians[k][d][i][j] +=
2384 (output_data.shape_values(first + n, k) *
2385 mapping_data
2387 [k][d][n][i][j]) +
2388 (output_data.shape_gradients[first + n][k][i] *
2389 mapping_data
2390 .jacobian_pushed_forward_grads[k][d][n][j]) +
2391 (output_data.shape_gradients[first + n][k][j] *
2392 mapping_data
2393 .jacobian_pushed_forward_grads[k][d][i][n]) -
2394 (output_data.shape_gradients[first + d][k][n] *
2395 mapping_data
2396 .jacobian_pushed_forward_grads[k][n][i][j]);
2397
2398 transformed_shape_hessians[k][d][i][j] -=
2399 (output_data.shape_values(first + d, k) *
2400 mapping_data
2402 [k][n][n][i][j]) +
2403 (output_data.shape_gradients[first + d][k][i] *
2404 mapping_data
2405 .jacobian_pushed_forward_grads[k][n][n][j]) +
2406 (output_data.shape_gradients[first + d][k][j] *
2407 mapping_data
2408 .jacobian_pushed_forward_grads[k][n][n][i]);
2409 for (unsigned int m = 0; m < spacedim; ++m)
2410 {
2411 transformed_shape_hessians[k][d][i][j] -=
2412 (mapping_data
2414 [m] *
2415 mapping_data
2417 [j] *
2418 output_data.shape_values(first + n, k)) +
2419 (mapping_data
2420 .jacobian_pushed_forward_grads[k][d][m]
2421 [j] *
2422 mapping_data
2423 .jacobian_pushed_forward_grads[k][m][i]
2424 [n] *
2425 output_data.shape_values(first + n, k));
2426
2427 transformed_shape_hessians[k][d][i][j] +=
2428 (mapping_data
2430 [m] *
2431 mapping_data
2433 [j] *
2434 output_data.shape_values(first + d, k)) +
2435 (mapping_data
2436 .jacobian_pushed_forward_grads[k][n][m]
2437 [j] *
2438 mapping_data
2439 .jacobian_pushed_forward_grads[k][m][i]
2440 [n] *
2441 output_data.shape_values(first + d, k));
2442 }
2443 }
2444
2445 for (unsigned int k = 0; k < n_q_points; ++k)
2446 for (unsigned int d = 0; d < dim; ++d)
2447 output_data.shape_hessians[first + d][k] =
2448 dof_sign * transformed_shape_hessians[k][d];
2449
2450 break;
2451 }
2452
2453 case mapping_nedelec:
2454 {
2455 for (unsigned int k = 0; k < n_q_points; ++k)
2456 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2457 fe_data.shape_grad_grads[dof_index][k + offset];
2458
2460 transformed_shape_hessians =
2461 make_array_view(fe_data.transformed_shape_hessians,
2462 offset,
2463 n_q_points);
2464 mapping.transform(
2465 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2466 offset,
2467 n_q_points),
2469 mapping_internal,
2470 transformed_shape_hessians);
2471
2472 for (unsigned int k = 0; k < n_q_points; ++k)
2473 for (unsigned int d = 0; d < spacedim; ++d)
2474 for (unsigned int n = 0; n < spacedim; ++n)
2475 for (unsigned int i = 0; i < spacedim; ++i)
2476 for (unsigned int j = 0; j < spacedim; ++j)
2477 {
2478 transformed_shape_hessians[k][d][i][j] -=
2479 (output_data.shape_values(first + n, k) *
2480 mapping_data
2482 [k][n][d][i][j]) +
2483 (output_data.shape_gradients[first + d][k][n] *
2484 mapping_data
2485 .jacobian_pushed_forward_grads[k][n][i][j]) +
2486 (output_data.shape_gradients[first + n][k][i] *
2487 mapping_data
2488 .jacobian_pushed_forward_grads[k][n][d][j]) +
2489 (output_data.shape_gradients[first + n][k][j] *
2490 mapping_data
2491 .jacobian_pushed_forward_grads[k][n][i][d]);
2492 }
2493
2494 for (unsigned int k = 0; k < n_q_points; ++k)
2495 for (unsigned int d = 0; d < dim; ++d)
2496 output_data.shape_hessians[first + d][k] =
2497 dof_sign * transformed_shape_hessians[k][d];
2498
2499 break;
2500 }
2501
2502 default:
2504 }
2505 }
2506
2507 // third derivatives are not implemented
2508 if (fe_data.update_each & update_3rd_derivatives)
2509 {
2511 }
2512 }
2513}
2514
2515
2516
2517template <int dim, int spacedim>
2520 const UpdateFlags flags) const
2521{
2523
2524 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2525 {
2526 const MappingKind mapping_kind = get_mapping_kind(i);
2527
2528 switch (mapping_kind)
2529 {
2530 case mapping_none:
2531 {
2532 if (flags & update_values)
2533 out |= update_values;
2534
2535 if (flags & update_gradients)
2538
2539 if (flags & update_hessians)
2543 break;
2544 }
2546 case mapping_piola:
2547 {
2548 if (flags & update_values)
2549 out |= update_values | update_piola;
2550
2551 if (flags & update_gradients)
2556
2557 if (flags & update_hessians)
2562
2563 break;
2564 }
2565
2566
2568 {
2569 if (flags & update_values)
2570 out |= update_values | update_piola;
2571
2572 if (flags & update_gradients)
2577
2578 if (flags & update_hessians)
2583
2584 break;
2585 }
2586
2587 case mapping_nedelec:
2588 case mapping_covariant:
2589 {
2590 if (flags & update_values)
2592
2593 if (flags & update_gradients)
2597
2598 if (flags & update_hessians)
2603
2604 break;
2605 }
2606
2607 default:
2608 {
2610 }
2611 }
2612 }
2613
2614 return out;
2615}
2616
2617#endif
2618// explicit instantiations
2619#include "fe/fe_poly_tensor.inst"
2620
2621
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:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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:352