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