Reference documentation for deal.II version GIT a3f17f8a20 2023-06-02 10:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mapping_fe_field.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
23 #include <deal.II/base/utilities.h>
24 
26 
27 #include <deal.II/fe/fe_q.h>
28 #include <deal.II/fe/fe_system.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 #include <deal.II/fe/mapping.h>
33 
35 
40 #include <deal.II/lac/la_vector.h>
47 #include <deal.II/lac/vector.h>
48 
50 
51 #include <fstream>
52 #include <memory>
53 #include <numeric>
54 
55 
56 
58 
59 
60 template <int dim, int spacedim, typename VectorType>
63  const ComponentMask & mask)
64  : unit_tangentials()
65  , n_shape_functions(fe.n_dofs_per_cell())
66  , mask(mask)
67  , local_dof_indices(fe.n_dofs_per_cell())
68  , local_dof_values(fe.n_dofs_per_cell())
69 {}
70 
71 
72 
73 template <int dim, int spacedim, typename VectorType>
74 std::size_t
76  const
77 {
78  Assert(false, ExcNotImplemented());
79  return 0;
80 }
81 
82 
83 
84 template <int dim, int spacedim, typename VectorType>
85 double &
87  const unsigned int qpoint,
88  const unsigned int shape_nr)
89 {
90  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
91  return shape_values[qpoint * n_shape_functions + shape_nr];
92 }
93 
94 
95 template <int dim, int spacedim, typename VectorType>
96 const Tensor<1, dim> &
98  const unsigned int qpoint,
99  const unsigned int shape_nr) const
100 {
101  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
102  shape_derivatives.size());
103  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
104 }
105 
106 
107 
108 template <int dim, int spacedim, typename VectorType>
111  const unsigned int qpoint,
112  const unsigned int shape_nr)
113 {
114  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
115  shape_derivatives.size());
116  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
117 }
118 
119 
120 template <int dim, int spacedim, typename VectorType>
121 const Tensor<2, dim> &
123  const unsigned int qpoint,
124  const unsigned int shape_nr) const
125 {
126  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
127  shape_second_derivatives.size());
128  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
129 }
130 
131 
132 
133 template <int dim, int spacedim, typename VectorType>
136  const unsigned int qpoint,
137  const unsigned int shape_nr)
138 {
139  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
140  shape_second_derivatives.size());
141  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
142 }
143 
144 
145 template <int dim, int spacedim, typename VectorType>
146 const Tensor<3, dim> &
148  const unsigned int qpoint,
149  const unsigned int shape_nr) const
150 {
151  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
152  shape_third_derivatives.size());
153  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
154 }
155 
156 
157 
158 template <int dim, int spacedim, typename VectorType>
161  const unsigned int qpoint,
162  const unsigned int shape_nr)
163 {
164  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
165  shape_third_derivatives.size());
166  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
167 }
168 
169 
170 template <int dim, int spacedim, typename VectorType>
171 const Tensor<4, dim> &
173  const unsigned int qpoint,
174  const unsigned int shape_nr) const
175 {
176  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
177  shape_fourth_derivatives.size());
178  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
179 }
180 
181 
182 
183 template <int dim, int spacedim, typename VectorType>
186  const unsigned int qpoint,
187  const unsigned int shape_nr)
188 {
189  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
190  shape_fourth_derivatives.size());
191  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
192 }
193 
194 
195 
196 template <int dim, int spacedim, typename VectorType>
199  const VectorType & euler_vector,
200  const ComponentMask & mask)
202  , uses_level_dofs(false)
204  , euler_dof_handler(&euler_dof_handler)
205  , fe_mask(mask.size() != 0u ?
206  mask :
208  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
209  true))
210  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
211  , fe_values(this->euler_dof_handler->get_fe(),
212  reference_cell.template get_nodal_type_quadrature<dim>(),
214 {
215  unsigned int size = 0;
216  for (unsigned int i = 0; i < fe_mask.size(); ++i)
217  {
218  if (fe_mask[i])
219  fe_to_real[i] = size++;
220  }
221  AssertDimension(size, spacedim);
222 }
223 
224 
225 
226 template <int dim, int spacedim, typename VectorType>
228  const DoFHandler<dim, spacedim> &euler_dof_handler,
229  const std::vector<VectorType> & euler_vector,
230  const ComponentMask & mask)
231  : reference_cell(euler_dof_handler.get_fe().reference_cell())
232  , uses_level_dofs(true)
233  , euler_dof_handler(&euler_dof_handler)
234  , fe_mask(mask.size() != 0u ?
235  mask :
237  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
238  true))
239  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
240  , fe_values(this->euler_dof_handler->get_fe(),
241  reference_cell.template get_nodal_type_quadrature<dim>(),
243 {
244  unsigned int size = 0;
245  for (unsigned int i = 0; i < fe_mask.size(); ++i)
246  {
247  if (fe_mask[i])
248  fe_to_real[i] = size++;
249  }
250  AssertDimension(size, spacedim);
251 
252  Assert(euler_dof_handler.has_level_dofs(),
253  ExcMessage("The underlying DoFHandler object did not call "
254  "distribute_mg_dofs(). In this case, the construction via "
255  "level vectors does not make sense."));
257  euler_dof_handler.get_triangulation().n_global_levels());
258  this->euler_vector.clear();
259  this->euler_vector.resize(euler_vector.size());
260  for (unsigned int i = 0; i < euler_vector.size(); ++i)
261  this->euler_vector[i] = &euler_vector[i];
262 }
263 
264 
265 
266 template <int dim, int spacedim, typename VectorType>
268  const DoFHandler<dim, spacedim> &euler_dof_handler,
269  const MGLevelObject<VectorType> &euler_vector,
270  const ComponentMask & mask)
271  : reference_cell(euler_dof_handler.get_fe().reference_cell())
272  , uses_level_dofs(true)
273  , euler_dof_handler(&euler_dof_handler)
274  , fe_mask(mask.size() != 0u ?
275  mask :
277  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
278  true))
279  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
280  , fe_values(this->euler_dof_handler->get_fe(),
281  reference_cell.template get_nodal_type_quadrature<dim>(),
283 {
284  unsigned int size = 0;
285  for (unsigned int i = 0; i < fe_mask.size(); ++i)
286  {
287  if (fe_mask[i])
288  fe_to_real[i] = size++;
289  }
290  AssertDimension(size, spacedim);
291 
292  Assert(euler_dof_handler.has_level_dofs(),
293  ExcMessage("The underlying DoFHandler object did not call "
294  "distribute_mg_dofs(). In this case, the construction via "
295  "level vectors does not make sense."));
296  AssertDimension(euler_vector.max_level() + 1,
297  euler_dof_handler.get_triangulation().n_global_levels());
298  this->euler_vector.clear();
299  this->euler_vector.resize(
300  euler_dof_handler.get_triangulation().n_global_levels());
301  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
302  ++i)
303  this->euler_vector[i] = &euler_vector[i];
304 }
305 
306 
307 
308 template <int dim, int spacedim, typename VectorType>
311  : reference_cell(mapping.reference_cell)
312  , uses_level_dofs(mapping.uses_level_dofs)
313  , euler_vector(mapping.euler_vector)
314  , euler_dof_handler(mapping.euler_dof_handler)
315  , fe_mask(mapping.fe_mask)
316  , fe_to_real(mapping.fe_to_real)
317  , fe_values(mapping.euler_dof_handler->get_fe(),
318  reference_cell.template get_nodal_type_quadrature<dim>(),
320 {}
321 
322 
323 
324 template <int dim, int spacedim, typename VectorType>
325 inline const double &
327  const unsigned int qpoint,
328  const unsigned int shape_nr) const
329 {
330  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
331  return shape_values[qpoint * n_shape_functions + shape_nr];
332 }
333 
334 
335 
336 template <int dim, int spacedim, typename VectorType>
337 bool
339 {
340  return false;
341 }
342 
343 
344 
345 template <int dim, int spacedim, typename VectorType>
346 bool
348  const ReferenceCell &reference_cell) const
349 {
351  ExcMessage("The dimension of your mapping (" +
352  Utilities::to_string(dim) +
353  ") and the reference cell cell_type (" +
355  " ) do not agree."));
356 
357  return this->reference_cell == reference_cell;
358 }
359 
360 
361 
362 template <int dim, int spacedim, typename VectorType>
363 boost::container::small_vector<Point<spacedim>,
366  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
367 {
368  // we transform our tria iterator into a dof iterator so we can access
369  // data not associated with triangulations
370  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
371  *cell, euler_dof_handler);
372 
373  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
374  AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
376  euler_dof_handler->get_fe().n_components());
377  if (uses_level_dofs)
378  {
379  AssertIndexRange(cell->level(), euler_vector.size());
380  AssertDimension(euler_vector[cell->level()]->size(),
381  euler_dof_handler->n_dofs(cell->level()));
382  }
383  else
384  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
385 
386  {
387  std::lock_guard<std::mutex> lock(fe_values_mutex);
388  fe_values.reinit(dof_cell);
389  }
390  const unsigned int dofs_per_cell =
391  euler_dof_handler->get_fe().n_dofs_per_cell();
392  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
393  if (uses_level_dofs)
394  dof_cell->get_mg_dof_indices(dof_indices);
395  else
396  dof_cell->get_dof_indices(dof_indices);
397 
398  const VectorType &vector =
399  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
400 
401  boost::container::small_vector<Point<spacedim>,
403  vertices(cell->n_vertices());
404  for (unsigned int i = 0; i < dofs_per_cell; ++i)
405  {
406  const unsigned int comp = fe_to_real
407  [euler_dof_handler->get_fe().system_to_component_index(i).first];
408  if (comp != numbers::invalid_unsigned_int)
409  {
410  typename VectorType::value_type value =
411  internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
412  if (euler_dof_handler->get_fe().is_primitive(i))
413  for (const unsigned int v : cell->vertex_indices())
414  vertices[v][comp] += fe_values.shape_value(i, v) * value;
415  else
416  Assert(false, ExcNotImplemented());
417  }
418  }
419 
420  return vertices;
421 }
422 
423 
424 
425 template <int dim, int spacedim, typename VectorType>
426 void
428  const std::vector<Point<dim>> & unit_points,
430 {
431  const auto fe = &euler_dof_handler->get_fe();
432  const unsigned int n_points = unit_points.size();
433 
434  for (unsigned int point = 0; point < n_points; ++point)
435  {
436  if (data.shape_values.size() != 0)
437  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
438  data.shape(point, i) = fe->shape_value(i, unit_points[point]);
439 
440  if (data.shape_derivatives.size() != 0)
441  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
442  data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
443 
444  if (data.shape_second_derivatives.size() != 0)
445  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
446  data.second_derivative(point, i) =
447  fe->shape_grad_grad(i, unit_points[point]);
448 
449  if (data.shape_third_derivatives.size() != 0)
450  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
451  data.third_derivative(point, i) =
452  fe->shape_3rd_derivative(i, unit_points[point]);
453 
454  if (data.shape_fourth_derivatives.size() != 0)
455  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
456  data.fourth_derivative(point, i) =
457  fe->shape_4th_derivative(i, unit_points[point]);
458  }
459 }
460 
461 
462 template <int dim, int spacedim, typename VectorType>
465  const UpdateFlags in) const
466 {
467  // add flags if the respective quantities are necessary to compute
468  // what we need. note that some flags appear in both conditions and
469  // in subsequent set operations. this leads to some circular
470  // logic. the only way to treat this is to iterate. since there are
471  // 5 if-clauses in the loop, it will take at most 4 iterations to
472  // converge. do them:
473  UpdateFlags out = in;
474  for (unsigned int i = 0; i < 5; ++i)
475  {
476  // The following is a little incorrect:
477  // If not applied on a face,
478  // update_boundary_forms does not
479  // make sense. On the other hand,
480  // it is necessary on a
481  // face. Currently,
482  // update_boundary_forms is simply
483  // ignored for the interior of a
484  // cell.
486  out |= update_boundary_forms;
487 
488  if (out &
492 
493  if (out &
498 
499  // The contravariant transformation is used in the Piola
500  // transformation, which requires the determinant of the Jacobi
501  // matrix of the transformation. Because we have no way of
502  // knowing here whether the finite element wants to use the
503  // contravariant or the Piola transforms, we add the volume elements
504  // to the list of flags to be updated for each cell.
506  out |= update_volume_elements;
507 
508  if (out & update_normal_vectors)
509  out |= update_volume_elements;
510  }
511 
512  return out;
513 }
514 
515 
516 
517 template <int dim, int spacedim, typename VectorType>
518 void
520  const UpdateFlags update_flags,
521  const Quadrature<dim> &q,
522  const unsigned int n_original_q_points,
523  InternalData & data) const
524 {
525  // store the flags in the internal data object so we can access them
526  // in fill_fe_*_values(). use the transitive hull of the required
527  // flags
528  data.update_each = requires_update_flags(update_flags);
529 
530  const unsigned int n_q_points = q.size();
531 
532  // see if we need the (transformation) shape function values
533  // and/or gradients and resize the necessary arrays
534  if (data.update_each & update_quadrature_points)
535  data.shape_values.resize(data.n_shape_functions * n_q_points);
536 
537  if (data.update_each &
541  data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
542 
543  if (data.update_each & update_covariant_transformation)
544  data.covariant.resize(n_original_q_points);
545 
546  if (data.update_each & update_contravariant_transformation)
547  data.contravariant.resize(n_original_q_points);
548 
549  if (data.update_each & update_volume_elements)
550  data.volume_elements.resize(n_original_q_points);
551 
552  if (data.update_each &
554  data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
555 
556  if (data.update_each & (update_jacobian_2nd_derivatives |
558  data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
559 
560  if (data.update_each & (update_jacobian_3rd_derivatives |
562  data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
563 
565 
566  // This (for face values and simplices) can be different for different calls,
567  // so always copy
568  data.quadrature_weights = q.get_weights();
569 }
570 
571 
572 template <int dim, int spacedim, typename VectorType>
573 void
575  const UpdateFlags update_flags,
576  const Quadrature<dim> &q,
577  const unsigned int n_original_q_points,
578  InternalData & data) const
579 {
580  compute_data(update_flags, q, n_original_q_points, data);
581 
582  if (dim > 1)
583  {
584  if (data.update_each & update_boundary_forms)
585  {
586  data.aux.resize(
587  dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
588 
589 
590  // TODO: only a single reference cell type possible...
591  const auto n_faces = reference_cell.n_faces();
592 
593  // Compute tangentials to the unit cell.
594  for (unsigned int i = 0; i < n_faces; ++i)
595  {
596  data.unit_tangentials[i].resize(n_original_q_points);
597  std::fill(
598  data.unit_tangentials[i].begin(),
599  data.unit_tangentials[i].end(),
600  reference_cell.template unit_tangential_vectors<dim>(i, 0));
601  if (dim > 2)
602  {
603  data.unit_tangentials[n_faces + i].resize(
604  n_original_q_points);
605  std::fill(
606  data.unit_tangentials[n_faces + i].begin(),
607  data.unit_tangentials[n_faces + i].end(),
608  reference_cell.template unit_tangential_vectors<dim>(i, 1));
609  }
610  }
611  }
612  }
613 }
614 
615 
616 template <int dim, int spacedim, typename VectorType>
617 typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
619  const UpdateFlags update_flags,
620  const Quadrature<dim> &quadrature) const
621 {
622  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
623  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
624  auto &data = dynamic_cast<InternalData &>(*data_ptr);
625  this->compute_data(update_flags, quadrature, quadrature.size(), data);
626 
627  return data_ptr;
628 }
629 
630 
631 
632 template <int dim, int spacedim, typename VectorType>
633 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
635  const UpdateFlags update_flags,
636  const hp::QCollection<dim - 1> &quadrature) const
637 {
638  AssertDimension(quadrature.size(), 1);
639 
640  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
641  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
642  auto & data = dynamic_cast<InternalData &>(*data_ptr);
643  const Quadrature<dim> q(
645  this->compute_face_data(update_flags, q, quadrature[0].size(), data);
646 
647  return data_ptr;
648 }
649 
650 
651 template <int dim, int spacedim, typename VectorType>
652 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
654  const UpdateFlags update_flags,
655  const Quadrature<dim - 1> &quadrature) const
656 {
657  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
658  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
659  auto & data = dynamic_cast<InternalData &>(*data_ptr);
660  const Quadrature<dim> q(
662  this->compute_face_data(update_flags, q, quadrature.size(), data);
663 
664  return data_ptr;
665 }
666 
667 
668 
669 namespace internal
670 {
671  namespace MappingFEFieldImplementation
672  {
673  namespace
674  {
681  template <int dim, int spacedim, typename VectorType>
682  void
683  maybe_compute_q_points(
684  const typename ::QProjector<dim>::DataSetDescriptor data_set,
685  const typename ::MappingFEField<dim, spacedim, VectorType>::
686  InternalData & data,
688  const ComponentMask & fe_mask,
689  const std::vector<unsigned int> & fe_to_real,
690  std::vector<Point<spacedim>> & quadrature_points)
691  {
692  const UpdateFlags update_flags = data.update_each;
693 
694  if (update_flags & update_quadrature_points)
695  {
696  for (unsigned int point = 0; point < quadrature_points.size();
697  ++point)
698  {
699  Point<spacedim> result;
700  const double * shape = &data.shape(point + data_set, 0);
701 
702  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
703  {
704  const unsigned int comp_k =
705  fe.system_to_component_index(k).first;
706  if (fe_mask[comp_k])
707  result[fe_to_real[comp_k]] +=
708  data.local_dof_values[k] * shape[k];
709  }
710 
711  quadrature_points[point] = result;
712  }
713  }
714  }
715 
723  template <int dim, int spacedim, typename VectorType>
724  void
725  maybe_update_Jacobians(
726  const typename ::QProjector<dim>::DataSetDescriptor data_set,
727  const typename ::MappingFEField<dim, spacedim, VectorType>::
728  InternalData & data,
730  const ComponentMask & fe_mask,
731  const std::vector<unsigned int> & fe_to_real)
732  {
733  const UpdateFlags update_flags = data.update_each;
734 
735  // then Jacobians
736  if (update_flags & update_contravariant_transformation)
737  {
738  const unsigned int n_q_points = data.contravariant.size();
739 
740  Assert(data.n_shape_functions > 0, ExcInternalError());
741 
742  for (unsigned int point = 0; point < n_q_points; ++point)
743  {
744  const Tensor<1, dim> *data_derv =
745  &data.derivative(point + data_set, 0);
746 
747  Tensor<1, dim> result[spacedim];
748 
749  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
750  {
751  const unsigned int comp_k =
752  fe.system_to_component_index(k).first;
753  if (fe_mask[comp_k])
754  result[fe_to_real[comp_k]] +=
755  data.local_dof_values[k] * data_derv[k];
756  }
757 
758  // write result into contravariant data
759  for (unsigned int i = 0; i < spacedim; ++i)
760  {
761  data.contravariant[point][i] = result[i];
762  }
763  }
764  }
765 
766  if (update_flags & update_covariant_transformation)
767  {
768  AssertDimension(data.covariant.size(), data.contravariant.size());
769  for (unsigned int point = 0; point < data.contravariant.size();
770  ++point)
771  data.covariant[point] =
772  (data.contravariant[point]).covariant_form();
773  }
774 
775  if (update_flags & update_volume_elements)
776  {
777  AssertDimension(data.contravariant.size(),
778  data.volume_elements.size());
779  for (unsigned int point = 0; point < data.contravariant.size();
780  ++point)
781  data.volume_elements[point] =
782  data.contravariant[point].determinant();
783  }
784  }
785 
792  template <int dim, int spacedim, typename VectorType>
793  void
795  const typename ::QProjector<dim>::DataSetDescriptor data_set,
796  const typename ::MappingFEField<dim, spacedim, VectorType>::
797  InternalData & data,
798  const FiniteElement<dim, spacedim> & fe,
799  const ComponentMask & fe_mask,
800  const std::vector<unsigned int> & fe_to_real,
801  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
802  {
803  const UpdateFlags update_flags = data.update_each;
804  if (update_flags & update_jacobian_grads)
805  {
806  const unsigned int n_q_points = jacobian_grads.size();
807 
808  for (unsigned int point = 0; point < n_q_points; ++point)
809  {
810  const Tensor<2, dim> *second =
811  &data.second_derivative(point + data_set, 0);
812 
814 
815  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
816  {
817  const unsigned int comp_k =
818  fe.system_to_component_index(k).first;
819  if (fe_mask[comp_k])
820  for (unsigned int j = 0; j < dim; ++j)
821  for (unsigned int l = 0; l < dim; ++l)
822  result[fe_to_real[comp_k]][j][l] +=
823  (second[k][j][l] * data.local_dof_values[k]);
824  }
825 
826  // never touch any data for j=dim in case dim<spacedim, so
827  // it will always be zero as it was initialized
828  for (unsigned int i = 0; i < spacedim; ++i)
829  for (unsigned int j = 0; j < dim; ++j)
830  for (unsigned int l = 0; l < dim; ++l)
831  jacobian_grads[point][i][j][l] = result[i][j][l];
832  }
833  }
834  }
835 
842  template <int dim, int spacedim, typename VectorType>
843  void
845  const typename ::QProjector<dim>::DataSetDescriptor data_set,
846  const typename ::MappingFEField<dim, spacedim, VectorType>::
847  InternalData & data,
849  const ComponentMask & fe_mask,
850  const std::vector<unsigned int> & fe_to_real,
851  std::vector<Tensor<3, spacedim>> & jacobian_pushed_forward_grads)
852  {
853  const UpdateFlags update_flags = data.update_each;
854  if (update_flags & update_jacobian_pushed_forward_grads)
855  {
856  const unsigned int n_q_points =
857  jacobian_pushed_forward_grads.size();
858 
859  double tmp[spacedim][spacedim][spacedim];
860  for (unsigned int point = 0; point < n_q_points; ++point)
861  {
862  const Tensor<2, dim> *second =
863  &data.second_derivative(point + data_set, 0);
864 
866 
867  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
868  {
869  const unsigned int comp_k =
870  fe.system_to_component_index(k).first;
871  if (fe_mask[comp_k])
872  for (unsigned int j = 0; j < dim; ++j)
873  for (unsigned int l = 0; l < dim; ++l)
874  result[fe_to_real[comp_k]][j][l] +=
875  (second[k][j][l] * data.local_dof_values[k]);
876  }
877 
878  // first push forward the j-components
879  for (unsigned int i = 0; i < spacedim; ++i)
880  for (unsigned int j = 0; j < spacedim; ++j)
881  for (unsigned int l = 0; l < dim; ++l)
882  {
883  tmp[i][j][l] =
884  result[i][0][l] * data.covariant[point][j][0];
885  for (unsigned int jr = 1; jr < dim; ++jr)
886  {
887  tmp[i][j][l] +=
888  result[i][jr][l] * data.covariant[point][j][jr];
889  }
890  }
891 
892  // now, pushing forward the l-components
893  for (unsigned int i = 0; i < spacedim; ++i)
894  for (unsigned int j = 0; j < spacedim; ++j)
895  for (unsigned int l = 0; l < spacedim; ++l)
896  {
897  jacobian_pushed_forward_grads[point][i][j][l] =
898  tmp[i][j][0] * data.covariant[point][l][0];
899  for (unsigned int lr = 1; lr < dim; ++lr)
900  {
901  jacobian_pushed_forward_grads[point][i][j][l] +=
902  tmp[i][j][lr] * data.covariant[point][l][lr];
903  }
904  }
905  }
906  }
907  }
908 
915  template <int dim, int spacedim, typename VectorType>
916  void
918  const typename ::QProjector<dim>::DataSetDescriptor data_set,
919  const typename ::MappingFEField<dim, spacedim, VectorType>::
920  InternalData & data,
921  const FiniteElement<dim, spacedim> & fe,
922  const ComponentMask & fe_mask,
923  const std::vector<unsigned int> & fe_to_real,
924  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
925  {
926  const UpdateFlags update_flags = data.update_each;
927  if (update_flags & update_jacobian_2nd_derivatives)
928  {
929  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
930 
931  for (unsigned int point = 0; point < n_q_points; ++point)
932  {
933  const Tensor<3, dim> *third =
934  &data.third_derivative(point + data_set, 0);
935 
937 
938  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
939  {
940  const unsigned int comp_k =
941  fe.system_to_component_index(k).first;
942  if (fe_mask[comp_k])
943  for (unsigned int j = 0; j < dim; ++j)
944  for (unsigned int l = 0; l < dim; ++l)
945  for (unsigned int m = 0; m < dim; ++m)
946  result[fe_to_real[comp_k]][j][l][m] +=
947  (third[k][j][l][m] * data.local_dof_values[k]);
948  }
949 
950  // never touch any data for j=dim in case dim<spacedim, so
951  // it will always be zero as it was initialized
952  for (unsigned int i = 0; i < spacedim; ++i)
953  for (unsigned int j = 0; j < dim; ++j)
954  for (unsigned int l = 0; l < dim; ++l)
955  for (unsigned int m = 0; m < dim; ++m)
956  jacobian_2nd_derivatives[point][i][j][l][m] =
957  result[i][j][l][m];
958  }
959  }
960  }
961 
969  template <int dim, int spacedim, typename VectorType>
970  void
972  const typename ::QProjector<dim>::DataSetDescriptor data_set,
973  const typename ::MappingFEField<dim, spacedim, VectorType>::
974  InternalData & data,
976  const ComponentMask & fe_mask,
977  const std::vector<unsigned int> & fe_to_real,
978  std::vector<Tensor<4, spacedim>>
979  &jacobian_pushed_forward_2nd_derivatives)
980  {
981  const UpdateFlags update_flags = data.update_each;
983  {
984  const unsigned int n_q_points =
985  jacobian_pushed_forward_2nd_derivatives.size();
986 
987  double tmp[spacedim][spacedim][spacedim][spacedim];
988  for (unsigned int point = 0; point < n_q_points; ++point)
989  {
990  const Tensor<3, dim> *third =
991  &data.third_derivative(point + data_set, 0);
992 
994 
995  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
996  {
997  const unsigned int comp_k =
998  fe.system_to_component_index(k).first;
999  if (fe_mask[comp_k])
1000  for (unsigned int j = 0; j < dim; ++j)
1001  for (unsigned int l = 0; l < dim; ++l)
1002  for (unsigned int m = 0; m < dim; ++m)
1003  result[fe_to_real[comp_k]][j][l][m] +=
1004  (third[k][j][l][m] * data.local_dof_values[k]);
1005  }
1006 
1007  // push forward the j-coordinate
1008  for (unsigned int i = 0; i < spacedim; ++i)
1009  for (unsigned int j = 0; j < spacedim; ++j)
1010  for (unsigned int l = 0; l < dim; ++l)
1011  for (unsigned int m = 0; m < dim; ++m)
1012  {
1013  jacobian_pushed_forward_2nd_derivatives
1014  [point][i][j][l][m] =
1015  result[i][0][l][m] * data.covariant[point][j][0];
1016  for (unsigned int jr = 1; jr < dim; ++jr)
1017  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1018  [l][m] +=
1019  result[i][jr][l][m] *
1020  data.covariant[point][j][jr];
1021  }
1022 
1023  // push forward the l-coordinate
1024  for (unsigned int i = 0; i < spacedim; ++i)
1025  for (unsigned int j = 0; j < spacedim; ++j)
1026  for (unsigned int l = 0; l < spacedim; ++l)
1027  for (unsigned int m = 0; m < dim; ++m)
1028  {
1029  tmp[i][j][l][m] =
1030  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1031  [0][m] *
1032  data.covariant[point][l][0];
1033  for (unsigned int lr = 1; lr < dim; ++lr)
1034  tmp[i][j][l][m] +=
1035  jacobian_pushed_forward_2nd_derivatives[point][i]
1036  [j][lr]
1037  [m] *
1038  data.covariant[point][l][lr];
1039  }
1040 
1041  // push forward the m-coordinate
1042  for (unsigned int i = 0; i < spacedim; ++i)
1043  for (unsigned int j = 0; j < spacedim; ++j)
1044  for (unsigned int l = 0; l < spacedim; ++l)
1045  for (unsigned int m = 0; m < spacedim; ++m)
1046  {
1047  jacobian_pushed_forward_2nd_derivatives
1048  [point][i][j][l][m] =
1049  tmp[i][j][l][0] * data.covariant[point][m][0];
1050  for (unsigned int mr = 1; mr < dim; ++mr)
1051  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1052  [l][m] +=
1053  tmp[i][j][l][mr] * data.covariant[point][m][mr];
1054  }
1055  }
1056  }
1057  }
1058 
1065  template <int dim, int spacedim, typename VectorType>
1066  void
1068  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1069  const typename ::MappingFEField<dim, spacedim, VectorType>::
1070  InternalData & data,
1071  const FiniteElement<dim, spacedim> & fe,
1072  const ComponentMask & fe_mask,
1073  const std::vector<unsigned int> & fe_to_real,
1074  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1075  {
1076  const UpdateFlags update_flags = data.update_each;
1077  if (update_flags & update_jacobian_3rd_derivatives)
1078  {
1079  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1080 
1081  for (unsigned int point = 0; point < n_q_points; ++point)
1082  {
1083  const Tensor<4, dim> *fourth =
1084  &data.fourth_derivative(point + data_set, 0);
1085 
1087 
1088  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1089  {
1090  const unsigned int comp_k =
1091  fe.system_to_component_index(k).first;
1092  if (fe_mask[comp_k])
1093  for (unsigned int j = 0; j < dim; ++j)
1094  for (unsigned int l = 0; l < dim; ++l)
1095  for (unsigned int m = 0; m < dim; ++m)
1096  for (unsigned int n = 0; n < dim; ++n)
1097  result[fe_to_real[comp_k]][j][l][m][n] +=
1098  (fourth[k][j][l][m][n] *
1099  data.local_dof_values[k]);
1100  }
1101 
1102  // never touch any data for j,l,m,n=dim in case
1103  // dim<spacedim, so it will always be zero as it was
1104  // initialized
1105  for (unsigned int i = 0; i < spacedim; ++i)
1106  for (unsigned int j = 0; j < dim; ++j)
1107  for (unsigned int l = 0; l < dim; ++l)
1108  for (unsigned int m = 0; m < dim; ++m)
1109  for (unsigned int n = 0; n < dim; ++n)
1110  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1111  result[i][j][l][m][n];
1112  }
1113  }
1114  }
1115 
1123  template <int dim, int spacedim, typename VectorType>
1124  void
1126  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1127  const typename ::MappingFEField<dim, spacedim, VectorType>::
1128  InternalData & data,
1129  const FiniteElement<dim, spacedim> &fe,
1130  const ComponentMask & fe_mask,
1131  const std::vector<unsigned int> & fe_to_real,
1132  std::vector<Tensor<5, spacedim>>
1133  &jacobian_pushed_forward_3rd_derivatives)
1134  {
1135  const UpdateFlags update_flags = data.update_each;
1137  {
1138  const unsigned int n_q_points =
1139  jacobian_pushed_forward_3rd_derivatives.size();
1140 
1141  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1142  for (unsigned int point = 0; point < n_q_points; ++point)
1143  {
1144  const Tensor<4, dim> *fourth =
1145  &data.fourth_derivative(point + data_set, 0);
1146 
1148 
1149  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1150  {
1151  const unsigned int comp_k =
1152  fe.system_to_component_index(k).first;
1153  if (fe_mask[comp_k])
1154  for (unsigned int j = 0; j < dim; ++j)
1155  for (unsigned int l = 0; l < dim; ++l)
1156  for (unsigned int m = 0; m < dim; ++m)
1157  for (unsigned int n = 0; n < dim; ++n)
1158  result[fe_to_real[comp_k]][j][l][m][n] +=
1159  (fourth[k][j][l][m][n] *
1160  data.local_dof_values[k]);
1161  }
1162 
1163  // push-forward the j-coordinate
1164  for (unsigned int i = 0; i < spacedim; ++i)
1165  for (unsigned int j = 0; j < spacedim; ++j)
1166  for (unsigned int l = 0; l < dim; ++l)
1167  for (unsigned int m = 0; m < dim; ++m)
1168  for (unsigned int n = 0; n < dim; ++n)
1169  {
1170  tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1171  data.covariant[point][j][0];
1172  for (unsigned int jr = 1; jr < dim; ++jr)
1173  tmp[i][j][l][m][n] +=
1174  result[i][jr][l][m][n] *
1175  data.covariant[point][j][jr];
1176  }
1177 
1178  // push-forward the l-coordinate
1179  for (unsigned int i = 0; i < spacedim; ++i)
1180  for (unsigned int j = 0; j < spacedim; ++j)
1181  for (unsigned int l = 0; l < spacedim; ++l)
1182  for (unsigned int m = 0; m < dim; ++m)
1183  for (unsigned int n = 0; n < dim; ++n)
1184  {
1185  jacobian_pushed_forward_3rd_derivatives
1186  [point][i][j][l][m][n] =
1187  tmp[i][j][0][m][n] *
1188  data.covariant[point][l][0];
1189  for (unsigned int lr = 1; lr < dim; ++lr)
1190  jacobian_pushed_forward_3rd_derivatives[point][i]
1191  [j][l][m]
1192  [n] +=
1193  tmp[i][j][lr][m][n] *
1194  data.covariant[point][l][lr];
1195  }
1196 
1197  // push-forward the m-coordinate
1198  for (unsigned int i = 0; i < spacedim; ++i)
1199  for (unsigned int j = 0; j < spacedim; ++j)
1200  for (unsigned int l = 0; l < spacedim; ++l)
1201  for (unsigned int m = 0; m < spacedim; ++m)
1202  for (unsigned int n = 0; n < dim; ++n)
1203  {
1204  tmp[i][j][l][m][n] =
1205  jacobian_pushed_forward_3rd_derivatives[point][i]
1206  [j][l][0]
1207  [n] *
1208  data.covariant[point][m][0];
1209  for (unsigned int mr = 1; mr < dim; ++mr)
1210  tmp[i][j][l][m][n] +=
1211  jacobian_pushed_forward_3rd_derivatives[point]
1212  [i][j][l]
1213  [mr][n] *
1214  data.covariant[point][m][mr];
1215  }
1216 
1217  // push-forward the n-coordinate
1218  for (unsigned int i = 0; i < spacedim; ++i)
1219  for (unsigned int j = 0; j < spacedim; ++j)
1220  for (unsigned int l = 0; l < spacedim; ++l)
1221  for (unsigned int m = 0; m < spacedim; ++m)
1222  for (unsigned int n = 0; n < spacedim; ++n)
1223  {
1224  jacobian_pushed_forward_3rd_derivatives
1225  [point][i][j][l][m][n] =
1226  tmp[i][j][l][m][0] *
1227  data.covariant[point][n][0];
1228  for (unsigned int nr = 1; nr < dim; ++nr)
1229  jacobian_pushed_forward_3rd_derivatives[point][i]
1230  [j][l][m]
1231  [n] +=
1232  tmp[i][j][l][m][nr] *
1233  data.covariant[point][n][nr];
1234  }
1235  }
1236  }
1237  }
1238 
1239 
1249  template <int dim, int spacedim, typename VectorType>
1250  void
1252  const ::Mapping<dim, spacedim> &mapping,
1253  const typename ::Triangulation<dim, spacedim>::cell_iterator
1254  & cell,
1255  const unsigned int face_no,
1256  const unsigned int subface_no,
1257  const typename QProjector<dim>::DataSetDescriptor data_set,
1258  const typename ::MappingFEField<dim, spacedim, VectorType>::
1259  InternalData &data,
1261  &output_data)
1262  {
1263  const UpdateFlags update_flags = data.update_each;
1264 
1265  if (update_flags & update_boundary_forms)
1266  {
1267  const unsigned int n_q_points = output_data.boundary_forms.size();
1268  if (update_flags & update_normal_vectors)
1269  AssertDimension(output_data.normal_vectors.size(), n_q_points);
1270  if (update_flags & update_JxW_values)
1271  AssertDimension(output_data.JxW_values.size(), n_q_points);
1272 
1273  // map the unit tangentials to the real cell. checking for d!=dim-1
1274  // eliminates compiler warnings regarding unsigned int expressions <
1275  // 0.
1276  for (unsigned int d = 0; d != dim - 1; ++d)
1277  {
1278  Assert(face_no + cell->n_faces() * d <
1279  data.unit_tangentials.size(),
1280  ExcInternalError());
1281  Assert(
1282  data.aux[d].size() <=
1283  data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1284  ExcInternalError());
1285 
1286  mapping.transform(
1288  data.unit_tangentials[face_no + cell->n_faces() * d]),
1290  data,
1291  make_array_view(data.aux[d]));
1292  }
1293 
1294  // if dim==spacedim, we can use the unit tangentials to compute the
1295  // boundary form by simply taking the cross product
1296  if (dim == spacedim)
1297  {
1298  for (unsigned int i = 0; i < n_q_points; ++i)
1299  switch (dim)
1300  {
1301  case 1:
1302  // in 1d, we don't have access to any of the data.aux
1303  // fields (because it has only dim-1 components), but we
1304  // can still compute the boundary form by simply looking
1305  // at the number of the face
1306  output_data.boundary_forms[i][0] =
1307  (face_no == 0 ? -1 : +1);
1308  break;
1309  case 2:
1310  output_data.boundary_forms[i] =
1311  cross_product_2d(data.aux[0][i]);
1312  break;
1313  case 3:
1314  output_data.boundary_forms[i] =
1315  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1316  break;
1317  default:
1318  Assert(false, ExcNotImplemented());
1319  }
1320  }
1321  else //(dim < spacedim)
1322  {
1323  // in the codim-one case, the boundary form results from the
1324  // cross product of all the face tangential vectors and the cell
1325  // normal vector
1326  //
1327  // to compute the cell normal, use the same method used in
1328  // fill_fe_values for cells above
1329  AssertDimension(data.contravariant.size(), n_q_points);
1330 
1331  for (unsigned int point = 0; point < n_q_points; ++point)
1332  {
1333  if (dim == 1)
1334  {
1335  // J is a tangent vector
1336  output_data.boundary_forms[point] =
1337  data.contravariant[point].transpose()[0];
1338  output_data.boundary_forms[point] /=
1339  (face_no == 0 ? -1. : +1.) *
1340  output_data.boundary_forms[point].norm();
1341  }
1342 
1343  if (dim == 2)
1344  {
1346  data.contravariant[point].transpose();
1347 
1348  Tensor<1, spacedim> cell_normal =
1349  cross_product_3d(DX_t[0], DX_t[1]);
1350  cell_normal /= cell_normal.norm();
1351 
1352  // then compute the face normal from the face tangent
1353  // and the cell normal:
1354  output_data.boundary_forms[point] =
1355  cross_product_3d(data.aux[0][point], cell_normal);
1356  }
1357  }
1358  }
1359 
1360  if (update_flags & (update_normal_vectors | update_JxW_values))
1361  for (unsigned int i = 0; i < output_data.boundary_forms.size();
1362  ++i)
1363  {
1364  if (update_flags & update_JxW_values)
1365  {
1366  output_data.JxW_values[i] =
1367  output_data.boundary_forms[i].norm() *
1368  data.quadrature_weights[i + data_set];
1369 
1370  if (subface_no != numbers::invalid_unsigned_int)
1371  {
1372  // TODO
1373  const double area_ratio =
1375  cell->subface_case(face_no), subface_no);
1376  output_data.JxW_values[i] *= area_ratio;
1377  }
1378  }
1379 
1380  if (update_flags & update_normal_vectors)
1381  output_data.normal_vectors[i] =
1382  Point<spacedim>(output_data.boundary_forms[i] /
1383  output_data.boundary_forms[i].norm());
1384  }
1385  }
1386  }
1387 
1394  template <int dim, int spacedim, typename VectorType>
1395  void
1397  const ::Mapping<dim, spacedim> &mapping,
1398  const typename ::Triangulation<dim, spacedim>::cell_iterator
1399  & cell,
1400  const unsigned int face_no,
1401  const unsigned int subface_no,
1402  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1403  const typename ::MappingFEField<dim, spacedim, VectorType>::
1404  InternalData & data,
1405  const FiniteElement<dim, spacedim> &fe,
1406  const ComponentMask & fe_mask,
1407  const std::vector<unsigned int> & fe_to_real,
1409  &output_data)
1410  {
1411  maybe_compute_q_points<dim, spacedim, VectorType>(
1412  data_set,
1413  data,
1414  fe,
1415  fe_mask,
1416  fe_to_real,
1417  output_data.quadrature_points);
1418 
1419  maybe_update_Jacobians<dim, spacedim, VectorType>(
1420  data_set, data, fe, fe_mask, fe_to_real);
1421 
1422  const UpdateFlags update_flags = data.update_each;
1423  const unsigned int n_q_points = data.contravariant.size();
1424 
1425  if (update_flags & update_jacobians)
1426  for (unsigned int point = 0; point < n_q_points; ++point)
1427  output_data.jacobians[point] = data.contravariant[point];
1428 
1429  if (update_flags & update_inverse_jacobians)
1430  for (unsigned int point = 0; point < n_q_points; ++point)
1431  output_data.inverse_jacobians[point] =
1432  data.covariant[point].transpose();
1433 
1434  maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1435  data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1436 
1437  maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1438  data_set,
1439  data,
1440  fe,
1441  fe_mask,
1442  fe_to_real,
1443  output_data.jacobian_pushed_forward_grads);
1444 
1445  maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1446  data_set,
1447  data,
1448  fe,
1449  fe_mask,
1450  fe_to_real,
1451  output_data.jacobian_2nd_derivatives);
1452 
1454  spacedim,
1455  VectorType>(
1456  data_set,
1457  data,
1458  fe,
1459  fe_mask,
1460  fe_to_real,
1462 
1463  maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1464  data_set,
1465  data,
1466  fe,
1467  fe_mask,
1468  fe_to_real,
1469  output_data.jacobian_3rd_derivatives);
1470 
1472  spacedim,
1473  VectorType>(
1474  data_set,
1475  data,
1476  fe,
1477  fe_mask,
1478  fe_to_real,
1480 
1481  maybe_compute_face_data<dim, spacedim, VectorType>(
1482  mapping, cell, face_no, subface_no, data_set, data, output_data);
1483  }
1484  } // namespace
1485  } // namespace MappingFEFieldImplementation
1486 } // namespace internal
1487 
1488 
1489 // Note that the CellSimilarity flag is modifiable, since MappingFEField can
1490 // need to recalculate data even when cells are similar.
1491 template <int dim, int spacedim, typename VectorType>
1494  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1496  const Quadrature<dim> & quadrature,
1497  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1499  &output_data) const
1500 {
1501  // convert data object to internal data for this class. fails with an
1502  // exception if that is not possible
1503  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1504  ExcInternalError());
1505  const InternalData &data = static_cast<const InternalData &>(internal_data);
1506 
1507  const unsigned int n_q_points = quadrature.size();
1508 
1509  update_internal_dofs(cell, data);
1510 
1511  internal::MappingFEFieldImplementation::
1512  maybe_compute_q_points<dim, spacedim, VectorType>(
1514  data,
1515  euler_dof_handler->get_fe(),
1516  fe_mask,
1517  fe_to_real,
1518  output_data.quadrature_points);
1519 
1520  internal::MappingFEFieldImplementation::
1521  maybe_update_Jacobians<dim, spacedim, VectorType>(
1523  data,
1524  euler_dof_handler->get_fe(),
1525  fe_mask,
1526  fe_to_real);
1527 
1528  const UpdateFlags update_flags = data.update_each;
1529  const std::vector<double> &weights = quadrature.get_weights();
1530 
1531  // Multiply quadrature weights by absolute value of Jacobian determinants or
1532  // the area element g=sqrt(DX^t DX) in case of codim > 0
1533 
1534  if (update_flags & (update_normal_vectors | update_JxW_values))
1535  {
1536  AssertDimension(output_data.JxW_values.size(), n_q_points);
1537 
1538  Assert(!(update_flags & update_normal_vectors) ||
1539  (output_data.normal_vectors.size() == n_q_points),
1540  ExcDimensionMismatch(output_data.normal_vectors.size(),
1541  n_q_points));
1542 
1543 
1544  for (unsigned int point = 0; point < n_q_points; ++point)
1545  {
1546  if (dim == spacedim)
1547  {
1548  const double det = data.volume_elements[point];
1549 
1550  // check for distorted cells.
1551 
1552  // TODO: this allows for anisotropies of up to 1e6 in 3d and
1553  // 1e12 in 2d. might want to find a finer
1554  // (dimension-independent) criterion
1555  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1556  cell->diameter() / std::sqrt(double(dim))),
1558  cell->center(), det, point)));
1559  output_data.JxW_values[point] = weights[point] * det;
1560  }
1561  // if dim==spacedim, then there is no cell normal to
1562  // compute. since this is for FEValues (and not FEFaceValues),
1563  // there are also no face normals to compute
1564  else // codim>0 case
1565  {
1566  Tensor<1, spacedim> DX_t[dim];
1567  for (unsigned int i = 0; i < spacedim; ++i)
1568  for (unsigned int j = 0; j < dim; ++j)
1569  DX_t[j][i] = data.contravariant[point][i][j];
1570 
1571  Tensor<2, dim> G; // First fundamental form
1572  for (unsigned int i = 0; i < dim; ++i)
1573  for (unsigned int j = 0; j < dim; ++j)
1574  G[i][j] = DX_t[i] * DX_t[j];
1575 
1576  output_data.JxW_values[point] =
1577  std::sqrt(determinant(G)) * weights[point];
1578 
1579  if (update_flags & update_normal_vectors)
1580  {
1581  Assert(spacedim - dim == 1,
1582  ExcMessage("There is no cell normal in codim 2."));
1583 
1584  if (dim == 1)
1585  output_data.normal_vectors[point] =
1586  cross_product_2d(-DX_t[0]);
1587  else
1588  {
1589  Assert(dim == 2, ExcInternalError());
1590 
1591  // dim-1==1 for the second argument, but this
1592  // avoids a compiler warning about array bounds:
1593  output_data.normal_vectors[point] =
1594  cross_product_3d(DX_t[0], DX_t[dim - 1]);
1595  }
1596 
1597  output_data.normal_vectors[point] /=
1598  output_data.normal_vectors[point].norm();
1599 
1600  if (cell->direction_flag() == false)
1601  output_data.normal_vectors[point] *= -1.;
1602  }
1603  } // codim>0 case
1604  }
1605  }
1606 
1607  // copy values from InternalData to vector given by reference
1608  if (update_flags & update_jacobians)
1609  {
1610  AssertDimension(output_data.jacobians.size(), n_q_points);
1611  for (unsigned int point = 0; point < n_q_points; ++point)
1612  output_data.jacobians[point] = data.contravariant[point];
1613  }
1614 
1615  // copy values from InternalData to vector given by reference
1616  if (update_flags & update_inverse_jacobians)
1617  {
1618  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1619  for (unsigned int point = 0; point < n_q_points; ++point)
1620  output_data.inverse_jacobians[point] =
1621  data.covariant[point].transpose();
1622  }
1623 
1624  // calculate derivatives of the Jacobians
1625  internal::MappingFEFieldImplementation::
1626  maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1628  data,
1629  euler_dof_handler->get_fe(),
1630  fe_mask,
1631  fe_to_real,
1632  output_data.jacobian_grads);
1633 
1634  // calculate derivatives of the Jacobians pushed forward to real cell
1635  // coordinates
1636  internal::MappingFEFieldImplementation::
1637  maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1639  data,
1640  euler_dof_handler->get_fe(),
1641  fe_mask,
1642  fe_to_real,
1643  output_data.jacobian_pushed_forward_grads);
1644 
1645  // calculate hessians of the Jacobians
1646  internal::MappingFEFieldImplementation::
1647  maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1649  data,
1650  euler_dof_handler->get_fe(),
1651  fe_mask,
1652  fe_to_real,
1653  output_data.jacobian_2nd_derivatives);
1654 
1655  // calculate hessians of the Jacobians pushed forward to real cell coordinates
1656  internal::MappingFEFieldImplementation::
1657  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1658  spacedim,
1659  VectorType>(
1661  data,
1662  euler_dof_handler->get_fe(),
1663  fe_mask,
1664  fe_to_real,
1666 
1667  // calculate gradients of the hessians of the Jacobians
1668  internal::MappingFEFieldImplementation::
1669  maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1671  data,
1672  euler_dof_handler->get_fe(),
1673  fe_mask,
1674  fe_to_real,
1675  output_data.jacobian_3rd_derivatives);
1676 
1677  // calculate gradients of the hessians of the Jacobians pushed forward to real
1678  // cell coordinates
1679  internal::MappingFEFieldImplementation::
1680  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1681  spacedim,
1682  VectorType>(
1684  data,
1685  euler_dof_handler->get_fe(),
1686  fe_mask,
1687  fe_to_real,
1689 
1691 }
1692 
1693 
1694 
1695 template <int dim, int spacedim, typename VectorType>
1696 void
1698  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1699  const unsigned int face_no,
1700  const hp::QCollection<dim - 1> & quadrature,
1701  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1703  &output_data) const
1704 {
1705  AssertDimension(quadrature.size(), 1);
1706 
1707  // convert data object to internal data for this class. fails with an
1708  // exception if that is not possible
1709  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1710  ExcInternalError());
1711  const InternalData &data = static_cast<const InternalData &>(internal_data);
1712 
1713  update_internal_dofs(cell, data);
1714 
1715  internal::MappingFEFieldImplementation::
1716  do_fill_fe_face_values<dim, spacedim, VectorType>(
1717  *this,
1718  cell,
1719  face_no,
1722  face_no,
1723  cell->face_orientation(face_no),
1724  cell->face_flip(face_no),
1725  cell->face_rotation(face_no),
1726  quadrature[0].size()),
1727  data,
1728  euler_dof_handler->get_fe(),
1729  fe_mask,
1730  fe_to_real,
1731  output_data);
1732 }
1733 
1734 
1735 template <int dim, int spacedim, typename VectorType>
1736 void
1738  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1739  const unsigned int face_no,
1740  const unsigned int subface_no,
1741  const Quadrature<dim - 1> & quadrature,
1742  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1744  &output_data) const
1745 {
1746  // convert data object to internal data for this class. fails with an
1747  // exception if that is not possible
1748  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1749  ExcInternalError());
1750  const InternalData &data = static_cast<const InternalData &>(internal_data);
1751 
1752  update_internal_dofs(cell, data);
1753 
1755  spacedim,
1756  VectorType>(
1757  *this,
1758  cell,
1759  face_no,
1762  face_no,
1763  subface_no,
1764  cell->face_orientation(face_no),
1765  cell->face_flip(face_no),
1766  cell->face_rotation(face_no),
1767  quadrature.size(),
1768  cell->subface_case(face_no)),
1769  data,
1770  euler_dof_handler->get_fe(),
1771  fe_mask,
1772  fe_to_real,
1773  output_data);
1774 }
1775 
1776 
1777 
1778 template <int dim, int spacedim, typename VectorType>
1779 void
1781  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1783  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1785  &output_data) const
1786 {
1787  AssertDimension(dim, spacedim);
1788  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1789  ExcInternalError());
1790  const InternalData &data = static_cast<const InternalData &>(internal_data);
1791 
1792  const unsigned int n_q_points = quadrature.size();
1793 
1794  update_internal_dofs(cell, data);
1795 
1796  internal::MappingFEFieldImplementation::
1797  maybe_compute_q_points<dim, spacedim, VectorType>(
1799  data,
1800  euler_dof_handler->get_fe(),
1801  fe_mask,
1802  fe_to_real,
1803  output_data.quadrature_points);
1804 
1805  internal::MappingFEFieldImplementation::
1806  maybe_update_Jacobians<dim, spacedim, VectorType>(
1808  data,
1809  euler_dof_handler->get_fe(),
1810  fe_mask,
1811  fe_to_real);
1812 
1813  const UpdateFlags update_flags = data.update_each;
1814  const std::vector<double> &weights = quadrature.get_weights();
1815 
1816  if (update_flags & (update_normal_vectors | update_JxW_values))
1817  {
1818  AssertDimension(output_data.JxW_values.size(), n_q_points);
1819 
1820  Assert(!(update_flags & update_normal_vectors) ||
1821  (output_data.normal_vectors.size() == n_q_points),
1822  ExcDimensionMismatch(output_data.normal_vectors.size(),
1823  n_q_points));
1824 
1825 
1826  for (unsigned int point = 0; point < n_q_points; ++point)
1827  {
1828  const double det = data.volume_elements[point];
1829 
1830  // check for distorted cells.
1831 
1832  // TODO: this allows for anisotropies of up to 1e6 in 3d and
1833  // 1e12 in 2d. might want to find a finer
1834  // (dimension-independent) criterion
1835  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1836  cell->diameter() / std::sqrt(double(dim))),
1838  cell->center(), det, point)));
1839 
1840  // The normals are n = J^{-T} * \hat{n} before normalizing.
1841  Tensor<1, spacedim> normal;
1842  for (unsigned int d = 0; d < spacedim; d++)
1843  normal[d] =
1844  data.covariant[point][d] * quadrature.normal_vector(point);
1845 
1846  output_data.JxW_values[point] = weights[point] * det * normal.norm();
1847 
1848  if ((update_flags & update_normal_vectors) != 0u)
1849  {
1850  normal /= normal.norm();
1851  output_data.normal_vectors[point] = normal;
1852  }
1853  }
1854 
1855  // copy values from InternalData to vector given by reference
1856  if (update_flags & update_jacobians)
1857  {
1858  AssertDimension(output_data.jacobians.size(), n_q_points);
1859  for (unsigned int point = 0; point < n_q_points; ++point)
1860  output_data.jacobians[point] = data.contravariant[point];
1861  }
1862 
1863  // copy values from InternalData to vector given by reference
1864  if (update_flags & update_inverse_jacobians)
1865  {
1866  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1867  for (unsigned int point = 0; point < n_q_points; ++point)
1868  output_data.inverse_jacobians[point] =
1869  data.covariant[point].transpose();
1870  }
1871 
1872  // calculate derivatives of the Jacobians
1873  internal::MappingFEFieldImplementation::
1874  maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1876  data,
1877  euler_dof_handler->get_fe(),
1878  fe_mask,
1879  fe_to_real,
1880  output_data.jacobian_grads);
1881 
1882  // calculate derivatives of the Jacobians pushed forward to real cell
1883  // coordinates
1884  internal::MappingFEFieldImplementation::
1885  maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1887  data,
1888  euler_dof_handler->get_fe(),
1889  fe_mask,
1890  fe_to_real,
1891  output_data.jacobian_pushed_forward_grads);
1892 
1893  // calculate hessians of the Jacobians
1894  internal::MappingFEFieldImplementation::
1895  maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1897  data,
1898  euler_dof_handler->get_fe(),
1899  fe_mask,
1900  fe_to_real,
1901  output_data.jacobian_2nd_derivatives);
1902 
1903  // calculate hessians of the Jacobians pushed forward to real cell
1904  // coordinates
1905  internal::MappingFEFieldImplementation::
1906  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1907  spacedim,
1908  VectorType>(
1910  data,
1911  euler_dof_handler->get_fe(),
1912  fe_mask,
1913  fe_to_real,
1915 
1916  // calculate gradients of the hessians of the Jacobians
1917  internal::MappingFEFieldImplementation::
1918  maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1920  data,
1921  euler_dof_handler->get_fe(),
1922  fe_mask,
1923  fe_to_real,
1924  output_data.jacobian_3rd_derivatives);
1925 
1926  // calculate gradients of the hessians of the Jacobians pushed forward to
1927  // real cell coordinates
1928  internal::MappingFEFieldImplementation::
1929  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1930  spacedim,
1931  VectorType>(
1933  data,
1934  euler_dof_handler->get_fe(),
1935  fe_mask,
1936  fe_to_real,
1938  }
1939 }
1940 
1941 namespace internal
1942 {
1943  namespace MappingFEFieldImplementation
1944  {
1945  namespace
1946  {
1947  template <int dim, int spacedim, int rank, typename VectorType>
1948  void
1950  const ArrayView<const Tensor<rank, dim>> & input,
1951  const MappingKind mapping_kind,
1952  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1953  const ArrayView<Tensor<rank, spacedim>> & output)
1954  {
1955  AssertDimension(input.size(), output.size());
1956  Assert((dynamic_cast<
1957  const typename ::
1958  MappingFEField<dim, spacedim, VectorType>::InternalData *>(
1959  &mapping_data) != nullptr),
1960  ExcInternalError());
1961  const typename ::MappingFEField<dim, spacedim, VectorType>::
1962  InternalData &data = static_cast<
1963  const typename ::MappingFEField<dim, spacedim, VectorType>::
1964  InternalData &>(mapping_data);
1965 
1966  switch (mapping_kind)
1967  {
1968  case mapping_contravariant:
1969  {
1970  Assert(
1971  data.update_each & update_contravariant_transformation,
1973  "update_contravariant_transformation"));
1974 
1975  for (unsigned int i = 0; i < output.size(); ++i)
1976  output[i] =
1977  apply_transformation(data.contravariant[i], input[i]);
1978 
1979  return;
1980  }
1981 
1982  case mapping_piola:
1983  {
1984  Assert(
1985  data.update_each & update_contravariant_transformation,
1987  "update_contravariant_transformation"));
1988  Assert(
1989  data.update_each & update_volume_elements,
1991  "update_volume_elements"));
1992  Assert(rank == 1, ExcMessage("Only for rank 1"));
1993  for (unsigned int i = 0; i < output.size(); ++i)
1994  {
1995  output[i] =
1996  apply_transformation(data.contravariant[i], input[i]);
1997  output[i] /= data.volume_elements[i];
1998  }
1999  return;
2000  }
2001 
2002 
2003  // We still allow this operation as in the
2004  // reference cell Derivatives are Tensor
2005  // rather than DerivativeForm
2006  case mapping_covariant:
2007  {
2008  Assert(
2009  data.update_each & update_contravariant_transformation,
2011  "update_contravariant_transformation"));
2012 
2013  for (unsigned int i = 0; i < output.size(); ++i)
2014  output[i] = apply_transformation(data.covariant[i], input[i]);
2015 
2016  return;
2017  }
2018 
2019  default:
2020  Assert(false, ExcNotImplemented());
2021  }
2022  }
2023 
2024 
2025  template <int dim, int spacedim, int rank, typename VectorType>
2026  void
2028  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
2029  const MappingKind mapping_kind,
2030  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2031  const ArrayView<Tensor<rank + 1, spacedim>> & output)
2032  {
2033  AssertDimension(input.size(), output.size());
2034  Assert((dynamic_cast<
2035  const typename ::
2036  MappingFEField<dim, spacedim, VectorType>::InternalData *>(
2037  &mapping_data) != nullptr),
2038  ExcInternalError());
2039  const typename ::MappingFEField<dim, spacedim, VectorType>::
2040  InternalData &data = static_cast<
2041  const typename ::MappingFEField<dim, spacedim, VectorType>::
2042  InternalData &>(mapping_data);
2043 
2044  switch (mapping_kind)
2045  {
2046  case mapping_covariant:
2047  {
2048  Assert(
2049  data.update_each & update_contravariant_transformation,
2051  "update_contravariant_transformation"));
2052 
2053  for (unsigned int i = 0; i < output.size(); ++i)
2054  output[i] = apply_transformation(data.covariant[i], input[i]);
2055 
2056  return;
2057  }
2058  default:
2059  Assert(false, ExcNotImplemented());
2060  }
2061  }
2062  } // namespace
2063  } // namespace MappingFEFieldImplementation
2064 } // namespace internal
2065 
2066 
2067 
2068 template <int dim, int spacedim, typename VectorType>
2069 void
2071  const ArrayView<const Tensor<1, dim>> & input,
2072  const MappingKind mapping_kind,
2073  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2074  const ArrayView<Tensor<1, spacedim>> & output) const
2075 {
2076  AssertDimension(input.size(), output.size());
2077 
2078  internal::MappingFEFieldImplementation::
2079  transform_fields<dim, spacedim, 1, VectorType>(input,
2080  mapping_kind,
2081  mapping_data,
2082  output);
2083 }
2084 
2085 
2086 
2087 template <int dim, int spacedim, typename VectorType>
2088 void
2090  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2091  const MappingKind mapping_kind,
2092  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2093  const ArrayView<Tensor<2, spacedim>> & output) const
2094 {
2095  AssertDimension(input.size(), output.size());
2096 
2097  internal::MappingFEFieldImplementation::
2098  transform_differential_forms<dim, spacedim, 1, VectorType>(input,
2099  mapping_kind,
2100  mapping_data,
2101  output);
2102 }
2103 
2104 
2105 
2106 template <int dim, int spacedim, typename VectorType>
2107 void
2109  const ArrayView<const Tensor<2, dim>> &input,
2110  const MappingKind,
2111  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2112  const ArrayView<Tensor<2, spacedim>> & output) const
2113 {
2114  (void)input;
2115  (void)output;
2116  (void)mapping_data;
2117  AssertDimension(input.size(), output.size());
2118 
2119  AssertThrow(false, ExcNotImplemented());
2120 }
2121 
2122 
2123 
2124 template <int dim, int spacedim, typename VectorType>
2125 void
2127  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2128  const MappingKind mapping_kind,
2129  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2130  const ArrayView<Tensor<3, spacedim>> & output) const
2131 {
2132  AssertDimension(input.size(), output.size());
2133  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2134  ExcInternalError());
2135  const InternalData &data = static_cast<const InternalData &>(mapping_data);
2136 
2137  switch (mapping_kind)
2138  {
2140  {
2141  Assert(data.update_each & update_contravariant_transformation,
2143  "update_covariant_transformation"));
2144 
2145  for (unsigned int q = 0; q < output.size(); ++q)
2146  for (unsigned int i = 0; i < spacedim; ++i)
2147  for (unsigned int j = 0; j < spacedim; ++j)
2148  for (unsigned int k = 0; k < spacedim; ++k)
2149  {
2150  output[q][i][j][k] = data.covariant[q][j][0] *
2151  data.covariant[q][k][0] *
2152  input[q][i][0][0];
2153  for (unsigned int J = 0; J < dim; ++J)
2154  {
2155  const unsigned int K0 = (0 == J) ? 1 : 0;
2156  for (unsigned int K = K0; K < dim; ++K)
2157  output[q][i][j][k] += data.covariant[q][j][J] *
2158  data.covariant[q][k][K] *
2159  input[q][i][J][K];
2160  }
2161  }
2162  return;
2163  }
2164 
2165  default:
2166  Assert(false, ExcNotImplemented());
2167  }
2168 }
2169 
2170 
2171 
2172 template <int dim, int spacedim, typename VectorType>
2173 void
2175  const ArrayView<const Tensor<3, dim>> &input,
2176  const MappingKind /*mapping_kind*/,
2177  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2178  const ArrayView<Tensor<3, spacedim>> & output) const
2179 {
2180  (void)input;
2181  (void)output;
2182  (void)mapping_data;
2183  AssertDimension(input.size(), output.size());
2184 
2185  AssertThrow(false, ExcNotImplemented());
2186 }
2187 
2188 
2189 
2190 template <int dim, int spacedim, typename VectorType>
2193  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2194  const Point<dim> & p) const
2195 {
2196  // Use the get_data function to create an InternalData with data vectors of
2197  // the right size and transformation shape values already computed at point
2198  // p.
2199  const Quadrature<dim> point_quadrature(p);
2200  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2201  get_data(update_quadrature_points | update_jacobians, point_quadrature));
2202  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2203  ExcInternalError());
2204 
2205  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2206 
2207  return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
2208 }
2209 
2210 
2211 template <int dim, int spacedim, typename VectorType>
2214  const InternalData &data) const
2215 {
2216  Point<spacedim> p_real;
2217 
2218  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2219  {
2220  unsigned int comp_i =
2221  euler_dof_handler->get_fe().system_to_component_index(i).first;
2222  if (fe_mask[comp_i])
2223  p_real[fe_to_real[comp_i]] +=
2224  data.local_dof_values[i] * data.shape(0, i);
2225  }
2226 
2227  return p_real;
2228 }
2229 
2230 
2231 
2232 template <int dim, int spacedim, typename VectorType>
2233 Point<dim>
2235  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2236  const Point<spacedim> & p) const
2237 {
2238  // first a Newton iteration based on the real mapping. It uses the center
2239  // point of the cell as a starting point
2240  Point<dim> initial_p_unit;
2241  try
2242  {
2243  initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
2244  .transform_real_to_unit_cell(cell, p);
2245  }
2246  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2247  {
2248  // mirror the conditions of the code below to determine if we need to
2249  // use an arbitrary starting point or if we just need to rethrow the
2250  // exception
2251  for (unsigned int d = 0; d < dim; ++d)
2252  initial_p_unit[d] = 0.5;
2253  }
2254 
2255  // TODO
2256  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2257 
2258  // for (unsigned int d=0; d<dim; ++d)
2259  // initial_p_unit[d] = 0.;
2260 
2261  const Quadrature<dim> point_quadrature(initial_p_unit);
2262 
2264  if (spacedim > dim)
2265  update_flags |= update_jacobian_grads;
2266  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2267  get_data(update_flags, point_quadrature));
2268  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2269  ExcInternalError());
2270 
2271  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2272 
2273  return do_transform_real_to_unit_cell(cell,
2274  p,
2275  initial_p_unit,
2276  dynamic_cast<InternalData &>(*mdata));
2277 }
2278 
2279 
2280 template <int dim, int spacedim, typename VectorType>
2281 Point<dim>
2283  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2284  const Point<spacedim> & p,
2285  const Point<dim> & initial_p_unit,
2286  InternalData & mdata) const
2287 {
2288  const unsigned int n_shapes = mdata.shape_values.size();
2289  (void)n_shapes;
2290  Assert(n_shapes != 0, ExcInternalError());
2291  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2292 
2293 
2294  // Newton iteration to solve
2295  // f(x)=p(x)-p=0
2296  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2297  // The start value was set to be the
2298  // linear approximation to the cell
2299  // The shape values and derivatives
2300  // of the mapping at this point are
2301  // previously computed.
2302  // f(x)
2303  Point<dim> p_unit = initial_p_unit;
2304  Point<dim> f;
2305  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
2307  Tensor<1, spacedim> p_minus_F = p - p_real;
2308  const double eps = 1.e-12 * cell->diameter();
2309  const unsigned int newton_iteration_limit = 20;
2310  unsigned int newton_iteration = 0;
2311  while (p_minus_F.norm_square() > eps * eps)
2312  {
2313  // f'(x)
2314  Point<spacedim> DF[dim];
2315  Tensor<2, dim> df;
2316  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2317  {
2318  const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2319  const unsigned int comp_k =
2320  euler_dof_handler->get_fe().system_to_component_index(k).first;
2321  if (fe_mask[comp_k])
2322  for (unsigned int j = 0; j < dim; ++j)
2323  DF[j][fe_to_real[comp_k]] +=
2324  mdata.local_dof_values[k] * grad_k[j];
2325  }
2326  for (unsigned int j = 0; j < dim; ++j)
2327  {
2328  f[j] = DF[j] * p_minus_F;
2329  for (unsigned int l = 0; l < dim; ++l)
2330  df[j][l] = -DF[j] * DF[l];
2331  }
2332  // Solve [f'(x)]d=f(x)
2333  const Tensor<1, dim> delta =
2334  invert(df) * static_cast<const Tensor<1, dim> &>(f);
2335  // do a line search
2336  double step_length = 1;
2337  do
2338  {
2339  // update of p_unit. The
2340  // spacedimth component of
2341  // transformed point is simply
2342  // ignored in codimension one
2343  // case. When this component is
2344  // not zero, then we are
2345  // projecting the point to the
2346  // surface or curve identified
2347  // by the cell.
2348  Point<dim> p_unit_trial = p_unit;
2349  for (unsigned int i = 0; i < dim; ++i)
2350  p_unit_trial[i] -= step_length * delta[i];
2351  // shape values and derivatives
2352  // at new p_unit point
2353  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
2354  mdata);
2355  // f(x)
2356  Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
2357  const Tensor<1, spacedim> f_trial = p - p_real_trial;
2358  // see if we are making progress with the current step length
2359  // and if not, reduce it by a factor of two and try again
2360  if (f_trial.norm() < p_minus_F.norm())
2361  {
2362  p_real = p_real_trial;
2363  p_unit = p_unit_trial;
2364  p_minus_F = f_trial;
2365  break;
2366  }
2367  else if (step_length > 0.05)
2368  step_length /= 2;
2369  else
2370  goto failure;
2371  }
2372  while (true);
2373  ++newton_iteration;
2374  if (newton_iteration > newton_iteration_limit)
2375  goto failure;
2376  }
2377  return p_unit;
2378  // if we get to the following label, then we have either run out
2379  // of Newton iterations, or the line search has not converged.
2380  // in either case, we need to give up, so throw an exception that
2381  // can then be caught
2382 failure:
2383  AssertThrow(false,
2385  // ...the compiler wants us to return something, though we can
2386  // of course never get here...
2387  return {};
2388 }
2389 
2390 
2391 template <int dim, int spacedim, typename VectorType>
2392 unsigned int
2394 {
2395  return euler_dof_handler->get_fe().degree;
2396 }
2397 
2398 
2399 
2400 template <int dim, int spacedim, typename VectorType>
2403 {
2404  return this->fe_mask;
2405 }
2406 
2407 
2408 template <int dim, int spacedim, typename VectorType>
2409 std::unique_ptr<Mapping<dim, spacedim>>
2411 {
2412  return std::make_unique<MappingFEField<dim, spacedim, VectorType>>(*this);
2413 }
2414 
2415 
2416 template <int dim, int spacedim, typename VectorType>
2417 void
2419  const typename Triangulation<dim, spacedim>::cell_iterator & cell,
2421  const
2422 {
2423  Assert(euler_dof_handler != nullptr,
2424  ExcMessage("euler_dof_handler is empty"));
2425 
2426  typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
2428  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2429  if (uses_level_dofs)
2430  {
2431  AssertIndexRange(cell->level(), euler_vector.size());
2432  AssertDimension(euler_vector[cell->level()]->size(),
2433  euler_dof_handler->n_dofs(cell->level()));
2434  }
2435  else
2436  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2437 
2438  if (uses_level_dofs)
2439  dof_cell->get_mg_dof_indices(data.local_dof_indices);
2440  else
2441  dof_cell->get_dof_indices(data.local_dof_indices);
2442 
2443  const VectorType &vector =
2444  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2445 
2446  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2447  data.local_dof_values[i] =
2449  data.local_dof_indices[i]);
2450 }
2451 
2452 // explicit instantiations
2453 #define SPLIT_INSTANTIATIONS_COUNT 2
2454 #ifndef SPLIT_INSTANTIATIONS_INDEX
2455 # define SPLIT_INSTANTIATIONS_INDEX 0
2456 #endif
2457 #include "mapping_fe_field.inst"
2458 
2459 
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:704
unsigned int size() const
DerivativeForm< 1, spacedim, dim, Number > transpose() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
InternalData(const FiniteElement< dim, spacedim > &fe, const ComponentMask &mask)
std::vector< double > quadrature_weights
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< std::vector< Tensor< 1, spacedim > > > aux
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
std::vector< Tensor< 3, dim > > shape_third_derivatives
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< types::global_dof_index > local_dof_indices
virtual std::size_t memory_consumption() const override
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< Tensor< 1, dim > > shape_derivatives
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< double > local_dof_values
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< double > shape_values
std::vector< Tensor< 2, dim > > shape_second_derivatives
std::vector< double > volume_elements
virtual void compute_shapes_virtual(const std::vector< Point< dim >> &unit_points, typename MappingFEField< dim, spacedim, VectorType >::InternalData &data) const
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
std::vector< unsigned int > fe_to_real
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
const ComponentMask fe_mask
void update_internal_dofs(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename MappingFEField< dim, spacedim, VectorType >::InternalData &data) const
void compute_face_data(const UpdateFlags update_flags, const Quadrature< dim > &q, const unsigned int n_original_q_points, InternalData &data) const
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
std::vector< SmartPointer< const VectorType, MappingFEField< dim, spacedim, VectorType > > > euler_vector
virtual bool preserves_vertex_locations() const override
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
virtual void fill_fe_immersed_surface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const NonMatching::ImmersedSurfaceQuadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
friend class MappingFEField
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
void compute_data(const UpdateFlags update_flags, const Quadrature< dim > &q, const unsigned int n_original_q_points, InternalData &data) const
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
ReferenceCell reference_cell
SmartPointer< const DoFHandler< dim, spacedim >, MappingFEField< dim, spacedim, VectorType > > euler_dof_handler
unsigned int get_degree() const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Point< spacedim > do_transform_unit_to_real_cell(const InternalData &mdata) const
Threads::Mutex fe_values_mutex
const bool uses_level_dofs
ComponentMask get_component_mask() const
Point< dim > do_transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_p_unit, InternalData &mdata) const
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
FEValues< dim, spacedim > fe_values
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Abstract base class for mapping classes.
Definition: mapping.h:317
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition: point.h:112
static DataSetDescriptor cell()
Definition: qprojector.h:394
static DataSetDescriptor subface(const ReferenceCell &reference_cell, 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)
const std::vector< Point< dim > > & get_points() const
const std::vector< double > & get_weights() const
unsigned int size() const
unsigned int n_faces() const
unsigned int get_dimension() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
numbers::NumberTraits< Number >::real_type norm() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
unsigned int size() const
Definition: collection.h:264
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
std::vector< Tensor< 1, spacedim > > normal_vectors
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< Tensor< 1, spacedim > > boundary_forms
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 3 > vertices[4]
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Point< 2 > second
Definition: grid_out.cc:4616
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInactiveCell()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_volume_elements
Determinant of the Jacobian.
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_jacobian_3rd_derivatives
@ update_values
Shape function values.
@ update_jacobian_grads
Gradient of volume element.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_quadrature_points
Transformed quadrature points.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
MappingKind
Definition: mapping.h:78
@ mapping_piola
Definition: mapping.h:113
@ mapping_covariant_gradient
Definition: mapping.h:99
@ mapping_covariant
Definition: mapping.h:88
@ mapping_contravariant
Definition: mapping.h:93
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:285
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:189
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:490
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:480
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim >> &output)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 4, spacedim >> &jacobian_pushed_forward_2nd_derivatives)
void maybe_update_jacobian_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim >> &unit_points, const std::vector< Polynomials::Polynomial< double >> &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 5, spacedim >> &jacobian_pushed_forward_3rd_derivatives)
void maybe_compute_face_data(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int n_q_points, const std::vector< double > &weights, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void transform_fields(const ArrayView< const Tensor< rank, dim >> &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim >> &output)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)