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