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