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