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