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