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