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