deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_poly_tensor.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2005 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#include <deal.II/base/array_view.h>
15#include <deal.II/base/derivative_form.h>
16#include <deal.II/base/polynomials_abf.h>
17#include <deal.II/base/polynomials_bdm.h>
18#include <deal.II/base/polynomials_bernardi_raugel.h>
19#include <deal.II/base/polynomials_nedelec.h>
20#include <deal.II/base/polynomials_raviart_thomas.h>
21#include <deal.II/base/polynomials_rt_bubbles.h>
22#include <deal.II/base/qprojector.h>
23
24#include <deal.II/fe/fe_poly_tensor.h>
25#include <deal.II/fe/fe_values.h>
26#include <deal.II/fe/mapping_cartesian.h>
27
28#include <deal.II/grid/cell_id.h>
29#include <deal.II/grid/tria.h>
30#include <deal.II/grid/tria_orientation.h>
31
32DEAL_II_NAMESPACE_OPEN
33
34namespace internal
35{
36 namespace FE_PolyTensor
37 {
38 namespace
39 {
40 template <int spacedim>
41 void
42 get_dof_sign_change_h_div(
43 const typename ::Triangulation<1, spacedim>::cell_iterator &,
44 const FiniteElement<1, spacedim> &,
45 const std::vector<MappingKind> &,
46 std::vector<double> &)
47 {
48 // Nothing to do in 1d.
49 }
50
51
52
53 // TODO: This function is not a consistent fix of the orientation issue
54 // like in 3d. It is rather kept not to break legacy behavior in 2d but
55 // should be replaced. See also the implementation of
56 // FE_RaviartThomas<dim>::initialize_quad_dof_index_permutation_and_sign_change()
57 // or other H(div) conforming elements such as FE_ABF<dim> and
58 // FE_BDM<dim>.
59 template <int spacedim>
60 void
61 get_dof_sign_change_h_div(
62 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
63 const FiniteElement<2, spacedim> &fe,
64 const std::vector<MappingKind> &mapping_kind,
65 std::vector<double> &face_sign)
66 {
67 const unsigned int dim = 2;
68 // const unsigned int spacedim = 2;
69
70 const CellId this_cell_id = cell->id();
71
72 for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
73 {
74 typename ::Triangulation<dim, spacedim>::face_iterator face =
75 cell->face(f);
76 if (!face->at_boundary())
77 {
78 const unsigned int nn = cell->neighbor_face_no(f);
79 const typename ::Triangulation<dim,
80 spacedim>::cell_iterator
81 neighbor_cell_at_face = cell->neighbor(f);
82 const CellId neighbor_cell_id = neighbor_cell_at_face->id();
83
84 // Only fix sign if the orientation is opposite and only do so
85 // on the face dofs on the cell with smaller cell_id.
86 if (((nn + f) % 2 == 0) && this_cell_id < neighbor_cell_id)
87 for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j)
88 {
89 const unsigned int cell_j = fe.face_to_cell_index(j, f);
90
91 Assert(f * fe.n_dofs_per_face(f) + j < face_sign.size(),
92 ExcInternalError());
93 Assert(mapping_kind.size() == 1 ||
94 cell_j < mapping_kind.size(),
95 ExcInternalError());
96
97 // TODO: This is probably only going to work for those
98 // elements for which all dofs are face dofs
99 if ((mapping_kind.size() > 1 ?
100 mapping_kind[cell_j] :
101 mapping_kind[0]) == mapping_raviart_thomas)
102 face_sign[f * fe.n_dofs_per_face(f) + j] = -1.0;
103 }
104 }
105 }
106 }
107
108
109
110 template <int spacedim>
111 void
112 get_dof_sign_change_h_div(
113 const typename ::Triangulation<3, spacedim>::cell_iterator
114 & /*cell*/,
115 const FiniteElement<3, spacedim> & /*fe*/,
116 const std::vector<MappingKind> & /*mapping_kind*/,
117 std::vector<double> & /*face_sign*/)
118 {
119 // Nothing to do. In 3d we take care of it through the
120 // adjust_quad_dof_sign_for_face_orientation_table
121 }
122
123 template <int spacedim>
124 void
125 get_dof_sign_change_nedelec(
126 const typename ::Triangulation<1, spacedim>::cell_iterator
127 & /*cell*/,
128 const FiniteElement<1, spacedim> & /*fe*/,
129 const std::vector<MappingKind> & /*mapping_kind*/,
130 std::vector<double> & /*line_dof_sign*/)
131 {
132 // nothing to do in 1d
133 }
134
135 template <int spacedim>
136 void
137 get_dof_sign_change_nedelec(
138 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
139 const FiniteElement<2, spacedim> &fe,
140 const std::vector<MappingKind> &mapping_kind,
141 std::vector<double> &line_dof_sign)
142 {
143 // The Nedelec finite elements in two spatial dimensions have two types
144 // of dofs: the line dofs and the quad dofs. The line dofs are
145 // associated with the edges. They are shared between the neighbouring
146 // cells and, consequently, need to be adjusted to compensate for a
147 // possible mismatch in the edge orientation. The quad dofs, they are
148 // associated with the interiors of the cells, are not shared between
149 // the cells in two spatial dimensions and need no adjustments.
150 //
151 // The Nedelec finite elements in two spatial dimensions have
152 // 2*(k+1)*(k+2) dofs per cell. All dofs are distributed between the
153 // line and quad dofs as the following:
154 //
155 // 4*(k+1) line dofs; (k+1) dofs per line.
156 // 2*k*(k+1) quad dofs.
157 //
158 // The dofs are indexed in the following order: first all line dofs,
159 // then all quad dofs.
160 //
161 // Here we adjust the line dofs. The sign of a line dof needs to be
162 // changed if the edge on which the dof resides points in the opposite
163 // direction.
164 const unsigned int k = fe.tensor_degree() - 1;
165
166 for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell; ++l)
167 if (cell->line_orientation(l) !=
168 numbers::default_geometric_orientation &&
169 mapping_kind[0] == mapping_nedelec)
170 {
171 if (k == 0)
172 {
173 // The lowest order element (k=0) is straightforward, because
174 // there is a single dof per edge, which needs to be flipped:
175 line_dof_sign[l] = -1.0;
176 }
177 else
178 {
179 // The case k > 0 is a bit more complicated. As we adjust
180 // only edge dofs in this function, we need to concern
181 // ourselves with the first 4*(k+1) entries in line_dof_sign
182 // vector ignoring the rest. There are (k+1) dofs per edge.
183 // Let us consider the local dof indices on one edge,
184 // local_line_dof = 0...k. The shape functions with even
185 // indices are asymmetric. The corresponding dofs need sign
186 // adjustment if the edge points in the opposite direction.
187 // The shape functions with odd indices are symmetric. The
188 // corresponding dofs need no sign adjustment even if the edge
189 // points in the opposite direction. In the current context
190 // the notion of symmetry of a shape function means that a
191 // shape function looks exactly the same if it is looked upon
192 // from the centers of the two neighbouring cells that share
193 // it.
194 for (unsigned int local_line_dof = 0;
195 local_line_dof < (k + 1);
196 local_line_dof++)
197 if (local_line_dof % 2 == 0)
198 line_dof_sign[local_line_dof + l * (k + 1)] = -1.0;
199 }
200 }
201 }
202
203 template <int spacedim>
204 void
205 get_dof_sign_change_nedelec(
206 const typename ::Triangulation<3, spacedim>::cell_iterator &cell,
207 const FiniteElement<3, spacedim> &fe,
208 const std::vector<MappingKind> &mapping_kind,
209 std::vector<double> &line_dof_sign)
210 {
211 // This function does half of the job - it adjusts the sign of the
212 // line (edge) dofs. In the three-dimensional space the quad (face) dofs
213 // need to be adjusted as well. The quad dofs are treated by
214 // FE_Nedelec<dim>::initialize_quad_dof_index_permutation_and_sign_change()
215 // in fe_nedelec.cc. The dofs associated with the interior of the cells,
216 // the hex dofs, need no adjustments as they are not shared between the
217 // neighboring cells.
218
219 const unsigned int k = fe.tensor_degree() - 1;
220 // The order of the Nedelec elements equals the tensor degree minus one,
221 // k = n - 1. In the three-dimensional space the Nedelec elements of the
222 // lowermost order, k = 0, have only 12 line (edge) dofs. The Nedelec
223 // elements of the higher orders, k > 0, have 3*(k+1)*(k+2)^2 dofs in
224 // total if dim=3. The dofs in a cell are distributed between lines
225 // (edges), quads (faces), and the hex (the interior of the cell) as the
226 // following:
227 //
228 // 12*(k+1) line dofs; (k+1) dofs per line.
229 // 2*6*k*(k+1) quad dofs; 2*k*(k+1) dofs per quad.
230 // 3*(k+2)^2*(k+1) hex dofs.
231 //
232 // The dofs are indexed in the following order: first all line dofs,
233 // then all quad dofs, and then all hex dofs.
234 //
235 // Here we adjust only the line (edge) dofs. The line dofs need only
236 // sign adjustment. That is, no permutation of the line dofs is needed.
237 for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
238 if (cell->line_orientation(l) !=
239 numbers::default_geometric_orientation &&
240 mapping_kind[0] == mapping_nedelec)
241 {
242 if (k == 0)
243 {
244 // The lowest order element (k=0) is straightforward, because
245 // there is a single dof per edge, which needs to be flipped:
246 line_dof_sign[l] = -1.0;
247 }
248 else
249 {
250 // The case k > 0 is a bit more complicated. As we adjust
251 // only edge dofs in this function, we need to concern
252 // ourselves with the first 12*(k+1) entries in line_dof_sign
253 // vector ignoring the rest. There are (k+1) dofs per edge.
254 // Let us consider the local dof indices on one edge,
255 // local_line_dof = 0...k. The shape functions with even
256 // indices are asymmetric. The corresponding dofs need sign
257 // adjustment if the edge points in the opposite direction.
258 // The shape functions with odd indices are symmetric. The
259 // corresponding dofs need no sign adjustment even if the edge
260 // points in the opposite direction. In the current context
261 // the notion of symmetry of a shape function means that a
262 // shape function looks exactly the same if it is looked upon
263 // from the centers of the two neighbouring cells that share
264 // it.
265 for (unsigned int local_line_dof = 0;
266 local_line_dof < (k + 1);
267 local_line_dof++)
268 if (local_line_dof % 2 == 0)
269 line_dof_sign[local_line_dof + l * (k + 1)] = -1.0;
270 }
271 }
272 }
273 } // namespace
274 } // namespace FE_PolyTensor
275} // namespace internal
276
277#ifndef DOXYGEN
278
279template <int dim, int spacedim>
280FE_PolyTensor<dim, spacedim>::FE_PolyTensor(
281 const TensorPolynomialsBase<dim> &polynomials,
282 const FiniteElementData<dim> &fe_data,
283 const std::vector<bool> &restriction_is_additive_flags,
284 const std::vector<ComponentMask> &nonzero_components)
285 : FiniteElement<dim, spacedim>(fe_data,
286 restriction_is_additive_flags,
287 nonzero_components)
288 , mapping_kind({MappingKind::mapping_none})
289 , poly_space(polynomials.clone())
290{
291 cached_point[0] = -1;
292 // Set up the table converting
293 // components to base
294 // components. Since we have only
295 // one base element, everything
296 // remains zero except the
297 // component in the base, which is
298 // the component itself
299 for (unsigned int comp = 0; comp < this->n_components(); ++comp)
300 this->component_to_base_table[comp].first.second = comp;
301
302 if (dim == 3)
303 {
304 adjust_quad_dof_sign_for_face_orientation_table.resize(
305 this->n_unique_2d_subobjects());
306
307 for (unsigned int f = 0; f < this->n_unique_2d_subobjects(); ++f)
308 {
309 adjust_quad_dof_sign_for_face_orientation_table[f] =
310 Table<2, bool>(this->n_dofs_per_quad(f),
311 this->reference_cell().n_face_orientations(f));
312 adjust_quad_dof_sign_for_face_orientation_table[f].fill(false);
313 }
314 }
315}
316
317
318
319template <int dim, int spacedim>
320FE_PolyTensor<dim, spacedim>::FE_PolyTensor(const FE_PolyTensor &fe)
321 : FiniteElement<dim, spacedim>(fe)
322 , mapping_kind(fe.mapping_kind)
323 , adjust_quad_dof_sign_for_face_orientation_table(
324 fe.adjust_quad_dof_sign_for_face_orientation_table)
325 , poly_space(fe.poly_space->clone())
326 , inverse_node_matrix(fe.inverse_node_matrix)
327{}
328
329
330
331template <int dim, int spacedim>
332bool
333FE_PolyTensor<dim, spacedim>::single_mapping_kind() const
334{
335 return mapping_kind.size() == 1;
336}
337
338
339template <int dim, int spacedim>
340bool
341FE_PolyTensor<dim, spacedim>::adjust_quad_dof_sign_for_face_orientation(
342 const unsigned int index,
343 const unsigned int face,
344 const types::geometric_orientation combined_orientation) const
345{
346 // do nothing in 1d and 2d
347 if (dim < 3)
348 return false;
349
350 // The exception are discontinuous
351 // elements for which there should be no
352 // face dofs anyway (i.e. dofs_per_quad==0
353 // in 3d), so we don't need the table, but
354 // the function should also not have been
355 // called
356 AssertIndexRange(index, this->n_dofs_per_quad(face));
357 Assert(adjust_quad_dof_sign_for_face_orientation_table
358 [this->n_unique_2d_subobjects() == 1 ? 0 : face]
359 .n_elements() ==
360 this->reference_cell().n_face_orientations(face) *
361 this->n_dofs_per_quad(face),
362 ExcInternalError());
363
364 return adjust_quad_dof_sign_for_face_orientation_table
365 [this->n_unique_2d_subobjects() == 1 ? 0 : face](index,
366 combined_orientation);
367}
368
369
370template <int dim, int spacedim>
371MappingKind
372FE_PolyTensor<dim, spacedim>::get_mapping_kind(const unsigned int i) const
373{
374 if (single_mapping_kind())
375 return mapping_kind[0];
376
377 AssertIndexRange(i, mapping_kind.size());
378 return mapping_kind[i];
379}
380
381
382
383template <int dim, int spacedim>
384double
385FE_PolyTensor<dim, spacedim>::shape_value(const unsigned int,
386 const Point<dim> &) const
387
388{
389 Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
390 return 0.;
391}
392
393
394
395template <int dim, int spacedim>
396double
397FE_PolyTensor<dim, spacedim>::shape_value_component(
398 const unsigned int i,
399 const Point<dim> &p,
400 const unsigned int component) const
401{
402 AssertIndexRange(i, this->n_dofs_per_cell());
403 AssertIndexRange(component, dim);
404
405 std::scoped_lock lock(cache_mutex);
406
407 if (cached_point != p || cached_values.empty())
408 {
409 cached_point = p;
410 cached_values.resize(poly_space->n());
411
412 std::vector<Tensor<4, dim>> dummy1;
413 std::vector<Tensor<5, dim>> dummy2;
414 poly_space->evaluate(
415 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
416 }
417
418 double s = 0;
419 if (inverse_node_matrix.n_cols() == 0)
420 return cached_values[i][component];
421 else
422 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
423 s += inverse_node_matrix(j, i) * cached_values[j][component];
424 return s;
425}
426
427
428
429template <int dim, int spacedim>
430Tensor<1, dim>
431FE_PolyTensor<dim, spacedim>::shape_grad(const unsigned int,
432 const Point<dim> &) const
433{
434 Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
435 return Tensor<1, dim>();
436}
437
438
439
440template <int dim, int spacedim>
441Tensor<1, dim>
442FE_PolyTensor<dim, spacedim>::shape_grad_component(
443 const unsigned int i,
444 const Point<dim> &p,
445 const unsigned int component) const
446{
447 AssertIndexRange(i, this->n_dofs_per_cell());
448 AssertIndexRange(component, dim);
449
450 std::scoped_lock lock(cache_mutex);
451
452 if (cached_point != p || cached_grads.empty())
453 {
454 cached_point = p;
455 cached_grads.resize(poly_space->n());
456
457 std::vector<Tensor<4, dim>> dummy1;
458 std::vector<Tensor<5, dim>> dummy2;
459 poly_space->evaluate(
460 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
461 }
462
463 Tensor<1, dim> s;
464 if (inverse_node_matrix.n_cols() == 0)
465 return cached_grads[i][component];
466 else
467 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
468 s += inverse_node_matrix(j, i) * cached_grads[j][component];
469
470 return s;
471}
472
473
474
475template <int dim, int spacedim>
476Tensor<2, dim>
477FE_PolyTensor<dim, spacedim>::shape_grad_grad(const unsigned int,
478 const Point<dim> &) const
479{
480 Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
481 return Tensor<2, dim>();
482}
483
484
485
486template <int dim, int spacedim>
487Tensor<2, dim>
488FE_PolyTensor<dim, spacedim>::shape_grad_grad_component(
489 const unsigned int i,
490 const Point<dim> &p,
491 const unsigned int component) const
492{
493 AssertIndexRange(i, this->n_dofs_per_cell());
494 AssertIndexRange(component, dim);
495
496 std::scoped_lock lock(cache_mutex);
497
498 if (cached_point != p || cached_grad_grads.empty())
499 {
500 cached_point = p;
501 cached_grad_grads.resize(poly_space->n());
502
503 std::vector<Tensor<4, dim>> dummy1;
504 std::vector<Tensor<5, dim>> dummy2;
505 poly_space->evaluate(
506 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
507 }
508
509 Tensor<2, dim> s;
510 if (inverse_node_matrix.n_cols() == 0)
511 return cached_grad_grads[i][component];
512 else
513 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
514 s += inverse_node_matrix(i, j) * cached_grad_grads[j][component];
515
516 return s;
517}
518
519
520//---------------------------------------------------------------------------
521// Fill data of FEValues
522//---------------------------------------------------------------------------
523
524template <int dim, int spacedim>
525void
526FE_PolyTensor<dim, spacedim>::fill_fe_values(
527 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
528 const CellSimilarity::Similarity /*cell_similarity*/,
529 const Quadrature<dim> &quadrature,
530 const Mapping<dim, spacedim> &mapping,
531 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
532 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
533 &mapping_data,
534 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
535 ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
536 spacedim>
537 &output_data) const
538{
539 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
540 ExcInternalError());
541
542 compute_fill(cell,
543 QProjector<dim>::DataSetDescriptor::cell(),
544 quadrature.size(),
545 mapping,
546 mapping_internal,
547 mapping_data,
548 static_cast<const InternalData &>(fe_internal),
549 output_data);
550}
551
552
553
554template <int dim, int spacedim>
555void
556FE_PolyTensor<dim, spacedim>::fill_fe_face_values(
557 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
558 const unsigned int face_no,
559 const hp::QCollection<dim - 1> &quadrature,
560 const Mapping<dim, spacedim> &mapping,
561 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
562 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
563 &mapping_data,
564 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
565 ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
566 spacedim>
567 &output_data) const
568{
569 AssertDimension(quadrature.size(), 1);
570 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
571 ExcInternalError());
572
573 compute_fill(cell,
574 QProjector<dim>::DataSetDescriptor::face(
575 this->reference_cell(),
576 face_no,
577 // TODO: fix the custom implementation of orientation
578 dim == 2 ? numbers::default_geometric_orientation :
579 cell->combined_face_orientation(face_no),
580 quadrature[0].size()),
581 quadrature[0].size(),
582 mapping,
583 mapping_internal,
584 mapping_data,
585 static_cast<const InternalData &>(fe_internal),
586 output_data);
587}
588
589
590
591template <int dim, int spacedim>
592void
593FE_PolyTensor<dim, spacedim>::fill_fe_subface_values(
594 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
595 const unsigned int face_no,
596 const unsigned int sub_no,
597 const Quadrature<dim - 1> &quadrature,
598 const Mapping<dim, spacedim> &mapping,
599 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
600 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
601 &mapping_data,
602 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
603 ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
604 spacedim>
605 &output_data) const
606{
607 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
608 ExcInternalError());
609
610 compute_fill(cell,
611 QProjector<dim>::DataSetDescriptor::subface(
612 this->reference_cell(),
613 face_no,
614 sub_no,
615 // TODO: fix the custom implementation of orientation
616 dim == 2 ? numbers::default_geometric_orientation :
617 cell->combined_face_orientation(face_no),
618 quadrature.size(),
619 cell->subface_case(face_no)),
620 quadrature.size(),
621 mapping,
622 mapping_internal,
623 mapping_data,
624 static_cast<const InternalData &>(fe_internal),
625 output_data);
626}
627
628
629
630template <int dim, int spacedim>
631void
632FE_PolyTensor<dim, spacedim>::compute_fill(
633 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
634 const typename QProjector<dim>::DataSetDescriptor &offset,
635 const unsigned int n_q_points,
636 const Mapping<dim, spacedim> &mapping,
637 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
638 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
639 &mapping_data,
640 const typename FE_PolyTensor<dim, spacedim>::InternalData &fe_data,
641 internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
642 &output_data) const
643{
644 Assert(!(fe_data.update_each & update_values) ||
645 fe_data.shape_values.size()[0] == this->n_dofs_per_cell(),
646 ExcDimensionMismatch(fe_data.shape_values.size()[0],
647 this->n_dofs_per_cell()));
648
649 // TODO: The dof_sign_change only affects Nedelec elements and is not the
650 // correct thing on complicated meshes for higher order Nedelec elements.
651 // Something similar to FE_Q should be done to permute dofs and to change the
652 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
653 // class) is preferable.
654 std::fill(fe_data.dof_sign_change.begin(),
655 fe_data.dof_sign_change.end(),
656 1.0);
657 if (fe_data.update_each & update_values)
658 internal::FE_PolyTensor::get_dof_sign_change_nedelec(
659 cell, *this, this->mapping_kind, fe_data.dof_sign_change);
660
661 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
662 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
663 // Also nothing in 3d since we take care of it by using the
664 // adjust_quad_dof_sign_for_face_orientation_table.
665 if (fe_data.update_each & update_values)
666 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
667 *this,
668 this->mapping_kind,
669 fe_data.dof_sign_change);
670
671 // What is the first dof_index on a quad?
672 const unsigned int first_quad_index = this->get_first_quad_index();
673 // How many dofs per quad and how many quad dofs do we have at all?
674 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
675 const unsigned int n_quad_dofs =
676 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
677
678 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
679 ++dof_index)
680 {
681 /*
682 * This assumes that the dofs are ordered by first vertices, lines, quads
683 * and volume dofs. Note that in 2d this always gives false.
684 */
685 const bool is_quad_dof =
686 (dim == 2 ? false :
687 (first_quad_index <= dof_index) &&
688 (dof_index < first_quad_index + n_quad_dofs));
689
690 // TODO: This hack is not pretty and it is only here to handle the 2d
691 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
692 // handled by the
693 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
694 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
695 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
696 // it is handles with a table. This array is allocated in
697 // fe_poly_tensor.h.
698 double dof_sign = 1.0;
699 // under some circumstances fe_data.dof_sign_change is not allocated
700 if (fe_data.update_each & update_values)
701 dof_sign = fe_data.dof_sign_change[dof_index];
702
703 if (is_quad_dof)
704 {
705 /*
706 * Find the face belonging to this dof_index. This is integer
707 * division.
708 */
709 unsigned int face_index_from_dof_index =
710 (dof_index - first_quad_index) / (n_dofs_per_quad);
711
712 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
713
714 // Correct the dof_sign if necessary
715 if (adjust_quad_dof_sign_for_face_orientation(
716 local_quad_dof_index,
717 face_index_from_dof_index,
718 cell->combined_face_orientation(face_index_from_dof_index)))
719 dof_sign = -1.0;
720 }
721
722 const MappingKind mapping_kind = get_mapping_kind(dof_index);
723
724 const unsigned int first =
725 output_data.shape_function_to_row_table
726 [dof_index * this->n_components() +
727 this->get_nonzero_components(dof_index).first_selected_component()];
728
729 if (fe_data.update_each & update_values)
730 {
731 switch (mapping_kind)
732 {
733 case mapping_none:
734 {
735 for (unsigned int k = 0; k < n_q_points; ++k)
736 for (unsigned int d = 0; d < dim; ++d)
737 output_data.shape_values(first + d, k) =
738 fe_data.shape_values[dof_index][k + offset][d];
739 break;
740 }
741
742 case mapping_covariant:
743 case mapping_contravariant:
744 {
745 const ArrayView<Tensor<1, spacedim>>
746 transformed_shape_values =
747 make_array_view(fe_data.transformed_shape_values,
748 offset,
749 n_q_points);
750 mapping.transform(make_array_view(fe_data.shape_values,
751 dof_index,
752 offset,
753 n_q_points),
754 mapping_kind,
755 mapping_internal,
756 transformed_shape_values);
757
758 for (unsigned int k = 0; k < n_q_points; ++k)
759 for (unsigned int d = 0; d < spacedim; ++d)
760 output_data.shape_values(first + d, k) =
761 transformed_shape_values[k][d];
762
763 break;
764 }
765 case mapping_raviart_thomas:
766 case mapping_piola:
767 {
768 const ArrayView<Tensor<1, spacedim>>
769 transformed_shape_values =
770 make_array_view(fe_data.transformed_shape_values,
771 offset,
772 n_q_points);
773 mapping.transform(make_array_view(fe_data.shape_values,
774 dof_index,
775 offset,
776 n_q_points),
777 mapping_piola,
778 mapping_internal,
779 transformed_shape_values);
780 for (unsigned int k = 0; k < n_q_points; ++k)
781 for (unsigned int d = 0; d < spacedim; ++d)
782 output_data.shape_values(first + d, k) =
783 dof_sign * transformed_shape_values[k][d];
784 break;
785 }
786
787 case mapping_nedelec:
788 {
789 const ArrayView<Tensor<1, spacedim>>
790 transformed_shape_values =
791 make_array_view(fe_data.transformed_shape_values,
792 offset,
793 n_q_points);
794 mapping.transform(make_array_view(fe_data.shape_values,
795 dof_index,
796 offset,
797 n_q_points),
798 mapping_covariant,
799 mapping_internal,
800 transformed_shape_values);
801
802 for (unsigned int k = 0; k < n_q_points; ++k)
803 for (unsigned int d = 0; d < spacedim; ++d)
804 output_data.shape_values(first + d, k) =
805 dof_sign * transformed_shape_values[k][d];
806
807 break;
808 }
809
810 default:
811 DEAL_II_NOT_IMPLEMENTED();
812 }
813 }
814
815 if (fe_data.update_each & update_gradients)
816 {
817 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
818 make_array_view(fe_data.transformed_shape_grads,
819 offset,
820 n_q_points);
821 switch (mapping_kind)
822 {
823 case mapping_none:
824 {
825 mapping.transform(make_array_view(fe_data.shape_grads,
826 dof_index,
827 offset,
828 n_q_points),
829 mapping_covariant,
830 mapping_internal,
831 transformed_shape_grads);
832 for (unsigned int k = 0; k < n_q_points; ++k)
833 for (unsigned int d = 0; d < spacedim; ++d)
834 output_data.shape_gradients[first + d][k] =
835 transformed_shape_grads[k][d];
836 break;
837 }
838
839 case mapping_covariant:
840 {
841 mapping.transform(make_array_view(fe_data.shape_grads,
842 dof_index,
843 offset,
844 n_q_points),
845 mapping_covariant_gradient,
846 mapping_internal,
847 transformed_shape_grads);
848
849 for (unsigned int k = 0; k < n_q_points; ++k)
850 for (unsigned int d = 0; d < spacedim; ++d)
851 for (unsigned int n = 0; n < spacedim; ++n)
852 transformed_shape_grads[k][d] -=
853 output_data.shape_values(first + n, k) *
854 mapping_data.jacobian_pushed_forward_grads[k][n][d];
855
856 for (unsigned int k = 0; k < n_q_points; ++k)
857 for (unsigned int d = 0; d < spacedim; ++d)
858 output_data.shape_gradients[first + d][k] =
859 transformed_shape_grads[k][d];
860 break;
861 }
862
863 case mapping_contravariant:
864 {
865 for (unsigned int k = 0; k < n_q_points; ++k)
866 if constexpr (dim == spacedim)
867 fe_data.untransformed_shape_grads[k + offset] =
868 fe_data.shape_grads[dof_index][k + offset];
869 else
870 for (unsigned int d = 0; d < dim; ++d)
871 fe_data.untransformed_shape_grads[k + offset][d] =
872 fe_data.shape_grads[dof_index][k + offset][d];
873 mapping.transform(
874 make_array_view(fe_data.untransformed_shape_grads,
875 offset,
876 n_q_points),
877 mapping_contravariant_gradient,
878 mapping_internal,
879 transformed_shape_grads);
880
881 for (unsigned int k = 0; k < n_q_points; ++k)
882 for (unsigned int d = 0; d < spacedim; ++d)
883 for (unsigned int n = 0; n < spacedim; ++n)
884 transformed_shape_grads[k][d] +=
885 output_data.shape_values(first + n, k) *
886 mapping_data.jacobian_pushed_forward_grads[k][d][n];
887
888 for (unsigned int k = 0; k < n_q_points; ++k)
889 for (unsigned int d = 0; d < spacedim; ++d)
890 output_data.shape_gradients[first + d][k] =
891 transformed_shape_grads[k][d];
892
893 break;
894 }
895
896 case mapping_raviart_thomas:
897 case mapping_piola:
898 {
899 for (unsigned int k = 0; k < n_q_points; ++k)
900 if constexpr (dim == spacedim)
901 fe_data.untransformed_shape_grads[k + offset] =
902 fe_data.shape_grads[dof_index][k + offset];
903 else
904 for (unsigned int d = 0; d < dim; ++d)
905 fe_data.untransformed_shape_grads[k + offset][d] =
906 fe_data.shape_grads[dof_index][k + offset][d];
907 mapping.transform(
908 make_array_view(fe_data.untransformed_shape_grads,
909 offset,
910 n_q_points),
911 mapping_piola_gradient,
912 mapping_internal,
913 transformed_shape_grads);
914
915 for (unsigned int k = 0; k < n_q_points; ++k)
916 for (unsigned int d = 0; d < spacedim; ++d)
917 for (unsigned int n = 0; n < spacedim; ++n)
918 transformed_shape_grads[k][d] +=
919 (output_data.shape_values(first + n, k) *
920 mapping_data
921 .jacobian_pushed_forward_grads[k][d][n]) -
922 (output_data.shape_values(first + d, k) *
923 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
924
925 for (unsigned int k = 0; k < n_q_points; ++k)
926 for (unsigned int d = 0; d < spacedim; ++d)
927 output_data.shape_gradients[first + d][k] =
928 dof_sign * transformed_shape_grads[k][d];
929
930 break;
931 }
932
933 case mapping_nedelec:
934 {
935 // treat the gradients of
936 // this particular shape
937 // function at all
938 // q-points. if Dv is the
939 // gradient of the shape
940 // function on the unit
941 // cell, then
942 // (J^-T)Dv(J^-1) is the
943 // value we want to have on
944 // the real cell.
945 for (unsigned int k = 0; k < n_q_points; ++k)
946 if constexpr (dim == spacedim)
947 fe_data.untransformed_shape_grads[k + offset] =
948 fe_data.shape_grads[dof_index][k + offset];
949 else
950 for (unsigned int d = 0; d < dim; ++d)
951 fe_data.untransformed_shape_grads[k + offset][d] =
952 fe_data.shape_grads[dof_index][k + offset][d];
953
954 mapping.transform(
955 make_array_view(fe_data.untransformed_shape_grads,
956 offset,
957 n_q_points),
958 mapping_covariant_gradient,
959 mapping_internal,
960 transformed_shape_grads);
961
962 for (unsigned int k = 0; k < n_q_points; ++k)
963 for (unsigned int d = 0; d < spacedim; ++d)
964 for (unsigned int n = 0; n < spacedim; ++n)
965 transformed_shape_grads[k][d] -=
966 output_data.shape_values(first + n, k) *
967 mapping_data.jacobian_pushed_forward_grads[k][n][d];
968
969 for (unsigned int k = 0; k < n_q_points; ++k)
970 for (unsigned int d = 0; d < spacedim; ++d)
971 output_data.shape_gradients[first + d][k] =
972 dof_sign * transformed_shape_grads[k][d];
973
974 break;
975 }
976
977 default:
978 DEAL_II_NOT_IMPLEMENTED();
979 }
980 }
981
982 if (fe_data.update_each & update_hessians)
983 {
984 switch (mapping_kind)
985 {
986 case mapping_none:
987 {
988 const ArrayView<Tensor<3, spacedim>>
989 transformed_shape_hessians =
990 make_array_view(fe_data.transformed_shape_hessians,
991 offset,
992 n_q_points);
993 mapping.transform(make_array_view(fe_data.shape_grad_grads,
994 dof_index,
995 offset,
996 n_q_points),
997 mapping_covariant_gradient,
998 mapping_internal,
999 transformed_shape_hessians);
1000
1001 for (unsigned int k = 0; k < n_q_points; ++k)
1002 for (unsigned int d = 0; d < spacedim; ++d)
1003 for (unsigned int n = 0; n < spacedim; ++n)
1004 transformed_shape_hessians[k][d] -=
1005 output_data.shape_gradients[first + d][k][n] *
1006 mapping_data.jacobian_pushed_forward_grads[k][n];
1007
1008 for (unsigned int k = 0; k < n_q_points; ++k)
1009 for (unsigned int d = 0; d < spacedim; ++d)
1010 output_data.shape_hessians[first + d][k] =
1011 transformed_shape_hessians[k][d];
1012
1013 break;
1014 }
1015 case mapping_covariant:
1016 {
1017 for (unsigned int k = 0; k < n_q_points; ++k)
1018 if constexpr (dim == spacedim)
1019 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1020 fe_data.shape_grad_grads[dof_index][k + offset];
1021 else
1022 for (unsigned int d = 0; d < dim; ++d)
1023 for (unsigned int e = 0; e < dim; ++e)
1024 fe_data
1025 .untransformed_shape_hessian_tensors[k + offset][d]
1026 [e] =
1027 fe_data
1028 .shape_grad_grads[dof_index][k + offset][d][e];
1029
1030 const ArrayView<Tensor<3, spacedim>>
1031 transformed_shape_hessians =
1032 make_array_view(fe_data.transformed_shape_hessians,
1033 offset,
1034 n_q_points);
1035 mapping.transform(
1036 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1037 offset,
1038 n_q_points),
1039 mapping_covariant_hessian,
1040 mapping_internal,
1041 transformed_shape_hessians);
1042
1043 for (unsigned int k = 0; k < n_q_points; ++k)
1044 for (unsigned int d = 0; d < spacedim; ++d)
1045 for (unsigned int n = 0; n < spacedim; ++n)
1046 for (unsigned int i = 0; i < spacedim; ++i)
1047 for (unsigned int j = 0; j < spacedim; ++j)
1048 {
1049 transformed_shape_hessians[k][d][i][j] -=
1050 (output_data.shape_values(first + n, k) *
1051 mapping_data
1052 .jacobian_pushed_forward_2nd_derivatives
1053 [k][n][d][i][j]) +
1054 (output_data.shape_gradients[first + d][k][n] *
1055 mapping_data
1056 .jacobian_pushed_forward_grads[k][n][i][j]) +
1057 (output_data.shape_gradients[first + n][k][i] *
1058 mapping_data
1059 .jacobian_pushed_forward_grads[k][n][d][j]) +
1060 (output_data.shape_gradients[first + n][k][j] *
1061 mapping_data
1062 .jacobian_pushed_forward_grads[k][n][i][d]);
1063 }
1064
1065 for (unsigned int k = 0; k < n_q_points; ++k)
1066 for (unsigned int d = 0; d < spacedim; ++d)
1067 output_data.shape_hessians[first + d][k] =
1068 transformed_shape_hessians[k][d];
1069
1070 break;
1071 }
1072
1073 case mapping_contravariant:
1074 {
1075 for (unsigned int k = 0; k < n_q_points; ++k)
1076 if constexpr (dim == spacedim)
1077 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1078 fe_data.shape_grad_grads[dof_index][k + offset];
1079 else
1080 for (unsigned int d = 0; d < dim; ++d)
1081 for (unsigned int e = 0; e < dim; ++e)
1082 fe_data
1083 .untransformed_shape_hessian_tensors[k + offset][d]
1084 [e] =
1085 fe_data
1086 .shape_grad_grads[dof_index][k + offset][d][e];
1087
1088 const ArrayView<Tensor<3, spacedim>>
1089 transformed_shape_hessians =
1090 make_array_view(fe_data.transformed_shape_hessians,
1091 offset,
1092 n_q_points);
1093 mapping.transform(
1094 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1095 offset,
1096 n_q_points),
1097 mapping_contravariant_hessian,
1098 mapping_internal,
1099 transformed_shape_hessians);
1100
1101 for (unsigned int k = 0; k < n_q_points; ++k)
1102 for (unsigned int d = 0; d < spacedim; ++d)
1103 for (unsigned int n = 0; n < spacedim; ++n)
1104 for (unsigned int i = 0; i < spacedim; ++i)
1105 for (unsigned int j = 0; j < spacedim; ++j)
1106 {
1107 transformed_shape_hessians[k][d][i][j] +=
1108 (output_data.shape_values(first + n, k) *
1109 mapping_data
1110 .jacobian_pushed_forward_2nd_derivatives
1111 [k][d][n][i][j]) +
1112 (output_data.shape_gradients[first + n][k][i] *
1113 mapping_data
1114 .jacobian_pushed_forward_grads[k][d][n][j]) +
1115 (output_data.shape_gradients[first + n][k][j] *
1116 mapping_data
1117 .jacobian_pushed_forward_grads[k][d][i][n]) -
1118 (output_data.shape_gradients[first + d][k][n] *
1119 mapping_data
1120 .jacobian_pushed_forward_grads[k][n][i][j]);
1121 for (unsigned int m = 0; m < spacedim; ++m)
1122 transformed_shape_hessians[k][d][i][j] -=
1123 (mapping_data
1124 .jacobian_pushed_forward_grads[k][d][i]
1125 [m] *
1126 mapping_data
1127 .jacobian_pushed_forward_grads[k][m][n]
1128 [j] *
1129 output_data.shape_values(first + n, k)) +
1130 (mapping_data
1131 .jacobian_pushed_forward_grads[k][d][m]
1132 [j] *
1133 mapping_data
1134 .jacobian_pushed_forward_grads[k][m][i]
1135 [n] *
1136 output_data.shape_values(first + n, k));
1137 }
1138
1139 for (unsigned int k = 0; k < n_q_points; ++k)
1140 for (unsigned int d = 0; d < spacedim; ++d)
1141 output_data.shape_hessians[first + d][k] =
1142 transformed_shape_hessians[k][d];
1143
1144 break;
1145 }
1146
1147 case mapping_raviart_thomas:
1148 case mapping_piola:
1149 {
1150 for (unsigned int k = 0; k < n_q_points; ++k)
1151 if constexpr (dim == spacedim)
1152 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1153 fe_data.shape_grad_grads[dof_index][k + offset];
1154 else
1155 for (unsigned int d = 0; d < dim; ++d)
1156 for (unsigned int e = 0; e < dim; ++e)
1157 fe_data
1158 .untransformed_shape_hessian_tensors[k + offset][d]
1159 [e] =
1160 fe_data
1161 .shape_grad_grads[dof_index][k + offset][d][e];
1162
1163 const ArrayView<Tensor<3, spacedim>>
1164 transformed_shape_hessians =
1165 make_array_view(fe_data.transformed_shape_hessians,
1166 offset,
1167 n_q_points);
1168 mapping.transform(
1169 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1170 offset,
1171 n_q_points),
1172 mapping_piola_hessian,
1173 mapping_internal,
1174 transformed_shape_hessians);
1175
1176 for (unsigned int k = 0; k < n_q_points; ++k)
1177 for (unsigned int d = 0; d < spacedim; ++d)
1178 for (unsigned int n = 0; n < spacedim; ++n)
1179 for (unsigned int i = 0; i < spacedim; ++i)
1180 for (unsigned int j = 0; j < spacedim; ++j)
1181 {
1182 transformed_shape_hessians[k][d][i][j] +=
1183 (output_data.shape_values(first + n, k) *
1184 mapping_data
1185 .jacobian_pushed_forward_2nd_derivatives
1186 [k][d][n][i][j]) +
1187 (output_data.shape_gradients[first + n][k][i] *
1188 mapping_data
1189 .jacobian_pushed_forward_grads[k][d][n][j]) +
1190 (output_data.shape_gradients[first + n][k][j] *
1191 mapping_data
1192 .jacobian_pushed_forward_grads[k][d][i][n]) -
1193 (output_data.shape_gradients[first + d][k][n] *
1194 mapping_data
1195 .jacobian_pushed_forward_grads[k][n][i][j]);
1196
1197 transformed_shape_hessians[k][d][i][j] -=
1198 (output_data.shape_values(first + d, k) *
1199 mapping_data
1200 .jacobian_pushed_forward_2nd_derivatives
1201 [k][n][n][i][j]) +
1202 (output_data.shape_gradients[first + d][k][i] *
1203 mapping_data
1204 .jacobian_pushed_forward_grads[k][n][n][j]) +
1205 (output_data.shape_gradients[first + d][k][j] *
1206 mapping_data
1207 .jacobian_pushed_forward_grads[k][n][n][i]);
1208
1209 for (unsigned int m = 0; m < spacedim; ++m)
1210 {
1211 transformed_shape_hessians[k][d][i][j] -=
1212 (mapping_data
1213 .jacobian_pushed_forward_grads[k][d][i]
1214 [m] *
1215 mapping_data
1216 .jacobian_pushed_forward_grads[k][m][n]
1217 [j] *
1218 output_data.shape_values(first + n, k)) +
1219 (mapping_data
1220 .jacobian_pushed_forward_grads[k][d][m]
1221 [j] *
1222 mapping_data
1223 .jacobian_pushed_forward_grads[k][m][i]
1224 [n] *
1225 output_data.shape_values(first + n, k));
1226
1227 transformed_shape_hessians[k][d][i][j] +=
1228 (mapping_data
1229 .jacobian_pushed_forward_grads[k][n][i]
1230 [m] *
1231 mapping_data
1232 .jacobian_pushed_forward_grads[k][m][n]
1233 [j] *
1234 output_data.shape_values(first + d, k)) +
1235 (mapping_data
1236 .jacobian_pushed_forward_grads[k][n][m]
1237 [j] *
1238 mapping_data
1239 .jacobian_pushed_forward_grads[k][m][i]
1240 [n] *
1241 output_data.shape_values(first + d, k));
1242 }
1243 }
1244
1245 for (unsigned int k = 0; k < n_q_points; ++k)
1246 for (unsigned int d = 0; d < spacedim; ++d)
1247 output_data.shape_hessians[first + d][k] =
1248 dof_sign * transformed_shape_hessians[k][d];
1249
1250 break;
1251 }
1252
1253 case mapping_nedelec:
1254 {
1255 for (unsigned int k = 0; k < n_q_points; ++k)
1256 if constexpr (dim == spacedim)
1257 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1258 fe_data.shape_grad_grads[dof_index][k + offset];
1259 else
1260 for (unsigned int d = 0; d < dim; ++d)
1261 for (unsigned int e = 0; e < dim; ++e)
1262 fe_data
1263 .untransformed_shape_hessian_tensors[k + offset][d]
1264 [e] =
1265 fe_data
1266 .shape_grad_grads[dof_index][k + offset][d][e];
1267
1268 const ArrayView<Tensor<3, spacedim>>
1269 transformed_shape_hessians =
1270 make_array_view(fe_data.transformed_shape_hessians,
1271 offset,
1272 n_q_points);
1273 mapping.transform(
1274 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1275 offset,
1276 n_q_points),
1277 mapping_covariant_hessian,
1278 mapping_internal,
1279 transformed_shape_hessians);
1280
1281 for (unsigned int k = 0; k < n_q_points; ++k)
1282 for (unsigned int d = 0; d < spacedim; ++d)
1283 for (unsigned int n = 0; n < spacedim; ++n)
1284 for (unsigned int i = 0; i < spacedim; ++i)
1285 for (unsigned int j = 0; j < spacedim; ++j)
1286 {
1287 transformed_shape_hessians[k][d][i][j] -=
1288 (output_data.shape_values(first + n, k) *
1289 mapping_data
1290 .jacobian_pushed_forward_2nd_derivatives
1291 [k][n][d][i][j]) +
1292 (output_data.shape_gradients[first + d][k][n] *
1293 mapping_data
1294 .jacobian_pushed_forward_grads[k][n][i][j]) +
1295 (output_data.shape_gradients[first + n][k][i] *
1296 mapping_data
1297 .jacobian_pushed_forward_grads[k][n][d][j]) +
1298 (output_data.shape_gradients[first + n][k][j] *
1299 mapping_data
1300 .jacobian_pushed_forward_grads[k][n][i][d]);
1301 }
1302
1303 for (unsigned int k = 0; k < n_q_points; ++k)
1304 for (unsigned int d = 0; d < spacedim; ++d)
1305 output_data.shape_hessians[first + d][k] =
1306 dof_sign * transformed_shape_hessians[k][d];
1307
1308 break;
1309 }
1310
1311 default:
1312 DEAL_II_NOT_IMPLEMENTED();
1313 }
1314 }
1315
1316 // third derivatives are not implemented
1317 if (fe_data.update_each & update_3rd_derivatives)
1318 {
1319 DEAL_II_NOT_IMPLEMENTED();
1320 }
1321 }
1322}
1323
1324
1325
1326template <int dim, int spacedim>
1327UpdateFlags
1328FE_PolyTensor<dim, spacedim>::requires_update_flags(
1329 const UpdateFlags flags) const
1330{
1331 UpdateFlags out = update_default;
1332
1333 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1334 {
1335 const MappingKind mapping_kind = get_mapping_kind(i);
1336
1337 switch (mapping_kind)
1338 {
1339 case mapping_none:
1340 {
1341 if (flags & update_values)
1342 out |= update_values;
1343
1344 if (flags & update_gradients)
1345 out |= update_gradients | update_values |
1346 update_jacobian_pushed_forward_grads;
1347
1348 if (flags & update_hessians)
1349 out |= update_hessians | update_values | update_gradients |
1350 update_jacobian_pushed_forward_grads |
1351 update_jacobian_pushed_forward_2nd_derivatives;
1352 break;
1353 }
1354 case mapping_raviart_thomas:
1355 case mapping_piola:
1356 {
1357 if (flags & update_values)
1358 out |= update_values | update_piola;
1359
1360 if (flags & update_gradients)
1361 out |= update_gradients | update_values | update_piola |
1362 update_jacobian_pushed_forward_grads |
1363 update_covariant_transformation |
1364 update_contravariant_transformation;
1365
1366 if (flags & update_hessians)
1367 out |= update_hessians | update_piola | update_values |
1368 update_gradients | update_jacobian_pushed_forward_grads |
1369 update_jacobian_pushed_forward_2nd_derivatives |
1370 update_covariant_transformation;
1371
1372 break;
1373 }
1374
1375
1376 case mapping_contravariant:
1377 {
1378 if (flags & update_values)
1379 out |= update_values | update_piola;
1380
1381 if (flags & update_gradients)
1382 out |= update_gradients | update_values |
1383 update_jacobian_pushed_forward_grads |
1384 update_covariant_transformation |
1385 update_contravariant_transformation;
1386
1387 if (flags & update_hessians)
1388 out |= update_hessians | update_piola | update_values |
1389 update_gradients | update_jacobian_pushed_forward_grads |
1390 update_jacobian_pushed_forward_2nd_derivatives |
1391 update_covariant_transformation;
1392
1393 break;
1394 }
1395
1396 case mapping_nedelec:
1397 case mapping_covariant:
1398 {
1399 if (flags & update_values)
1400 out |= update_values | update_covariant_transformation;
1401
1402 if (flags & update_gradients)
1403 out |= update_gradients | update_values |
1404 update_jacobian_pushed_forward_grads |
1405 update_covariant_transformation;
1406
1407 if (flags & update_hessians)
1408 out |= update_hessians | update_values | update_gradients |
1409 update_jacobian_pushed_forward_grads |
1410 update_jacobian_pushed_forward_2nd_derivatives |
1411 update_covariant_transformation;
1412
1413 break;
1414 }
1415
1416 default:
1417 {
1418 DEAL_II_NOT_IMPLEMENTED();
1419 }
1420 }
1421 }
1422
1423 return out;
1424}
1425
1426#endif
1427// explicit instantiations
1428#include "fe/fe_poly_tensor.inst"
1429
1430
1431DEAL_II_NAMESPACE_CLOSE
array_view.h
make_array_view
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:1002
cell_id.h
FE_PolyTensor::InternalData
Definition fe_poly_tensor.h:514
FE_PolyTensor::InternalData::transformed_shape_hessians
std::vector< Tensor< 3, spacedim > > transformed_shape_hessians
Definition fe_poly_tensor.h:545
FE_PolyTensor::InternalData::shape_grad_grads
Table< 2, DerivativeForm< 2, dim, spacedim > > shape_grad_grads
Definition fe_poly_tensor.h:534
FE_PolyTensor::InternalData::shape_grads
Table< 2, DerivativeForm< 1, dim, spacedim > > shape_grads
Definition fe_poly_tensor.h:527
FE_PolyTensor::InternalData::transformed_shape_grads
std::vector< Tensor< 2, spacedim > > transformed_shape_grads
Definition fe_poly_tensor.h:542
FE_PolyTensor::InternalData::untransformed_shape_hessian_tensors
std::vector< Tensor< 3, dim > > untransformed_shape_hessian_tensors
Definition fe_poly_tensor.h:546
FE_PolyTensor::InternalData::untransformed_shape_grads
std::vector< Tensor< 2, dim > > untransformed_shape_grads
Definition fe_poly_tensor.h:543
FE_PolyTensor::InternalData::shape_values
Table< 2, Tensor< 1, dim > > shape_values
Definition fe_poly_tensor.h:520
FE_PolyTensor::InternalData::transformed_shape_values
std::vector< Tensor< 1, spacedim > > transformed_shape_values
Definition fe_poly_tensor.h:540
FE_PolyTensor::InternalData::dof_sign_change
std::vector< double > dof_sign_change
Definition fe_poly_tensor.h:539
FE_PolyTensor
Definition fe_poly_tensor.h:140
FE_PolyTensor::requires_update_flags
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
FE_PolyTensor::adjust_quad_dof_sign_for_face_orientation
bool adjust_quad_dof_sign_for_face_orientation(const unsigned int index, const unsigned int face_no, const types::geometric_orientation combined_orientation) const
FE_PolyTensor::shape_grad
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
FE_PolyTensor::shape_grad_grad_component
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
FE_PolyTensor::shape_value_component
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
FE_PolyTensor::shape_grad_grad
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
FE_PolyTensor::compute_fill
void compute_fill(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename QProjector< dim >::DataSetDescriptor &offset, const unsigned int n_q_points, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FE_PolyTensor< dim, spacedim >::InternalData &fe_internal, internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
FE_PolyTensor::fill_fe_face_values
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
FE_PolyTensor::shape_value
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
FE_PolyTensor::fill_fe_values
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
FE_PolyTensor::fill_fe_subface_values
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
FE_PolyTensor::get_mapping_kind
MappingKind get_mapping_kind(const unsigned int i) const
FE_PolyTensor::FE_PolyTensor
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)
FE_PolyTensor::single_mapping_kind
bool single_mapping_kind() const
FE_PolyTensor::shape_grad_component
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
FiniteElementData::n_dofs_per_face
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
FiniteElementData::tensor_degree
unsigned int tensor_degree() const
FiniteElement::face_to_cell_index
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
Mapping
Abstract base class for mapping classes.
Definition mapping.h:318
Mapping::transform
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
Point
Definition point.h:111
QProjector::DataSetDescriptor
Class storing the offset index into a Quadrature rule created by project_to_all_faces() or project_to...
Definition qprojector.h:204
QProjector
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:68
Quadrature::size
unsigned int size() const
hp::Collection::size
unsigned int size() const
Definition collection.h:314
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Definition mapping_related_data.h:113
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_grads
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Definition mapping_related_data.h:101
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
derivative_form.h
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exception_macros.h:705
fe_values.h
fe_poly_tensor.h
first
Point< 2 > first
Definition grid_out.cc:4639
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:62
update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_2nd_derivatives
Definition fe_update_flags.h:201
update_contravariant_transformation
@ update_contravariant_transformation
Contravariant transformation.
Definition fe_update_flags.h:170
update_jacobian_pushed_forward_grads
@ update_jacobian_pushed_forward_grads
Definition fe_update_flags.h:192
update_hessians
@ update_hessians
Second derivatives of shape functions.
Definition fe_update_flags.h:85
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:73
update_3rd_derivatives
@ update_3rd_derivatives
Third derivatives of shape functions.
Definition fe_update_flags.h:91
update_covariant_transformation
@ update_covariant_transformation
Covariant transformation.
Definition fe_update_flags.h:163
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:79
update_default
@ update_default
No update.
Definition fe_update_flags.h:64
update_piola
@ update_piola
Values needed for Piola transform.
Definition fe_update_flags.h:219
MappingKind
MappingKind
Definition mapping.h:79
mapping_piola
@ mapping_piola
Definition mapping.h:114
mapping_covariant_gradient
@ mapping_covariant_gradient
Definition mapping.h:100
mapping_covariant
@ mapping_covariant
Definition mapping.h:89
mapping_nedelec
@ mapping_nedelec
Definition mapping.h:129
mapping_contravariant
@ mapping_contravariant
Definition mapping.h:94
mapping_raviart_thomas
@ mapping_raviart_thomas
Definition mapping.h:134
mapping_contravariant_hessian
@ mapping_contravariant_hessian
Definition mapping.h:156
mapping_covariant_hessian
@ mapping_covariant_hessian
Definition mapping.h:150
mapping_contravariant_gradient
@ mapping_contravariant_gradient
Definition mapping.h:106
mapping_none
@ mapping_none
Definition mapping.h:84
mapping_piola_gradient
@ mapping_piola_gradient
Definition mapping.h:120
mapping_piola_hessian
@ mapping_piola_hessian
Definition mapping.h:162
mapping_cartesian.h
size
std::size_t size
Definition mpi.cc:733
GridGenerator::reference_cell
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell< dim > &reference_cell)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
PlasticityLab::component
*  *  *  ScaleZFunction< dim, Number, components >::ScaleZFunction *  component(component)
Definition ALE_Finite_Strain_Plasticity.h:8885
PlasticityLab::const
*  *  *  *  std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters   const
Definition ALE_Finite_Strain_Plasticity.h:9082
internal::combined_face_orientation
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
Definition tria_orientation.h:31
numbers::default_geometric_orientation
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:342
types::geometric_orientation
std::uint8_t geometric_orientation
Definition types.h:38
polynomials_abf.h
polynomials_bdm.h
polynomials_bernardi_raugel.h
polynomials_nedelec.h
polynomials_raviart_thomas.h
polynomials_rt_bubbles.h
tria_orientation.h
Generated by doxygen 1.9.8