Reference documentation for deal.II version Git e3a3ec7800 2020-08-07 14:08:19 +0200
\(\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\}}\)
fe_poly_tensor.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
26 
28 #include <deal.II/fe/fe_values.h>
30 
31 #include <deal.II/grid/tria.h>
32 
34 
35 namespace internal
36 {
37  namespace FE_PolyTensor
38  {
39  namespace
40  {
41  //---------------------------------------------------------------------------
42  // Utility method, which is used to determine the change of sign for
43  // the DoFs on the faces of the given cell.
44  //---------------------------------------------------------------------------
45 
52  template <int spacedim>
53  void
54  get_face_sign_change_rt(const ::Triangulation<1>::cell_iterator &,
56  const std::vector<MappingKind> &,
57  std::vector<double> &)
58  {
59  // nothing to do in 1d
60  }
61 
62 
63 
64  // template<int spacedim>
65  void
66  get_face_sign_change_rt(
67  const ::Triangulation<2>::cell_iterator &cell,
68  const FiniteElement<2, 2> & fe,
69  const std::vector<MappingKind> & mapping_kind,
70  std::vector<double> & face_sign)
71  {
72  const unsigned int dim = 2;
73  const unsigned int spacedim = 2;
74 
75  for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2;
76  f < GeometryInfo<dim>::faces_per_cell;
77  ++f)
78  {
80  cell->face(f);
81  if (!face->at_boundary())
82  {
83  const unsigned int nn = cell->neighbor_face_no(f);
84 
86  for (unsigned int j = 0; j < fe.n_dofs_per_face(); ++j)
87  {
88  const unsigned int cell_j = fe.face_to_cell_index(j, f);
89 
90  Assert(f * fe.n_dofs_per_face() + j < face_sign.size(),
92  Assert(mapping_kind.size() == 1 ||
93  cell_j < mapping_kind.size(),
95 
96  // TODO: This is probably only going to work for those
97  // elements for which all dofs are face dofs
98  if ((mapping_kind.size() > 1 ?
99  mapping_kind[cell_j] :
100  mapping_kind[0]) == mapping_raviart_thomas)
101  face_sign[f * fe.n_dofs_per_face() + j] = -1.0;
102  }
103  }
104  }
105  }
106 
107 
108 
109  template <int spacedim>
110  void
111  get_face_sign_change_rt(
112  const ::Triangulation<3>::cell_iterator & /*cell*/,
113  const FiniteElement<3, spacedim> & /*fe*/,
114  const std::vector<MappingKind> & /*mapping_kind*/,
115  std::vector<double> & /*face_sign*/)
116  {
117  // TODO: think about what it would take here
118  }
119 
120 
121 
122  template <int spacedim>
123  void
124  get_face_sign_change_nedelec(
125  const ::Triangulation<1>::cell_iterator & /*cell*/,
126  const FiniteElement<1, spacedim> & /*fe*/,
127  const std::vector<MappingKind> & /*mapping_kind*/,
128  std::vector<double> & /*face_sign*/)
129  {
130  // nothing to do in 1d
131  }
132 
133 
134 
135  template <int spacedim>
136  void
137  get_face_sign_change_nedelec(
138  const ::Triangulation<2>::cell_iterator & /*cell*/,
139  const FiniteElement<2, spacedim> & /*fe*/,
140  const std::vector<MappingKind> & /*mapping_kind*/,
141  std::vector<double> & /*face_sign*/)
142  {
143  // TODO: think about what it would take here
144  }
145 
146 
147  template <int spacedim>
148  void
149  get_face_sign_change_nedelec(
150  const ::Triangulation<3>::cell_iterator &cell,
151  const FiniteElement<3, spacedim> & /*fe*/,
152  const std::vector<MappingKind> &mapping_kind,
153  std::vector<double> & face_sign)
154  {
155  const unsigned int dim = 3;
156  // TODO: This is probably only going to work for those elements for
157  // which all dofs are face dofs
158  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
159  if (!(cell->line_orientation(l)) &&
160  mapping_kind[0] == mapping_nedelec)
161  face_sign[l] = -1.0;
162  }
163  } // namespace
164  } // namespace FE_PolyTensor
165 } // namespace internal
166 
167 
168 
169 template <int dim, int spacedim>
171  const TensorPolynomialsBase<dim> &polynomials,
172  const FiniteElementData<dim> & fe_data,
173  const std::vector<bool> & restriction_is_additive_flags,
174  const std::vector<ComponentMask> &nonzero_components)
175  : FiniteElement<dim, spacedim>(fe_data,
176  restriction_is_additive_flags,
177  nonzero_components)
178  , mapping_kind({MappingKind::mapping_none})
179  , poly_space(polynomials.clone())
180 {
181  cached_point(0) = -1;
182  // Set up the table converting
183  // components to base
184  // components. Since we have only
185  // one base element, everything
186  // remains zero except the
187  // component in the base, which is
188  // the component itself
189  for (unsigned int comp = 0; comp < this->n_components(); ++comp)
190  this->component_to_base_table[comp].first.second = comp;
191 }
192 
193 
194 
195 template <int dim, int spacedim>
197  : FiniteElement<dim, spacedim>(fe)
199  , poly_space(fe.poly_space->clone())
201 {}
202 
203 
204 
205 template <int dim, int spacedim>
206 bool
208 {
209  return mapping_kind.size() == 1;
210 }
211 
212 
213 
214 template <int dim, int spacedim>
217 {
218  if (single_mapping_kind())
219  return mapping_kind[0];
220 
221  AssertIndexRange(i, mapping_kind.size());
222  return mapping_kind[i];
223 }
224 
225 
226 
227 template <int dim, int spacedim>
228 double
230  const Point<dim> &) const
231 
232 {
234  return 0.;
235 }
236 
237 
238 
239 template <int dim, int spacedim>
240 double
242  const unsigned int i,
243  const Point<dim> & p,
244  const unsigned int component) const
245 {
246  AssertIndexRange(i, this->n_dofs_per_cell());
247  AssertIndexRange(component, dim);
248 
249  std::lock_guard<std::mutex> lock(cache_mutex);
250 
251  if (cached_point != p || cached_values.size() == 0)
252  {
253  cached_point = p;
254  cached_values.resize(poly_space->n());
255 
256  std::vector<Tensor<4, dim>> dummy1;
257  std::vector<Tensor<5, dim>> dummy2;
258  poly_space->evaluate(
259  p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
260  }
261 
262  double s = 0;
263  if (inverse_node_matrix.n_cols() == 0)
264  return cached_values[i][component];
265  else
266  for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
267  s += inverse_node_matrix(j, i) * cached_values[j][component];
268  return s;
269 }
270 
271 
272 
273 template <int dim, int spacedim>
276  const Point<dim> &) const
277 {
279  return Tensor<1, dim>();
280 }
281 
282 
283 
284 template <int dim, int spacedim>
287  const unsigned int i,
288  const Point<dim> & p,
289  const unsigned int component) const
290 {
291  AssertIndexRange(i, this->n_dofs_per_cell());
292  AssertIndexRange(component, dim);
293 
294  std::lock_guard<std::mutex> lock(cache_mutex);
295 
296  if (cached_point != p || cached_grads.size() == 0)
297  {
298  cached_point = p;
299  cached_grads.resize(poly_space->n());
300 
301  std::vector<Tensor<4, dim>> dummy1;
302  std::vector<Tensor<5, dim>> dummy2;
303  poly_space->evaluate(
304  p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
305  }
306 
307  Tensor<1, dim> s;
308  if (inverse_node_matrix.n_cols() == 0)
309  return cached_grads[i][component];
310  else
311  for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
312  s += inverse_node_matrix(j, i) * cached_grads[j][component];
313 
314  return s;
315 }
316 
317 
318 
319 template <int dim, int spacedim>
322  const Point<dim> &) const
323 {
325  return Tensor<2, dim>();
326 }
327 
328 
329 
330 template <int dim, int spacedim>
333  const unsigned int i,
334  const Point<dim> & p,
335  const unsigned int component) const
336 {
337  AssertIndexRange(i, this->n_dofs_per_cell());
338  AssertIndexRange(component, dim);
339 
340  std::lock_guard<std::mutex> lock(cache_mutex);
341 
342  if (cached_point != p || cached_grad_grads.size() == 0)
343  {
344  cached_point = p;
345  cached_grad_grads.resize(poly_space->n());
346 
347  std::vector<Tensor<4, dim>> dummy1;
348  std::vector<Tensor<5, dim>> dummy2;
349  poly_space->evaluate(
350  p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
351  }
352 
353  Tensor<2, dim> s;
354  if (inverse_node_matrix.n_cols() == 0)
355  return cached_grad_grads[i][component];
356  else
357  for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
358  s += inverse_node_matrix(i, j) * cached_grad_grads[j][component];
359 
360  return s;
361 }
362 
363 
364 //---------------------------------------------------------------------------
365 // Fill data of FEValues
366 //---------------------------------------------------------------------------
367 
368 template <int dim, int spacedim>
369 void
371  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
372  const CellSimilarity::Similarity cell_similarity,
373  const Quadrature<dim> & quadrature,
374  const Mapping<dim, spacedim> & mapping,
375  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
376  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
377  spacedim>
378  & mapping_data,
379  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
381  spacedim>
382  &output_data) const
383 {
384  // convert data object to internal
385  // data for this class. fails with
386  // an exception if that is not
387  // possible
388  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
389  ExcInternalError());
390  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
391 
392  const unsigned int n_q_points = quadrature.size();
393 
394  Assert(!(fe_data.update_each & update_values) ||
395  fe_data.shape_values.size()[0] == this->n_dofs_per_cell(),
396  ExcDimensionMismatch(fe_data.shape_values.size()[0],
397  this->n_dofs_per_cell()));
398  Assert(!(fe_data.update_each & update_values) ||
399  fe_data.shape_values.size()[1] == n_q_points,
400  ExcDimensionMismatch(fe_data.shape_values.size()[1], n_q_points));
401 
402  // Create table with sign changes, due to the special structure of the RT
403  // elements.
404  // TODO: Preliminary hack to demonstrate the overall principle!
405 
406  // Compute eventual sign changes depending on the neighborhood
407  // between two faces.
408  std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0);
409 
410  internal::FE_PolyTensor::get_face_sign_change_rt(cell,
411  *this,
412  this->mapping_kind,
413  fe_data.sign_change);
414 
415  internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
416  *this,
417  this->mapping_kind,
418  fe_data.sign_change);
419 
420 
421  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
422  {
424 
425  const unsigned int first =
426  output_data.shape_function_to_row_table[i * this->n_components() +
427  this->get_nonzero_components(i)
429 
430  // update the shape function values as necessary
431  //
432  // we only need to do this if the current cell is not a translation of
433  // the previous one; or, even if it is a translation, if we use mappings
434  // other than the standard mappings that require us to recompute values
435  // and derivatives because of possible sign changes
436  if (fe_data.update_each & update_values &&
437  ((cell_similarity != CellSimilarity::translation) ||
438  ((mapping_kind == mapping_piola) ||
439  (mapping_kind == mapping_raviart_thomas) ||
440  (mapping_kind == mapping_nedelec))))
441  {
442  switch (mapping_kind)
443  {
444  case mapping_none:
445  {
446  for (unsigned int k = 0; k < n_q_points; ++k)
447  for (unsigned int d = 0; d < dim; ++d)
448  output_data.shape_values(first + d, k) =
449  fe_data.shape_values[i][k][d];
450  break;
451  }
452 
453  case mapping_covariant:
455  {
456  mapping.transform(make_array_view(fe_data.shape_values, i),
457  mapping_kind,
458  mapping_internal,
460  fe_data.transformed_shape_values));
461 
462  for (unsigned int k = 0; k < n_q_points; ++k)
463  for (unsigned int d = 0; d < dim; ++d)
464  output_data.shape_values(first + d, k) =
465  fe_data.transformed_shape_values[k][d];
466 
467  break;
468  }
469 
471  case mapping_piola:
472  {
473  mapping.transform(make_array_view(fe_data.shape_values, i),
475  mapping_internal,
477  fe_data.transformed_shape_values));
478  for (unsigned int k = 0; k < n_q_points; ++k)
479  for (unsigned int d = 0; d < dim; ++d)
480  output_data.shape_values(first + d, k) =
481  fe_data.sign_change[i] *
482  fe_data.transformed_shape_values[k][d];
483  break;
484  }
485 
486  case mapping_nedelec:
487  {
488  mapping.transform(make_array_view(fe_data.shape_values, i),
490  mapping_internal,
492  fe_data.transformed_shape_values));
493 
494  for (unsigned int k = 0; k < n_q_points; ++k)
495  for (unsigned int d = 0; d < dim; ++d)
496  output_data.shape_values(first + d, k) =
497  fe_data.sign_change[i] *
498  fe_data.transformed_shape_values[k][d];
499 
500  break;
501  }
502 
503  default:
504  Assert(false, ExcNotImplemented());
505  }
506  }
507 
508  // update gradients. apply the same logic as above
509  if (fe_data.update_each & update_gradients &&
510  ((cell_similarity != CellSimilarity::translation) ||
511  ((mapping_kind == mapping_piola) ||
512  (mapping_kind == mapping_raviart_thomas) ||
513  (mapping_kind == mapping_nedelec))))
514 
515  {
516  switch (mapping_kind)
517  {
518  case mapping_none:
519  {
520  mapping.transform(make_array_view(fe_data.shape_grads, i),
522  mapping_internal,
524  fe_data.transformed_shape_grads));
525  for (unsigned int k = 0; k < n_q_points; ++k)
526  for (unsigned int d = 0; d < dim; ++d)
527  output_data.shape_gradients[first + d][k] =
528  fe_data.transformed_shape_grads[k][d];
529  break;
530  }
531  case mapping_covariant:
532  {
533  mapping.transform(make_array_view(fe_data.shape_grads, i),
535  mapping_internal,
537  fe_data.transformed_shape_grads));
538 
539  for (unsigned int k = 0; k < n_q_points; ++k)
540  for (unsigned int d = 0; d < spacedim; ++d)
541  for (unsigned int n = 0; n < spacedim; ++n)
542  fe_data.transformed_shape_grads[k][d] -=
543  output_data.shape_values(first + n, k) *
544  mapping_data.jacobian_pushed_forward_grads[k][n][d];
545 
546  for (unsigned int k = 0; k < n_q_points; ++k)
547  for (unsigned int d = 0; d < dim; ++d)
548  output_data.shape_gradients[first + d][k] =
549  fe_data.transformed_shape_grads[k][d];
550 
551  break;
552  }
554  {
555  for (unsigned int k = 0; k < n_q_points; ++k)
556  fe_data.untransformed_shape_grads[k] =
557  fe_data.shape_grads[i][k];
558  mapping.transform(
559  make_array_view(fe_data.untransformed_shape_grads),
561  mapping_internal,
562  make_array_view(fe_data.transformed_shape_grads));
563 
564  for (unsigned int k = 0; k < n_q_points; ++k)
565  for (unsigned int d = 0; d < spacedim; ++d)
566  for (unsigned int n = 0; n < spacedim; ++n)
567  fe_data.transformed_shape_grads[k][d] +=
568  output_data.shape_values(first + n, k) *
569  mapping_data.jacobian_pushed_forward_grads[k][d][n];
570 
571 
572  for (unsigned int k = 0; k < n_q_points; ++k)
573  for (unsigned int d = 0; d < dim; ++d)
574  output_data.shape_gradients[first + d][k] =
575  fe_data.transformed_shape_grads[k][d];
576 
577  break;
578  }
580  case mapping_piola:
581  {
582  for (unsigned int k = 0; k < n_q_points; ++k)
583  fe_data.untransformed_shape_grads[k] =
584  fe_data.shape_grads[i][k];
585  mapping.transform(
586  make_array_view(fe_data.untransformed_shape_grads),
588  mapping_internal,
589  make_array_view(fe_data.transformed_shape_grads));
590 
591  for (unsigned int k = 0; k < n_q_points; ++k)
592  for (unsigned int d = 0; d < spacedim; ++d)
593  for (unsigned int n = 0; n < spacedim; ++n)
594  fe_data.transformed_shape_grads[k][d] +=
595  (output_data.shape_values(first + n, k) *
596  mapping_data
597  .jacobian_pushed_forward_grads[k][d][n]) -
598  (output_data.shape_values(first + d, k) *
599  mapping_data.jacobian_pushed_forward_grads[k][n][n]);
600 
601  for (unsigned int k = 0; k < n_q_points; ++k)
602  for (unsigned int d = 0; d < dim; ++d)
603  output_data.shape_gradients[first + d][k] =
604  fe_data.sign_change[i] *
605  fe_data.transformed_shape_grads[k][d];
606 
607  break;
608  }
609 
610  case mapping_nedelec:
611  {
612  // treat the gradients of
613  // this particular shape
614  // function at all
615  // q-points. if Dv is the
616  // gradient of the shape
617  // function on the unit
618  // cell, then
619  // (J^-T)Dv(J^-1) is the
620  // value we want to have on
621  // the real cell.
622  for (unsigned int k = 0; k < n_q_points; ++k)
623  fe_data.untransformed_shape_grads[k] =
624  fe_data.shape_grads[i][k];
625 
626  mapping.transform(
627  make_array_view(fe_data.untransformed_shape_grads),
629  mapping_internal,
630  make_array_view(fe_data.transformed_shape_grads));
631 
632  for (unsigned int k = 0; k < n_q_points; ++k)
633  for (unsigned int d = 0; d < spacedim; ++d)
634  for (unsigned int n = 0; n < spacedim; ++n)
635  fe_data.transformed_shape_grads[k][d] -=
636  output_data.shape_values(first + n, k) *
637  mapping_data.jacobian_pushed_forward_grads[k][n][d];
638 
639  for (unsigned int k = 0; k < n_q_points; ++k)
640  for (unsigned int d = 0; d < dim; ++d)
641  output_data.shape_gradients[first + d][k] =
642  fe_data.sign_change[i] *
643  fe_data.transformed_shape_grads[k][d];
644 
645  break;
646  }
647 
648  default:
649  Assert(false, ExcNotImplemented());
650  }
651  }
652 
653  // update hessians. apply the same logic as above
654  if (fe_data.update_each & update_hessians &&
655  ((cell_similarity != CellSimilarity::translation) ||
656  ((mapping_kind == mapping_piola) ||
657  (mapping_kind == mapping_raviart_thomas) ||
658  (mapping_kind == mapping_nedelec))))
659 
660  {
661  switch (mapping_kind)
662  {
663  case mapping_none:
664  {
665  mapping.transform(
666  make_array_view(fe_data.shape_grad_grads, i),
668  mapping_internal,
669  make_array_view(fe_data.transformed_shape_hessians));
670 
671  for (unsigned int k = 0; k < n_q_points; ++k)
672  for (unsigned int d = 0; d < spacedim; ++d)
673  for (unsigned int n = 0; n < spacedim; ++n)
674  fe_data.transformed_shape_hessians[k][d] -=
675  output_data.shape_gradients[first + d][k][n] *
676  mapping_data.jacobian_pushed_forward_grads[k][n];
677 
678  for (unsigned int k = 0; k < n_q_points; ++k)
679  for (unsigned int d = 0; d < dim; ++d)
680  output_data.shape_hessians[first + d][k] =
681  fe_data.transformed_shape_hessians[k][d];
682 
683  break;
684  }
685  case mapping_covariant:
686  {
687  for (unsigned int k = 0; k < n_q_points; ++k)
688  fe_data.untransformed_shape_hessian_tensors[k] =
689  fe_data.shape_grad_grads[i][k];
690 
691  mapping.transform(
693  fe_data.untransformed_shape_hessian_tensors),
695  mapping_internal,
696  make_array_view(fe_data.transformed_shape_hessians));
697 
698  for (unsigned int k = 0; k < n_q_points; ++k)
699  for (unsigned int d = 0; d < spacedim; ++d)
700  for (unsigned int n = 0; n < spacedim; ++n)
701  for (unsigned int i = 0; i < spacedim; ++i)
702  for (unsigned int j = 0; j < spacedim; ++j)
703  {
704  fe_data.transformed_shape_hessians[k][d][i][j] -=
705  (output_data.shape_values(first + n, k) *
706  mapping_data
707  .jacobian_pushed_forward_2nd_derivatives
708  [k][n][d][i][j]) +
709  (output_data.shape_gradients[first + d][k][n] *
710  mapping_data
711  .jacobian_pushed_forward_grads[k][n][i][j]) +
712  (output_data.shape_gradients[first + n][k][i] *
713  mapping_data
714  .jacobian_pushed_forward_grads[k][n][d][j]) +
715  (output_data.shape_gradients[first + n][k][j] *
716  mapping_data
717  .jacobian_pushed_forward_grads[k][n][i][d]);
718  }
719 
720  for (unsigned int k = 0; k < n_q_points; ++k)
721  for (unsigned int d = 0; d < dim; ++d)
722  output_data.shape_hessians[first + d][k] =
723  fe_data.transformed_shape_hessians[k][d];
724 
725  break;
726  }
728  {
729  for (unsigned int k = 0; k < n_q_points; ++k)
730  fe_data.untransformed_shape_hessian_tensors[k] =
731  fe_data.shape_grad_grads[i][k];
732 
733  mapping.transform(
735  fe_data.untransformed_shape_hessian_tensors),
737  mapping_internal,
738  make_array_view(fe_data.transformed_shape_hessians));
739 
740  for (unsigned int k = 0; k < n_q_points; ++k)
741  for (unsigned int d = 0; d < spacedim; ++d)
742  for (unsigned int n = 0; n < spacedim; ++n)
743  for (unsigned int i = 0; i < spacedim; ++i)
744  for (unsigned int j = 0; j < spacedim; ++j)
745  {
746  fe_data.transformed_shape_hessians[k][d][i][j] +=
747  (output_data.shape_values(first + n, k) *
748  mapping_data
749  .jacobian_pushed_forward_2nd_derivatives
750  [k][d][n][i][j]) +
751  (output_data.shape_gradients[first + n][k][i] *
752  mapping_data
753  .jacobian_pushed_forward_grads[k][d][n][j]) +
754  (output_data.shape_gradients[first + n][k][j] *
755  mapping_data
756  .jacobian_pushed_forward_grads[k][d][i][n]) -
757  (output_data.shape_gradients[first + d][k][n] *
758  mapping_data
759  .jacobian_pushed_forward_grads[k][n][i][j]);
760  for (unsigned int m = 0; m < spacedim; ++m)
761  fe_data
762  .transformed_shape_hessians[k][d][i][j] -=
763  (mapping_data
764  .jacobian_pushed_forward_grads[k][d][i]
765  [m] *
766  mapping_data
767  .jacobian_pushed_forward_grads[k][m][n]
768  [j] *
769  output_data.shape_values(first + n, k)) +
770  (mapping_data
771  .jacobian_pushed_forward_grads[k][d][m]
772  [j] *
773  mapping_data
774  .jacobian_pushed_forward_grads[k][m][i]
775  [n] *
776  output_data.shape_values(first + n, k));
777  }
778 
779  for (unsigned int k = 0; k < n_q_points; ++k)
780  for (unsigned int d = 0; d < dim; ++d)
781  output_data.shape_hessians[first + d][k] =
782  fe_data.transformed_shape_hessians[k][d];
783 
784  break;
785  }
787  case mapping_piola:
788  {
789  for (unsigned int k = 0; k < n_q_points; ++k)
790  fe_data.untransformed_shape_hessian_tensors[k] =
791  fe_data.shape_grad_grads[i][k];
792 
793  mapping.transform(
795  fe_data.untransformed_shape_hessian_tensors),
797  mapping_internal,
798  make_array_view(fe_data.transformed_shape_hessians));
799 
800  for (unsigned int k = 0; k < n_q_points; ++k)
801  for (unsigned int d = 0; d < spacedim; ++d)
802  for (unsigned int n = 0; n < spacedim; ++n)
803  for (unsigned int i = 0; i < spacedim; ++i)
804  for (unsigned int j = 0; j < spacedim; ++j)
805  {
806  fe_data.transformed_shape_hessians[k][d][i][j] +=
807  (output_data.shape_values(first + n, k) *
808  mapping_data
809  .jacobian_pushed_forward_2nd_derivatives
810  [k][d][n][i][j]) +
811  (output_data.shape_gradients[first + n][k][i] *
812  mapping_data
813  .jacobian_pushed_forward_grads[k][d][n][j]) +
814  (output_data.shape_gradients[first + n][k][j] *
815  mapping_data
816  .jacobian_pushed_forward_grads[k][d][i][n]) -
817  (output_data.shape_gradients[first + d][k][n] *
818  mapping_data
819  .jacobian_pushed_forward_grads[k][n][i][j]);
820 
821  fe_data.transformed_shape_hessians[k][d][i][j] -=
822  (output_data.shape_values(first + d, k) *
823  mapping_data
824  .jacobian_pushed_forward_2nd_derivatives
825  [k][n][n][i][j]) +
826  (output_data.shape_gradients[first + d][k][i] *
827  mapping_data
828  .jacobian_pushed_forward_grads[k][n][n][j]) +
829  (output_data.shape_gradients[first + d][k][j] *
830  mapping_data
831  .jacobian_pushed_forward_grads[k][n][n][i]);
832 
833  for (unsigned int m = 0; m < spacedim; ++m)
834  {
835  fe_data
836  .transformed_shape_hessians[k][d][i][j] -=
837  (mapping_data
838  .jacobian_pushed_forward_grads[k][d][i]
839  [m] *
840  mapping_data
841  .jacobian_pushed_forward_grads[k][m][n]
842  [j] *
843  output_data.shape_values(first + n, k)) +
844  (mapping_data
845  .jacobian_pushed_forward_grads[k][d][m]
846  [j] *
847  mapping_data
848  .jacobian_pushed_forward_grads[k][m][i]
849  [n] *
850  output_data.shape_values(first + n, k));
851 
852  fe_data
853  .transformed_shape_hessians[k][d][i][j] +=
854  (mapping_data
855  .jacobian_pushed_forward_grads[k][n][i]
856  [m] *
857  mapping_data
858  .jacobian_pushed_forward_grads[k][m][n]
859  [j] *
860  output_data.shape_values(first + d, k)) +
861  (mapping_data
862  .jacobian_pushed_forward_grads[k][n][m]
863  [j] *
864  mapping_data
865  .jacobian_pushed_forward_grads[k][m][i]
866  [n] *
867  output_data.shape_values(first + d, k));
868  }
869  }
870 
871  for (unsigned int k = 0; k < n_q_points; ++k)
872  for (unsigned int d = 0; d < dim; ++d)
873  output_data.shape_hessians[first + d][k] =
874  fe_data.sign_change[i] *
875  fe_data.transformed_shape_hessians[k][d];
876 
877  break;
878  }
879 
880  case mapping_nedelec:
881  {
882  for (unsigned int k = 0; k < n_q_points; ++k)
883  fe_data.untransformed_shape_hessian_tensors[k] =
884  fe_data.shape_grad_grads[i][k];
885 
886  mapping.transform(
888  fe_data.untransformed_shape_hessian_tensors),
890  mapping_internal,
891  make_array_view(fe_data.transformed_shape_hessians));
892 
893  for (unsigned int k = 0; k < n_q_points; ++k)
894  for (unsigned int d = 0; d < spacedim; ++d)
895  for (unsigned int n = 0; n < spacedim; ++n)
896  for (unsigned int i = 0; i < spacedim; ++i)
897  for (unsigned int j = 0; j < spacedim; ++j)
898  {
899  fe_data.transformed_shape_hessians[k][d][i][j] -=
900  (output_data.shape_values(first + n, k) *
901  mapping_data
902  .jacobian_pushed_forward_2nd_derivatives
903  [k][n][d][i][j]) +
904  (output_data.shape_gradients[first + d][k][n] *
905  mapping_data
906  .jacobian_pushed_forward_grads[k][n][i][j]) +
907  (output_data.shape_gradients[first + n][k][i] *
908  mapping_data
909  .jacobian_pushed_forward_grads[k][n][d][j]) +
910  (output_data.shape_gradients[first + n][k][j] *
911  mapping_data
912  .jacobian_pushed_forward_grads[k][n][i][d]);
913  }
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_hessians[first + d][k] =
918  fe_data.sign_change[i] *
919  fe_data.transformed_shape_hessians[k][d];
920 
921  break;
922  }
923 
924  default:
925  Assert(false, ExcNotImplemented());
926  }
927  }
928 
929  // third derivatives are not implemented
930  if (fe_data.update_each & update_3rd_derivatives &&
931  ((cell_similarity != CellSimilarity::translation) ||
932  ((mapping_kind == mapping_piola) ||
933  (mapping_kind == mapping_raviart_thomas) ||
934  (mapping_kind == mapping_nedelec))))
935  {
936  Assert(false, ExcNotImplemented())
937  }
938  }
939 }
940 
941 
942 
943 template <int dim, int spacedim>
944 void
946  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
947  const unsigned int face_no,
948  const Quadrature<dim - 1> & quadrature,
949  const Mapping<dim, spacedim> & mapping,
950  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
951  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
952  spacedim>
953  & mapping_data,
954  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
956  spacedim>
957  &output_data) const
958 {
959  // convert data object to internal
960  // data for this class. fails with
961  // an exception if that is not
962  // possible
963  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
964  ExcInternalError());
965  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
966 
967  const unsigned int n_q_points = quadrature.size();
968  // offset determines which data set
969  // to take (all data sets for all
970  // faces are stored contiguously)
971 
972  const auto offset =
974  face_no,
975  cell->face_orientation(face_no),
976  cell->face_flip(face_no),
977  cell->face_rotation(face_no),
978  n_q_points);
979 
980  // TODO: Size assertions
981 
982  // Create table with sign changes, due to the special structure of the RT
983  // elements.
984  // TODO: Preliminary hack to demonstrate the overall prinicple!
985 
986  // Compute eventual sign changes depending
987  // on the neighborhood between two faces.
988  std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0);
989 
990  internal::FE_PolyTensor::get_face_sign_change_rt(cell,
991  *this,
992  this->mapping_kind,
993  fe_data.sign_change);
994 
995  internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
996  *this,
997  this->mapping_kind,
998  fe_data.sign_change);
999 
1000  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1001  {
1003 
1004  const unsigned int first =
1005  output_data.shape_function_to_row_table[i * this->n_components() +
1006  this->get_nonzero_components(i)
1008 
1009  if (fe_data.update_each & update_values)
1010  {
1011  switch (mapping_kind)
1012  {
1013  case mapping_none:
1014  {
1015  for (unsigned int k = 0; k < n_q_points; ++k)
1016  for (unsigned int d = 0; d < dim; ++d)
1017  output_data.shape_values(first + d, k) =
1018  fe_data.shape_values[i][k + offset][d];
1019  break;
1020  }
1021 
1022  case mapping_covariant:
1023  case mapping_contravariant:
1024  {
1026  transformed_shape_values =
1027  make_array_view(fe_data.transformed_shape_values,
1028  offset,
1029  n_q_points);
1030  mapping.transform(make_array_view(fe_data.shape_values,
1031  i,
1032  offset,
1033  n_q_points),
1034  mapping_kind,
1035  mapping_internal,
1036  transformed_shape_values);
1037 
1038  for (unsigned int k = 0; k < n_q_points; ++k)
1039  for (unsigned int d = 0; d < dim; ++d)
1040  output_data.shape_values(first + d, k) =
1041  transformed_shape_values[k][d];
1042 
1043  break;
1044  }
1046  case mapping_piola:
1047  {
1049  transformed_shape_values =
1050  make_array_view(fe_data.transformed_shape_values,
1051  offset,
1052  n_q_points);
1053  mapping.transform(make_array_view(fe_data.shape_values,
1054  i,
1055  offset,
1056  n_q_points),
1057  mapping_piola,
1058  mapping_internal,
1059  transformed_shape_values);
1060  for (unsigned int k = 0; k < n_q_points; ++k)
1061  for (unsigned int d = 0; d < dim; ++d)
1062  output_data.shape_values(first + d, k) =
1063  fe_data.sign_change[i] * transformed_shape_values[k][d];
1064  break;
1065  }
1066 
1067  case mapping_nedelec:
1068  {
1070  transformed_shape_values =
1071  make_array_view(fe_data.transformed_shape_values,
1072  offset,
1073  n_q_points);
1074  mapping.transform(make_array_view(fe_data.shape_values,
1075  i,
1076  offset,
1077  n_q_points),
1079  mapping_internal,
1080  transformed_shape_values);
1081 
1082  for (unsigned int k = 0; k < n_q_points; ++k)
1083  for (unsigned int d = 0; d < dim; ++d)
1084  output_data.shape_values(first + d, k) =
1085  fe_data.sign_change[i] * transformed_shape_values[k][d];
1086 
1087  break;
1088  }
1089 
1090  default:
1091  Assert(false, ExcNotImplemented());
1092  }
1093  }
1094 
1095  if (fe_data.update_each & update_gradients)
1096  {
1097  switch (mapping_kind)
1098  {
1099  case mapping_none:
1100  {
1101  const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1102  make_array_view(fe_data.transformed_shape_grads,
1103  offset,
1104  n_q_points);
1105  mapping.transform(
1106  make_array_view(fe_data.shape_grads, i, offset, n_q_points),
1108  mapping_internal,
1109  transformed_shape_grads);
1110  for (unsigned int k = 0; k < n_q_points; ++k)
1111  for (unsigned int d = 0; d < dim; ++d)
1112  output_data.shape_gradients[first + d][k] =
1113  transformed_shape_grads[k][d];
1114  break;
1115  }
1116 
1117  case mapping_covariant:
1118  {
1119  const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1120  make_array_view(fe_data.transformed_shape_grads,
1121  offset,
1122  n_q_points);
1123  mapping.transform(
1124  make_array_view(fe_data.shape_grads, i, offset, n_q_points),
1126  mapping_internal,
1127  transformed_shape_grads);
1128 
1129  for (unsigned int k = 0; k < n_q_points; ++k)
1130  for (unsigned int d = 0; d < spacedim; ++d)
1131  for (unsigned int n = 0; n < spacedim; ++n)
1132  transformed_shape_grads[k][d] -=
1133  output_data.shape_values(first + n, k) *
1134  mapping_data.jacobian_pushed_forward_grads[k][n][d];
1135 
1136  for (unsigned int k = 0; k < n_q_points; ++k)
1137  for (unsigned int d = 0; d < dim; ++d)
1138  output_data.shape_gradients[first + d][k] =
1139  transformed_shape_grads[k][d];
1140  break;
1141  }
1142 
1143  case mapping_contravariant:
1144  {
1145  const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1146  make_array_view(fe_data.transformed_shape_grads,
1147  offset,
1148  n_q_points);
1149  for (unsigned int k = 0; k < n_q_points; ++k)
1150  fe_data.untransformed_shape_grads[k + offset] =
1151  fe_data.shape_grads[i][k + offset];
1152  mapping.transform(
1153  make_array_view(fe_data.untransformed_shape_grads,
1154  offset,
1155  n_q_points),
1157  mapping_internal,
1158  transformed_shape_grads);
1159 
1160  for (unsigned int k = 0; k < n_q_points; ++k)
1161  for (unsigned int d = 0; d < spacedim; ++d)
1162  for (unsigned int n = 0; n < spacedim; ++n)
1163  transformed_shape_grads[k][d] +=
1164  output_data.shape_values(first + n, k) *
1165  mapping_data.jacobian_pushed_forward_grads[k][d][n];
1166 
1167  for (unsigned int k = 0; k < n_q_points; ++k)
1168  for (unsigned int d = 0; d < dim; ++d)
1169  output_data.shape_gradients[first + d][k] =
1170  transformed_shape_grads[k][d];
1171 
1172  break;
1173  }
1174 
1176  case mapping_piola:
1177  {
1178  const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1179  make_array_view(fe_data.transformed_shape_grads,
1180  offset,
1181  n_q_points);
1182  for (unsigned int k = 0; k < n_q_points; ++k)
1183  fe_data.untransformed_shape_grads[k + offset] =
1184  fe_data.shape_grads[i][k + offset];
1185  mapping.transform(
1186  make_array_view(fe_data.untransformed_shape_grads,
1187  offset,
1188  n_q_points),
1190  mapping_internal,
1191  transformed_shape_grads);
1192 
1193  for (unsigned int k = 0; k < n_q_points; ++k)
1194  for (unsigned int d = 0; d < spacedim; ++d)
1195  for (unsigned int n = 0; n < spacedim; ++n)
1196  transformed_shape_grads[k][d] +=
1197  (output_data.shape_values(first + n, k) *
1198  mapping_data
1199  .jacobian_pushed_forward_grads[k][d][n]) -
1200  (output_data.shape_values(first + d, k) *
1201  mapping_data.jacobian_pushed_forward_grads[k][n][n]);
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_gradients[first + d][k] =
1206  fe_data.sign_change[i] * transformed_shape_grads[k][d];
1207 
1208  break;
1209  }
1210 
1211  case mapping_nedelec:
1212  {
1213  // treat the gradients of
1214  // this particular shape
1215  // function at all
1216  // q-points. if Dv is the
1217  // gradient of the shape
1218  // function on the unit
1219  // cell, then
1220  // (J^-T)Dv(J^-1) is the
1221  // value we want to have on
1222  // the real cell.
1223  for (unsigned int k = 0; k < n_q_points; ++k)
1224  fe_data.untransformed_shape_grads[k + offset] =
1225  fe_data.shape_grads[i][k + offset];
1226 
1227  const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1228  make_array_view(fe_data.transformed_shape_grads,
1229  offset,
1230  n_q_points);
1231  mapping.transform(
1232  make_array_view(fe_data.untransformed_shape_grads,
1233  offset,
1234  n_q_points),
1236  mapping_internal,
1237  transformed_shape_grads);
1238 
1239  for (unsigned int k = 0; k < n_q_points; ++k)
1240  for (unsigned int d = 0; d < spacedim; ++d)
1241  for (unsigned int n = 0; n < spacedim; ++n)
1242  transformed_shape_grads[k][d] -=
1243  output_data.shape_values(first + n, k) *
1244  mapping_data.jacobian_pushed_forward_grads[k][n][d];
1245 
1246  for (unsigned int k = 0; k < n_q_points; ++k)
1247  for (unsigned int d = 0; d < dim; ++d)
1248  output_data.shape_gradients[first + d][k] =
1249  fe_data.sign_change[i] * transformed_shape_grads[k][d];
1250 
1251  break;
1252  }
1253 
1254  default:
1255  Assert(false, ExcNotImplemented());
1256  }
1257  }
1258 
1259  if (fe_data.update_each & update_hessians)
1260  {
1261  switch (mapping_kind)
1262  {
1263  case mapping_none:
1264  {
1266  transformed_shape_hessians =
1267  make_array_view(fe_data.transformed_shape_hessians,
1268  offset,
1269  n_q_points);
1270  mapping.transform(make_array_view(fe_data.shape_grad_grads,
1271  i,
1272  offset,
1273  n_q_points),
1275  mapping_internal,
1276  transformed_shape_hessians);
1277 
1278  for (unsigned int k = 0; k < n_q_points; ++k)
1279  for (unsigned int d = 0; d < spacedim; ++d)
1280  for (unsigned int n = 0; n < spacedim; ++n)
1281  transformed_shape_hessians[k][d] -=
1282  output_data.shape_gradients[first + d][k][n] *
1283  mapping_data.jacobian_pushed_forward_grads[k][n];
1284 
1285  for (unsigned int k = 0; k < n_q_points; ++k)
1286  for (unsigned int d = 0; d < dim; ++d)
1287  output_data.shape_hessians[first + d][k] =
1288  transformed_shape_hessians[k][d];
1289 
1290  break;
1291  }
1292  case mapping_covariant:
1293  {
1294  for (unsigned int k = 0; k < n_q_points; ++k)
1295  fe_data.untransformed_shape_hessian_tensors[k + offset] =
1296  fe_data.shape_grad_grads[i][k + offset];
1297 
1299  transformed_shape_hessians =
1300  make_array_view(fe_data.transformed_shape_hessians,
1301  offset,
1302  n_q_points);
1303  mapping.transform(
1304  make_array_view(fe_data.untransformed_shape_hessian_tensors,
1305  offset,
1306  n_q_points),
1308  mapping_internal,
1309  transformed_shape_hessians);
1310 
1311  for (unsigned int k = 0; k < n_q_points; ++k)
1312  for (unsigned int d = 0; d < spacedim; ++d)
1313  for (unsigned int n = 0; n < spacedim; ++n)
1314  for (unsigned int i = 0; i < spacedim; ++i)
1315  for (unsigned int j = 0; j < spacedim; ++j)
1316  {
1317  transformed_shape_hessians[k][d][i][j] -=
1318  (output_data.shape_values(first + n, k) *
1319  mapping_data
1320  .jacobian_pushed_forward_2nd_derivatives
1321  [k][n][d][i][j]) +
1322  (output_data.shape_gradients[first + d][k][n] *
1323  mapping_data
1324  .jacobian_pushed_forward_grads[k][n][i][j]) +
1325  (output_data.shape_gradients[first + n][k][i] *
1326  mapping_data
1327  .jacobian_pushed_forward_grads[k][n][d][j]) +
1328  (output_data.shape_gradients[first + n][k][j] *
1329  mapping_data
1330  .jacobian_pushed_forward_grads[k][n][i][d]);
1331  }
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_hessians[first + d][k] =
1336  transformed_shape_hessians[k][d];
1337 
1338  break;
1339  }
1340 
1341  case mapping_contravariant:
1342  {
1343  for (unsigned int k = 0; k < n_q_points; ++k)
1344  fe_data.untransformed_shape_hessian_tensors[k + offset] =
1345  fe_data.shape_grad_grads[i][k + offset];
1346 
1348  transformed_shape_hessians =
1349  make_array_view(fe_data.transformed_shape_hessians,
1350  offset,
1351  n_q_points);
1352  mapping.transform(
1353  make_array_view(fe_data.untransformed_shape_hessian_tensors,
1354  offset,
1355  n_q_points),
1357  mapping_internal,
1358  transformed_shape_hessians);
1359 
1360  for (unsigned int k = 0; k < n_q_points; ++k)
1361  for (unsigned int d = 0; d < spacedim; ++d)
1362  for (unsigned int n = 0; n < spacedim; ++n)
1363  for (unsigned int i = 0; i < spacedim; ++i)
1364  for (unsigned int j = 0; j < spacedim; ++j)
1365  {
1366  transformed_shape_hessians[k][d][i][j] +=
1367  (output_data.shape_values(first + n, k) *
1368  mapping_data
1369  .jacobian_pushed_forward_2nd_derivatives
1370  [k][d][n][i][j]) +
1371  (output_data.shape_gradients[first + n][k][i] *
1372  mapping_data
1373  .jacobian_pushed_forward_grads[k][d][n][j]) +
1374  (output_data.shape_gradients[first + n][k][j] *
1375  mapping_data
1376  .jacobian_pushed_forward_grads[k][d][i][n]) -
1377  (output_data.shape_gradients[first + d][k][n] *
1378  mapping_data
1379  .jacobian_pushed_forward_grads[k][n][i][j]);
1380  for (unsigned int m = 0; m < spacedim; ++m)
1381  transformed_shape_hessians[k][d][i][j] -=
1382  (mapping_data
1383  .jacobian_pushed_forward_grads[k][d][i]
1384  [m] *
1385  mapping_data
1386  .jacobian_pushed_forward_grads[k][m][n]
1387  [j] *
1388  output_data.shape_values(first + n, k)) +
1389  (mapping_data
1390  .jacobian_pushed_forward_grads[k][d][m]
1391  [j] *
1392  mapping_data
1393  .jacobian_pushed_forward_grads[k][m][i]
1394  [n] *
1395  output_data.shape_values(first + n, k));
1396  }
1397 
1398  for (unsigned int k = 0; k < n_q_points; ++k)
1399  for (unsigned int d = 0; d < dim; ++d)
1400  output_data.shape_hessians[first + d][k] =
1401  transformed_shape_hessians[k][d];
1402 
1403  break;
1404  }
1405 
1407  case mapping_piola:
1408  {
1409  for (unsigned int k = 0; k < n_q_points; ++k)
1410  fe_data.untransformed_shape_hessian_tensors[k + offset] =
1411  fe_data.shape_grad_grads[i][k + offset];
1412 
1414  transformed_shape_hessians =
1415  make_array_view(fe_data.transformed_shape_hessians,
1416  offset,
1417  n_q_points);
1418  mapping.transform(
1419  make_array_view(fe_data.untransformed_shape_hessian_tensors,
1420  offset,
1421  n_q_points),
1423  mapping_internal,
1424  transformed_shape_hessians);
1425 
1426  for (unsigned int k = 0; k < n_q_points; ++k)
1427  for (unsigned int d = 0; d < spacedim; ++d)
1428  for (unsigned int n = 0; n < spacedim; ++n)
1429  for (unsigned int i = 0; i < spacedim; ++i)
1430  for (unsigned int j = 0; j < spacedim; ++j)
1431  {
1432  transformed_shape_hessians[k][d][i][j] +=
1433  (output_data.shape_values(first + n, k) *
1434  mapping_data
1435  .jacobian_pushed_forward_2nd_derivatives
1436  [k][d][n][i][j]) +
1437  (output_data.shape_gradients[first + n][k][i] *
1438  mapping_data
1439  .jacobian_pushed_forward_grads[k][d][n][j]) +
1440  (output_data.shape_gradients[first + n][k][j] *
1441  mapping_data
1442  .jacobian_pushed_forward_grads[k][d][i][n]) -
1443  (output_data.shape_gradients[first + d][k][n] *
1444  mapping_data
1445  .jacobian_pushed_forward_grads[k][n][i][j]);
1446 
1447  transformed_shape_hessians[k][d][i][j] -=
1448  (output_data.shape_values(first + d, k) *
1449  mapping_data
1450  .jacobian_pushed_forward_2nd_derivatives
1451  [k][n][n][i][j]) +
1452  (output_data.shape_gradients[first + d][k][i] *
1453  mapping_data
1454  .jacobian_pushed_forward_grads[k][n][n][j]) +
1455  (output_data.shape_gradients[first + d][k][j] *
1456  mapping_data
1457  .jacobian_pushed_forward_grads[k][n][n][i]);
1458 
1459  for (unsigned int m = 0; m < spacedim; ++m)
1460  {
1461  transformed_shape_hessians[k][d][i][j] -=
1462  (mapping_data
1463  .jacobian_pushed_forward_grads[k][d][i]
1464  [m] *
1465  mapping_data
1466  .jacobian_pushed_forward_grads[k][m][n]
1467  [j] *
1468  output_data.shape_values(first + n, k)) +
1469  (mapping_data
1470  .jacobian_pushed_forward_grads[k][d][m]
1471  [j] *
1472  mapping_data
1473  .jacobian_pushed_forward_grads[k][m][i]
1474  [n] *
1475  output_data.shape_values(first + n, k));
1476 
1477  transformed_shape_hessians[k][d][i][j] +=
1478  (mapping_data
1479  .jacobian_pushed_forward_grads[k][n][i]
1480  [m] *
1481  mapping_data
1482  .jacobian_pushed_forward_grads[k][m][n]
1483  [j] *
1484  output_data.shape_values(first + d, k)) +
1485  (mapping_data
1486  .jacobian_pushed_forward_grads[k][n][m]
1487  [j] *
1488  mapping_data
1489  .jacobian_pushed_forward_grads[k][m][i]
1490  [n] *
1491  output_data.shape_values(first + d, k));
1492  }
1493  }
1494 
1495  for (unsigned int k = 0; k < n_q_points; ++k)
1496  for (unsigned int d = 0; d < dim; ++d)
1497  output_data.shape_hessians[first + d][k] =
1498  fe_data.sign_change[i] *
1499  transformed_shape_hessians[k][d];
1500 
1501  break;
1502  }
1503 
1504  case mapping_nedelec:
1505  {
1506  for (unsigned int k = 0; k < n_q_points; ++k)
1507  fe_data.untransformed_shape_hessian_tensors[k + offset] =
1508  fe_data.shape_grad_grads[i][k + offset];
1509 
1511  transformed_shape_hessians =
1512  make_array_view(fe_data.transformed_shape_hessians,
1513  offset,
1514  n_q_points);
1515  mapping.transform(
1516  make_array_view(fe_data.untransformed_shape_hessian_tensors,
1517  offset,
1518  n_q_points),
1520  mapping_internal,
1521  transformed_shape_hessians);
1522 
1523  for (unsigned int k = 0; k < n_q_points; ++k)
1524  for (unsigned int d = 0; d < spacedim; ++d)
1525  for (unsigned int n = 0; n < spacedim; ++n)
1526  for (unsigned int i = 0; i < spacedim; ++i)
1527  for (unsigned int j = 0; j < spacedim; ++j)
1528  {
1529  transformed_shape_hessians[k][d][i][j] -=
1530  (output_data.shape_values(first + n, k) *
1531  mapping_data
1532  .jacobian_pushed_forward_2nd_derivatives
1533  [k][n][d][i][j]) +
1534  (output_data.shape_gradients[first + d][k][n] *
1535  mapping_data
1536  .jacobian_pushed_forward_grads[k][n][i][j]) +
1537  (output_data.shape_gradients[first + n][k][i] *
1538  mapping_data
1539  .jacobian_pushed_forward_grads[k][n][d][j]) +
1540  (output_data.shape_gradients[first + n][k][j] *
1541  mapping_data
1542  .jacobian_pushed_forward_grads[k][n][i][d]);
1543  }
1544 
1545  for (unsigned int k = 0; k < n_q_points; ++k)
1546  for (unsigned int d = 0; d < dim; ++d)
1547  output_data.shape_hessians[first + d][k] =
1548  fe_data.sign_change[i] *
1549  transformed_shape_hessians[k][d];
1550 
1551  break;
1552  }
1553 
1554  default:
1555  Assert(false, ExcNotImplemented());
1556  }
1557  }
1558 
1559  // third derivatives are not implemented
1560  if (fe_data.update_each & update_3rd_derivatives)
1561  {
1562  Assert(false, ExcNotImplemented())
1563  }
1564  }
1565 }
1566 
1567 
1568 
1569 template <int dim, int spacedim>
1570 void
1572  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1573  const unsigned int face_no,
1574  const unsigned int sub_no,
1575  const Quadrature<dim - 1> & quadrature,
1576  const Mapping<dim, spacedim> & mapping,
1577  const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1578  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
1579  spacedim>
1580  & mapping_data,
1581  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1583  spacedim>
1584  &output_data) const
1585 {
1586  // convert data object to internal
1587  // data for this class. fails with
1588  // an exception if that is not
1589  // possible
1590  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1591  ExcInternalError());
1592  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1593 
1594  const unsigned int n_q_points = quadrature.size();
1595 
1596  // offset determines which data set
1597  // to take (all data sets for all
1598  // sub-faces are stored contiguously)
1599  const auto offset =
1601  face_no,
1602  sub_no,
1603  cell->face_orientation(face_no),
1604  cell->face_flip(face_no),
1605  cell->face_rotation(face_no),
1606  n_q_points,
1607  cell->subface_case(face_no));
1608 
1609  // Assert(mapping_kind == independent
1610  // || ( mapping_kind == independent_on_cartesian
1611  // && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0),
1612  // ExcNotImplemented());
1613  // TODO: Size assertions
1614 
1615  // TODO: Sign change for the face DoFs!
1616 
1617  // Compute eventual sign changes depending
1618  // on the neighborhood between two faces.
1619  std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0);
1620 
1621  internal::FE_PolyTensor::get_face_sign_change_rt(cell,
1622  *this,
1623  this->mapping_kind,
1624  fe_data.sign_change);
1625 
1626  internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
1627  *this,
1628  this->mapping_kind,
1629  fe_data.sign_change);
1630 
1631  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
1632  {
1634 
1635  const unsigned int first =
1636  output_data.shape_function_to_row_table[i * this->n_components() +
1637  this->get_nonzero_components(i)
1639 
1640  if (fe_data.update_each & update_values)
1641  {
1642  switch (mapping_kind)
1643  {
1644  case mapping_none:
1645  {
1646  for (unsigned int k = 0; k < n_q_points; ++k)
1647  for (unsigned int d = 0; d < dim; ++d)
1648  output_data.shape_values(first + d, k) =
1649  fe_data.shape_values[i][k + offset][d];
1650  break;
1651  }
1652 
1653  case mapping_covariant:
1654  case mapping_contravariant:
1655  {
1657  transformed_shape_values =
1658  make_array_view(fe_data.transformed_shape_values,
1659  offset,
1660  n_q_points);
1661  mapping.transform(make_array_view(fe_data.shape_values,
1662  i,
1663  offset,
1664  n_q_points),
1665  mapping_kind,
1666  mapping_internal,
1667  transformed_shape_values);
1668 
1669  for (unsigned int k = 0; k < n_q_points; ++k)
1670  for (unsigned int d = 0; d < dim; ++d)
1671  output_data.shape_values(first + d, k) =
1672  transformed_shape_values[k][d];
1673 
1674  break;
1675  }
1676 
1678  case mapping_piola:
1679  {
1681  transformed_shape_values =
1682  make_array_view(fe_data.transformed_shape_values,
1683  offset,
1684  n_q_points);
1685 
1686  mapping.transform(make_array_view(fe_data.shape_values,
1687  i,
1688  offset,
1689  n_q_points),
1690  mapping_piola,
1691  mapping_internal,
1692  transformed_shape_values);
1693  for (unsigned int k = 0; k < n_q_points; ++k)
1694  for (unsigned int d = 0; d < dim; ++d)
1695  output_data.shape_values(first + d, k) =
1696  fe_data.sign_change[i] * transformed_shape_values[k][d];
1697  break;
1698  }
1699 
1700  case mapping_nedelec:
1701  {
1703  transformed_shape_values =
1704  make_array_view(fe_data.transformed_shape_values,
1705  offset,
1706  n_q_points);
1707 
1708  mapping.transform(make_array_view(fe_data.shape_values,
1709  i,
1710  offset,
1711  n_q_points),
1713  mapping_internal,
1714  transformed_shape_values);
1715 
1716  for (unsigned int k = 0; k < n_q_points; ++k)
1717  for (unsigned int d = 0; d < dim; ++d)
1718  output_data.shape_values(first + d, k) =
1719  fe_data.sign_change[i] * transformed_shape_values[k][d];
1720 
1721  break;
1722  }
1723 
1724  default:
1725  Assert(false, ExcNotImplemented());
1726  }
1727  }
1728 
1729  if (fe_data.update_each & update_gradients)
1730  {
1731  const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1732  make_array_view(fe_data.transformed_shape_grads,
1733  offset,
1734  n_q_points);
1735  switch (mapping_kind)
1736  {
1737  case mapping_none:
1738  {
1739  mapping.transform(
1740  make_array_view(fe_data.shape_grads, i, offset, n_q_points),
1742  mapping_internal,
1743  transformed_shape_grads);
1744  for (unsigned int k = 0; k < n_q_points; ++k)
1745  for (unsigned int d = 0; d < dim; ++d)
1746  output_data.shape_gradients[first + d][k] =
1747  transformed_shape_grads[k][d];
1748  break;
1749  }
1750 
1751  case mapping_covariant:
1752  {
1753  mapping.transform(
1754  make_array_view(fe_data.shape_grads, i, offset, n_q_points),
1756  mapping_internal,
1757  transformed_shape_grads);
1758 
1759  for (unsigned int k = 0; k < n_q_points; ++k)
1760  for (unsigned int d = 0; d < spacedim; ++d)
1761  for (unsigned int n = 0; n < spacedim; ++n)
1762  transformed_shape_grads[k][d] -=
1763  output_data.shape_values(first + n, k) *
1764  mapping_data.jacobian_pushed_forward_grads[k][n][d];
1765 
1766  for (unsigned int k = 0; k < n_q_points; ++k)
1767  for (unsigned int d = 0; d < dim; ++d)
1768  output_data.shape_gradients[first + d][k] =
1769  transformed_shape_grads[k][d];
1770 
1771  break;
1772  }
1773 
1774  case mapping_contravariant:
1775  {
1776  for (unsigned int k = 0; k < n_q_points; ++k)
1777  fe_data.untransformed_shape_grads[k + offset] =
1778  fe_data.shape_grads[i][k + offset];
1779 
1780  mapping.transform(
1781  make_array_view(fe_data.untransformed_shape_grads,
1782  offset,
1783  n_q_points),
1785  mapping_internal,
1786  transformed_shape_grads);
1787 
1788  for (unsigned int k = 0; k < n_q_points; ++k)
1789  for (unsigned int d = 0; d < spacedim; ++d)
1790  for (unsigned int n = 0; n < spacedim; ++n)
1791  transformed_shape_grads[k][d] +=
1792  output_data.shape_values(first + n, k) *
1793  mapping_data.jacobian_pushed_forward_grads[k][d][n];
1794 
1795  for (unsigned int k = 0; k < n_q_points; ++k)
1796  for (unsigned int d = 0; d < dim; ++d)
1797  output_data.shape_gradients[first + d][k] =
1798  transformed_shape_grads[k][d];
1799 
1800  break;
1801  }
1802 
1804  case mapping_piola:
1805  {
1806  for (unsigned int k = 0; k < n_q_points; ++k)
1807  fe_data.untransformed_shape_grads[k + offset] =
1808  fe_data.shape_grads[i][k + offset];
1809 
1810  mapping.transform(
1811  make_array_view(fe_data.untransformed_shape_grads,
1812  offset,
1813  n_q_points),
1815  mapping_internal,
1816  transformed_shape_grads);
1817 
1818  for (unsigned int k = 0; k < n_q_points; ++k)
1819  for (unsigned int d = 0; d < spacedim; ++d)
1820  for (unsigned int n = 0; n < spacedim; ++n)
1821  transformed_shape_grads[k][d] +=
1822  (output_data.shape_values(first + n, k) *
1823  mapping_data
1824  .jacobian_pushed_forward_grads[k][d][n]) -
1825  (output_data.shape_values(first + d, k) *
1826  mapping_data.jacobian_pushed_forward_grads[k][n][n]);
1827 
1828  for (unsigned int k = 0; k < n_q_points; ++k)
1829  for (unsigned int d = 0; d < dim; ++d)
1830  output_data.shape_gradients[first + d][k] =
1831  fe_data.sign_change[i] * transformed_shape_grads[k][d];
1832 
1833  break;
1834  }
1835 
1836  case mapping_nedelec:
1837  {
1838  // this particular shape
1839  // function at all
1840  // q-points. if Dv is the
1841  // gradient of the shape
1842  // function on the unit
1843  // cell, then
1844  // (J^-T)Dv(J^-1) is the
1845  // value we want to have on
1846  // the real cell.
1847  for (unsigned int k = 0; k < n_q_points; ++k)
1848  fe_data.untransformed_shape_grads[k + offset] =
1849  fe_data.shape_grads[i][k + offset];
1850 
1851  mapping.transform(
1852  make_array_view(fe_data.untransformed_shape_grads,
1853  offset,
1854  n_q_points),
1856  mapping_internal,
1857  transformed_shape_grads);
1858 
1859  for (unsigned int k = 0; k < n_q_points; ++k)
1860  for (unsigned int d = 0; d < spacedim; ++d)
1861  for (unsigned int n = 0; n < spacedim; ++n)
1862  transformed_shape_grads[k][d] -=
1863  output_data.shape_values(first + n, k) *
1864  mapping_data.jacobian_pushed_forward_grads[k][n][d];
1865 
1866  for (unsigned int k = 0; k < n_q_points; ++k)
1867  for (unsigned int d = 0; d < dim; ++d)
1868  output_data.shape_gradients[first + d][k] =
1869  fe_data.sign_change[i] * transformed_shape_grads[k][d];
1870 
1871  break;
1872  }
1873 
1874  default:
1875  Assert(false, ExcNotImplemented());
1876  }
1877  }
1878 
1879  if (fe_data.update_each & update_hessians)
1880  {
1881  switch (mapping_kind)
1882  {
1883  case mapping_none:
1884  {
1886  transformed_shape_hessians =
1887  make_array_view(fe_data.transformed_shape_hessians,
1888  offset,
1889  n_q_points);
1890  mapping.transform(make_array_view(fe_data.shape_grad_grads,
1891  i,
1892  offset,
1893  n_q_points),
1895  mapping_internal,
1896  transformed_shape_hessians);
1897 
1898  for (unsigned int k = 0; k < n_q_points; ++k)
1899  for (unsigned int d = 0; d < spacedim; ++d)
1900  for (unsigned int n = 0; n < spacedim; ++n)
1901  transformed_shape_hessians[k][d] -=
1902  output_data.shape_gradients[first + d][k][n] *
1903  mapping_data.jacobian_pushed_forward_grads[k][n];
1904 
1905  for (unsigned int k = 0; k < n_q_points; ++k)
1906  for (unsigned int d = 0; d < dim; ++d)
1907  output_data.shape_hessians[first + d][k] =
1908  transformed_shape_hessians[k][d];
1909 
1910  break;
1911  }
1912  case mapping_covariant:
1913  {
1914  for (unsigned int k = 0; k < n_q_points; ++k)
1915  fe_data.untransformed_shape_hessian_tensors[k + offset] =
1916  fe_data.shape_grad_grads[i][k + offset];
1917 
1919  transformed_shape_hessians =
1920  make_array_view(fe_data.transformed_shape_hessians,
1921  offset,
1922  n_q_points);
1923  mapping.transform(
1924  make_array_view(fe_data.untransformed_shape_hessian_tensors,
1925  offset,
1926  n_q_points),
1928  mapping_internal,
1929  transformed_shape_hessians);
1930 
1931  for (unsigned int k = 0; k < n_q_points; ++k)
1932  for (unsigned int d = 0; d < spacedim; ++d)
1933  for (unsigned int n = 0; n < spacedim; ++n)
1934  for (unsigned int i = 0; i < spacedim; ++i)
1935  for (unsigned int j = 0; j < spacedim; ++j)
1936  {
1937  transformed_shape_hessians[k][d][i][j] -=
1938  (output_data.shape_values(first + n, k) *
1939  mapping_data
1940  .jacobian_pushed_forward_2nd_derivatives
1941  [k][n][d][i][j]) +
1942  (output_data.shape_gradients[first + d][k][n] *
1943  mapping_data
1944  .jacobian_pushed_forward_grads[k][n][i][j]) +
1945  (output_data.shape_gradients[first + n][k][i] *
1946  mapping_data
1947  .jacobian_pushed_forward_grads[k][n][d][j]) +
1948  (output_data.shape_gradients[first + n][k][j] *
1949  mapping_data
1950  .jacobian_pushed_forward_grads[k][n][i][d]);
1951  }
1952 
1953  for (unsigned int k = 0; k < n_q_points; ++k)
1954  for (unsigned int d = 0; d < dim; ++d)
1955  output_data.shape_hessians[first + d][k] =
1956  transformed_shape_hessians[k][d];
1957 
1958  break;
1959  }
1960 
1961  case mapping_contravariant:
1962  {
1963  for (unsigned int k = 0; k < n_q_points; ++k)
1964  fe_data.untransformed_shape_hessian_tensors[k + offset] =
1965  fe_data.shape_grad_grads[i][k + offset];
1966 
1968  transformed_shape_hessians =
1969  make_array_view(fe_data.transformed_shape_hessians,
1970  offset,
1971  n_q_points);
1972  mapping.transform(
1973  make_array_view(fe_data.untransformed_shape_hessian_tensors,
1974  offset,
1975  n_q_points),
1977  mapping_internal,
1978  transformed_shape_hessians);
1979 
1980  for (unsigned int k = 0; k < n_q_points; ++k)
1981  for (unsigned int d = 0; d < spacedim; ++d)
1982  for (unsigned int n = 0; n < spacedim; ++n)
1983  for (unsigned int i = 0; i < spacedim; ++i)
1984  for (unsigned int j = 0; j < spacedim; ++j)
1985  {
1986  transformed_shape_hessians[k][d][i][j] +=
1987  (output_data.shape_values(first + n, k) *
1988  mapping_data
1989  .jacobian_pushed_forward_2nd_derivatives
1990  [k][d][n][i][j]) +
1991  (output_data.shape_gradients[first + n][k][i] *
1992  mapping_data
1993  .jacobian_pushed_forward_grads[k][d][n][j]) +
1994  (output_data.shape_gradients[first + n][k][j] *
1995  mapping_data
1996  .jacobian_pushed_forward_grads[k][d][i][n]) -
1997  (output_data.shape_gradients[first + d][k][n] *
1998  mapping_data
1999  .jacobian_pushed_forward_grads[k][n][i][j]);
2000  for (unsigned int m = 0; m < spacedim; ++m)
2001  transformed_shape_hessians[k][d][i][j] -=
2002  (mapping_data
2003  .jacobian_pushed_forward_grads[k][d][i]
2004  [m] *
2005  mapping_data
2006  .jacobian_pushed_forward_grads[k][m][n]
2007  [j] *
2008  output_data.shape_values(first + n, k)) +
2009  (mapping_data
2010  .jacobian_pushed_forward_grads[k][d][m]
2011  [j] *
2012  mapping_data
2013  .jacobian_pushed_forward_grads[k][m][i]
2014  [n] *
2015  output_data.shape_values(first + n, k));
2016  }
2017 
2018  for (unsigned int k = 0; k < n_q_points; ++k)
2019  for (unsigned int d = 0; d < dim; ++d)
2020  output_data.shape_hessians[first + d][k] =
2021  transformed_shape_hessians[k][d];
2022 
2023  break;
2024  }
2025 
2027  case mapping_piola:
2028  {
2029  for (unsigned int k = 0; k < n_q_points; ++k)
2030  fe_data.untransformed_shape_hessian_tensors[k + offset] =
2031  fe_data.shape_grad_grads[i][k + offset];
2032 
2034  transformed_shape_hessians =
2035  make_array_view(fe_data.transformed_shape_hessians,
2036  offset,
2037  n_q_points);
2038  mapping.transform(
2039  make_array_view(fe_data.untransformed_shape_hessian_tensors,
2040  offset,
2041  n_q_points),
2043  mapping_internal,
2044  transformed_shape_hessians);
2045 
2046  for (unsigned int k = 0; k < n_q_points; ++k)
2047  for (unsigned int d = 0; d < spacedim; ++d)
2048  for (unsigned int n = 0; n < spacedim; ++n)
2049  for (unsigned int i = 0; i < spacedim; ++i)
2050  for (unsigned int j = 0; j < spacedim; ++j)
2051  {
2052  transformed_shape_hessians[k][d][i][j] +=
2053  (output_data.shape_values(first + n, k) *
2054  mapping_data
2055  .jacobian_pushed_forward_2nd_derivatives
2056  [k][d][n][i][j]) +
2057  (output_data.shape_gradients[first + n][k][i] *
2058  mapping_data
2059  .jacobian_pushed_forward_grads[k][d][n][j]) +
2060  (output_data.shape_gradients[first + n][k][j] *
2061  mapping_data
2062  .jacobian_pushed_forward_grads[k][d][i][n]) -
2063  (output_data.shape_gradients[first + d][k][n] *
2064  mapping_data
2065  .jacobian_pushed_forward_grads[k][n][i][j]);
2066 
2067  transformed_shape_hessians[k][d][i][j] -=
2068  (output_data.shape_values(first + d, k) *
2069  mapping_data
2070  .jacobian_pushed_forward_2nd_derivatives
2071  [k][n][n][i][j]) +
2072  (output_data.shape_gradients[first + d][k][i] *
2073  mapping_data
2074  .jacobian_pushed_forward_grads[k][n][n][j]) +
2075  (output_data.shape_gradients[first + d][k][j] *
2076  mapping_data
2077  .jacobian_pushed_forward_grads[k][n][n][i]);
2078  for (unsigned int m = 0; m < spacedim; ++m)
2079  {
2080  transformed_shape_hessians[k][d][i][j] -=
2081  (mapping_data
2082  .jacobian_pushed_forward_grads[k][d][i]
2083  [m] *
2084  mapping_data
2085  .jacobian_pushed_forward_grads[k][m][n]
2086  [j] *
2087  output_data.shape_values(first + n, k)) +
2088  (mapping_data
2089  .jacobian_pushed_forward_grads[k][d][m]
2090  [j] *
2091  mapping_data
2092  .jacobian_pushed_forward_grads[k][m][i]
2093  [n] *
2094  output_data.shape_values(first + n, k));
2095 
2096  transformed_shape_hessians[k][d][i][j] +=
2097  (mapping_data
2098  .jacobian_pushed_forward_grads[k][n][i]
2099  [m] *
2100  mapping_data
2101  .jacobian_pushed_forward_grads[k][m][n]
2102  [j] *
2103  output_data.shape_values(first + d, k)) +
2104  (mapping_data
2105  .jacobian_pushed_forward_grads[k][n][m]
2106  [j] *
2107  mapping_data
2108  .jacobian_pushed_forward_grads[k][m][i]
2109  [n] *
2110  output_data.shape_values(first + d, k));
2111  }
2112  }
2113 
2114  for (unsigned int k = 0; k < n_q_points; ++k)
2115  for (unsigned int d = 0; d < dim; ++d)
2116  output_data.shape_hessians[first + d][k] =
2117  fe_data.sign_change[i] *
2118  transformed_shape_hessians[k][d];
2119 
2120  break;
2121  }
2122 
2123  case mapping_nedelec:
2124  {
2125  for (unsigned int k = 0; k < n_q_points; ++k)
2126  fe_data.untransformed_shape_hessian_tensors[k + offset] =
2127  fe_data.shape_grad_grads[i][k + offset];
2128 
2130  transformed_shape_hessians =
2131  make_array_view(fe_data.transformed_shape_hessians,
2132  offset,
2133  n_q_points);
2134  mapping.transform(
2135  make_array_view(fe_data.untransformed_shape_hessian_tensors,
2136  offset,
2137  n_q_points),
2139  mapping_internal,
2140  transformed_shape_hessians);
2141 
2142  for (unsigned int k = 0; k < n_q_points; ++k)
2143  for (unsigned int d = 0; d < spacedim; ++d)
2144  for (unsigned int n = 0; n < spacedim; ++n)
2145  for (unsigned int i = 0; i < spacedim; ++i)
2146  for (unsigned int j = 0; j < spacedim; ++j)
2147  {
2148  transformed_shape_hessians[k][d][i][j] -=
2149  (output_data.shape_values(first + n, k) *
2150  mapping_data
2151  .jacobian_pushed_forward_2nd_derivatives
2152  [k][n][d][i][j]) +
2153  (output_data.shape_gradients[first + d][k][n] *
2154  mapping_data
2155  .jacobian_pushed_forward_grads[k][n][i][j]) +
2156  (output_data.shape_gradients[first + n][k][i] *
2157  mapping_data
2158  .jacobian_pushed_forward_grads[k][n][d][j]) +
2159  (output_data.shape_gradients[first + n][k][j] *
2160  mapping_data
2161  .jacobian_pushed_forward_grads[k][n][i][d]);
2162  }
2163 
2164  for (unsigned int k = 0; k < n_q_points; ++k)
2165  for (unsigned int d = 0; d < dim; ++d)
2166  output_data.shape_hessians[first + d][k] =
2167  fe_data.sign_change[i] *
2168  transformed_shape_hessians[k][d];
2169 
2170  break;
2171  }
2172 
2173  default:
2174  Assert(false, ExcNotImplemented());
2175  }
2176  }
2177 
2178  // third derivatives are not implemented
2179  if (fe_data.update_each & update_3rd_derivatives)
2180  {
2181  Assert(false, ExcNotImplemented())
2182  }
2183  }
2184 }
2185 
2186 
2187 
2188 template <int dim, int spacedim>
2191  const UpdateFlags flags) const
2192 {
2194 
2195  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2196  {
2198 
2199  switch (mapping_kind)
2200  {
2201  case mapping_none:
2202  {
2203  if (flags & update_values)
2204  out |= update_values;
2205 
2206  if (flags & update_gradients)
2207  out |= update_gradients | update_values |
2209 
2210  if (flags & update_hessians)
2211  out |= update_hessians | update_values | update_gradients |
2214  break;
2215  }
2217  case mapping_piola:
2218  {
2219  if (flags & update_values)
2220  out |= update_values | update_piola;
2221 
2222  if (flags & update_gradients)
2223  out |= update_gradients | update_values | update_piola |
2227 
2228  if (flags & update_hessians)
2229  out |= update_hessians | update_piola | update_values |
2230  update_gradients | update_jacobian_pushed_forward_grads |
2233 
2234  break;
2235  }
2236 
2237 
2238  case mapping_contravariant:
2239  {
2240  if (flags & update_values)
2241  out |= update_values | update_piola;
2242 
2243  if (flags & update_gradients)
2244  out |= update_gradients | update_values |
2248 
2249  if (flags & update_hessians)
2250  out |= update_hessians | update_piola | update_values |
2251  update_gradients | update_jacobian_pushed_forward_grads |
2254 
2255  break;
2256  }
2257 
2258  case mapping_nedelec:
2259  case mapping_covariant:
2260  {
2261  if (flags & update_values)
2262  out |= update_values | update_covariant_transformation;
2263 
2264  if (flags & update_gradients)
2265  out |= update_gradients | update_values |
2268 
2269  if (flags & update_hessians)
2270  out |= update_hessians | update_values | update_gradients |
2274 
2275  break;
2276  }
2277 
2278  default:
2279  {
2280  Assert(false, ExcNotImplemented());
2281  }
2282  }
2283  }
2284 
2285  return out;
2286 }
2287 
2288 
2289 // explicit instantiations
2290 #include "fe_poly_tensor.inst"
2291 
2292 
Shape function values.
Contravariant transformation.
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
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
std::vector< Tensor< 2, dim > > cached_grads
MappingKind get_mapping_kind(const unsigned int i) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
bool single_mapping_kind() const
std::vector< Tensor< 1, dim > > cached_values
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition: fe.h:2562
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const =0
FullMatrix< double > inverse_node_matrix
MappingKind
Definition: mapping.h:62
ReferenceCell::Type reference_cell_type() const
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
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
No update.
Third derivatives of shape functions.
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
#define Assert(cond, exc)
Definition: exceptions.h:1411
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:301
std::mutex cache_mutex
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:3243
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:665
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
Point< dim > cached_point
std::vector< Tensor< 3, dim > > cached_grad_grads
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Definition: fe.cc:552
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
Second derivatives of shape functions.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int size() const
Point< 2 > first
Definition: grid_out.cc:4338
const std::unique_ptr< const TensorPolynomialsBase< dim > > poly_space
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
unsigned int n_components() const
unsigned int n_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
Shape function gradients.
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static ::ExceptionBase & ExcNotImplemented()
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_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(const TensorPolynomialsBase< dim > &polynomials, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Values needed for Piola transform.
Covariant transformation.
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
std::vector< MappingKind > mapping_kind
static DataSetDescriptor face(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:1184