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