deal.II version 9.8.0
\(\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 < dim; ++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 < dim; ++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 < dim; ++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 < dim; ++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 < dim; ++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 fe_data.untransformed_shape_grads[k + offset] =
867 fe_data.shape_grads[dof_index][k + offset];
868 mapping.transform(
869 make_array_view(fe_data.untransformed_shape_grads,
870 offset,
871 n_q_points),
872 mapping_contravariant_gradient,
873 mapping_internal,
874 transformed_shape_grads);
875
876 for (unsigned int k = 0; k < n_q_points; ++k)
877 for (unsigned int d = 0; d < spacedim; ++d)
878 for (unsigned int n = 0; n < spacedim; ++n)
879 transformed_shape_grads[k][d] +=
880 output_data.shape_values(first + n, k) *
881 mapping_data.jacobian_pushed_forward_grads[k][d][n];
882
883 for (unsigned int k = 0; k < n_q_points; ++k)
884 for (unsigned int d = 0; d < dim; ++d)
885 output_data.shape_gradients[first + d][k] =
886 transformed_shape_grads[k][d];
887
888 break;
889 }
890
891 case mapping_raviart_thomas:
892 case mapping_piola:
893 {
894 for (unsigned int k = 0; k < n_q_points; ++k)
895 fe_data.untransformed_shape_grads[k + offset] =
896 fe_data.shape_grads[dof_index][k + offset];
897 mapping.transform(
898 make_array_view(fe_data.untransformed_shape_grads,
899 offset,
900 n_q_points),
901 mapping_piola_gradient,
902 mapping_internal,
903 transformed_shape_grads);
904
905 for (unsigned int k = 0; k < n_q_points; ++k)
906 for (unsigned int d = 0; d < spacedim; ++d)
907 for (unsigned int n = 0; n < spacedim; ++n)
908 transformed_shape_grads[k][d] +=
909 (output_data.shape_values(first + n, k) *
910 mapping_data
911 .jacobian_pushed_forward_grads[k][d][n]) -
912 (output_data.shape_values(first + d, k) *
913 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
914
915 for (unsigned int k = 0; k < n_q_points; ++k)
916 for (unsigned int d = 0; d < dim; ++d)
917 output_data.shape_gradients[first + d][k] =
918 dof_sign * transformed_shape_grads[k][d];
919
920 break;
921 }
922
923 case mapping_nedelec:
924 {
925 // treat the gradients of
926 // this particular shape
927 // function at all
928 // q-points. if Dv is the
929 // gradient of the shape
930 // function on the unit
931 // cell, then
932 // (J^-T)Dv(J^-1) is the
933 // value we want to have on
934 // the real cell.
935 for (unsigned int k = 0; k < n_q_points; ++k)
936 fe_data.untransformed_shape_grads[k + offset] =
937 fe_data.shape_grads[dof_index][k + offset];
938
939 mapping.transform(
940 make_array_view(fe_data.untransformed_shape_grads,
941 offset,
942 n_q_points),
943 mapping_covariant_gradient,
944 mapping_internal,
945 transformed_shape_grads);
946
947 for (unsigned int k = 0; k < n_q_points; ++k)
948 for (unsigned int d = 0; d < spacedim; ++d)
949 for (unsigned int n = 0; n < spacedim; ++n)
950 transformed_shape_grads[k][d] -=
951 output_data.shape_values(first + n, k) *
952 mapping_data.jacobian_pushed_forward_grads[k][n][d];
953
954 for (unsigned int k = 0; k < n_q_points; ++k)
955 for (unsigned int d = 0; d < dim; ++d)
956 output_data.shape_gradients[first + d][k] =
957 dof_sign * transformed_shape_grads[k][d];
958
959 break;
960 }
961
962 default:
963 DEAL_II_NOT_IMPLEMENTED();
964 }
965 }
966
967 if (fe_data.update_each & update_hessians)
968 {
969 switch (mapping_kind)
970 {
971 case mapping_none:
972 {
973 const ArrayView<Tensor<3, spacedim>>
974 transformed_shape_hessians =
975 make_array_view(fe_data.transformed_shape_hessians,
976 offset,
977 n_q_points);
978 mapping.transform(make_array_view(fe_data.shape_grad_grads,
979 dof_index,
980 offset,
981 n_q_points),
982 mapping_covariant_gradient,
983 mapping_internal,
984 transformed_shape_hessians);
985
986 for (unsigned int k = 0; k < n_q_points; ++k)
987 for (unsigned int d = 0; d < spacedim; ++d)
988 for (unsigned int n = 0; n < spacedim; ++n)
989 transformed_shape_hessians[k][d] -=
990 output_data.shape_gradients[first + d][k][n] *
991 mapping_data.jacobian_pushed_forward_grads[k][n];
992
993 for (unsigned int k = 0; k < n_q_points; ++k)
994 for (unsigned int d = 0; d < dim; ++d)
995 output_data.shape_hessians[first + d][k] =
996 transformed_shape_hessians[k][d];
997
998 break;
999 }
1000 case mapping_covariant:
1001 {
1002 for (unsigned int k = 0; k < n_q_points; ++k)
1003 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1004 fe_data.shape_grad_grads[dof_index][k + offset];
1005
1006 const ArrayView<Tensor<3, spacedim>>
1007 transformed_shape_hessians =
1008 make_array_view(fe_data.transformed_shape_hessians,
1009 offset,
1010 n_q_points);
1011 mapping.transform(
1012 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1013 offset,
1014 n_q_points),
1015 mapping_covariant_hessian,
1016 mapping_internal,
1017 transformed_shape_hessians);
1018
1019 for (unsigned int k = 0; k < n_q_points; ++k)
1020 for (unsigned int d = 0; d < spacedim; ++d)
1021 for (unsigned int n = 0; n < spacedim; ++n)
1022 for (unsigned int i = 0; i < spacedim; ++i)
1023 for (unsigned int j = 0; j < spacedim; ++j)
1024 {
1025 transformed_shape_hessians[k][d][i][j] -=
1026 (output_data.shape_values(first + n, k) *
1027 mapping_data
1028 .jacobian_pushed_forward_2nd_derivatives
1029 [k][n][d][i][j]) +
1030 (output_data.shape_gradients[first + d][k][n] *
1031 mapping_data
1032 .jacobian_pushed_forward_grads[k][n][i][j]) +
1033 (output_data.shape_gradients[first + n][k][i] *
1034 mapping_data
1035 .jacobian_pushed_forward_grads[k][n][d][j]) +
1036 (output_data.shape_gradients[first + n][k][j] *
1037 mapping_data
1038 .jacobian_pushed_forward_grads[k][n][i][d]);
1039 }
1040
1041 for (unsigned int k = 0; k < n_q_points; ++k)
1042 for (unsigned int d = 0; d < dim; ++d)
1043 output_data.shape_hessians[first + d][k] =
1044 transformed_shape_hessians[k][d];
1045
1046 break;
1047 }
1048
1049 case mapping_contravariant:
1050 {
1051 for (unsigned int k = 0; k < n_q_points; ++k)
1052 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1053 fe_data.shape_grad_grads[dof_index][k + offset];
1054
1055 const ArrayView<Tensor<3, spacedim>>
1056 transformed_shape_hessians =
1057 make_array_view(fe_data.transformed_shape_hessians,
1058 offset,
1059 n_q_points);
1060 mapping.transform(
1061 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1062 offset,
1063 n_q_points),
1064 mapping_contravariant_hessian,
1065 mapping_internal,
1066 transformed_shape_hessians);
1067
1068 for (unsigned int k = 0; k < n_q_points; ++k)
1069 for (unsigned int d = 0; d < spacedim; ++d)
1070 for (unsigned int n = 0; n < spacedim; ++n)
1071 for (unsigned int i = 0; i < spacedim; ++i)
1072 for (unsigned int j = 0; j < spacedim; ++j)
1073 {
1074 transformed_shape_hessians[k][d][i][j] +=
1075 (output_data.shape_values(first + n, k) *
1076 mapping_data
1077 .jacobian_pushed_forward_2nd_derivatives
1078 [k][d][n][i][j]) +
1079 (output_data.shape_gradients[first + n][k][i] *
1080 mapping_data
1081 .jacobian_pushed_forward_grads[k][d][n][j]) +
1082 (output_data.shape_gradients[first + n][k][j] *
1083 mapping_data
1084 .jacobian_pushed_forward_grads[k][d][i][n]) -
1085 (output_data.shape_gradients[first + d][k][n] *
1086 mapping_data
1087 .jacobian_pushed_forward_grads[k][n][i][j]);
1088 for (unsigned int m = 0; m < spacedim; ++m)
1089 transformed_shape_hessians[k][d][i][j] -=
1090 (mapping_data
1091 .jacobian_pushed_forward_grads[k][d][i]
1092 [m] *
1093 mapping_data
1094 .jacobian_pushed_forward_grads[k][m][n]
1095 [j] *
1096 output_data.shape_values(first + n, k)) +
1097 (mapping_data
1098 .jacobian_pushed_forward_grads[k][d][m]
1099 [j] *
1100 mapping_data
1101 .jacobian_pushed_forward_grads[k][m][i]
1102 [n] *
1103 output_data.shape_values(first + n, k));
1104 }
1105
1106 for (unsigned int k = 0; k < n_q_points; ++k)
1107 for (unsigned int d = 0; d < dim; ++d)
1108 output_data.shape_hessians[first + d][k] =
1109 transformed_shape_hessians[k][d];
1110
1111 break;
1112 }
1113
1114 case mapping_raviart_thomas:
1115 case mapping_piola:
1116 {
1117 for (unsigned int k = 0; k < n_q_points; ++k)
1118 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1119 fe_data.shape_grad_grads[dof_index][k + offset];
1120
1121 const ArrayView<Tensor<3, spacedim>>
1122 transformed_shape_hessians =
1123 make_array_view(fe_data.transformed_shape_hessians,
1124 offset,
1125 n_q_points);
1126 mapping.transform(
1127 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1128 offset,
1129 n_q_points),
1130 mapping_piola_hessian,
1131 mapping_internal,
1132 transformed_shape_hessians);
1133
1134 for (unsigned int k = 0; k < n_q_points; ++k)
1135 for (unsigned int d = 0; d < spacedim; ++d)
1136 for (unsigned int n = 0; n < spacedim; ++n)
1137 for (unsigned int i = 0; i < spacedim; ++i)
1138 for (unsigned int j = 0; j < spacedim; ++j)
1139 {
1140 transformed_shape_hessians[k][d][i][j] +=
1141 (output_data.shape_values(first + n, k) *
1142 mapping_data
1143 .jacobian_pushed_forward_2nd_derivatives
1144 [k][d][n][i][j]) +
1145 (output_data.shape_gradients[first + n][k][i] *
1146 mapping_data
1147 .jacobian_pushed_forward_grads[k][d][n][j]) +
1148 (output_data.shape_gradients[first + n][k][j] *
1149 mapping_data
1150 .jacobian_pushed_forward_grads[k][d][i][n]) -
1151 (output_data.shape_gradients[first + d][k][n] *
1152 mapping_data
1153 .jacobian_pushed_forward_grads[k][n][i][j]);
1154
1155 transformed_shape_hessians[k][d][i][j] -=
1156 (output_data.shape_values(first + d, k) *
1157 mapping_data
1158 .jacobian_pushed_forward_2nd_derivatives
1159 [k][n][n][i][j]) +
1160 (output_data.shape_gradients[first + d][k][i] *
1161 mapping_data
1162 .jacobian_pushed_forward_grads[k][n][n][j]) +
1163 (output_data.shape_gradients[first + d][k][j] *
1164 mapping_data
1165 .jacobian_pushed_forward_grads[k][n][n][i]);
1166
1167 for (unsigned int m = 0; m < spacedim; ++m)
1168 {
1169 transformed_shape_hessians[k][d][i][j] -=
1170 (mapping_data
1171 .jacobian_pushed_forward_grads[k][d][i]
1172 [m] *
1173 mapping_data
1174 .jacobian_pushed_forward_grads[k][m][n]
1175 [j] *
1176 output_data.shape_values(first + n, k)) +
1177 (mapping_data
1178 .jacobian_pushed_forward_grads[k][d][m]
1179 [j] *
1180 mapping_data
1181 .jacobian_pushed_forward_grads[k][m][i]
1182 [n] *
1183 output_data.shape_values(first + n, k));
1184
1185 transformed_shape_hessians[k][d][i][j] +=
1186 (mapping_data
1187 .jacobian_pushed_forward_grads[k][n][i]
1188 [m] *
1189 mapping_data
1190 .jacobian_pushed_forward_grads[k][m][n]
1191 [j] *
1192 output_data.shape_values(first + d, k)) +
1193 (mapping_data
1194 .jacobian_pushed_forward_grads[k][n][m]
1195 [j] *
1196 mapping_data
1197 .jacobian_pushed_forward_grads[k][m][i]
1198 [n] *
1199 output_data.shape_values(first + d, k));
1200 }
1201 }
1202
1203 for (unsigned int k = 0; k < n_q_points; ++k)
1204 for (unsigned int d = 0; d < dim; ++d)
1205 output_data.shape_hessians[first + d][k] =
1206 dof_sign * transformed_shape_hessians[k][d];
1207
1208 break;
1209 }
1210
1211 case mapping_nedelec:
1212 {
1213 for (unsigned int k = 0; k < n_q_points; ++k)
1214 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1215 fe_data.shape_grad_grads[dof_index][k + offset];
1216
1217 const ArrayView<Tensor<3, spacedim>>
1218 transformed_shape_hessians =
1219 make_array_view(fe_data.transformed_shape_hessians,
1220 offset,
1221 n_q_points);
1222 mapping.transform(
1223 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1224 offset,
1225 n_q_points),
1226 mapping_covariant_hessian,
1227 mapping_internal,
1228 transformed_shape_hessians);
1229
1230 for (unsigned int k = 0; k < n_q_points; ++k)
1231 for (unsigned int d = 0; d < spacedim; ++d)
1232 for (unsigned int n = 0; n < spacedim; ++n)
1233 for (unsigned int i = 0; i < spacedim; ++i)
1234 for (unsigned int j = 0; j < spacedim; ++j)
1235 {
1236 transformed_shape_hessians[k][d][i][j] -=
1237 (output_data.shape_values(first + n, k) *
1238 mapping_data
1239 .jacobian_pushed_forward_2nd_derivatives
1240 [k][n][d][i][j]) +
1241 (output_data.shape_gradients[first + d][k][n] *
1242 mapping_data
1243 .jacobian_pushed_forward_grads[k][n][i][j]) +
1244 (output_data.shape_gradients[first + n][k][i] *
1245 mapping_data
1246 .jacobian_pushed_forward_grads[k][n][d][j]) +
1247 (output_data.shape_gradients[first + n][k][j] *
1248 mapping_data
1249 .jacobian_pushed_forward_grads[k][n][i][d]);
1250 }
1251
1252 for (unsigned int k = 0; k < n_q_points; ++k)
1253 for (unsigned int d = 0; d < dim; ++d)
1254 output_data.shape_hessians[first + d][k] =
1255 dof_sign * transformed_shape_hessians[k][d];
1256
1257 break;
1258 }
1259
1260 default:
1261 DEAL_II_NOT_IMPLEMENTED();
1262 }
1263 }
1264
1265 // third derivatives are not implemented
1266 if (fe_data.update_each & update_3rd_derivatives)
1267 {
1268 DEAL_II_NOT_IMPLEMENTED();
1269 }
1270 }
1271}
1272
1273
1274
1275template <int dim, int spacedim>
1276UpdateFlags
1277FE_PolyTensor<dim, spacedim>::requires_update_flags(
1278 const UpdateFlags flags) const
1279{
1280 UpdateFlags out = update_default;
1281
1282 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1283 {
1284 const MappingKind mapping_kind = get_mapping_kind(i);
1285
1286 switch (mapping_kind)
1287 {
1288 case mapping_none:
1289 {
1290 if (flags & update_values)
1291 out |= update_values;
1292
1293 if (flags & update_gradients)
1294 out |= update_gradients | update_values |
1295 update_jacobian_pushed_forward_grads;
1296
1297 if (flags & update_hessians)
1298 out |= update_hessians | update_values | update_gradients |
1299 update_jacobian_pushed_forward_grads |
1300 update_jacobian_pushed_forward_2nd_derivatives;
1301 break;
1302 }
1303 case mapping_raviart_thomas:
1304 case mapping_piola:
1305 {
1306 if (flags & update_values)
1307 out |= update_values | update_piola;
1308
1309 if (flags & update_gradients)
1310 out |= update_gradients | update_values | update_piola |
1311 update_jacobian_pushed_forward_grads |
1312 update_covariant_transformation |
1313 update_contravariant_transformation;
1314
1315 if (flags & update_hessians)
1316 out |= update_hessians | update_piola | update_values |
1317 update_gradients | update_jacobian_pushed_forward_grads |
1318 update_jacobian_pushed_forward_2nd_derivatives |
1319 update_covariant_transformation;
1320
1321 break;
1322 }
1323
1324
1325 case mapping_contravariant:
1326 {
1327 if (flags & update_values)
1328 out |= update_values | update_piola;
1329
1330 if (flags & update_gradients)
1331 out |= update_gradients | update_values |
1332 update_jacobian_pushed_forward_grads |
1333 update_covariant_transformation |
1334 update_contravariant_transformation;
1335
1336 if (flags & update_hessians)
1337 out |= update_hessians | update_piola | update_values |
1338 update_gradients | update_jacobian_pushed_forward_grads |
1339 update_jacobian_pushed_forward_2nd_derivatives |
1340 update_covariant_transformation;
1341
1342 break;
1343 }
1344
1345 case mapping_nedelec:
1346 case mapping_covariant:
1347 {
1348 if (flags & update_values)
1349 out |= update_values | update_covariant_transformation;
1350
1351 if (flags & update_gradients)
1352 out |= update_gradients | update_values |
1353 update_jacobian_pushed_forward_grads |
1354 update_covariant_transformation;
1355
1356 if (flags & update_hessians)
1357 out |= update_hessians | update_values | update_gradients |
1358 update_jacobian_pushed_forward_grads |
1359 update_jacobian_pushed_forward_2nd_derivatives |
1360 update_covariant_transformation;
1361
1362 break;
1363 }
1364
1365 default:
1366 {
1367 DEAL_II_NOT_IMPLEMENTED();
1368 }
1369 }
1370 }
1371
1372 return out;
1373}
1374
1375#endif
1376// explicit instantiations
1377#include "fe/fe_poly_tensor.inst"
1378
1379
1380DEAL_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
CellAccessor::subface_case
internal::SubfaceCase< dim > subface_case(const unsigned int face_no) const
FE_PolyTensor::InternalData::transformed_shape_hessians
std::vector< Tensor< 3, spacedim > > transformed_shape_hessians
Definition fe_poly_tensor.h:521
FE_PolyTensor::InternalData::shape_grad_grads
Table< 2, DerivativeForm< 2, dim, spacedim > > shape_grad_grads
Definition fe_poly_tensor.h:510
FE_PolyTensor::InternalData::shape_grads
Table< 2, DerivativeForm< 1, dim, spacedim > > shape_grads
Definition fe_poly_tensor.h:503
FE_PolyTensor::InternalData::transformed_shape_grads
std::vector< Tensor< 2, spacedim > > transformed_shape_grads
Definition fe_poly_tensor.h:518
FE_PolyTensor::InternalData::untransformed_shape_hessian_tensors
std::vector< Tensor< 3, dim > > untransformed_shape_hessian_tensors
Definition fe_poly_tensor.h:522
FE_PolyTensor::InternalData::untransformed_shape_grads
std::vector< Tensor< 2, dim > > untransformed_shape_grads
Definition fe_poly_tensor.h:519
FE_PolyTensor::InternalData::shape_values
Table< 2, Tensor< 1, dim > > shape_values
Definition fe_poly_tensor.h:496
FE_PolyTensor::InternalData::transformed_shape_values
std::vector< Tensor< 1, spacedim > > transformed_shape_values
Definition fe_poly_tensor.h:516
FE_PolyTensor::InternalData::dof_sign_change
std::vector< double > dof_sign_change
Definition fe_poly_tensor.h:515
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
FiniteElement::InternalDataBase
friend class InternalDataBase
Definition fe.h:3217
MappingQ::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 override
Definition mapping_q.cc:1435
Quadrature::size
unsigned int size() const
hp::Collection::size
unsigned int size() const
Definition collection.h:314
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
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
FiniteElement::ExcFENotPrimitive
static ::ExceptionBase & ExcFENotPrimitive()
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Triangulation::cell_iterator
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition tria.h:1624
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
TriaAccessor< dim, dim, dim >::combined_face_orientation
types::geometric_orientation combined_face_orientation(const unsigned int face) const
mapping_cartesian.h
mapping
MappingQ< dim, spacedim > StaticMappingQ1< dim, spacedim >::mapping
Definition mapping_q1.h:101
GridGenerator::reference_cell
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell< dim > &reference_cell)
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:8901
internal::FE_PolyTensor
Definition fe_poly_tensor.cc:37
numbers::default_geometric_orientation
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:350
types::geometric_orientation
std::uint8_t geometric_orientation
Definition types.h:38
n_q_points
*  const unsigned int n_q_points
Definition nonlinear-heat_transfer_with_AD_NOX.h:1386
polynomials_abf.h
polynomials_bdm.h
polynomials_bernardi_raugel.h
polynomials_nedelec.h
polynomials_raviart_thomas.h
polynomials_rt_bubbles.h
GeometryInfo::faces_per_cell
static constexpr unsigned int faces_per_cell
Definition geometry_info.h:1984
tria_orientation.h
Generated by doxygen 1.16.1