Reference documentation for deal.II version GIT d9d8a449a2 2022-08-17 08:45: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 #include <deal.II/fe/mapping_q1.h>
34 
36 
41 #include <deal.II/lac/la_vector.h>
48 #include <deal.II/lac/vector.h>
49 
51 
52 #include <fstream>
53 #include <memory>
54 #include <numeric>
55 
56 
57 
59 
60 
61 template <int dim, int spacedim, typename VectorType>
64  const ComponentMask & mask)
65  : unit_tangentials()
66  , n_shape_functions(fe.n_dofs_per_cell())
67  , mask(mask)
68  , local_dof_indices(fe.n_dofs_per_cell())
69  , local_dof_values(fe.n_dofs_per_cell())
70 {}
71 
72 
73 
74 template <int dim, int spacedim, typename VectorType>
75 std::size_t
77  const
78 {
79  Assert(false, ExcNotImplemented());
80  return 0;
81 }
82 
83 
84 
85 template <int dim, int spacedim, typename VectorType>
86 double &
88  const unsigned int qpoint,
89  const unsigned int shape_nr)
90 {
91  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
92  return shape_values[qpoint * n_shape_functions + shape_nr];
93 }
94 
95 
96 template <int dim, int spacedim, typename VectorType>
97 const Tensor<1, dim> &
99  const unsigned int qpoint,
100  const unsigned int shape_nr) const
101 {
102  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
103  shape_derivatives.size());
104  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
105 }
106 
107 
108 
109 template <int dim, int spacedim, typename VectorType>
112  const unsigned int qpoint,
113  const unsigned int shape_nr)
114 {
115  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
116  shape_derivatives.size());
117  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
118 }
119 
120 
121 template <int dim, int spacedim, typename VectorType>
122 const Tensor<2, dim> &
124  const unsigned int qpoint,
125  const unsigned int shape_nr) const
126 {
127  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
128  shape_second_derivatives.size());
129  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
130 }
131 
132 
133 
134 template <int dim, int spacedim, typename VectorType>
137  const unsigned int qpoint,
138  const unsigned int shape_nr)
139 {
140  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
141  shape_second_derivatives.size());
142  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
143 }
144 
145 
146 template <int dim, int spacedim, typename VectorType>
147 const Tensor<3, dim> &
149  const unsigned int qpoint,
150  const unsigned int shape_nr) const
151 {
152  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
153  shape_third_derivatives.size());
154  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
155 }
156 
157 
158 
159 template <int dim, int spacedim, typename VectorType>
162  const unsigned int qpoint,
163  const unsigned int shape_nr)
164 {
165  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
166  shape_third_derivatives.size());
167  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
168 }
169 
170 
171 template <int dim, int spacedim, typename VectorType>
172 const Tensor<4, dim> &
174  const unsigned int qpoint,
175  const unsigned int shape_nr) const
176 {
177  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
178  shape_fourth_derivatives.size());
179  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
180 }
181 
182 
183 
184 template <int dim, int spacedim, typename VectorType>
187  const unsigned int qpoint,
188  const unsigned int shape_nr)
189 {
190  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
191  shape_fourth_derivatives.size());
192  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
193 }
194 
195 
196 
197 template <int dim, int spacedim, typename VectorType>
200  const VectorType & euler_vector,
201  const ComponentMask & mask)
203  , uses_level_dofs(false)
205  , euler_dof_handler(&euler_dof_handler)
206  , fe_mask(mask.size() != 0u ?
207  mask :
209  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
210  true))
211  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
212  , fe_values(this->euler_dof_handler->get_fe(),
213  reference_cell.template get_nodal_type_quadrature<dim>(),
215 {
216  unsigned int size = 0;
217  for (unsigned int i = 0; i < fe_mask.size(); ++i)
218  {
219  if (fe_mask[i])
220  fe_to_real[i] = size++;
221  }
222  AssertDimension(size, spacedim);
223 }
224 
225 
226 
227 template <int dim, int spacedim, typename VectorType>
229  const DoFHandler<dim, spacedim> &euler_dof_handler,
230  const std::vector<VectorType> & euler_vector,
231  const ComponentMask & mask)
232  : reference_cell(euler_dof_handler.get_fe().reference_cell())
233  , uses_level_dofs(true)
234  , euler_dof_handler(&euler_dof_handler)
235  , fe_mask(mask.size() != 0u ?
236  mask :
238  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
239  true))
240  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
241  , fe_values(this->euler_dof_handler->get_fe(),
242  reference_cell.template get_nodal_type_quadrature<dim>(),
244 {
245  unsigned int size = 0;
246  for (unsigned int i = 0; i < fe_mask.size(); ++i)
247  {
248  if (fe_mask[i])
249  fe_to_real[i] = size++;
250  }
251  AssertDimension(size, spacedim);
252 
253  Assert(euler_dof_handler.has_level_dofs(),
254  ExcMessage("The underlying DoFHandler object did not call "
255  "distribute_mg_dofs(). In this case, the construction via "
256  "level vectors does not make sense."));
258  euler_dof_handler.get_triangulation().n_global_levels());
259  this->euler_vector.clear();
260  this->euler_vector.resize(euler_vector.size());
261  for (unsigned int i = 0; i < euler_vector.size(); ++i)
262  this->euler_vector[i] = &euler_vector[i];
263 }
264 
265 
266 
267 template <int dim, int spacedim, typename VectorType>
269  const DoFHandler<dim, spacedim> &euler_dof_handler,
270  const MGLevelObject<VectorType> &euler_vector,
271  const ComponentMask & mask)
272  : reference_cell(euler_dof_handler.get_fe().reference_cell())
273  , uses_level_dofs(true)
274  , euler_dof_handler(&euler_dof_handler)
275  , fe_mask(mask.size() != 0u ?
276  mask :
278  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
279  true))
280  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
281  , fe_values(this->euler_dof_handler->get_fe(),
282  reference_cell.template get_nodal_type_quadrature<dim>(),
284 {
285  unsigned int size = 0;
286  for (unsigned int i = 0; i < fe_mask.size(); ++i)
287  {
288  if (fe_mask[i])
289  fe_to_real[i] = size++;
290  }
291  AssertDimension(size, spacedim);
292 
293  Assert(euler_dof_handler.has_level_dofs(),
294  ExcMessage("The underlying DoFHandler object did not call "
295  "distribute_mg_dofs(). In this case, the construction via "
296  "level vectors does not make sense."));
297  AssertDimension(euler_vector.max_level() + 1,
298  euler_dof_handler.get_triangulation().n_global_levels());
299  this->euler_vector.clear();
300  this->euler_vector.resize(
301  euler_dof_handler.get_triangulation().n_global_levels());
302  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
303  ++i)
304  this->euler_vector[i] = &euler_vector[i];
305 }
306 
307 
308 
309 template <int dim, int spacedim, typename VectorType>
312  : reference_cell(mapping.reference_cell)
313  , uses_level_dofs(mapping.uses_level_dofs)
314  , euler_vector(mapping.euler_vector)
315  , euler_dof_handler(mapping.euler_dof_handler)
316  , fe_mask(mapping.fe_mask)
317  , fe_to_real(mapping.fe_to_real)
318  , fe_values(mapping.euler_dof_handler->get_fe(),
319  reference_cell.template get_nodal_type_quadrature<dim>(),
321 {}
322 
323 
324 
325 template <int dim, int spacedim, typename VectorType>
326 inline const double &
328  const unsigned int qpoint,
329  const unsigned int shape_nr) const
330 {
331  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
332  return shape_values[qpoint * n_shape_functions + shape_nr];
333 }
334 
335 
336 
337 template <int dim, int spacedim, typename VectorType>
338 bool
340 {
341  return false;
342 }
343 
344 
345 
346 template <int dim, int spacedim, typename VectorType>
347 bool
349  const ReferenceCell &reference_cell) const
350 {
352  ExcMessage("The dimension of your mapping (" +
353  Utilities::to_string(dim) +
354  ") and the reference cell cell_type (" +
356  " ) do not agree."));
357 
358  return this->reference_cell == reference_cell;
359 }
360 
361 
362 
363 template <int dim, int spacedim, typename VectorType>
364 boost::container::small_vector<Point<spacedim>,
367  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
368 {
369  // we transform our tria iterator into a dof iterator so we can access
370  // data not associated with triangulations
371  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
372  *cell, euler_dof_handler);
373 
374  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
375  AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
377  euler_dof_handler->get_fe().n_components());
378  if (uses_level_dofs)
379  {
380  AssertIndexRange(cell->level(), euler_vector.size());
381  AssertDimension(euler_vector[cell->level()]->size(),
382  euler_dof_handler->n_dofs(cell->level()));
383  }
384  else
385  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
386 
387  {
388  std::lock_guard<std::mutex> lock(fe_values_mutex);
389  fe_values.reinit(dof_cell);
390  }
391  const unsigned int dofs_per_cell =
392  euler_dof_handler->get_fe().n_dofs_per_cell();
393  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
394  if (uses_level_dofs)
395  dof_cell->get_mg_dof_indices(dof_indices);
396  else
397  dof_cell->get_dof_indices(dof_indices);
398 
399  const VectorType &vector =
400  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
401 
402  boost::container::small_vector<Point<spacedim>,
404  vertices(cell->n_vertices());
405  for (unsigned int i = 0; i < dofs_per_cell; ++i)
406  {
407  const unsigned int comp = fe_to_real
408  [euler_dof_handler->get_fe().system_to_component_index(i).first];
409  if (comp != numbers::invalid_unsigned_int)
410  {
411  typename VectorType::value_type value =
412  internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
413  if (euler_dof_handler->get_fe().is_primitive(i))
414  for (const unsigned int v : cell->vertex_indices())
415  vertices[v][comp] += fe_values.shape_value(i, v) * value;
416  else
417  Assert(false, ExcNotImplemented());
418  }
419  }
420 
421  return vertices;
422 }
423 
424 
425 
426 template <int dim, int spacedim, typename VectorType>
427 void
429  const std::vector<Point<dim>> & unit_points,
431 {
432  const auto fe = &euler_dof_handler->get_fe();
433  const unsigned int n_points = unit_points.size();
434 
435  for (unsigned int point = 0; point < n_points; ++point)
436  {
437  if (data.shape_values.size() != 0)
438  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
439  data.shape(point, i) = fe->shape_value(i, unit_points[point]);
440 
441  if (data.shape_derivatives.size() != 0)
442  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
443  data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
444 
445  if (data.shape_second_derivatives.size() != 0)
446  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
447  data.second_derivative(point, i) =
448  fe->shape_grad_grad(i, unit_points[point]);
449 
450  if (data.shape_third_derivatives.size() != 0)
451  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
452  data.third_derivative(point, i) =
453  fe->shape_3rd_derivative(i, unit_points[point]);
454 
455  if (data.shape_fourth_derivatives.size() != 0)
456  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
457  data.fourth_derivative(point, i) =
458  fe->shape_4th_derivative(i, unit_points[point]);
459  }
460 }
461 
462 
463 template <int dim, int spacedim, typename VectorType>
466  const UpdateFlags in) const
467 {
468  // add flags if the respective quantities are necessary to compute
469  // what we need. note that some flags appear in both conditions and
470  // in subsequent set operations. this leads to some circular
471  // logic. the only way to treat this is to iterate. since there are
472  // 5 if-clauses in the loop, it will take at most 4 iterations to
473  // converge. do them:
474  UpdateFlags out = in;
475  for (unsigned int i = 0; i < 5; ++i)
476  {
477  // The following is a little incorrect:
478  // If not applied on a face,
479  // update_boundary_forms does not
480  // make sense. On the other hand,
481  // it is necessary on a
482  // face. Currently,
483  // update_boundary_forms is simply
484  // ignored for the interior of a
485  // cell.
487  out |= update_boundary_forms;
488 
493 
494  if (out &
499 
500  // The contravariant transformation
501  // is a Piola transformation, which
502  // requires the determinant of the
503  // Jacobi matrix of the transformation.
504  // Therefore these values have to be
505  // updated for each cell.
507  out |= update_JxW_values;
508 
509  if (out & update_normal_vectors)
510  out |= update_JxW_values;
511  }
512 
513  return out;
514 }
515 
516 
517 
518 template <int dim, int spacedim, typename VectorType>
519 void
521  const UpdateFlags update_flags,
522  const Quadrature<dim> &q,
523  const unsigned int n_original_q_points,
524  InternalData & data) const
525 {
526  // store the flags in the internal data object so we can access them
527  // in fill_fe_*_values(). use the transitive hull of the required
528  // flags
529  data.update_each = requires_update_flags(update_flags);
530 
531  const unsigned int n_q_points = q.size();
532 
533  // see if we need the (transformation) shape function values
534  // and/or gradients and resize the necessary arrays
535  if (data.update_each & update_quadrature_points)
536  data.shape_values.resize(data.n_shape_functions * n_q_points);
537 
538  if (data.update_each &
542  data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
543 
544  if (data.update_each & update_covariant_transformation)
545  data.covariant.resize(n_original_q_points);
546 
547  if (data.update_each & update_contravariant_transformation)
548  data.contravariant.resize(n_original_q_points);
549 
550  if (data.update_each & update_volume_elements)
551  data.volume_elements.resize(n_original_q_points);
552 
553  if (data.update_each &
555  data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
556 
557  if (data.update_each & (update_jacobian_2nd_derivatives |
559  data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
560 
561  if (data.update_each & (update_jacobian_3rd_derivatives |
563  data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
564 
566 
567  // This (for face values and simplices) can be different for different calls,
568  // so always copy
569  data.quadrature_weights = q.get_weights();
570 }
571 
572 
573 template <int dim, int spacedim, typename VectorType>
574 void
576  const UpdateFlags update_flags,
577  const Quadrature<dim> &q,
578  const unsigned int n_original_q_points,
579  InternalData & data) const
580 {
581  compute_data(update_flags, q, n_original_q_points, data);
582 
583  if (dim > 1)
584  {
585  if (data.update_each & update_boundary_forms)
586  {
587  data.aux.resize(
588  dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
589 
590 
591  // TODO: only a single reference cell type possible...
592  const auto n_faces = reference_cell.n_faces();
593 
594  // Compute tangentials to the unit cell.
595  for (unsigned int i = 0; i < n_faces; ++i)
596  {
597  data.unit_tangentials[i].resize(n_original_q_points);
598  std::fill(
599  data.unit_tangentials[i].begin(),
600  data.unit_tangentials[i].end(),
601  reference_cell.template unit_tangential_vectors<dim>(i, 0));
602  if (dim > 2)
603  {
604  data.unit_tangentials[n_faces + i].resize(
605  n_original_q_points);
606  std::fill(
607  data.unit_tangentials[n_faces + i].begin(),
608  data.unit_tangentials[n_faces + i].end(),
609  reference_cell.template unit_tangential_vectors<dim>(i, 1));
610  }
611  }
612  }
613  }
614 }
615 
616 
617 template <int dim, int spacedim, typename VectorType>
618 typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
620  const UpdateFlags update_flags,
621  const Quadrature<dim> &quadrature) const
622 {
623  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
624  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
625  auto &data = dynamic_cast<InternalData &>(*data_ptr);
626  this->compute_data(update_flags, quadrature, quadrature.size(), data);
627 
628  return data_ptr;
629 }
630 
631 
632 
633 template <int dim, int spacedim, typename VectorType>
634 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
636  const UpdateFlags update_flags,
637  const hp::QCollection<dim - 1> &quadrature) const
638 {
639  AssertDimension(quadrature.size(), 1);
640 
641  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
642  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
643  auto & data = dynamic_cast<InternalData &>(*data_ptr);
644  const Quadrature<dim> q(
646  this->compute_face_data(update_flags, q, quadrature[0].size(), data);
647 
648  return data_ptr;
649 }
650 
651 
652 template <int dim, int spacedim, typename VectorType>
653 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
655  const UpdateFlags update_flags,
656  const Quadrature<dim - 1> &quadrature) const
657 {
658  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
659  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
660  auto & data = dynamic_cast<InternalData &>(*data_ptr);
661  const Quadrature<dim> q(
663  this->compute_face_data(update_flags, q, quadrature.size(), data);
664 
665  return data_ptr;
666 }
667 
668 
669 
670 namespace internal
671 {
672  namespace MappingFEFieldImplementation
673  {
674  namespace
675  {
682  template <int dim, int spacedim, typename VectorType>
683  void
685  const typename ::QProjector<dim>::DataSetDescriptor data_set,
686  const typename ::MappingFEField<dim, spacedim, VectorType>::
687  InternalData & data,
689  const ComponentMask & fe_mask,
690  const std::vector<unsigned int> & fe_to_real,
691  std::vector<Point<spacedim>> & quadrature_points)
692  {
693  const UpdateFlags update_flags = data.update_each;
694 
695  if (update_flags & update_quadrature_points)
696  {
697  for (unsigned int point = 0; point < quadrature_points.size();
698  ++point)
699  {
700  Point<spacedim> result;
701  const double * shape = &data.shape(point + data_set, 0);
702 
703  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
704  {
705  const unsigned int comp_k =
706  fe.system_to_component_index(k).first;
707  if (fe_mask[comp_k])
708  result[fe_to_real[comp_k]] +=
709  data.local_dof_values[k] * shape[k];
710  }
711 
712  quadrature_points[point] = result;
713  }
714  }
715  }
716 
724  template <int dim, int spacedim, typename VectorType>
725  void
727  const typename ::QProjector<dim>::DataSetDescriptor data_set,
728  const typename ::MappingFEField<dim, spacedim, VectorType>::
729  InternalData & data,
731  const ComponentMask & fe_mask,
732  const std::vector<unsigned int> & fe_to_real)
733  {
734  const UpdateFlags update_flags = data.update_each;
735 
736  // then Jacobians
737  if (update_flags & update_contravariant_transformation)
738  {
739  const unsigned int n_q_points = data.contravariant.size();
740 
741  Assert(data.n_shape_functions > 0, ExcInternalError());
742 
743  for (unsigned int point = 0; point < n_q_points; ++point)
744  {
745  const Tensor<1, dim> *data_derv =
746  &data.derivative(point + data_set, 0);
747 
748  Tensor<1, dim> result[spacedim];
749 
750  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
751  {
752  const unsigned int comp_k =
753  fe.system_to_component_index(k).first;
754  if (fe_mask[comp_k])
755  result[fe_to_real[comp_k]] +=
756  data.local_dof_values[k] * data_derv[k];
757  }
758 
759  // write result into contravariant data
760  for (unsigned int i = 0; i < spacedim; ++i)
761  {
762  data.contravariant[point][i] = result[i];
763  }
764  }
765  }
766 
767  if (update_flags & update_covariant_transformation)
768  {
769  AssertDimension(data.covariant.size(), data.contravariant.size());
770  for (unsigned int point = 0; point < data.contravariant.size();
771  ++point)
772  data.covariant[point] =
773  (data.contravariant[point]).covariant_form();
774  }
775 
776  if (update_flags & update_volume_elements)
777  {
778  AssertDimension(data.contravariant.size(),
779  data.volume_elements.size());
780  for (unsigned int point = 0; point < data.contravariant.size();
781  ++point)
782  data.volume_elements[point] =
783  data.contravariant[point].determinant();
784  }
785  }
786 
793  template <int dim, int spacedim, typename VectorType>
794  void
796  const typename ::QProjector<dim>::DataSetDescriptor data_set,
797  const typename ::MappingFEField<dim, spacedim, VectorType>::
798  InternalData & data,
799  const FiniteElement<dim, spacedim> & fe,
800  const ComponentMask & fe_mask,
801  const std::vector<unsigned int> & fe_to_real,
802  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
803  {
804  const UpdateFlags update_flags = data.update_each;
805  if (update_flags & update_jacobian_grads)
806  {
807  const unsigned int n_q_points = jacobian_grads.size();
808 
809  for (unsigned int point = 0; point < n_q_points; ++point)
810  {
811  const Tensor<2, dim> *second =
812  &data.second_derivative(point + data_set, 0);
813 
815 
816  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
817  {
818  const unsigned int comp_k =
819  fe.system_to_component_index(k).first;
820  if (fe_mask[comp_k])
821  for (unsigned int j = 0; j < dim; ++j)
822  for (unsigned int l = 0; l < dim; ++l)
823  result[fe_to_real[comp_k]][j][l] +=
824  (second[k][j][l] * data.local_dof_values[k]);
825  }
826 
827  // never touch any data for j=dim in case dim<spacedim, so
828  // it will always be zero as it was initialized
829  for (unsigned int i = 0; i < spacedim; ++i)
830  for (unsigned int j = 0; j < dim; ++j)
831  for (unsigned int l = 0; l < dim; ++l)
832  jacobian_grads[point][i][j][l] = result[i][j][l];
833  }
834  }
835  }
836 
843  template <int dim, int spacedim, typename VectorType>
844  void
846  const typename ::QProjector<dim>::DataSetDescriptor data_set,
847  const typename ::MappingFEField<dim, spacedim, VectorType>::
848  InternalData & data,
850  const ComponentMask & fe_mask,
851  const std::vector<unsigned int> & fe_to_real,
852  std::vector<Tensor<3, spacedim>> & jacobian_pushed_forward_grads)
853  {
854  const UpdateFlags update_flags = data.update_each;
855  if (update_flags & update_jacobian_pushed_forward_grads)
856  {
857  const unsigned int n_q_points =
858  jacobian_pushed_forward_grads.size();
859 
860  double tmp[spacedim][spacedim][spacedim];
861  for (unsigned int point = 0; point < n_q_points; ++point)
862  {
863  const Tensor<2, dim> *second =
864  &data.second_derivative(point + data_set, 0);
865 
867 
868  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
869  {
870  const unsigned int comp_k =
871  fe.system_to_component_index(k).first;
872  if (fe_mask[comp_k])
873  for (unsigned int j = 0; j < dim; ++j)
874  for (unsigned int l = 0; l < dim; ++l)
875  result[fe_to_real[comp_k]][j][l] +=
876  (second[k][j][l] * data.local_dof_values[k]);
877  }
878 
879  // first push forward the j-components
880  for (unsigned int i = 0; i < spacedim; ++i)
881  for (unsigned int j = 0; j < spacedim; ++j)
882  for (unsigned int l = 0; l < dim; ++l)
883  {
884  tmp[i][j][l] =
885  result[i][0][l] * data.covariant[point][j][0];
886  for (unsigned int jr = 1; jr < dim; ++jr)
887  {
888  tmp[i][j][l] +=
889  result[i][jr][l] * data.covariant[point][j][jr];
890  }
891  }
892 
893  // now, pushing forward the l-components
894  for (unsigned int i = 0; i < spacedim; ++i)
895  for (unsigned int j = 0; j < spacedim; ++j)
896  for (unsigned int l = 0; l < spacedim; ++l)
897  {
898  jacobian_pushed_forward_grads[point][i][j][l] =
899  tmp[i][j][0] * data.covariant[point][l][0];
900  for (unsigned int lr = 1; lr < dim; ++lr)
901  {
902  jacobian_pushed_forward_grads[point][i][j][l] +=
903  tmp[i][j][lr] * data.covariant[point][l][lr];
904  }
905  }
906  }
907  }
908  }
909 
916  template <int dim, int spacedim, typename VectorType>
917  void
919  const typename ::QProjector<dim>::DataSetDescriptor data_set,
920  const typename ::MappingFEField<dim, spacedim, VectorType>::
921  InternalData & data,
922  const FiniteElement<dim, spacedim> & fe,
923  const ComponentMask & fe_mask,
924  const std::vector<unsigned int> & fe_to_real,
925  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
926  {
927  const UpdateFlags update_flags = data.update_each;
928  if (update_flags & update_jacobian_2nd_derivatives)
929  {
930  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
931 
932  for (unsigned int point = 0; point < n_q_points; ++point)
933  {
934  const Tensor<3, dim> *third =
935  &data.third_derivative(point + data_set, 0);
936 
938 
939  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
940  {
941  const unsigned int comp_k =
942  fe.system_to_component_index(k).first;
943  if (fe_mask[comp_k])
944  for (unsigned int j = 0; j < dim; ++j)
945  for (unsigned int l = 0; l < dim; ++l)
946  for (unsigned int m = 0; m < dim; ++m)
947  result[fe_to_real[comp_k]][j][l][m] +=
948  (third[k][j][l][m] * data.local_dof_values[k]);
949  }
950 
951  // never touch any data for j=dim in case dim<spacedim, so
952  // it will always be zero as it was initialized
953  for (unsigned int i = 0; i < spacedim; ++i)
954  for (unsigned int j = 0; j < dim; ++j)
955  for (unsigned int l = 0; l < dim; ++l)
956  for (unsigned int m = 0; m < dim; ++m)
957  jacobian_2nd_derivatives[point][i][j][l][m] =
958  result[i][j][l][m];
959  }
960  }
961  }
962 
970  template <int dim, int spacedim, typename VectorType>
971  void
973  const typename ::QProjector<dim>::DataSetDescriptor data_set,
974  const typename ::MappingFEField<dim, spacedim, VectorType>::
975  InternalData & data,
977  const ComponentMask & fe_mask,
978  const std::vector<unsigned int> & fe_to_real,
979  std::vector<Tensor<4, spacedim>>
980  &jacobian_pushed_forward_2nd_derivatives)
981  {
982  const UpdateFlags update_flags = data.update_each;
984  {
985  const unsigned int n_q_points =
986  jacobian_pushed_forward_2nd_derivatives.size();
987 
988  double tmp[spacedim][spacedim][spacedim][spacedim];
989  for (unsigned int point = 0; point < n_q_points; ++point)
990  {
991  const Tensor<3, dim> *third =
992  &data.third_derivative(point + data_set, 0);
993 
995 
996  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
997  {
998  const unsigned int comp_k =
999  fe.system_to_component_index(k).first;
1000  if (fe_mask[comp_k])
1001  for (unsigned int j = 0; j < dim; ++j)
1002  for (unsigned int l = 0; l < dim; ++l)
1003  for (unsigned int m = 0; m < dim; ++m)
1004  result[fe_to_real[comp_k]][j][l][m] +=
1005  (third[k][j][l][m] * data.local_dof_values[k]);
1006  }
1007 
1008  // push forward the j-coordinate
1009  for (unsigned int i = 0; i < spacedim; ++i)
1010  for (unsigned int j = 0; j < spacedim; ++j)
1011  for (unsigned int l = 0; l < dim; ++l)
1012  for (unsigned int m = 0; m < dim; ++m)
1013  {
1014  jacobian_pushed_forward_2nd_derivatives
1015  [point][i][j][l][m] =
1016  result[i][0][l][m] * data.covariant[point][j][0];
1017  for (unsigned int jr = 1; jr < dim; ++jr)
1018  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1019  [l][m] +=
1020  result[i][jr][l][m] *
1021  data.covariant[point][j][jr];
1022  }
1023 
1024  // push forward the l-coordinate
1025  for (unsigned int i = 0; i < spacedim; ++i)
1026  for (unsigned int j = 0; j < spacedim; ++j)
1027  for (unsigned int l = 0; l < spacedim; ++l)
1028  for (unsigned int m = 0; m < dim; ++m)
1029  {
1030  tmp[i][j][l][m] =
1031  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1032  [0][m] *
1033  data.covariant[point][l][0];
1034  for (unsigned int lr = 1; lr < dim; ++lr)
1035  tmp[i][j][l][m] +=
1036  jacobian_pushed_forward_2nd_derivatives[point][i]
1037  [j][lr]
1038  [m] *
1039  data.covariant[point][l][lr];
1040  }
1041 
1042  // push forward the m-coordinate
1043  for (unsigned int i = 0; i < spacedim; ++i)
1044  for (unsigned int j = 0; j < spacedim; ++j)
1045  for (unsigned int l = 0; l < spacedim; ++l)
1046  for (unsigned int m = 0; m < spacedim; ++m)
1047  {
1048  jacobian_pushed_forward_2nd_derivatives
1049  [point][i][j][l][m] =
1050  tmp[i][j][l][0] * data.covariant[point][m][0];
1051  for (unsigned int mr = 1; mr < dim; ++mr)
1052  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1053  [l][m] +=
1054  tmp[i][j][l][mr] * data.covariant[point][m][mr];
1055  }
1056  }
1057  }
1058  }
1059 
1066  template <int dim, int spacedim, typename VectorType>
1067  void
1069  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1070  const typename ::MappingFEField<dim, spacedim, VectorType>::
1071  InternalData & data,
1072  const FiniteElement<dim, spacedim> & fe,
1073  const ComponentMask & fe_mask,
1074  const std::vector<unsigned int> & fe_to_real,
1075  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1076  {
1077  const UpdateFlags update_flags = data.update_each;
1078  if (update_flags & update_jacobian_3rd_derivatives)
1079  {
1080  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1081 
1082  for (unsigned int point = 0; point < n_q_points; ++point)
1083  {
1084  const Tensor<4, dim> *fourth =
1085  &data.fourth_derivative(point + data_set, 0);
1086 
1088 
1089  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1090  {
1091  const unsigned int comp_k =
1092  fe.system_to_component_index(k).first;
1093  if (fe_mask[comp_k])
1094  for (unsigned int j = 0; j < dim; ++j)
1095  for (unsigned int l = 0; l < dim; ++l)
1096  for (unsigned int m = 0; m < dim; ++m)
1097  for (unsigned int n = 0; n < dim; ++n)
1098  result[fe_to_real[comp_k]][j][l][m][n] +=
1099  (fourth[k][j][l][m][n] *
1100  data.local_dof_values[k]);
1101  }
1102 
1103  // never touch any data for j,l,m,n=dim in case
1104  // dim<spacedim, so it will always be zero as it was
1105  // initialized
1106  for (unsigned int i = 0; i < spacedim; ++i)
1107  for (unsigned int j = 0; j < dim; ++j)
1108  for (unsigned int l = 0; l < dim; ++l)
1109  for (unsigned int m = 0; m < dim; ++m)
1110  for (unsigned int n = 0; n < dim; ++n)
1111  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1112  result[i][j][l][m][n];
1113  }
1114  }
1115  }
1116 
1124  template <int dim, int spacedim, typename VectorType>
1125  void
1127  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1128  const typename ::MappingFEField<dim, spacedim, VectorType>::
1129  InternalData & data,
1130  const FiniteElement<dim, spacedim> &fe,
1131  const ComponentMask & fe_mask,
1132  const std::vector<unsigned int> & fe_to_real,
1133  std::vector<Tensor<5, spacedim>>
1134  &jacobian_pushed_forward_3rd_derivatives)
1135  {
1136  const UpdateFlags update_flags = data.update_each;
1138  {
1139  const unsigned int n_q_points =
1140  jacobian_pushed_forward_3rd_derivatives.size();
1141 
1142  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1143  for (unsigned int point = 0; point < n_q_points; ++point)
1144  {
1145  const Tensor<4, dim> *fourth =
1146  &data.fourth_derivative(point + data_set, 0);
1147 
1149 
1150  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1151  {
1152  const unsigned int comp_k =
1153  fe.system_to_component_index(k).first;
1154  if (fe_mask[comp_k])
1155  for (unsigned int j = 0; j < dim; ++j)
1156  for (unsigned int l = 0; l < dim; ++l)
1157  for (unsigned int m = 0; m < dim; ++m)
1158  for (unsigned int n = 0; n < dim; ++n)
1159  result[fe_to_real[comp_k]][j][l][m][n] +=
1160  (fourth[k][j][l][m][n] *
1161  data.local_dof_values[k]);
1162  }
1163 
1164  // push-forward the j-coordinate
1165  for (unsigned int i = 0; i < spacedim; ++i)
1166  for (unsigned int j = 0; j < spacedim; ++j)
1167  for (unsigned int l = 0; l < dim; ++l)
1168  for (unsigned int m = 0; m < dim; ++m)
1169  for (unsigned int n = 0; n < dim; ++n)
1170  {
1171  tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1172  data.covariant[point][j][0];
1173  for (unsigned int jr = 1; jr < dim; ++jr)
1174  tmp[i][j][l][m][n] +=
1175  result[i][jr][l][m][n] *
1176  data.covariant[point][j][jr];
1177  }
1178 
1179  // push-forward the l-coordinate
1180  for (unsigned int i = 0; i < spacedim; ++i)
1181  for (unsigned int j = 0; j < spacedim; ++j)
1182  for (unsigned int l = 0; l < spacedim; ++l)
1183  for (unsigned int m = 0; m < dim; ++m)
1184  for (unsigned int n = 0; n < dim; ++n)
1185  {
1186  jacobian_pushed_forward_3rd_derivatives
1187  [point][i][j][l][m][n] =
1188  tmp[i][j][0][m][n] *
1189  data.covariant[point][l][0];
1190  for (unsigned int lr = 1; lr < dim; ++lr)
1191  jacobian_pushed_forward_3rd_derivatives[point][i]
1192  [j][l][m]
1193  [n] +=
1194  tmp[i][j][lr][m][n] *
1195  data.covariant[point][l][lr];
1196  }
1197 
1198  // push-forward the m-coordinate
1199  for (unsigned int i = 0; i < spacedim; ++i)
1200  for (unsigned int j = 0; j < spacedim; ++j)
1201  for (unsigned int l = 0; l < spacedim; ++l)
1202  for (unsigned int m = 0; m < spacedim; ++m)
1203  for (unsigned int n = 0; n < dim; ++n)
1204  {
1205  tmp[i][j][l][m][n] =
1206  jacobian_pushed_forward_3rd_derivatives[point][i]
1207  [j][l][0]
1208  [n] *
1209  data.covariant[point][m][0];
1210  for (unsigned int mr = 1; mr < dim; ++mr)
1211  tmp[i][j][l][m][n] +=
1212  jacobian_pushed_forward_3rd_derivatives[point]
1213  [i][j][l]
1214  [mr][n] *
1215  data.covariant[point][m][mr];
1216  }
1217 
1218  // push-forward the n-coordinate
1219  for (unsigned int i = 0; i < spacedim; ++i)
1220  for (unsigned int j = 0; j < spacedim; ++j)
1221  for (unsigned int l = 0; l < spacedim; ++l)
1222  for (unsigned int m = 0; m < spacedim; ++m)
1223  for (unsigned int n = 0; n < spacedim; ++n)
1224  {
1225  jacobian_pushed_forward_3rd_derivatives
1226  [point][i][j][l][m][n] =
1227  tmp[i][j][l][m][0] *
1228  data.covariant[point][n][0];
1229  for (unsigned int nr = 1; nr < dim; ++nr)
1230  jacobian_pushed_forward_3rd_derivatives[point][i]
1231  [j][l][m]
1232  [n] +=
1233  tmp[i][j][l][m][nr] *
1234  data.covariant[point][n][nr];
1235  }
1236  }
1237  }
1238  }
1239 
1240 
1250  template <int dim, int spacedim, typename VectorType>
1251  void
1253  const ::Mapping<dim, spacedim> &mapping,
1254  const typename ::Triangulation<dim, spacedim>::cell_iterator
1255  & cell,
1256  const unsigned int face_no,
1257  const unsigned int subface_no,
1258  const typename QProjector<dim>::DataSetDescriptor data_set,
1259  const typename ::MappingFEField<dim, spacedim, VectorType>::
1260  InternalData &data,
1262  &output_data)
1263  {
1264  const UpdateFlags update_flags = data.update_each;
1265 
1266  if (update_flags & update_boundary_forms)
1267  {
1268  const unsigned int n_q_points = output_data.boundary_forms.size();
1269  if (update_flags & update_normal_vectors)
1270  AssertDimension(output_data.normal_vectors.size(), n_q_points);
1271  if (update_flags & update_JxW_values)
1272  AssertDimension(output_data.JxW_values.size(), n_q_points);
1273 
1274  // map the unit tangentials to the real cell. checking for d!=dim-1
1275  // eliminates compiler warnings regarding unsigned int expressions <
1276  // 0.
1277  for (unsigned int d = 0; d != dim - 1; ++d)
1278  {
1279  Assert(face_no + cell->n_faces() * d <
1280  data.unit_tangentials.size(),
1281  ExcInternalError());
1282  Assert(
1283  data.aux[d].size() <=
1284  data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1285  ExcInternalError());
1286 
1287  mapping.transform(
1289  data.unit_tangentials[face_no + cell->n_faces() * d]),
1291  data,
1292  make_array_view(data.aux[d]));
1293  }
1294 
1295  // if dim==spacedim, we can use the unit tangentials to compute the
1296  // boundary form by simply taking the cross product
1297  if (dim == spacedim)
1298  {
1299  for (unsigned int i = 0; i < n_q_points; ++i)
1300  switch (dim)
1301  {
1302  case 1:
1303  // in 1d, we don't have access to any of the data.aux
1304  // fields (because it has only dim-1 components), but we
1305  // can still compute the boundary form by simply looking
1306  // at the number of the face
1307  output_data.boundary_forms[i][0] =
1308  (face_no == 0 ? -1 : +1);
1309  break;
1310  case 2:
1311  output_data.boundary_forms[i] =
1312  cross_product_2d(data.aux[0][i]);
1313  break;
1314  case 3:
1315  output_data.boundary_forms[i] =
1316  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1317  break;
1318  default:
1319  Assert(false, ExcNotImplemented());
1320  }
1321  }
1322  else //(dim < spacedim)
1323  {
1324  // in the codim-one case, the boundary form results from the
1325  // cross product of all the face tangential vectors and the cell
1326  // normal vector
1327  //
1328  // to compute the cell normal, use the same method used in
1329  // fill_fe_values for cells above
1330  AssertDimension(data.contravariant.size(), n_q_points);
1331 
1332  for (unsigned int point = 0; point < n_q_points; ++point)
1333  {
1334  if (dim == 1)
1335  {
1336  // J is a tangent vector
1337  output_data.boundary_forms[point] =
1338  data.contravariant[point].transpose()[0];
1339  output_data.boundary_forms[point] /=
1340  (face_no == 0 ? -1. : +1.) *
1341  output_data.boundary_forms[point].norm();
1342  }
1343 
1344  if (dim == 2)
1345  {
1347  data.contravariant[point].transpose();
1348 
1349  Tensor<1, spacedim> cell_normal =
1350  cross_product_3d(DX_t[0], DX_t[1]);
1351  cell_normal /= cell_normal.norm();
1352 
1353  // then compute the face normal from the face tangent
1354  // and the cell normal:
1355  output_data.boundary_forms[point] =
1356  cross_product_3d(data.aux[0][point], cell_normal);
1357  }
1358  }
1359  }
1360 
1361  if (update_flags & (update_normal_vectors | update_JxW_values))
1362  for (unsigned int i = 0; i < output_data.boundary_forms.size();
1363  ++i)
1364  {
1365  if (update_flags & update_JxW_values)
1366  {
1367  output_data.JxW_values[i] =
1368  output_data.boundary_forms[i].norm() *
1369  data.quadrature_weights[i + data_set];
1370 
1371  if (subface_no != numbers::invalid_unsigned_int)
1372  {
1373  // TODO
1374  const double area_ratio =
1376  cell->subface_case(face_no), subface_no);
1377  output_data.JxW_values[i] *= area_ratio;
1378  }
1379  }
1380 
1381  if (update_flags & update_normal_vectors)
1382  output_data.normal_vectors[i] =
1383  Point<spacedim>(output_data.boundary_forms[i] /
1384  output_data.boundary_forms[i].norm());
1385  }
1386 
1387  if (update_flags & update_jacobians)
1388  for (unsigned int point = 0; point < n_q_points; ++point)
1389  output_data.jacobians[point] = data.contravariant[point];
1390 
1391  if (update_flags & update_inverse_jacobians)
1392  for (unsigned int point = 0; point < n_q_points; ++point)
1393  output_data.inverse_jacobians[point] =
1394  data.covariant[point].transpose();
1395  }
1396  }
1397 
1404  template <int dim, int spacedim, typename VectorType>
1405  void
1407  const ::Mapping<dim, spacedim> &mapping,
1408  const typename ::Triangulation<dim, spacedim>::cell_iterator
1409  & cell,
1410  const unsigned int face_no,
1411  const unsigned int subface_no,
1412  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1413  const typename ::MappingFEField<dim, spacedim, VectorType>::
1414  InternalData & data,
1415  const FiniteElement<dim, spacedim> &fe,
1416  const ComponentMask & fe_mask,
1417  const std::vector<unsigned int> & fe_to_real,
1419  &output_data)
1420  {
1421  maybe_compute_q_points<dim, spacedim, VectorType>(
1422  data_set,
1423  data,
1424  fe,
1425  fe_mask,
1426  fe_to_real,
1427  output_data.quadrature_points);
1428 
1429  maybe_update_Jacobians<dim, spacedim, VectorType>(
1430  data_set, data, fe, fe_mask, fe_to_real);
1431 
1432  maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1433  data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1434 
1435  maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1436  data_set,
1437  data,
1438  fe,
1439  fe_mask,
1440  fe_to_real,
1441  output_data.jacobian_pushed_forward_grads);
1442 
1443  maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1444  data_set,
1445  data,
1446  fe,
1447  fe_mask,
1448  fe_to_real,
1449  output_data.jacobian_2nd_derivatives);
1450 
1452  spacedim,
1453  VectorType>(
1454  data_set,
1455  data,
1456  fe,
1457  fe_mask,
1458  fe_to_real,
1460 
1461  maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1462  data_set,
1463  data,
1464  fe,
1465  fe_mask,
1466  fe_to_real,
1467  output_data.jacobian_3rd_derivatives);
1468 
1470  spacedim,
1471  VectorType>(
1472  data_set,
1473  data,
1474  fe,
1475  fe_mask,
1476  fe_to_real,
1478 
1479  maybe_compute_face_data<dim, spacedim, VectorType>(
1480  mapping, cell, face_no, subface_no, data_set, data, output_data);
1481  }
1482  } // namespace
1483  } // namespace MappingFEFieldImplementation
1484 } // namespace internal
1485 
1486 
1487 // Note that the CellSimilarity flag is modifiable, since MappingFEField can
1488 // need to recalculate data even when cells are similar.
1489 template <int dim, int spacedim, typename VectorType>
1492  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1494  const Quadrature<dim> & quadrature,
1495  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1497  &output_data) const
1498 {
1499  // convert data object to internal data for this class. fails with an
1500  // exception if that is not possible
1501  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1502  ExcInternalError());
1503  const InternalData &data = static_cast<const InternalData &>(internal_data);
1504 
1505  const unsigned int n_q_points = quadrature.size();
1506 
1507  update_internal_dofs(cell, data);
1508 
1509  internal::MappingFEFieldImplementation::
1510  maybe_compute_q_points<dim, spacedim, VectorType>(
1512  data,
1513  euler_dof_handler->get_fe(),
1514  fe_mask,
1515  fe_to_real,
1516  output_data.quadrature_points);
1517 
1518  internal::MappingFEFieldImplementation::
1519  maybe_update_Jacobians<dim, spacedim, VectorType>(
1521  data,
1522  euler_dof_handler->get_fe(),
1523  fe_mask,
1524  fe_to_real);
1525 
1526  const UpdateFlags update_flags = data.update_each;
1527  const std::vector<double> &weights = quadrature.get_weights();
1528 
1529  // Multiply quadrature weights by absolute value of Jacobian determinants or
1530  // the area element g=sqrt(DX^t DX) in case of codim > 0
1531 
1532  if (update_flags & (update_normal_vectors | update_JxW_values))
1533  {
1534  AssertDimension(output_data.JxW_values.size(), n_q_points);
1535 
1536  Assert(!(update_flags & update_normal_vectors) ||
1537  (output_data.normal_vectors.size() == n_q_points),
1538  ExcDimensionMismatch(output_data.normal_vectors.size(),
1539  n_q_points));
1540 
1541 
1542  for (unsigned int point = 0; point < n_q_points; ++point)
1543  {
1544  if (dim == spacedim)
1545  {
1546  const double det = data.contravariant[point].determinant();
1547 
1548  // check for distorted cells.
1549 
1550  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1551  // 1e12 in 2D. might want to find a finer
1552  // (dimension-independent) criterion
1553  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1554  cell->diameter() / std::sqrt(double(dim))),
1556  cell->center(), det, point)));
1557  output_data.JxW_values[point] = weights[point] * det;
1558  }
1559  // if dim==spacedim, then there is no cell normal to
1560  // compute. since this is for FEValues (and not FEFaceValues),
1561  // there are also no face normals to compute
1562  else // codim>0 case
1563  {
1564  Tensor<1, spacedim> DX_t[dim];
1565  for (unsigned int i = 0; i < spacedim; ++i)
1566  for (unsigned int j = 0; j < dim; ++j)
1567  DX_t[j][i] = data.contravariant[point][i][j];
1568 
1569  Tensor<2, dim> G; // First fundamental form
1570  for (unsigned int i = 0; i < dim; ++i)
1571  for (unsigned int j = 0; j < dim; ++j)
1572  G[i][j] = DX_t[i] * DX_t[j];
1573 
1574  output_data.JxW_values[point] =
1575  std::sqrt(determinant(G)) * weights[point];
1576 
1577  if (update_flags & update_normal_vectors)
1578  {
1579  Assert(spacedim - dim == 1,
1580  ExcMessage("There is no cell normal in codim 2."));
1581 
1582  if (dim == 1)
1583  output_data.normal_vectors[point] =
1584  cross_product_2d(-DX_t[0]);
1585  else
1586  {
1587  Assert(dim == 2, ExcInternalError());
1588 
1589  // dim-1==1 for the second argument, but this
1590  // avoids a compiler warning about array bounds:
1591  output_data.normal_vectors[point] =
1592  cross_product_3d(DX_t[0], DX_t[dim - 1]);
1593  }
1594 
1595  output_data.normal_vectors[point] /=
1596  output_data.normal_vectors[point].norm();
1597 
1598  if (cell->direction_flag() == false)
1599  output_data.normal_vectors[point] *= -1.;
1600  }
1601  } // codim>0 case
1602  }
1603  }
1604 
1605  // copy values from InternalData to vector given by reference
1606  if (update_flags & update_jacobians)
1607  {
1608  AssertDimension(output_data.jacobians.size(), n_q_points);
1609  for (unsigned int point = 0; point < n_q_points; ++point)
1610  output_data.jacobians[point] = data.contravariant[point];
1611  }
1612 
1613  // copy values from InternalData to vector given by reference
1614  if (update_flags & update_inverse_jacobians)
1615  {
1616  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1617  for (unsigned int point = 0; point < n_q_points; ++point)
1618  output_data.inverse_jacobians[point] =
1619  data.covariant[point].transpose();
1620  }
1621 
1622  // calculate derivatives of the Jacobians
1623  internal::MappingFEFieldImplementation::
1624  maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1626  data,
1627  euler_dof_handler->get_fe(),
1628  fe_mask,
1629  fe_to_real,
1630  output_data.jacobian_grads);
1631 
1632  // calculate derivatives of the Jacobians pushed forward to real cell
1633  // coordinates
1634  internal::MappingFEFieldImplementation::
1635  maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1637  data,
1638  euler_dof_handler->get_fe(),
1639  fe_mask,
1640  fe_to_real,
1641  output_data.jacobian_pushed_forward_grads);
1642 
1643  // calculate hessians of the Jacobians
1644  internal::MappingFEFieldImplementation::
1645  maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1647  data,
1648  euler_dof_handler->get_fe(),
1649  fe_mask,
1650  fe_to_real,
1651  output_data.jacobian_2nd_derivatives);
1652 
1653  // calculate hessians of the Jacobians pushed forward to real cell coordinates
1654  internal::MappingFEFieldImplementation::
1655  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1656  spacedim,
1657  VectorType>(
1659  data,
1660  euler_dof_handler->get_fe(),
1661  fe_mask,
1662  fe_to_real,
1664 
1665  // calculate gradients of the hessians of the Jacobians
1666  internal::MappingFEFieldImplementation::
1667  maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1669  data,
1670  euler_dof_handler->get_fe(),
1671  fe_mask,
1672  fe_to_real,
1673  output_data.jacobian_3rd_derivatives);
1674 
1675  // calculate gradients of the hessians of the Jacobians pushed forward to real
1676  // cell coordinates
1677  internal::MappingFEFieldImplementation::
1678  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1679  spacedim,
1680  VectorType>(
1682  data,
1683  euler_dof_handler->get_fe(),
1684  fe_mask,
1685  fe_to_real,
1687 
1689 }
1690 
1691 
1692 
1693 template <int dim, int spacedim, typename VectorType>
1694 void
1696  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1697  const unsigned int face_no,
1698  const hp::QCollection<dim - 1> & quadrature,
1699  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1701  &output_data) const
1702 {
1703  AssertDimension(quadrature.size(), 1);
1704 
1705  // convert data object to internal data for this class. fails with an
1706  // exception if that is not possible
1707  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1708  ExcInternalError());
1709  const InternalData &data = static_cast<const InternalData &>(internal_data);
1710 
1711  update_internal_dofs(cell, data);
1712 
1713  internal::MappingFEFieldImplementation::
1714  do_fill_fe_face_values<dim, spacedim, VectorType>(
1715  *this,
1716  cell,
1717  face_no,
1720  face_no,
1721  cell->face_orientation(face_no),
1722  cell->face_flip(face_no),
1723  cell->face_rotation(face_no),
1724  quadrature[0].size()),
1725  data,
1726  euler_dof_handler->get_fe(),
1727  fe_mask,
1728  fe_to_real,
1729  output_data);
1730 }
1731 
1732 
1733 template <int dim, int spacedim, typename VectorType>
1734 void
1736  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1737  const unsigned int face_no,
1738  const unsigned int subface_no,
1739  const Quadrature<dim - 1> & quadrature,
1740  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1742  &output_data) const
1743 {
1744  // convert data object to internal data for this class. fails with an
1745  // exception if that is not possible
1746  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1747  ExcInternalError());
1748  const InternalData &data = static_cast<const InternalData &>(internal_data);
1749 
1750  update_internal_dofs(cell, data);
1751 
1753  spacedim,
1754  VectorType>(
1755  *this,
1756  cell,
1757  face_no,
1760  face_no,
1761  subface_no,
1762  cell->face_orientation(face_no),
1763  cell->face_flip(face_no),
1764  cell->face_rotation(face_no),
1765  quadrature.size(),
1766  cell->subface_case(face_no)),
1767  data,
1768  euler_dof_handler->get_fe(),
1769  fe_mask,
1770  fe_to_real,
1771  output_data);
1772 }
1773 
1774 
1775 namespace internal
1776 {
1777  namespace MappingFEFieldImplementation
1778  {
1779  namespace
1780  {
1781  template <int dim, int spacedim, int rank, typename VectorType>
1782  void
1784  const ArrayView<const Tensor<rank, dim>> & input,
1785  const MappingKind mapping_kind,
1786  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1787  const ArrayView<Tensor<rank, spacedim>> & output)
1788  {
1789  AssertDimension(input.size(), output.size());
1790  Assert((dynamic_cast<
1791  const typename ::
1792  MappingFEField<dim, spacedim, VectorType>::InternalData *>(
1793  &mapping_data) != nullptr),
1794  ExcInternalError());
1795  const typename ::MappingFEField<dim, spacedim, VectorType>::
1796  InternalData &data = static_cast<
1797  const typename ::MappingFEField<dim, spacedim, VectorType>::
1798  InternalData &>(mapping_data);
1799 
1800  switch (mapping_kind)
1801  {
1802  case mapping_contravariant:
1803  {
1804  Assert(
1805  data.update_each & update_contravariant_transformation,
1807  "update_contravariant_transformation"));
1808 
1809  for (unsigned int i = 0; i < output.size(); ++i)
1810  output[i] =
1811  apply_transformation(data.contravariant[i], input[i]);
1812 
1813  return;
1814  }
1815 
1816  case mapping_piola:
1817  {
1818  Assert(
1819  data.update_each & update_contravariant_transformation,
1821  "update_contravariant_transformation"));
1822  Assert(
1823  data.update_each & update_volume_elements,
1825  "update_volume_elements"));
1826  Assert(rank == 1, ExcMessage("Only for rank 1"));
1827  for (unsigned int i = 0; i < output.size(); ++i)
1828  {
1829  output[i] =
1830  apply_transformation(data.contravariant[i], input[i]);
1831  output[i] /= data.volume_elements[i];
1832  }
1833  return;
1834  }
1835 
1836 
1837  // We still allow this operation as in the
1838  // reference cell Derivatives are Tensor
1839  // rather than DerivativeForm
1840  case mapping_covariant:
1841  {
1842  Assert(
1843  data.update_each & update_contravariant_transformation,
1845  "update_contravariant_transformation"));
1846 
1847  for (unsigned int i = 0; i < output.size(); ++i)
1848  output[i] = apply_transformation(data.covariant[i], input[i]);
1849 
1850  return;
1851  }
1852 
1853  default:
1854  Assert(false, ExcNotImplemented());
1855  }
1856  }
1857 
1858 
1859  template <int dim, int spacedim, int rank, typename VectorType>
1860  void
1862  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
1863  const MappingKind mapping_kind,
1864  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1865  const ArrayView<Tensor<rank + 1, spacedim>> & output)
1866  {
1867  AssertDimension(input.size(), output.size());
1868  Assert((dynamic_cast<
1869  const typename ::
1870  MappingFEField<dim, spacedim, VectorType>::InternalData *>(
1871  &mapping_data) != nullptr),
1872  ExcInternalError());
1873  const typename ::MappingFEField<dim, spacedim, VectorType>::
1874  InternalData &data = static_cast<
1875  const typename ::MappingFEField<dim, spacedim, VectorType>::
1876  InternalData &>(mapping_data);
1877 
1878  switch (mapping_kind)
1879  {
1880  case mapping_covariant:
1881  {
1882  Assert(
1883  data.update_each & update_contravariant_transformation,
1885  "update_contravariant_transformation"));
1886 
1887  for (unsigned int i = 0; i < output.size(); ++i)
1888  output[i] = apply_transformation(data.covariant[i], input[i]);
1889 
1890  return;
1891  }
1892  default:
1893  Assert(false, ExcNotImplemented());
1894  }
1895  }
1896  } // namespace
1897  } // namespace MappingFEFieldImplementation
1898 } // namespace internal
1899 
1900 
1901 
1902 template <int dim, int spacedim, typename VectorType>
1903 void
1905  const ArrayView<const Tensor<1, dim>> & input,
1906  const MappingKind mapping_kind,
1907  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1908  const ArrayView<Tensor<1, spacedim>> & output) const
1909 {
1910  AssertDimension(input.size(), output.size());
1911 
1912  internal::MappingFEFieldImplementation::
1913  transform_fields<dim, spacedim, 1, VectorType>(input,
1914  mapping_kind,
1915  mapping_data,
1916  output);
1917 }
1918 
1919 
1920 
1921 template <int dim, int spacedim, typename VectorType>
1922 void
1924  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1925  const MappingKind mapping_kind,
1926  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1927  const ArrayView<Tensor<2, spacedim>> & output) const
1928 {
1929  AssertDimension(input.size(), output.size());
1930 
1931  internal::MappingFEFieldImplementation::
1932  transform_differential_forms<dim, spacedim, 1, VectorType>(input,
1933  mapping_kind,
1934  mapping_data,
1935  output);
1936 }
1937 
1938 
1939 
1940 template <int dim, int spacedim, typename VectorType>
1941 void
1943  const ArrayView<const Tensor<2, dim>> &input,
1944  const MappingKind,
1945  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1946  const ArrayView<Tensor<2, spacedim>> & output) const
1947 {
1948  (void)input;
1949  (void)output;
1950  (void)mapping_data;
1951  AssertDimension(input.size(), output.size());
1952 
1953  AssertThrow(false, ExcNotImplemented());
1954 }
1955 
1956 
1957 
1958 template <int dim, int spacedim, typename VectorType>
1959 void
1961  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1962  const MappingKind mapping_kind,
1963  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1964  const ArrayView<Tensor<3, spacedim>> & output) const
1965 {
1966  AssertDimension(input.size(), output.size());
1967  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1968  ExcInternalError());
1969  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1970 
1971  switch (mapping_kind)
1972  {
1974  {
1975  Assert(data.update_each & update_contravariant_transformation,
1977  "update_covariant_transformation"));
1978 
1979  for (unsigned int q = 0; q < output.size(); ++q)
1980  for (unsigned int i = 0; i < spacedim; ++i)
1981  for (unsigned int j = 0; j < spacedim; ++j)
1982  for (unsigned int k = 0; k < spacedim; ++k)
1983  {
1984  output[q][i][j][k] = data.covariant[q][j][0] *
1985  data.covariant[q][k][0] *
1986  input[q][i][0][0];
1987  for (unsigned int J = 0; J < dim; ++J)
1988  {
1989  const unsigned int K0 = (0 == J) ? 1 : 0;
1990  for (unsigned int K = K0; K < dim; ++K)
1991  output[q][i][j][k] += data.covariant[q][j][J] *
1992  data.covariant[q][k][K] *
1993  input[q][i][J][K];
1994  }
1995  }
1996  return;
1997  }
1998 
1999  default:
2000  Assert(false, ExcNotImplemented());
2001  }
2002 }
2003 
2004 
2005 
2006 template <int dim, int spacedim, typename VectorType>
2007 void
2009  const ArrayView<const Tensor<3, dim>> &input,
2010  const MappingKind /*mapping_kind*/,
2011  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2012  const ArrayView<Tensor<3, spacedim>> & output) const
2013 {
2014  (void)input;
2015  (void)output;
2016  (void)mapping_data;
2017  AssertDimension(input.size(), output.size());
2018 
2019  AssertThrow(false, ExcNotImplemented());
2020 }
2021 
2022 
2023 
2024 template <int dim, int spacedim, typename VectorType>
2027  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2028  const Point<dim> & p) const
2029 {
2030  // Use the get_data function to create an InternalData with data vectors of
2031  // the right size and transformation shape values already computed at point
2032  // p.
2033  const Quadrature<dim> point_quadrature(p);
2034  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2035  get_data(update_quadrature_points | update_jacobians, point_quadrature));
2036  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2037  ExcInternalError());
2038 
2039  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2040 
2041  return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
2042 }
2043 
2044 
2045 template <int dim, int spacedim, typename VectorType>
2048  const InternalData &data) const
2049 {
2050  Point<spacedim> p_real;
2051 
2052  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2053  {
2054  unsigned int comp_i =
2055  euler_dof_handler->get_fe().system_to_component_index(i).first;
2056  if (fe_mask[comp_i])
2057  p_real[fe_to_real[comp_i]] +=
2058  data.local_dof_values[i] * data.shape(0, i);
2059  }
2060 
2061  return p_real;
2062 }
2063 
2064 
2065 
2066 template <int dim, int spacedim, typename VectorType>
2067 Point<dim>
2069  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2070  const Point<spacedim> & p) const
2071 {
2072  // first a Newton iteration based on the real mapping. It uses the center
2073  // point of the cell as a starting point
2074  Point<dim> initial_p_unit;
2075  try
2076  {
2077  initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
2078  .transform_real_to_unit_cell(cell, p);
2079  }
2080  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2081  {
2082  // mirror the conditions of the code below to determine if we need to
2083  // use an arbitrary starting point or if we just need to rethrow the
2084  // exception
2085  for (unsigned int d = 0; d < dim; ++d)
2086  initial_p_unit[d] = 0.5;
2087  }
2088 
2089  // TODO
2090  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2091 
2092  // for (unsigned int d=0; d<dim; ++d)
2093  // initial_p_unit[d] = 0.;
2094 
2095  const Quadrature<dim> point_quadrature(initial_p_unit);
2096 
2098  if (spacedim > dim)
2099  update_flags |= update_jacobian_grads;
2100  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2101  get_data(update_flags, point_quadrature));
2102  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2103  ExcInternalError());
2104 
2105  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2106 
2107  return do_transform_real_to_unit_cell(cell,
2108  p,
2109  initial_p_unit,
2110  dynamic_cast<InternalData &>(*mdata));
2111 }
2112 
2113 
2114 template <int dim, int spacedim, typename VectorType>
2115 Point<dim>
2117  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2118  const Point<spacedim> & p,
2119  const Point<dim> & initial_p_unit,
2120  InternalData & mdata) const
2121 {
2122  const unsigned int n_shapes = mdata.shape_values.size();
2123  (void)n_shapes;
2124  Assert(n_shapes != 0, ExcInternalError());
2125  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2126 
2127 
2128  // Newton iteration to solve
2129  // f(x)=p(x)-p=0
2130  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2131  // The start value was set to be the
2132  // linear approximation to the cell
2133  // The shape values and derivatives
2134  // of the mapping at this point are
2135  // previously computed.
2136  // f(x)
2137  Point<dim> p_unit = initial_p_unit;
2138  Point<dim> f;
2139  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
2141  Tensor<1, spacedim> p_minus_F = p - p_real;
2142  const double eps = 1.e-12 * cell->diameter();
2143  const unsigned int newton_iteration_limit = 20;
2144  unsigned int newton_iteration = 0;
2145  while (p_minus_F.norm_square() > eps * eps)
2146  {
2147  // f'(x)
2148  Point<spacedim> DF[dim];
2149  Tensor<2, dim> df;
2150  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2151  {
2152  const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2153  const unsigned int comp_k =
2154  euler_dof_handler->get_fe().system_to_component_index(k).first;
2155  if (fe_mask[comp_k])
2156  for (unsigned int j = 0; j < dim; ++j)
2157  DF[j][fe_to_real[comp_k]] +=
2158  mdata.local_dof_values[k] * grad_k[j];
2159  }
2160  for (unsigned int j = 0; j < dim; ++j)
2161  {
2162  f[j] = DF[j] * p_minus_F;
2163  for (unsigned int l = 0; l < dim; ++l)
2164  df[j][l] = -DF[j] * DF[l];
2165  }
2166  // Solve [f'(x)]d=f(x)
2167  const Tensor<1, dim> delta =
2168  invert(df) * static_cast<const Tensor<1, dim> &>(f);
2169  // do a line search
2170  double step_length = 1;
2171  do
2172  {
2173  // update of p_unit. The
2174  // spacedimth component of
2175  // transformed point is simply
2176  // ignored in codimension one
2177  // case. When this component is
2178  // not zero, then we are
2179  // projecting the point to the
2180  // surface or curve identified
2181  // by the cell.
2182  Point<dim> p_unit_trial = p_unit;
2183  for (unsigned int i = 0; i < dim; ++i)
2184  p_unit_trial[i] -= step_length * delta[i];
2185  // shape values and derivatives
2186  // at new p_unit point
2187  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
2188  mdata);
2189  // f(x)
2190  Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
2191  const Tensor<1, spacedim> f_trial = p - p_real_trial;
2192  // see if we are making progress with the current step length
2193  // and if not, reduce it by a factor of two and try again
2194  if (f_trial.norm() < p_minus_F.norm())
2195  {
2196  p_real = p_real_trial;
2197  p_unit = p_unit_trial;
2198  p_minus_F = f_trial;
2199  break;
2200  }
2201  else if (step_length > 0.05)
2202  step_length /= 2;
2203  else
2204  goto failure;
2205  }
2206  while (true);
2207  ++newton_iteration;
2208  if (newton_iteration > newton_iteration_limit)
2209  goto failure;
2210  }
2211  return p_unit;
2212  // if we get to the following label, then we have either run out
2213  // of Newton iterations, or the line search has not converged.
2214  // in either case, we need to give up, so throw an exception that
2215  // can then be caught
2216 failure:
2217  AssertThrow(false,
2219  // ...the compiler wants us to return something, though we can
2220  // of course never get here...
2221  return {};
2222 }
2223 
2224 
2225 template <int dim, int spacedim, typename VectorType>
2226 unsigned int
2228 {
2229  return euler_dof_handler->get_fe().degree;
2230 }
2231 
2232 
2233 
2234 template <int dim, int spacedim, typename VectorType>
2237 {
2238  return this->fe_mask;
2239 }
2240 
2241 
2242 template <int dim, int spacedim, typename VectorType>
2243 std::unique_ptr<Mapping<dim, spacedim>>
2245 {
2246  return std::make_unique<MappingFEField<dim, spacedim, VectorType>>(*this);
2247 }
2248 
2249 
2250 template <int dim, int spacedim, typename VectorType>
2251 void
2253  const typename Triangulation<dim, spacedim>::cell_iterator & cell,
2255  const
2256 {
2257  Assert(euler_dof_handler != nullptr,
2258  ExcMessage("euler_dof_handler is empty"));
2259 
2260  typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
2262  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2263  if (uses_level_dofs)
2264  {
2265  AssertIndexRange(cell->level(), euler_vector.size());
2266  AssertDimension(euler_vector[cell->level()]->size(),
2267  euler_dof_handler->n_dofs(cell->level()));
2268  }
2269  else
2270  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2271 
2272  if (uses_level_dofs)
2273  dof_cell->get_mg_dof_indices(data.local_dof_indices);
2274  else
2275  dof_cell->get_dof_indices(data.local_dof_indices);
2276 
2277  const VectorType &vector =
2278  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2279 
2280  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2281  data.local_dof_values[i] =
2283  data.local_dof_indices[i]);
2284 }
2285 
2286 // explicit instantiations
2287 #define SPLIT_INSTANTIATIONS_COUNT 2
2288 #ifndef SPLIT_INSTANTIATIONS_INDEX
2289 # define SPLIT_INSTANTIATIONS_INDEX 0
2290 #endif
2291 #include "mapping_fe_field.inst"
2292 
2293 
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:699
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
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:311
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
Definition: point.h:111
static DataSetDescriptor cell()
Definition: qprojector.h:361
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< Point< spacedim > > quadrature_points
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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:4606
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInactiveCell()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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:72
@ mapping_piola
Definition: mapping.h:107
@ mapping_covariant_gradient
Definition: mapping.h:93
@ mapping_covariant
Definition: mapping.h:82
@ mapping_contravariant
Definition: mapping.h:87
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:260
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:190
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:493
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:482
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_pushed_forward_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 5, spacedim >> &jacobian_pushed_forward_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 4, spacedim >> &jacobian_pushed_forward_2nd_derivatives)
void maybe_compute_q_points(const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Point< spacedim >> &quadrature_points)
void maybe_update_Jacobians(const CellSimilarity::Similarity cell_similarity, const typename ::QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data)
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, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim >> &jacobian_grads)
void maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim >> &jacobian_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim >> &jacobian_pushed_forward_grads)
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)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim >> &jacobian_2nd_derivatives)
static const unsigned int invalid_unsigned_int
Definition: types.h:201
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 Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)