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