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