Reference documentation for deal.II version Git 08727cc441 2020-07-02 15:45:42 -0400
\(\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 - 2020 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, typename DoFHandlerType>
64  const ComponentMask & mask)
65  : unit_tangentials()
66  , n_shape_functions(fe.dofs_per_cell)
67  , mask(mask)
68  , local_dof_indices(fe.dofs_per_cell)
69  , local_dof_values(fe.dofs_per_cell)
70 {}
71 
72 
73 
74 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
75 std::size_t
78 {
79  Assert(false, ExcNotImplemented());
80  return 0;
81 }
82 
83 
84 
85 template <int dim, int spacedim, typename DoFHandlerType, 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 DoFHandlerType, typename VectorType>
97 const Tensor<1, dim> &
99  derivative(const unsigned int qpoint, const unsigned int shape_nr) const
100 {
101  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
102  shape_derivatives.size());
103  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
104 }
105 
106 
107 
108 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
111  derivative(const unsigned int qpoint, 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 DoFHandlerType, typename VectorType>
120 const Tensor<2, dim> &
122  second_derivative(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 DoFHandlerType, typename VectorType>
135  second_derivative(const unsigned int qpoint, const unsigned int shape_nr)
136 {
137  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
138  shape_second_derivatives.size());
139  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
140 }
141 
142 
143 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
144 const Tensor<3, dim> &
146  third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
147 {
148  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
149  shape_third_derivatives.size());
150  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
151 }
152 
153 
154 
155 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
158  third_derivative(const unsigned int qpoint, const unsigned int shape_nr)
159 {
160  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
161  shape_third_derivatives.size());
162  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
163 }
164 
165 
166 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
167 const Tensor<4, dim> &
169  fourth_derivative(const unsigned int qpoint,
170  const unsigned int shape_nr) const
171 {
172  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
173  shape_fourth_derivatives.size());
174  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
175 }
176 
177 
178 
179 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
182  fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr)
183 {
184  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
185  shape_fourth_derivatives.size());
186  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
187 }
188 
189 
190 
191 namespace
192 {
193  // Construct a quadrature formula containing the vertices of the reference
194  // cell in dimension dim (with invalid weights)
195  template <int dim>
197  get_vertex_quadrature()
198  {
199  static Quadrature<dim> quad;
200  if (quad.size() == 0)
201  {
202  std::vector<Point<dim>> points(GeometryInfo<dim>::vertices_per_cell);
203  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
205  quad = Quadrature<dim>(points);
206  }
207  return quad;
208  }
209 } // namespace
210 
211 
212 
213 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
215  const DoFHandlerType &euler_dof_handler,
216  const VectorType & euler_vector,
217  const ComponentMask & mask)
218  : uses_level_dofs(false)
219  , euler_vector({&euler_vector})
221  , fe_mask(mask.size() ?
222  mask :
224  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
225  true))
227  , fe_values(this->euler_dof_handler->get_fe(),
228  get_vertex_quadrature<dim>(),
230 {
231  unsigned int size = 0;
232  for (unsigned int i = 0; i < fe_mask.size(); ++i)
233  {
234  if (fe_mask[i])
235  fe_to_real[i] = size++;
236  }
237  AssertDimension(size, spacedim);
238 }
239 
240 
241 
242 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
244  const DoFHandlerType & euler_dof_handler,
245  const std::vector<VectorType> &euler_vector,
246  const ComponentMask & mask)
247  : uses_level_dofs(true)
248  , euler_dof_handler(&euler_dof_handler)
249  , fe_mask(mask.size() ?
250  mask :
252  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
253  true))
255  , fe_values(this->euler_dof_handler->get_fe(),
256  get_vertex_quadrature<dim>(),
258 {
259  unsigned int size = 0;
260  for (unsigned int i = 0; i < fe_mask.size(); ++i)
261  {
262  if (fe_mask[i])
263  fe_to_real[i] = size++;
264  }
265  AssertDimension(size, spacedim);
266 
267  Assert(euler_dof_handler.has_level_dofs(),
268  ExcMessage("The underlying DoFHandler object did not call "
269  "distribute_mg_dofs(). In this case, the construction via "
270  "level vectors does not make sense."));
271  AssertDimension(euler_vector.size(),
272  euler_dof_handler.get_triangulation().n_global_levels());
273  this->euler_vector.clear();
274  this->euler_vector.resize(euler_vector.size());
275  for (unsigned int i = 0; i < euler_vector.size(); ++i)
276  this->euler_vector[i] = &euler_vector[i];
277 }
278 
279 
280 
281 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
283  const DoFHandlerType & euler_dof_handler,
285  const ComponentMask & mask)
286  : uses_level_dofs(true)
287  , euler_dof_handler(&euler_dof_handler)
288  , fe_mask(mask.size() ?
289  mask :
291  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
292  true))
294  , fe_values(this->euler_dof_handler->get_fe(),
295  get_vertex_quadrature<dim>(),
297 {
298  unsigned int size = 0;
299  for (unsigned int i = 0; i < fe_mask.size(); ++i)
300  {
301  if (fe_mask[i])
302  fe_to_real[i] = size++;
303  }
304  AssertDimension(size, spacedim);
305 
306  Assert(euler_dof_handler.has_level_dofs(),
307  ExcMessage("The underlying DoFHandler object did not call "
308  "distribute_mg_dofs(). In this case, the construction via "
309  "level vectors does not make sense."));
310  AssertDimension(euler_vector.max_level() + 1,
311  euler_dof_handler.get_triangulation().n_global_levels());
312  this->euler_vector.clear();
313  this->euler_vector.resize(
314  euler_dof_handler.get_triangulation().n_global_levels());
315  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
316  ++i)
317  this->euler_vector[i] = &euler_vector[i];
318 }
319 
320 
321 
322 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
326  , euler_vector(mapping.euler_vector)
328  , fe_mask(mapping.fe_mask)
329  , fe_to_real(mapping.fe_to_real)
330  , fe_values(mapping.euler_dof_handler->get_fe(),
331  get_vertex_quadrature<dim>(),
333 {}
334 
335 
336 
337 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
338 inline const double &
340  const unsigned int qpoint,
341  const unsigned int shape_nr) const
342 {
343  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
344  return shape_values[qpoint * n_shape_functions + shape_nr];
345 }
346 
347 
348 
349 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
350 bool
353 {
354  return false;
355 }
356 
357 
358 
359 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
360 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
362  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
363 {
364  // we transform our tria iterator into a dof iterator so we can access
365  // data not associated with triangulations
366  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
367  *cell, euler_dof_handler);
368 
369  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
371  fe_values.n_quadrature_points);
373  euler_dof_handler->get_fe().n_components());
374  if (uses_level_dofs)
375  {
376  AssertIndexRange(cell->level(), euler_vector.size());
377  AssertDimension(euler_vector[cell->level()]->size(),
378  euler_dof_handler->n_dofs(cell->level()));
379  }
380  else
381  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
382 
383  {
384  std::lock_guard<std::mutex> lock(fe_values_mutex);
385  fe_values.reinit(dof_cell);
386  }
387  const unsigned int dofs_per_cell = euler_dof_handler->get_fe().dofs_per_cell;
388  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
389  if (uses_level_dofs)
390  dof_cell->get_mg_dof_indices(dof_indices);
391  else
392  dof_cell->get_dof_indices(dof_indices);
393 
394  const VectorType &vector =
395  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
396 
397  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
398  for (unsigned int i = 0; i < dofs_per_cell; ++i)
399  {
400  const unsigned int comp = fe_to_real
401  [euler_dof_handler->get_fe().system_to_component_index(i).first];
402  if (comp != numbers::invalid_unsigned_int)
403  {
404  typename VectorType::value_type value =
405  internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
406  if (euler_dof_handler->get_fe().is_primitive(i))
407  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
408  vertices[v][comp] += fe_values.shape_value(i, v) * value;
409  else
410  Assert(false, ExcNotImplemented());
411  }
412  }
413 
414  return vertices;
415 }
416 
417 
418 
419 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
420 void
423  const std::vector<Point<dim>> &unit_points,
425  InternalData &data) const
426 {
427  const auto fe = &euler_dof_handler->get_fe();
428  const unsigned int n_points = unit_points.size();
429 
430  for (unsigned int point = 0; point < n_points; ++point)
431  {
432  if (data.shape_values.size() != 0)
433  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
434  data.shape(point, i) = fe->shape_value(i, unit_points[point]);
435 
436  if (data.shape_derivatives.size() != 0)
437  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
438  data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
439 
440  if (data.shape_second_derivatives.size() != 0)
441  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
442  data.second_derivative(point, i) =
443  fe->shape_grad_grad(i, unit_points[point]);
444 
445  if (data.shape_third_derivatives.size() != 0)
446  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
447  data.third_derivative(point, i) =
448  fe->shape_3rd_derivative(i, unit_points[point]);
449 
450  if (data.shape_fourth_derivatives.size() != 0)
451  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
452  data.fourth_derivative(point, i) =
453  fe->shape_4th_derivative(i, unit_points[point]);
454  }
455 }
456 
457 
458 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
462 {
463  // add flags if the respective quantities are necessary to compute
464  // what we need. note that some flags appear in both conditions and
465  // in subsequent set operations. this leads to some circular
466  // logic. the only way to treat this is to iterate. since there are
467  // 5 if-clauses in the loop, it will take at most 4 iterations to
468  // converge. do them:
469  UpdateFlags out = in;
470  for (unsigned int i = 0; i < 5; ++i)
471  {
472  // The following is a little incorrect:
473  // If not applied on a face,
474  // update_boundary_forms does not
475  // make sense. On the other hand,
476  // it is necessary on a
477  // face. Currently,
478  // update_boundary_forms is simply
479  // ignored for the interior of a
480  // cell.
482  out |= update_boundary_forms;
483 
488 
489  if (out &
494 
495  // The contravariant transformation
496  // is a Piola transformation, which
497  // requires the determinant of the
498  // Jacobi matrix of the transformation.
499  // Therefore these values have to be
500  // updated for each cell.
502  out |= update_JxW_values;
503 
504  if (out & update_normal_vectors)
505  out |= update_JxW_values;
506  }
507 
508  return out;
509 }
510 
511 
512 
513 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
514 void
516  const UpdateFlags update_flags,
517  const Quadrature<dim> &q,
518  const unsigned int n_original_q_points,
519  InternalData & data) const
520 {
521  // store the flags in the internal data object so we can access them
522  // in fill_fe_*_values(). use the transitive hull of the required
523  // flags
524  data.update_each = requires_update_flags(update_flags);
525 
526  const unsigned int n_q_points = q.size();
527 
528  // see if we need the (transformation) shape function values
529  // and/or gradients and resize the necessary arrays
531  data.shape_values.resize(data.n_shape_functions * n_q_points);
532 
533  if (data.update_each &
537  data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
538 
540  data.covariant.resize(n_original_q_points);
541 
543  data.contravariant.resize(n_original_q_points);
544 
546  data.volume_elements.resize(n_original_q_points);
547 
548  if (data.update_each &
550  data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
551 
554  data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
555 
558  data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
559 
561 }
562 
563 
564 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
565 void
567  const UpdateFlags update_flags,
568  const Quadrature<dim> &q,
569  const unsigned int n_original_q_points,
570  InternalData & data) const
571 {
572  compute_data(update_flags, q, n_original_q_points, data);
573 
574  if (dim > 1)
575  {
577  {
578  data.aux.resize(
579  dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
580 
581  // Compute tangentials to the unit cell.
582  for (const unsigned int i : GeometryInfo<dim>::face_indices())
583  {
584  data.unit_tangentials[i].resize(n_original_q_points);
585  std::fill(data.unit_tangentials[i].begin(),
586  data.unit_tangentials[i].end(),
588  if (dim > 2)
589  {
591  .resize(n_original_q_points);
592  std::fill(
594  .begin(),
596  .end(),
598  }
599  }
600  }
601  }
602 }
603 
604 
605 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
606 typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
608  const UpdateFlags update_flags,
609  const Quadrature<dim> &quadrature) const
610 {
611  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
612  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
613  auto &data = dynamic_cast<InternalData &>(*data_ptr);
614  this->compute_data(update_flags, quadrature, quadrature.size(), data);
615 
616  return data_ptr;
617 }
618 
619 
620 
621 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
622 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
624  const UpdateFlags update_flags,
625  const Quadrature<dim - 1> &quadrature) const
626 {
627  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
628  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
629  auto & data = dynamic_cast<InternalData &>(*data_ptr);
631  this->compute_face_data(update_flags, q, quadrature.size(), data);
632 
633  return data_ptr;
634 }
635 
636 
637 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
638 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
640  const UpdateFlags update_flags,
641  const Quadrature<dim - 1> &quadrature) const
642 {
643  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
644  std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
645  auto & data = dynamic_cast<InternalData &>(*data_ptr);
647  this->compute_face_data(update_flags, q, quadrature.size(), data);
648 
649  return data_ptr;
650 }
651 
652 
653 
654 namespace internal
655 {
656  namespace MappingFEFieldImplementation
657  {
658  namespace
659  {
666  template <int dim,
667  int spacedim,
668  typename VectorType,
669  typename DoFHandlerType>
670  void
671  maybe_compute_q_points(
672  const typename ::QProjector<dim>::DataSetDescriptor data_set,
673  const typename ::
675  InternalData & data,
677  const ComponentMask & fe_mask,
678  const std::vector<unsigned int> & fe_to_real,
679  std::vector<Point<spacedim>> & quadrature_points)
680  {
681  const UpdateFlags update_flags = data.update_each;
682 
683  if (update_flags & update_quadrature_points)
684  {
685  for (unsigned int point = 0; point < quadrature_points.size();
686  ++point)
687  {
688  Point<spacedim> result;
689  const double * shape = &data.shape(point + data_set, 0);
690 
691  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
692  {
693  const unsigned int comp_k =
694  fe.system_to_component_index(k).first;
695  if (fe_mask[comp_k])
696  result[fe_to_real[comp_k]] +=
697  data.local_dof_values[k] * shape[k];
698  }
699 
700  quadrature_points[point] = result;
701  }
702  }
703  }
704 
712  template <int dim,
713  int spacedim,
714  typename VectorType,
715  typename DoFHandlerType>
716  void
717  maybe_update_Jacobians(
718  const typename ::QProjector<dim>::DataSetDescriptor data_set,
719  const typename ::
721  InternalData & data,
723  const ComponentMask & fe_mask,
724  const std::vector<unsigned int> & fe_to_real)
725  {
726  const UpdateFlags update_flags = data.update_each;
727 
728  // then Jacobians
729  if (update_flags & update_contravariant_transformation)
730  {
731  const unsigned int n_q_points = data.contravariant.size();
732 
733  Assert(data.n_shape_functions > 0, ExcInternalError());
734 
735  for (unsigned int point = 0; point < n_q_points; ++point)
736  {
737  const Tensor<1, dim> *data_derv =
738  &data.derivative(point + data_set, 0);
739 
740  Tensor<1, dim> result[spacedim];
741 
742  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
743  {
744  const unsigned int comp_k =
745  fe.system_to_component_index(k).first;
746  if (fe_mask[comp_k])
747  result[fe_to_real[comp_k]] +=
748  data.local_dof_values[k] * data_derv[k];
749  }
750 
751  // write result into contravariant data
752  for (unsigned int i = 0; i < spacedim; ++i)
753  {
754  data.contravariant[point][i] = result[i];
755  }
756  }
757  }
758 
759  if (update_flags & update_covariant_transformation)
760  {
761  AssertDimension(data.covariant.size(), data.contravariant.size());
762  for (unsigned int point = 0; point < data.contravariant.size();
763  ++point)
764  data.covariant[point] =
765  (data.contravariant[point]).covariant_form();
766  }
767 
768  if (update_flags & update_volume_elements)
769  {
770  AssertDimension(data.contravariant.size(),
771  data.volume_elements.size());
772  for (unsigned int point = 0; point < data.contravariant.size();
773  ++point)
774  data.volume_elements[point] =
775  data.contravariant[point].determinant();
776  }
777  }
778 
785  template <int dim,
786  int spacedim,
787  typename VectorType,
788  typename DoFHandlerType>
789  void
790  maybe_update_jacobian_grads(
791  const typename ::QProjector<dim>::DataSetDescriptor data_set,
792  const typename ::
794  InternalData & data,
795  const FiniteElement<dim, spacedim> & fe,
796  const ComponentMask & fe_mask,
797  const std::vector<unsigned int> & fe_to_real,
798  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
799  {
800  const UpdateFlags update_flags = data.update_each;
801  if (update_flags & update_jacobian_grads)
802  {
803  const unsigned int n_q_points = jacobian_grads.size();
804 
805  for (unsigned int point = 0; point < n_q_points; ++point)
806  {
807  const Tensor<2, dim> *second =
808  &data.second_derivative(point + data_set, 0);
809 
811 
812  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
813  {
814  const unsigned int comp_k =
815  fe.system_to_component_index(k).first;
816  if (fe_mask[comp_k])
817  for (unsigned int j = 0; j < dim; ++j)
818  for (unsigned int l = 0; l < dim; ++l)
819  result[fe_to_real[comp_k]][j][l] +=
820  (second[k][j][l] * data.local_dof_values[k]);
821  }
822 
823  // never touch any data for j=dim in case dim<spacedim, so
824  // it will always be zero as it was initialized
825  for (unsigned int i = 0; i < spacedim; ++i)
826  for (unsigned int j = 0; j < dim; ++j)
827  for (unsigned int l = 0; l < dim; ++l)
828  jacobian_grads[point][i][j][l] = result[i][j][l];
829  }
830  }
831  }
832 
839  template <int dim,
840  int spacedim,
841  typename VectorType,
842  typename DoFHandlerType>
843  void
844  maybe_update_jacobian_pushed_forward_grads(
845  const typename ::QProjector<dim>::DataSetDescriptor data_set,
846  const typename ::
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,
917  int spacedim,
918  typename VectorType,
919  typename DoFHandlerType>
920  void
921  maybe_update_jacobian_2nd_derivatives(
922  const typename ::QProjector<dim>::DataSetDescriptor data_set,
923  const typename ::
925  InternalData & data,
926  const FiniteElement<dim, spacedim> & fe,
927  const ComponentMask & fe_mask,
928  const std::vector<unsigned int> & fe_to_real,
929  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
930  {
931  const UpdateFlags update_flags = data.update_each;
932  if (update_flags & update_jacobian_2nd_derivatives)
933  {
934  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
935 
936  for (unsigned int point = 0; point < n_q_points; ++point)
937  {
938  const Tensor<3, dim> *third =
939  &data.third_derivative(point + data_set, 0);
940 
942 
943  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
944  {
945  const unsigned int comp_k =
946  fe.system_to_component_index(k).first;
947  if (fe_mask[comp_k])
948  for (unsigned int j = 0; j < dim; ++j)
949  for (unsigned int l = 0; l < dim; ++l)
950  for (unsigned int m = 0; m < dim; ++m)
951  result[fe_to_real[comp_k]][j][l][m] +=
952  (third[k][j][l][m] * data.local_dof_values[k]);
953  }
954 
955  // never touch any data for j=dim in case dim<spacedim, so
956  // it will always be zero as it was initialized
957  for (unsigned int i = 0; i < spacedim; ++i)
958  for (unsigned int j = 0; j < dim; ++j)
959  for (unsigned int l = 0; l < dim; ++l)
960  for (unsigned int m = 0; m < dim; ++m)
961  jacobian_2nd_derivatives[point][i][j][l][m] =
962  result[i][j][l][m];
963  }
964  }
965  }
966 
974  template <int dim,
975  int spacedim,
976  typename VectorType,
977  typename DoFHandlerType>
978  void
979  maybe_update_jacobian_pushed_forward_2nd_derivatives(
980  const typename ::QProjector<dim>::DataSetDescriptor data_set,
981  const typename ::
983  InternalData & data,
985  const ComponentMask & fe_mask,
986  const std::vector<unsigned int> & fe_to_real,
987  std::vector<Tensor<4, spacedim>>
988  &jacobian_pushed_forward_2nd_derivatives)
989  {
990  const UpdateFlags update_flags = data.update_each;
992  {
993  const unsigned int n_q_points =
994  jacobian_pushed_forward_2nd_derivatives.size();
995 
996  double tmp[spacedim][spacedim][spacedim][spacedim];
997  for (unsigned int point = 0; point < n_q_points; ++point)
998  {
999  const Tensor<3, dim> *third =
1000  &data.third_derivative(point + data_set, 0);
1001 
1003 
1004  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1005  {
1006  const unsigned int comp_k =
1007  fe.system_to_component_index(k).first;
1008  if (fe_mask[comp_k])
1009  for (unsigned int j = 0; j < dim; ++j)
1010  for (unsigned int l = 0; l < dim; ++l)
1011  for (unsigned int m = 0; m < dim; ++m)
1012  result[fe_to_real[comp_k]][j][l][m] +=
1013  (third[k][j][l][m] * data.local_dof_values[k]);
1014  }
1015 
1016  // push forward the j-coordinate
1017  for (unsigned int i = 0; i < spacedim; ++i)
1018  for (unsigned int j = 0; j < spacedim; ++j)
1019  for (unsigned int l = 0; l < dim; ++l)
1020  for (unsigned int m = 0; m < dim; ++m)
1021  {
1022  jacobian_pushed_forward_2nd_derivatives
1023  [point][i][j][l][m] =
1024  result[i][0][l][m] * data.covariant[point][j][0];
1025  for (unsigned int jr = 1; jr < dim; ++jr)
1026  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1027  [l][m] +=
1028  result[i][jr][l][m] *
1029  data.covariant[point][j][jr];
1030  }
1031 
1032  // push forward the l-coordinate
1033  for (unsigned int i = 0; i < spacedim; ++i)
1034  for (unsigned int j = 0; j < spacedim; ++j)
1035  for (unsigned int l = 0; l < spacedim; ++l)
1036  for (unsigned int m = 0; m < dim; ++m)
1037  {
1038  tmp[i][j][l][m] =
1039  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1040  [0][m] *
1041  data.covariant[point][l][0];
1042  for (unsigned int lr = 1; lr < dim; ++lr)
1043  tmp[i][j][l][m] +=
1044  jacobian_pushed_forward_2nd_derivatives[point][i]
1045  [j][lr]
1046  [m] *
1047  data.covariant[point][l][lr];
1048  }
1049 
1050  // push forward the m-coordinate
1051  for (unsigned int i = 0; i < spacedim; ++i)
1052  for (unsigned int j = 0; j < spacedim; ++j)
1053  for (unsigned int l = 0; l < spacedim; ++l)
1054  for (unsigned int m = 0; m < spacedim; ++m)
1055  {
1056  jacobian_pushed_forward_2nd_derivatives
1057  [point][i][j][l][m] =
1058  tmp[i][j][l][0] * data.covariant[point][m][0];
1059  for (unsigned int mr = 1; mr < dim; ++mr)
1060  jacobian_pushed_forward_2nd_derivatives[point][i][j]
1061  [l][m] +=
1062  tmp[i][j][l][mr] * data.covariant[point][m][mr];
1063  }
1064  }
1065  }
1066  }
1067 
1074  template <int dim,
1075  int spacedim,
1076  typename VectorType,
1077  typename DoFHandlerType>
1078  void
1079  maybe_update_jacobian_3rd_derivatives(
1080  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1081  const typename ::
1083  InternalData & data,
1084  const FiniteElement<dim, spacedim> & fe,
1085  const ComponentMask & fe_mask,
1086  const std::vector<unsigned int> & fe_to_real,
1087  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1088  {
1089  const UpdateFlags update_flags = data.update_each;
1090  if (update_flags & update_jacobian_3rd_derivatives)
1091  {
1092  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1093 
1094  for (unsigned int point = 0; point < n_q_points; ++point)
1095  {
1096  const Tensor<4, dim> *fourth =
1097  &data.fourth_derivative(point + data_set, 0);
1098 
1100 
1101  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1102  {
1103  const unsigned int comp_k =
1104  fe.system_to_component_index(k).first;
1105  if (fe_mask[comp_k])
1106  for (unsigned int j = 0; j < dim; ++j)
1107  for (unsigned int l = 0; l < dim; ++l)
1108  for (unsigned int m = 0; m < dim; ++m)
1109  for (unsigned int n = 0; n < dim; ++n)
1110  result[fe_to_real[comp_k]][j][l][m][n] +=
1111  (fourth[k][j][l][m][n] *
1112  data.local_dof_values[k]);
1113  }
1114 
1115  // never touch any data for j,l,m,n=dim in case
1116  // dim<spacedim, so it will always be zero as it was
1117  // initialized
1118  for (unsigned int i = 0; i < spacedim; ++i)
1119  for (unsigned int j = 0; j < dim; ++j)
1120  for (unsigned int l = 0; l < dim; ++l)
1121  for (unsigned int m = 0; m < dim; ++m)
1122  for (unsigned int n = 0; n < dim; ++n)
1123  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1124  result[i][j][l][m][n];
1125  }
1126  }
1127  }
1128 
1136  template <int dim,
1137  int spacedim,
1138  typename VectorType,
1139  typename DoFHandlerType>
1140  void
1141  maybe_update_jacobian_pushed_forward_3rd_derivatives(
1142  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1143  const typename ::
1145  InternalData & data,
1146  const FiniteElement<dim, spacedim> &fe,
1147  const ComponentMask & fe_mask,
1148  const std::vector<unsigned int> & fe_to_real,
1149  std::vector<Tensor<5, spacedim>>
1150  &jacobian_pushed_forward_3rd_derivatives)
1151  {
1152  const UpdateFlags update_flags = data.update_each;
1154  {
1155  const unsigned int n_q_points =
1156  jacobian_pushed_forward_3rd_derivatives.size();
1157 
1158  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1159  for (unsigned int point = 0; point < n_q_points; ++point)
1160  {
1161  const Tensor<4, dim> *fourth =
1162  &data.fourth_derivative(point + data_set, 0);
1163 
1165 
1166  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1167  {
1168  const unsigned int comp_k =
1169  fe.system_to_component_index(k).first;
1170  if (fe_mask[comp_k])
1171  for (unsigned int j = 0; j < dim; ++j)
1172  for (unsigned int l = 0; l < dim; ++l)
1173  for (unsigned int m = 0; m < dim; ++m)
1174  for (unsigned int n = 0; n < dim; ++n)
1175  result[fe_to_real[comp_k]][j][l][m][n] +=
1176  (fourth[k][j][l][m][n] *
1177  data.local_dof_values[k]);
1178  }
1179 
1180  // push-forward the j-coordinate
1181  for (unsigned int i = 0; i < spacedim; ++i)
1182  for (unsigned int j = 0; j < spacedim; ++j)
1183  for (unsigned int l = 0; l < dim; ++l)
1184  for (unsigned int m = 0; m < dim; ++m)
1185  for (unsigned int n = 0; n < dim; ++n)
1186  {
1187  tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1188  data.covariant[point][j][0];
1189  for (unsigned int jr = 1; jr < dim; ++jr)
1190  tmp[i][j][l][m][n] +=
1191  result[i][jr][l][m][n] *
1192  data.covariant[point][j][jr];
1193  }
1194 
1195  // push-forward the l-coordinate
1196  for (unsigned int i = 0; i < spacedim; ++i)
1197  for (unsigned int j = 0; j < spacedim; ++j)
1198  for (unsigned int l = 0; l < spacedim; ++l)
1199  for (unsigned int m = 0; m < dim; ++m)
1200  for (unsigned int n = 0; n < dim; ++n)
1201  {
1202  jacobian_pushed_forward_3rd_derivatives
1203  [point][i][j][l][m][n] =
1204  tmp[i][j][0][m][n] *
1205  data.covariant[point][l][0];
1206  for (unsigned int lr = 1; lr < dim; ++lr)
1207  jacobian_pushed_forward_3rd_derivatives[point][i]
1208  [j][l][m]
1209  [n] +=
1210  tmp[i][j][lr][m][n] *
1211  data.covariant[point][l][lr];
1212  }
1213 
1214  // push-forward the m-coordinate
1215  for (unsigned int i = 0; i < spacedim; ++i)
1216  for (unsigned int j = 0; j < spacedim; ++j)
1217  for (unsigned int l = 0; l < spacedim; ++l)
1218  for (unsigned int m = 0; m < spacedim; ++m)
1219  for (unsigned int n = 0; n < dim; ++n)
1220  {
1221  tmp[i][j][l][m][n] =
1222  jacobian_pushed_forward_3rd_derivatives[point][i]
1223  [j][l][0]
1224  [n] *
1225  data.covariant[point][m][0];
1226  for (unsigned int mr = 1; mr < dim; ++mr)
1227  tmp[i][j][l][m][n] +=
1228  jacobian_pushed_forward_3rd_derivatives[point]
1229  [i][j][l]
1230  [mr][n] *
1231  data.covariant[point][m][mr];
1232  }
1233 
1234  // push-forward the n-coordinate
1235  for (unsigned int i = 0; i < spacedim; ++i)
1236  for (unsigned int j = 0; j < spacedim; ++j)
1237  for (unsigned int l = 0; l < spacedim; ++l)
1238  for (unsigned int m = 0; m < spacedim; ++m)
1239  for (unsigned int n = 0; n < spacedim; ++n)
1240  {
1241  jacobian_pushed_forward_3rd_derivatives
1242  [point][i][j][l][m][n] =
1243  tmp[i][j][l][m][0] *
1244  data.covariant[point][n][0];
1245  for (unsigned int nr = 1; nr < dim; ++nr)
1246  jacobian_pushed_forward_3rd_derivatives[point][i]
1247  [j][l][m]
1248  [n] +=
1249  tmp[i][j][l][m][nr] *
1250  data.covariant[point][n][nr];
1251  }
1252  }
1253  }
1254  }
1255 
1256 
1266  template <int dim,
1267  int spacedim,
1268  typename VectorType,
1269  typename DoFHandlerType>
1270  void
1271  maybe_compute_face_data(
1272  const ::Mapping<dim, spacedim> &mapping,
1273  const typename ::Triangulation<dim, spacedim>::cell_iterator
1274  & cell,
1275  const unsigned int face_no,
1276  const unsigned int subface_no,
1277  const std::vector<double> &weights,
1278  const typename ::
1280  InternalData &data,
1282  &output_data)
1283  {
1284  const UpdateFlags update_flags = data.update_each;
1285 
1286  if (update_flags & update_boundary_forms)
1287  {
1288  const unsigned int n_q_points = output_data.boundary_forms.size();
1289  if (update_flags & update_normal_vectors)
1290  AssertDimension(output_data.normal_vectors.size(), n_q_points);
1291  if (update_flags & update_JxW_values)
1292  AssertDimension(output_data.JxW_values.size(), n_q_points);
1293 
1294  // map the unit tangentials to the real cell. checking for d!=dim-1
1295  // eliminates compiler warnings regarding unsigned int expressions <
1296  // 0.
1297  for (unsigned int d = 0; d != dim - 1; ++d)
1298  {
1300  data.unit_tangentials.size(),
1301  ExcInternalError());
1302  Assert(
1303  data.aux[d].size() <=
1304  data
1305  .unit_tangentials[face_no +
1307  .size(),
1308  ExcInternalError());
1309 
1310  mapping.transform(
1312  data
1313  .unit_tangentials[face_no +
1316  data,
1317  make_array_view(data.aux[d]));
1318  }
1319 
1320  // if dim==spacedim, we can use the unit tangentials to compute the
1321  // boundary form by simply taking the cross product
1322  if (dim == spacedim)
1323  {
1324  for (unsigned int i = 0; i < n_q_points; ++i)
1325  switch (dim)
1326  {
1327  case 1:
1328  // in 1d, we don't have access to any of the data.aux
1329  // fields (because it has only dim-1 components), but we
1330  // can still compute the boundary form by simply looking
1331  // at the number of the face
1332  output_data.boundary_forms[i][0] =
1333  (face_no == 0 ? -1 : +1);
1334  break;
1335  case 2:
1336  output_data.boundary_forms[i] =
1337  cross_product_2d(data.aux[0][i]);
1338  break;
1339  case 3:
1340  output_data.boundary_forms[i] =
1341  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1342  break;
1343  default:
1344  Assert(false, ExcNotImplemented());
1345  }
1346  }
1347  else //(dim < spacedim)
1348  {
1349  // in the codim-one case, the boundary form results from the
1350  // cross product of all the face tangential vectors and the cell
1351  // normal vector
1352  //
1353  // to compute the cell normal, use the same method used in
1354  // fill_fe_values for cells above
1355  AssertDimension(data.contravariant.size(), n_q_points);
1356 
1357  for (unsigned int point = 0; point < n_q_points; ++point)
1358  {
1359  if (dim == 1)
1360  {
1361  // J is a tangent vector
1362  output_data.boundary_forms[point] =
1363  data.contravariant[point].transpose()[0];
1364  output_data.boundary_forms[point] /=
1365  (face_no == 0 ? -1. : +1.) *
1366  output_data.boundary_forms[point].norm();
1367  }
1368 
1369  if (dim == 2)
1370  {
1372  data.contravariant[point].transpose();
1373 
1374  Tensor<1, spacedim> cell_normal =
1375  cross_product_3d(DX_t[0], DX_t[1]);
1376  cell_normal /= cell_normal.norm();
1377 
1378  // then compute the face normal from the face tangent
1379  // and the cell normal:
1380  output_data.boundary_forms[point] =
1381  cross_product_3d(data.aux[0][point], cell_normal);
1382  }
1383  }
1384  }
1385 
1386  if (update_flags & (update_normal_vectors | update_JxW_values))
1387  for (unsigned int i = 0; i < output_data.boundary_forms.size();
1388  ++i)
1389  {
1390  if (update_flags & update_JxW_values)
1391  {
1392  output_data.JxW_values[i] =
1393  output_data.boundary_forms[i].norm() * weights[i];
1394 
1395  if (subface_no != numbers::invalid_unsigned_int)
1396  {
1397  const double area_ratio =
1399  cell->subface_case(face_no), subface_no);
1400  output_data.JxW_values[i] *= area_ratio;
1401  }
1402  }
1403 
1404  if (update_flags & update_normal_vectors)
1405  output_data.normal_vectors[i] =
1406  Point<spacedim>(output_data.boundary_forms[i] /
1407  output_data.boundary_forms[i].norm());
1408  }
1409 
1410  if (update_flags & update_jacobians)
1411  for (unsigned int point = 0; point < n_q_points; ++point)
1412  output_data.jacobians[point] = data.contravariant[point];
1413 
1414  if (update_flags & update_inverse_jacobians)
1415  for (unsigned int point = 0; point < n_q_points; ++point)
1416  output_data.inverse_jacobians[point] =
1417  data.covariant[point].transpose();
1418  }
1419  }
1420 
1427  template <int dim,
1428  int spacedim,
1429  typename VectorType,
1430  typename DoFHandlerType>
1431  void
1432  do_fill_fe_face_values(
1433  const ::Mapping<dim, spacedim> &mapping,
1434  const typename ::Triangulation<dim, spacedim>::cell_iterator
1435  & cell,
1436  const unsigned int face_no,
1437  const unsigned int subface_no,
1438  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1439  const Quadrature<dim - 1> & quadrature,
1440  const typename ::
1442  InternalData & data,
1443  const FiniteElement<dim, spacedim> &fe,
1444  const ComponentMask & fe_mask,
1445  const std::vector<unsigned int> & fe_to_real,
1447  &output_data)
1448  {
1449  maybe_compute_q_points<dim, spacedim, VectorType, DoFHandlerType>(
1450  data_set,
1451  data,
1452  fe,
1453  fe_mask,
1454  fe_to_real,
1455  output_data.quadrature_points);
1456 
1457  maybe_update_Jacobians<dim, spacedim, VectorType, DoFHandlerType>(
1458  data_set, data, fe, fe_mask, fe_to_real);
1459 
1460  maybe_update_jacobian_grads<dim, spacedim, VectorType, DoFHandlerType>(
1461  data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1462 
1463  maybe_update_jacobian_pushed_forward_grads<dim,
1464  spacedim,
1465  VectorType,
1466  DoFHandlerType>(
1467  data_set,
1468  data,
1469  fe,
1470  fe_mask,
1471  fe_to_real,
1472  output_data.jacobian_pushed_forward_grads);
1473 
1474  maybe_update_jacobian_2nd_derivatives<dim,
1475  spacedim,
1476  VectorType,
1477  DoFHandlerType>(
1478  data_set,
1479  data,
1480  fe,
1481  fe_mask,
1482  fe_to_real,
1483  output_data.jacobian_2nd_derivatives);
1484 
1485  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1486  spacedim,
1487  VectorType,
1488  DoFHandlerType>(
1489  data_set,
1490  data,
1491  fe,
1492  fe_mask,
1493  fe_to_real,
1495 
1496  maybe_update_jacobian_3rd_derivatives<dim,
1497  spacedim,
1498  VectorType,
1499  DoFHandlerType>(
1500  data_set,
1501  data,
1502  fe,
1503  fe_mask,
1504  fe_to_real,
1505  output_data.jacobian_3rd_derivatives);
1506 
1507  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1508  spacedim,
1509  VectorType,
1510  DoFHandlerType>(
1511  data_set,
1512  data,
1513  fe,
1514  fe_mask,
1515  fe_to_real,
1517 
1518  maybe_compute_face_data<dim, spacedim, VectorType, DoFHandlerType>(
1519  mapping,
1520  cell,
1521  face_no,
1522  subface_no,
1523  quadrature.get_weights(),
1524  data,
1525  output_data);
1526  }
1527  } // namespace
1528  } // namespace MappingFEFieldImplementation
1529 } // namespace internal
1530 
1531 
1532 // Note that the CellSimilarity flag is modifiable, since MappingFEField can
1533 // need to recalculate data even when cells are similar.
1534 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1537  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1539  const Quadrature<dim> & quadrature,
1540  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1542  &output_data) const
1543 {
1544  // convert data object to internal data for this class. fails with an
1545  // exception if that is not possible
1546  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1547  ExcInternalError());
1548  const InternalData &data = static_cast<const InternalData &>(internal_data);
1549 
1550  const unsigned int n_q_points = quadrature.size();
1551 
1552  update_internal_dofs(cell, data);
1553 
1554  internal::MappingFEFieldImplementation::
1555  maybe_compute_q_points<dim, spacedim, VectorType, DoFHandlerType>(
1557  data,
1558  euler_dof_handler->get_fe(),
1559  fe_mask,
1560  fe_to_real,
1561  output_data.quadrature_points);
1562 
1563  internal::MappingFEFieldImplementation::
1564  maybe_update_Jacobians<dim, spacedim, VectorType, DoFHandlerType>(
1566  data,
1567  euler_dof_handler->get_fe(),
1568  fe_mask,
1569  fe_to_real);
1570 
1571  const UpdateFlags update_flags = data.update_each;
1572  const std::vector<double> &weights = quadrature.get_weights();
1573 
1574  // Multiply quadrature weights by absolute value of Jacobian determinants or
1575  // the area element g=sqrt(DX^t DX) in case of codim > 0
1576 
1577  if (update_flags & (update_normal_vectors | update_JxW_values))
1578  {
1579  AssertDimension(output_data.JxW_values.size(), n_q_points);
1580 
1581  Assert(!(update_flags & update_normal_vectors) ||
1582  (output_data.normal_vectors.size() == n_q_points),
1583  ExcDimensionMismatch(output_data.normal_vectors.size(),
1584  n_q_points));
1585 
1586 
1587  for (unsigned int point = 0; point < n_q_points; ++point)
1588  {
1589  if (dim == spacedim)
1590  {
1591  const double det = data.contravariant[point].determinant();
1592 
1593  // check for distorted cells.
1594 
1595  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1596  // 1e12 in 2D. might want to find a finer
1597  // (dimension-independent) criterion
1598  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1599  cell->diameter() / std::sqrt(double(dim))),
1601  cell->center(), det, point)));
1602  output_data.JxW_values[point] = weights[point] * det;
1603  }
1604  // if dim==spacedim, then there is no cell normal to
1605  // compute. since this is for FEValues (and not FEFaceValues),
1606  // there are also no face normals to compute
1607  else // codim>0 case
1608  {
1609  Tensor<1, spacedim> DX_t[dim];
1610  for (unsigned int i = 0; i < spacedim; ++i)
1611  for (unsigned int j = 0; j < dim; ++j)
1612  DX_t[j][i] = data.contravariant[point][i][j];
1613 
1614  Tensor<2, dim> G; // First fundamental form
1615  for (unsigned int i = 0; i < dim; ++i)
1616  for (unsigned int j = 0; j < dim; ++j)
1617  G[i][j] = DX_t[i] * DX_t[j];
1618 
1619  output_data.JxW_values[point] =
1620  std::sqrt(determinant(G)) * weights[point];
1621 
1622  if (update_flags & update_normal_vectors)
1623  {
1624  Assert(spacedim - dim == 1,
1625  ExcMessage("There is no cell normal in codim 2."));
1626 
1627  if (dim == 1)
1628  output_data.normal_vectors[point] =
1629  cross_product_2d(-DX_t[0]);
1630  else // dim == 2
1631  output_data.normal_vectors[point] =
1632  cross_product_3d(DX_t[0], DX_t[1]);
1633 
1634  output_data.normal_vectors[point] /=
1635  output_data.normal_vectors[point].norm();
1636 
1637  if (cell->direction_flag() == false)
1638  output_data.normal_vectors[point] *= -1.;
1639  }
1640  } // codim>0 case
1641  }
1642  }
1643 
1644  // copy values from InternalData to vector given by reference
1645  if (update_flags & update_jacobians)
1646  {
1647  AssertDimension(output_data.jacobians.size(), n_q_points);
1648  for (unsigned int point = 0; point < n_q_points; ++point)
1649  output_data.jacobians[point] = data.contravariant[point];
1650  }
1651 
1652  // copy values from InternalData to vector given by reference
1653  if (update_flags & update_inverse_jacobians)
1654  {
1655  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1656  for (unsigned int point = 0; point < n_q_points; ++point)
1657  output_data.inverse_jacobians[point] =
1658  data.covariant[point].transpose();
1659  }
1660 
1661  // calculate derivatives of the Jacobians
1662  internal::MappingFEFieldImplementation::
1663  maybe_update_jacobian_grads<dim, spacedim, VectorType, DoFHandlerType>(
1665  data,
1666  euler_dof_handler->get_fe(),
1667  fe_mask,
1668  fe_to_real,
1669  output_data.jacobian_grads);
1670 
1671  // calculate derivatives of the Jacobians pushed forward to real cell
1672  // coordinates
1673  internal::MappingFEFieldImplementation::
1674  maybe_update_jacobian_pushed_forward_grads<dim,
1675  spacedim,
1676  VectorType,
1677  DoFHandlerType>(
1679  data,
1680  euler_dof_handler->get_fe(),
1681  fe_mask,
1682  fe_to_real,
1683  output_data.jacobian_pushed_forward_grads);
1684 
1685  // calculate hessians of the Jacobians
1686  internal::MappingFEFieldImplementation::maybe_update_jacobian_2nd_derivatives<
1687  dim,
1688  spacedim,
1689  VectorType,
1690  DoFHandlerType>(QProjector<dim>::DataSetDescriptor::cell(),
1691  data,
1692  euler_dof_handler->get_fe(),
1693  fe_mask,
1694  fe_to_real,
1695  output_data.jacobian_2nd_derivatives);
1696 
1697  // calculate hessians of the Jacobians pushed forward to real cell coordinates
1698  internal::MappingFEFieldImplementation::
1699  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1700  spacedim,
1701  VectorType,
1702  DoFHandlerType>(
1704  data,
1705  euler_dof_handler->get_fe(),
1706  fe_mask,
1707  fe_to_real,
1709 
1710  // calculate gradients of the hessians of the Jacobians
1711  internal::MappingFEFieldImplementation::maybe_update_jacobian_3rd_derivatives<
1712  dim,
1713  spacedim,
1714  VectorType,
1715  DoFHandlerType>(QProjector<dim>::DataSetDescriptor::cell(),
1716  data,
1717  euler_dof_handler->get_fe(),
1718  fe_mask,
1719  fe_to_real,
1720  output_data.jacobian_3rd_derivatives);
1721 
1722  // calculate gradients of the hessians of the Jacobians pushed forward to real
1723  // cell coordinates
1724  internal::MappingFEFieldImplementation::
1725  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1726  spacedim,
1727  VectorType,
1728  DoFHandlerType>(
1730  data,
1731  euler_dof_handler->get_fe(),
1732  fe_mask,
1733  fe_to_real,
1735 
1737 }
1738 
1739 
1740 
1741 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1742 void
1744  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1745  const unsigned int face_no,
1746  const Quadrature<dim - 1> & quadrature,
1747  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1749  &output_data) const
1750 {
1751  // convert data object to internal data for this class. fails with an
1752  // exception if that is not possible
1753  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1754  ExcInternalError());
1755  const InternalData &data = static_cast<const InternalData &>(internal_data);
1756 
1757  update_internal_dofs(cell, data);
1758 
1759  internal::MappingFEFieldImplementation::
1760  do_fill_fe_face_values<dim, spacedim, VectorType, DoFHandlerType>(
1761  *this,
1762  cell,
1763  face_no,
1766  cell->face_orientation(face_no),
1767  cell->face_flip(face_no),
1768  cell->face_rotation(face_no),
1769  quadrature.size()),
1770  quadrature,
1771  data,
1772  euler_dof_handler->get_fe(),
1773  fe_mask,
1774  fe_to_real,
1775  output_data);
1776 }
1777 
1778 
1779 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1780 void
1783  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1784  const unsigned int face_no,
1785  const unsigned int subface_no,
1786  const Quadrature<dim - 1> & quadrature,
1787  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1789  &output_data) const
1790 {
1791  // convert data object to internal data for this class. fails with an
1792  // exception if that is not possible
1793  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1794  ExcInternalError());
1795  const InternalData &data = static_cast<const InternalData &>(internal_data);
1796 
1797  update_internal_dofs(cell, data);
1798 
1799  internal::MappingFEFieldImplementation::
1800  do_fill_fe_face_values<dim, spacedim, VectorType, DoFHandlerType>(
1801  *this,
1802  cell,
1803  face_no,
1806  subface_no,
1807  cell->face_orientation(
1808  face_no),
1809  cell->face_flip(face_no),
1810  cell->face_rotation(face_no),
1811  quadrature.size(),
1812  cell->subface_case(face_no)),
1813  quadrature,
1814  data,
1815  euler_dof_handler->get_fe(),
1816  fe_mask,
1817  fe_to_real,
1818  output_data);
1819 }
1820 
1821 
1822 namespace internal
1823 {
1824  namespace MappingFEFieldImplementation
1825  {
1826  namespace
1827  {
1828  template <int dim,
1829  int spacedim,
1830  int rank,
1831  typename VectorType,
1832  typename DoFHandlerType>
1833  void
1834  transform_fields(
1835  const ArrayView<const Tensor<rank, dim>> & input,
1836  const MappingKind mapping_kind,
1837  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1838  const ArrayView<Tensor<rank, spacedim>> & output)
1839  {
1840  AssertDimension(input.size(), output.size());
1841  Assert((dynamic_cast<
1842  const typename ::
1843  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1844  InternalData *>(&mapping_data) != nullptr),
1845  ExcInternalError());
1846  const typename ::MappingFEField<dim,
1847  spacedim,
1848  VectorType,
1849  DoFHandlerType>::InternalData
1850  &data = static_cast<
1851  const typename ::
1852  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1853  InternalData &>(mapping_data);
1854 
1855  switch (mapping_kind)
1856  {
1857  case mapping_contravariant:
1858  {
1859  Assert(
1860  data.update_each & update_contravariant_transformation,
1862  "update_contravariant_transformation"));
1863 
1864  for (unsigned int i = 0; i < output.size(); ++i)
1865  output[i] =
1866  apply_transformation(data.contravariant[i], input[i]);
1867 
1868  return;
1869  }
1870 
1871  case mapping_piola:
1872  {
1873  Assert(
1874  data.update_each & update_contravariant_transformation,
1876  "update_contravariant_transformation"));
1877  Assert(
1878  data.update_each & update_volume_elements,
1880  "update_volume_elements"));
1881  Assert(rank == 1, ExcMessage("Only for rank 1"));
1882  for (unsigned int i = 0; i < output.size(); ++i)
1883  {
1884  output[i] =
1885  apply_transformation(data.contravariant[i], input[i]);
1886  output[i] /= data.volume_elements[i];
1887  }
1888  return;
1889  }
1890 
1891 
1892  // We still allow this operation as in the
1893  // reference cell Derivatives are Tensor
1894  // rather than DerivativeForm
1895  case mapping_covariant:
1896  {
1897  Assert(
1898  data.update_each & update_contravariant_transformation,
1900  "update_contravariant_transformation"));
1901 
1902  for (unsigned int i = 0; i < output.size(); ++i)
1903  output[i] = apply_transformation(data.covariant[i], input[i]);
1904 
1905  return;
1906  }
1907 
1908  default:
1909  Assert(false, ExcNotImplemented());
1910  }
1911  }
1912 
1913 
1914  template <int dim,
1915  int spacedim,
1916  int rank,
1917  typename VectorType,
1918  typename DoFHandlerType>
1919  void
1920  transform_differential_forms(
1921  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
1922  const MappingKind mapping_kind,
1923  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1924  const ArrayView<Tensor<rank + 1, spacedim>> & output)
1925  {
1926  AssertDimension(input.size(), output.size());
1927  Assert((dynamic_cast<
1928  const typename ::
1929  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1930  InternalData *>(&mapping_data) != nullptr),
1931  ExcInternalError());
1932  const typename ::MappingFEField<dim,
1933  spacedim,
1934  VectorType,
1935  DoFHandlerType>::InternalData
1936  &data = static_cast<
1937  const typename ::
1938  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1939  InternalData &>(mapping_data);
1940 
1941  switch (mapping_kind)
1942  {
1943  case mapping_covariant:
1944  {
1945  Assert(
1946  data.update_each & update_contravariant_transformation,
1948  "update_contravariant_transformation"));
1949 
1950  for (unsigned int i = 0; i < output.size(); ++i)
1951  output[i] = apply_transformation(data.covariant[i], input[i]);
1952 
1953  return;
1954  }
1955  default:
1956  Assert(false, ExcNotImplemented());
1957  }
1958  }
1959  } // namespace
1960  } // namespace MappingFEFieldImplementation
1961 } // namespace internal
1962 
1963 
1964 
1965 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1966 void
1968  const ArrayView<const Tensor<1, dim>> & input,
1969  const MappingKind mapping_kind,
1970  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1971  const ArrayView<Tensor<1, spacedim>> & output) const
1972 {
1973  AssertDimension(input.size(), output.size());
1974 
1975  internal::MappingFEFieldImplementation::
1976  transform_fields<dim, spacedim, 1, VectorType, DoFHandlerType>(input,
1977  mapping_kind,
1978  mapping_data,
1979  output);
1980 }
1981 
1982 
1983 
1984 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1985 void
1987  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1988  const MappingKind mapping_kind,
1989  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1990  const ArrayView<Tensor<2, spacedim>> & output) const
1991 {
1992  AssertDimension(input.size(), output.size());
1993 
1994  internal::MappingFEFieldImplementation::
1995  transform_differential_forms<dim, spacedim, 1, VectorType, DoFHandlerType>(
1996  input, mapping_kind, mapping_data, output);
1997 }
1998 
1999 
2000 
2001 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2002 void
2004  const ArrayView<const Tensor<2, dim>> &input,
2005  const MappingKind,
2006  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2007  const ArrayView<Tensor<2, spacedim>> & output) const
2008 {
2009  (void)input;
2010  (void)output;
2011  (void)mapping_data;
2012  AssertDimension(input.size(), output.size());
2013 
2014  AssertThrow(false, ExcNotImplemented());
2015 }
2016 
2017 
2018 
2019 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2020 void
2022  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2023  const MappingKind mapping_kind,
2024  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2025  const ArrayView<Tensor<3, spacedim>> & output) const
2026 {
2027  AssertDimension(input.size(), output.size());
2028  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2029  ExcInternalError());
2030  const InternalData &data = static_cast<const InternalData &>(mapping_data);
2031 
2032  switch (mapping_kind)
2033  {
2035  {
2038  "update_covariant_transformation"));
2039 
2040  for (unsigned int q = 0; q < output.size(); ++q)
2041  for (unsigned int i = 0; i < spacedim; ++i)
2042  for (unsigned int j = 0; j < spacedim; ++j)
2043  for (unsigned int k = 0; k < spacedim; ++k)
2044  {
2045  output[q][i][j][k] = data.covariant[q][j][0] *
2046  data.covariant[q][k][0] *
2047  input[q][i][0][0];
2048  for (unsigned int J = 0; J < dim; ++J)
2049  {
2050  const unsigned int K0 = (0 == J) ? 1 : 0;
2051  for (unsigned int K = K0; K < dim; ++K)
2052  output[q][i][j][k] += data.covariant[q][j][J] *
2053  data.covariant[q][k][K] *
2054  input[q][i][J][K];
2055  }
2056  }
2057  return;
2058  }
2059 
2060  default:
2061  Assert(false, ExcNotImplemented());
2062  }
2063 }
2064 
2065 
2066 
2067 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2068 void
2070  const ArrayView<const Tensor<3, dim>> &input,
2071  const MappingKind /*mapping_kind*/,
2072  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2073  const ArrayView<Tensor<3, spacedim>> & output) const
2074 {
2075  (void)input;
2076  (void)output;
2077  (void)mapping_data;
2078  AssertDimension(input.size(), output.size());
2079 
2080  AssertThrow(false, ExcNotImplemented());
2081 }
2082 
2083 
2084 
2085 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2089  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2090  const Point<dim> & p) const
2091 {
2092  // Use the get_data function to create an InternalData with data vectors of
2093  // the right size and transformation shape values already computed at point
2094  // p.
2095  const Quadrature<dim> point_quadrature(p);
2096  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2097  get_data(update_quadrature_points | update_jacobians, point_quadrature));
2098  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2099  ExcInternalError());
2100 
2101  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2102 
2103  return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
2104 }
2105 
2106 
2107 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2111 {
2112  Point<spacedim> p_real;
2113 
2114  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2115  {
2116  unsigned int comp_i =
2117  euler_dof_handler->get_fe().system_to_component_index(i).first;
2118  if (fe_mask[comp_i])
2119  p_real[fe_to_real[comp_i]] +=
2120  data.local_dof_values[i] * data.shape(0, i);
2121  }
2122 
2123  return p_real;
2124 }
2125 
2126 
2127 
2128 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2129 Point<dim>
2132  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2133  const Point<spacedim> & p) const
2134 {
2135  // first a Newton iteration based on the real mapping. It uses the center
2136  // point of the cell as a starting point
2137  Point<dim> initial_p_unit;
2138  try
2139  {
2140  initial_p_unit =
2141  StaticMappingQ1<dim, spacedim>::mapping.transform_real_to_unit_cell(
2142  cell, p);
2143  }
2144  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2145  {
2146  // mirror the conditions of the code below to determine if we need to
2147  // use an arbitrary starting point or if we just need to rethrow the
2148  // exception
2149  for (unsigned int d = 0; d < dim; ++d)
2150  initial_p_unit[d] = 0.5;
2151  }
2152 
2153  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2154 
2155  // for (unsigned int d=0; d<dim; ++d)
2156  // initial_p_unit[d] = 0.;
2157 
2158  const Quadrature<dim> point_quadrature(initial_p_unit);
2159 
2161  if (spacedim > dim)
2162  update_flags |= update_jacobian_grads;
2163  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2164  get_data(update_flags, point_quadrature));
2165  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2166  ExcInternalError());
2167 
2168  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2169 
2170  return do_transform_real_to_unit_cell(cell,
2171  p,
2172  initial_p_unit,
2173  dynamic_cast<InternalData &>(*mdata));
2174 }
2175 
2176 
2177 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2178 Point<dim>
2181  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2182  const Point<spacedim> & p,
2183  const Point<dim> & initial_p_unit,
2184  InternalData & mdata) const
2185 {
2186  const unsigned int n_shapes = mdata.shape_values.size();
2187  (void)n_shapes;
2188  Assert(n_shapes != 0, ExcInternalError());
2189  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2190 
2191 
2192  // Newton iteration to solve
2193  // f(x)=p(x)-p=0
2194  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2195  // The start value was set to be the
2196  // linear approximation to the cell
2197  // The shape values and derivatives
2198  // of the mapping at this point are
2199  // previously computed.
2200  // f(x)
2201  Point<dim> p_unit = initial_p_unit;
2202  Point<dim> f;
2203  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
2205  Tensor<1, spacedim> p_minus_F = p - p_real;
2206  const double eps = 1.e-12 * cell->diameter();
2207  const unsigned int newton_iteration_limit = 20;
2208  unsigned int newton_iteration = 0;
2209  while (p_minus_F.norm_square() > eps * eps)
2210  {
2211  // f'(x)
2212  Point<spacedim> DF[dim];
2213  Tensor<2, dim> df;
2214  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2215  {
2216  const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2217  const unsigned int comp_k =
2218  euler_dof_handler->get_fe().system_to_component_index(k).first;
2219  if (fe_mask[comp_k])
2220  for (unsigned int j = 0; j < dim; ++j)
2221  DF[j][fe_to_real[comp_k]] +=
2222  mdata.local_dof_values[k] * grad_k[j];
2223  }
2224  for (unsigned int j = 0; j < dim; ++j)
2225  {
2226  f[j] = DF[j] * p_minus_F;
2227  for (unsigned int l = 0; l < dim; ++l)
2228  df[j][l] = -DF[j] * DF[l];
2229  }
2230  // Solve [f'(x)]d=f(x)
2231  const Tensor<1, dim> delta =
2232  invert(df) * static_cast<const Tensor<1, dim> &>(f);
2233  // do a line search
2234  double step_length = 1;
2235  do
2236  {
2237  // update of p_unit. The
2238  // spacedimth component of
2239  // transformed point is simply
2240  // ignored in codimension one
2241  // case. When this component is
2242  // not zero, then we are
2243  // projecting the point to the
2244  // surface or curve identified
2245  // by the cell.
2246  Point<dim> p_unit_trial = p_unit;
2247  for (unsigned int i = 0; i < dim; ++i)
2248  p_unit_trial[i] -= step_length * delta[i];
2249  // shape values and derivatives
2250  // at new p_unit point
2251  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
2252  mdata);
2253  // f(x)
2254  Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
2255  const Tensor<1, spacedim> f_trial = p - p_real_trial;
2256  // see if we are making progress with the current step length
2257  // and if not, reduce it by a factor of two and try again
2258  if (f_trial.norm() < p_minus_F.norm())
2259  {
2260  p_real = p_real_trial;
2261  p_unit = p_unit_trial;
2262  p_minus_F = f_trial;
2263  break;
2264  }
2265  else if (step_length > 0.05)
2266  step_length /= 2;
2267  else
2268  goto failure;
2269  }
2270  while (true);
2271  ++newton_iteration;
2272  if (newton_iteration > newton_iteration_limit)
2273  goto failure;
2274  }
2275  return p_unit;
2276  // if we get to the following label, then we have either run out
2277  // of Newton iterations, or the line search has not converged.
2278  // in either case, we need to give up, so throw an exception that
2279  // can then be caught
2280 failure:
2281  AssertThrow(false,
2283  // ...the compiler wants us to return something, though we can
2284  // of course never get here...
2285  return Point<dim>();
2286 }
2287 
2288 
2289 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2290 unsigned int
2292 {
2293  return euler_dof_handler->get_fe().degree;
2294 }
2295 
2296 
2297 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2300  const
2301 {
2302  return this->fe_mask;
2303 }
2304 
2305 
2306 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2307 std::unique_ptr<Mapping<dim, spacedim>>
2309 {
2310  return std::make_unique<
2312 }
2313 
2314 
2315 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2316 void
2318  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2320  InternalData &data) const
2321 {
2322  Assert(euler_dof_handler != nullptr,
2323  ExcMessage("euler_dof_handler is empty"));
2324 
2325  typename DoFHandlerType::cell_iterator dof_cell(*cell, euler_dof_handler);
2326  Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2327  if (uses_level_dofs)
2328  {
2329  AssertIndexRange(cell->level(), euler_vector.size());
2330  AssertDimension(euler_vector[cell->level()]->size(),
2331  euler_dof_handler->n_dofs(cell->level()));
2332  }
2333  else
2334  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2335 
2336  if (uses_level_dofs)
2337  dof_cell->get_mg_dof_indices(data.local_dof_indices);
2338  else
2339  dof_cell->get_dof_indices(data.local_dof_indices);
2340 
2341  const VectorType &vector =
2342  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2343 
2344  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2345  data.local_dof_values[i] =
2347  data.local_dof_indices[i]);
2348 }
2349 
2350 // explicit instantiations
2351 #define SPLIT_INSTANTIATIONS_COUNT 2
2352 #ifndef SPLIT_INSTANTIATIONS_INDEX
2353 # define SPLIT_INSTANTIATIONS_INDEX 0
2354 #endif
2355 #include "mapping_fe_field.inst"
2356 
2357 
Transformed quadrature weights.
unsigned int max_level() const
Shape function values.
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Tensor< 1, spacedim, Number > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number > &grad_F, const Tensor< 1, dim, Number > &d_x)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
void compute_face_data(const UpdateFlags update_flags, const Quadrature< dim > &q, const unsigned int n_original_q_points, InternalData &data) const
Contravariant transformation.
std::mutex fe_values_mutex
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=StaticMappingQ1< dim, spacedim >::mapping, const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:444
const std::vector< Point< dim > > & get_points() const
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
const std::vector< double > & get_weights() const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
virtual void compute_shapes_virtual(const std::vector< Point< dim >> &unit_points, typename MappingFEField< dim, spacedim, VectorType, DoFHandlerType >::InternalData &data) const
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
std::vector< unsigned int > fe_to_real
Determinant of the Jacobian.
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
static Point< dim > unit_cell_vertex(const unsigned int vertex)
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
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
const ComponentMask fe_mask
Point< 2 > second
Definition: grid_out.cc:4353
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:342
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Point< spacedim > do_transform_unit_to_real_cell(const InternalData &mdata) const
unsigned int min_level() const
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2430
unsigned int get_degree() const
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< 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
std::vector< Tensor< 3, dim > > shape_third_derivatives
std::vector< Tensor< 1, dim > > shape_derivatives
static ::ExceptionBase & ExcMessage(std::string arg1)
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< Tensor< 2, dim > > shape_second_derivatives
#define Assert(cond, exc)
Definition: exceptions.h:1403
UpdateFlags
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:301
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:664
std::vector< double > volume_elements
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
unsigned int size() const
Point< 3 > vertices[4]
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
Gradient of volume element.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
std::vector< double > local_dof_values
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
unsigned int size() const
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< double > shape_values
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< std::vector< Tensor< 1, spacedim > > > aux
friend class MappingFEField
std::vector< Point< spacedim > > quadrature_points
static Point< dim > project_to_unit_cell(const Point< dim > &p)
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
void compute_data(const UpdateFlags update_flags, const Quadrature< dim > &q, const unsigned int n_original_q_points, InternalData &data) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3062
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
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
const bool uses_level_dofs
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
Normal vectors.
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
SmartPointer< const DoFHandlerType, MappingFEField< dim, spacedim, VectorType, DoFHandlerType > > euler_dof_handler
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, DoFHandlerType > > > euler_vector
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
static ::ExceptionBase & ExcNotImplemented()
ComponentMask get_component_mask() const
void update_internal_dofs(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename MappingFEField< dim, spacedim, VectorType, DoFHandlerType >::InternalData &data) const
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:371
virtual bool preserves_vertex_locations() const override
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
numbers::NumberTraits< Number >::real_type norm() const
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
static ::ExceptionBase & ExcInactiveCell()
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2405
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
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
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
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Covariant transformation.
std::vector< Tensor< 1, spacedim > > normal_vectors
virtual std::size_t memory_consumption() const override
InternalData(const FiniteElement< dim, spacedim > &fe, const ComponentMask &mask)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
UpdateFlags update_each
Definition: mapping.h:625
static ::ExceptionBase & ExcInternalError()
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: quadrature.cc:1139