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