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