Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mapping_fe_field.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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 
16 #include <deal.II/base/array_view.h>
17 #include <deal.II/base/memory_consumption.h>
18 #include <deal.II/base/polynomial.h>
19 #include <deal.II/base/qprojector.h>
20 #include <deal.II/base/quadrature.h>
21 #include <deal.II/base/quadrature_lib.h>
22 #include <deal.II/base/std_cxx14/memory.h>
23 #include <deal.II/base/tensor_product_polynomials.h>
24 #include <deal.II/base/utilities.h>
25 
26 #include <deal.II/dofs/dof_accessor.h>
27 
28 #include <deal.II/fe/fe_q.h>
29 #include <deal.II/fe/fe_system.h>
30 #include <deal.II/fe/fe_tools.h>
31 #include <deal.II/fe/fe_values.h>
32 #include <deal.II/fe/mapping.h>
33 #include <deal.II/fe/mapping_fe_field.h>
34 #include <deal.II/fe/mapping_q1.h>
35 
36 #include <deal.II/grid/tria_iterator.h>
37 
38 #include <deal.II/lac/block_vector.h>
39 #include <deal.II/lac/full_matrix.h>
40 #include <deal.II/lac/la_parallel_block_vector.h>
41 #include <deal.II/lac/la_parallel_vector.h>
42 #include <deal.II/lac/la_vector.h>
43 #include <deal.II/lac/petsc_block_vector.h>
44 #include <deal.II/lac/petsc_vector.h>
45 #include <deal.II/lac/trilinos_epetra_vector.h>
46 #include <deal.II/lac/trilinos_parallel_block_vector.h>
47 #include <deal.II/lac/trilinos_tpetra_vector.h>
48 #include <deal.II/lac/trilinos_vector.h>
49 #include <deal.II/lac/vector.h>
50 
51 #include <deal.II/numerics/vector_tools.h>
52 
53 #include <fstream>
54 #include <memory>
55 #include <numeric>
56 
57 
58 
59 DEAL_II_NAMESPACE_OPEN
60 
61 
62 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
65  const ComponentMask & mask)
66  : unit_tangentials()
67  , n_shape_functions(fe.dofs_per_cell)
68  , mask(mask)
69  , local_dof_indices(fe.dofs_per_cell)
70  , local_dof_values(fe.dofs_per_cell)
71 {}
72 
73 
74 
75 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
76 std::size_t
79 {
80  Assert(false, ExcNotImplemented());
81  return 0;
82 }
83 
84 
85 
86 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
87 double &
89  const unsigned int qpoint,
90  const unsigned int shape_nr)
91 {
92  Assert(qpoint * n_shape_functions + shape_nr < shape_values.size(),
93  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
94  0,
95  shape_values.size()));
96  return shape_values[qpoint * n_shape_functions + shape_nr];
97 }
98 
99 
100 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
101 const Tensor<1, dim> &
103  derivative(const unsigned int qpoint, const unsigned int shape_nr) const
104 {
105  Assert(qpoint * n_shape_functions + shape_nr < shape_derivatives.size(),
106  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
107  0,
108  shape_derivatives.size()));
109  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
110 }
111 
112 
113 
114 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
117  derivative(const unsigned int qpoint, const unsigned int shape_nr)
118 {
119  Assert(qpoint * n_shape_functions + shape_nr < shape_derivatives.size(),
120  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
121  0,
122  shape_derivatives.size()));
123  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
124 }
125 
126 
127 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
128 const Tensor<2, dim> &
130  second_derivative(const unsigned int qpoint,
131  const unsigned int shape_nr) const
132 {
133  Assert(qpoint * n_shape_functions + shape_nr <
135  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
136  0,
137  shape_second_derivatives.size()));
138  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
139 }
140 
141 
142 
143 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
146  second_derivative(const unsigned int qpoint, const unsigned int shape_nr)
147 {
148  Assert(qpoint * n_shape_functions + shape_nr <
150  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
151  0,
152  shape_second_derivatives.size()));
153  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
154 }
155 
156 
157 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
158 const Tensor<3, dim> &
160  third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
161 {
162  Assert(qpoint * n_shape_functions + shape_nr < shape_third_derivatives.size(),
163  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
164  0,
165  shape_third_derivatives.size()));
166  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
167 }
168 
169 
170 
171 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
174  third_derivative(const unsigned int qpoint, const unsigned int shape_nr)
175 {
176  Assert(qpoint * n_shape_functions + shape_nr < shape_third_derivatives.size(),
177  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
178  0,
179  shape_third_derivatives.size()));
180  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
181 }
182 
183 
184 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
185 const Tensor<4, dim> &
187  fourth_derivative(const unsigned int qpoint,
188  const unsigned int shape_nr) const
189 {
190  Assert(qpoint * n_shape_functions + shape_nr <
192  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
193  0,
194  shape_fourth_derivatives.size()));
195  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
196 }
197 
198 
199 
200 template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
203  fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr)
204 {
205  Assert(qpoint * n_shape_functions + shape_nr <
207  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
208  0,
209  shape_fourth_derivatives.size()));
210  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
211 }
212 
213 
214 
215 namespace
216 {
217  // Construct a quadrature formula containing the vertices of the reference
218  // cell in dimension dim (with invalid weights)
219  template <int dim>
221  get_vertex_quadrature()
222  {
223  static Quadrature<dim> quad;
224  if (quad.size() == 0)
225  {
226  std::vector<Point<dim>> points(GeometryInfo<dim>::vertices_per_cell);
227  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
229  quad = Quadrature<dim>(points);
230  }
231  return quad;
232  }
233 } // namespace
234 
235 
236 
237 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
239  const DoFHandlerType &euler_dof_handler,
240  const VectorType & euler_vector,
241  const ComponentMask & mask)
242  : uses_level_dofs(false)
243  , euler_vector({&euler_vector})
245  , fe_mask(mask.size() ?
246  mask :
248  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
249  true))
251  , fe_values(this->euler_dof_handler->get_fe(),
252  get_vertex_quadrature<dim>(),
254 {
255  unsigned int size = 0;
256  for (unsigned int i = 0; i < fe_mask.size(); ++i)
257  {
258  if (fe_mask[i])
259  fe_to_real[i] = size++;
260  }
261  AssertDimension(size, spacedim);
262 }
263 
264 
265 
266 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
268  const DoFHandlerType & euler_dof_handler,
269  const std::vector<VectorType> &euler_vector,
270  const ComponentMask & mask)
271  : uses_level_dofs(true)
272  , euler_dof_handler(&euler_dof_handler)
273  , fe_mask(mask.size() ?
274  mask :
276  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
277  true))
278  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
279  , fe_values(this->euler_dof_handler->get_fe(),
280  get_vertex_quadrature<dim>(),
282 {
283  unsigned int size = 0;
284  for (unsigned int i = 0; i < fe_mask.size(); ++i)
285  {
286  if (fe_mask[i])
287  fe_to_real[i] = size++;
288  }
289  AssertDimension(size, spacedim);
290 
291  Assert(euler_dof_handler.has_level_dofs(),
292  ExcMessage("The underlying DoFHandler object did not call "
293  "distribute_mg_dofs(). In this case, the construction via "
294  "level vectors does not make sense."));
295  AssertDimension(euler_vector.size(),
296  euler_dof_handler.get_triangulation().n_global_levels());
297  this->euler_vector.clear();
298  this->euler_vector.resize(euler_vector.size());
299  for (unsigned int i = 0; i < euler_vector.size(); ++i)
300  this->euler_vector[i] = &euler_vector[i];
301 }
302 
303 
304 
305 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
307  const DoFHandlerType & euler_dof_handler,
309  const ComponentMask & mask)
310  : uses_level_dofs(true)
311  , euler_dof_handler(&euler_dof_handler)
312  , fe_mask(mask.size() ?
313  mask :
315  euler_dof_handler.get_fe().get_nonzero_components(0).size(),
316  true))
317  , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
318  , fe_values(this->euler_dof_handler->get_fe(),
319  get_vertex_quadrature<dim>(),
321 {
322  unsigned int size = 0;
323  for (unsigned int i = 0; i < fe_mask.size(); ++i)
324  {
325  if (fe_mask[i])
326  fe_to_real[i] = size++;
327  }
328  AssertDimension(size, spacedim);
329 
330  Assert(euler_dof_handler.has_level_dofs(),
331  ExcMessage("The underlying DoFHandler object did not call "
332  "distribute_mg_dofs(). In this case, the construction via "
333  "level vectors does not make sense."));
334  AssertDimension(euler_vector.max_level() + 1,
335  euler_dof_handler.get_triangulation().n_global_levels());
336  this->euler_vector.clear();
337  this->euler_vector.resize(
338  euler_dof_handler.get_triangulation().n_global_levels());
339  for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
340  ++i)
341  this->euler_vector[i] = &euler_vector[i];
342 }
343 
344 
345 
346 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
350  , euler_vector(mapping.euler_vector)
352  , fe_mask(mapping.fe_mask)
353  , fe_to_real(mapping.fe_to_real)
354  , fe_values(mapping.euler_dof_handler->get_fe(),
355  get_vertex_quadrature<dim>(),
357 {}
358 
359 
360 
361 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
362 inline const double &
364  const unsigned int qpoint,
365  const unsigned int shape_nr) const
366 {
367  Assert(qpoint * n_shape_functions + shape_nr < shape_values.size(),
368  ExcIndexRange(qpoint * n_shape_functions + shape_nr,
369  0,
370  shape_values.size()));
371  return shape_values[qpoint * n_shape_functions + shape_nr];
372 }
373 
374 
375 
376 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
377 bool
380 {
381  return false;
382 }
383 
384 
385 
386 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
387 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
389  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
390 {
391  // we transform our tria iterator into a dof iterator so we can access
392  // data not associated with triangulations
393  const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
394  *cell, euler_dof_handler);
395 
396  Assert(uses_level_dofs || dof_cell->active() == true, ExcInactiveCell());
398  fe_values.n_quadrature_points);
400  euler_dof_handler->get_fe().n_components());
401  if (uses_level_dofs)
402  {
403  AssertIndexRange(cell->level(), euler_vector.size());
404  AssertDimension(euler_vector[cell->level()]->size(),
405  euler_dof_handler->n_dofs(cell->level()));
406  }
407  else
408  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
409 
410  {
411  std::lock_guard<std::mutex> lock(fe_values_mutex);
412  fe_values.reinit(dof_cell);
413  }
414  const unsigned int dofs_per_cell = euler_dof_handler->get_fe().dofs_per_cell;
415  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
416  if (uses_level_dofs)
417  dof_cell->get_mg_dof_indices(dof_indices);
418  else
419  dof_cell->get_dof_indices(dof_indices);
420 
421  const VectorType &vector =
422  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
423 
424  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
425  for (unsigned int i = 0; i < dofs_per_cell; ++i)
426  {
427  const unsigned int comp = fe_to_real
428  [euler_dof_handler->get_fe().system_to_component_index(i).first];
429  if (comp != numbers::invalid_unsigned_int)
430  {
431  typename VectorType::value_type value =
432  internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
433  if (euler_dof_handler->get_fe().is_primitive(i))
434  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
435  ++v)
436  vertices[v][comp] += fe_values.shape_value(i, v) * value;
437  else
438  Assert(false, ExcNotImplemented());
439  }
440  }
441 
442  return vertices;
443 }
444 
445 
446 
447 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
448 void
451  const std::vector<Point<dim>> &unit_points,
453  InternalData &data) const
454 {
455  const auto fe = &euler_dof_handler->get_fe();
456  const unsigned int n_points = unit_points.size();
457 
458  for (unsigned int point = 0; point < n_points; ++point)
459  {
460  if (data.shape_values.size() != 0)
461  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
462  data.shape(point, i) = fe->shape_value(i, unit_points[point]);
463 
464  if (data.shape_derivatives.size() != 0)
465  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
466  data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
467 
468  if (data.shape_second_derivatives.size() != 0)
469  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
470  data.second_derivative(point, i) =
471  fe->shape_grad_grad(i, unit_points[point]);
472 
473  if (data.shape_third_derivatives.size() != 0)
474  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
475  data.third_derivative(point, i) =
476  fe->shape_3rd_derivative(i, unit_points[point]);
477 
478  if (data.shape_fourth_derivatives.size() != 0)
479  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
480  data.fourth_derivative(point, i) =
481  fe->shape_4th_derivative(i, unit_points[point]);
482  }
483 }
484 
485 
486 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
490 {
491  // add flags if the respective quantities are necessary to compute
492  // what we need. note that some flags appear in both conditions and
493  // in subsequent set operations. this leads to some circular
494  // logic. the only way to treat this is to iterate. since there are
495  // 5 if-clauses in the loop, it will take at most 4 iterations to
496  // converge. do them:
497  UpdateFlags out = in;
498  for (unsigned int i = 0; i < 5; ++i)
499  {
500  // The following is a little incorrect:
501  // If not applied on a face,
502  // update_boundary_forms does not
503  // make sense. On the other hand,
504  // it is necessary on a
505  // face. Currently,
506  // update_boundary_forms is simply
507  // ignored for the interior of a
508  // cell.
510  out |= update_boundary_forms;
511 
516 
517  if (out &
522 
523  // The contravariant transformation
524  // is a Piola transformation, which
525  // requires the determinant of the
526  // Jacobi matrix of the transformation.
527  // Therefore these values have to be
528  // updated for each cell.
530  out |= update_JxW_values;
531 
532  if (out & update_normal_vectors)
533  out |= update_JxW_values;
534  }
535 
536  return out;
537 }
538 
539 
540 
541 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
542 void
544  const UpdateFlags update_flags,
545  const Quadrature<dim> &q,
546  const unsigned int n_original_q_points,
547  InternalData & data) const
548 {
549  // store the flags in the internal data object so we can access them
550  // in fill_fe_*_values(). use the transitive hull of the required
551  // flags
552  data.update_each = requires_update_flags(update_flags);
553 
554  const unsigned int n_q_points = q.size();
555 
556  // see if we need the (transformation) shape function values
557  // and/or gradients and resize the necessary arrays
559  data.shape_values.resize(data.n_shape_functions * n_q_points);
560 
561  if (data.update_each &
565  data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
566 
568  data.covariant.resize(n_original_q_points);
569 
571  data.contravariant.resize(n_original_q_points);
572 
574  data.volume_elements.resize(n_original_q_points);
575 
576  if (data.update_each &
578  data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
579 
582  data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
583 
586  data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
587 
589 }
590 
591 
592 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
593 void
595  const UpdateFlags update_flags,
596  const Quadrature<dim> &q,
597  const unsigned int n_original_q_points,
598  InternalData & data) const
599 {
600  compute_data(update_flags, q, n_original_q_points, data);
601 
602  if (dim > 1)
603  {
605  {
606  data.aux.resize(
607  dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
608 
609  // Compute tangentials to the unit cell.
610  for (unsigned int i = 0; i < data.unit_tangentials.size(); ++i)
611  data.unit_tangentials[i].resize(n_original_q_points);
612 
613  if (dim == 2)
614  {
615  // ensure a counterclockwise
616  // orientation of tangentials
617  static const int tangential_orientation[4] = {-1, 1, 1, -1};
618  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
619  ++i)
620  {
621  Tensor<1, dim> tang;
622  tang[1 - i / 2] = tangential_orientation[i];
623  std::fill(data.unit_tangentials[i].begin(),
624  data.unit_tangentials[i].end(),
625  tang);
626  }
627  }
628  else if (dim == 3)
629  {
630  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
631  ++i)
632  {
633  Tensor<1, dim> tang1, tang2;
634 
635  const unsigned int nd =
637 
638  // first tangential
639  // vector in direction
640  // of the (nd+1)%3 axis
641  // and inverted in case
642  // of unit inward normal
643  tang1[(nd + 1) % dim] =
645  // second tangential
646  // vector in direction
647  // of the (nd+2)%3 axis
648  tang2[(nd + 2) % dim] = 1.;
649 
650  // same unit tangents
651  // for all quadrature
652  // points on this face
653  std::fill(data.unit_tangentials[i].begin(),
654  data.unit_tangentials[i].end(),
655  tang1);
656  std::fill(
658  .begin(),
660  .end(),
661  tang2);
662  }
663  }
664  }
665  }
666 }
667 
668 
669 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
670 typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
672  const UpdateFlags update_flags,
673  const Quadrature<dim> &quadrature) const
674 {
675  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
676  std_cxx14::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
677  auto &data = dynamic_cast<InternalData &>(*data_ptr);
678  this->compute_data(update_flags, quadrature, quadrature.size(), data);
679 
680  return data_ptr;
681 }
682 
683 
684 
685 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
686 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
688  const UpdateFlags update_flags,
689  const Quadrature<dim - 1> &quadrature) const
690 {
691  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
692  std_cxx14::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
693  auto & data = dynamic_cast<InternalData &>(*data_ptr);
695  this->compute_face_data(update_flags, q, quadrature.size(), data);
696 
697  return data_ptr;
698 }
699 
700 
701 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
702 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
704  const UpdateFlags update_flags,
705  const Quadrature<dim - 1> &quadrature) const
706 {
707  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
708  std_cxx14::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
709  auto & data = dynamic_cast<InternalData &>(*data_ptr);
711  this->compute_face_data(update_flags, q, quadrature.size(), data);
712 
713  return data_ptr;
714 }
715 
716 
717 
718 namespace internal
719 {
720  namespace MappingFEFieldImplementation
721  {
722  namespace
723  {
730  template <int dim,
731  int spacedim,
732  typename VectorType,
733  typename DoFHandlerType>
734  void
735  maybe_compute_q_points(
736  const typename ::QProjector<dim>::DataSetDescriptor data_set,
737  const typename ::
739  InternalData & data,
741  const ComponentMask & fe_mask,
742  const std::vector<unsigned int> & fe_to_real,
743  std::vector<Point<spacedim>> & quadrature_points)
744  {
745  const UpdateFlags update_flags = data.update_each;
746 
747  if (update_flags & update_quadrature_points)
748  {
749  for (unsigned int point = 0; point < quadrature_points.size();
750  ++point)
751  {
752  Point<spacedim> result;
753  const double * shape = &data.shape(point + data_set, 0);
754 
755  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
756  {
757  unsigned int comp_k = fe.system_to_component_index(k).first;
758  if (fe_mask[comp_k])
759  result[fe_to_real[comp_k]] +=
760  data.local_dof_values[k] * shape[k];
761  }
762 
763  quadrature_points[point] = result;
764  }
765  }
766  }
767 
775  template <int dim,
776  int spacedim,
777  typename VectorType,
778  typename DoFHandlerType>
779  void
780  maybe_update_Jacobians(
781  const CellSimilarity::Similarity cell_similarity,
782  const typename ::QProjector<dim>::DataSetDescriptor data_set,
783  const typename ::
785  InternalData & data,
787  const ComponentMask & fe_mask,
788  const std::vector<unsigned int> & fe_to_real)
789  {
790  const UpdateFlags update_flags = data.update_each;
791 
792  // then Jacobians
793  if (update_flags & update_contravariant_transformation)
794  {
795  // if the current cell is just a translation of the previous one, no
796  // need to recompute jacobians...
797  if (cell_similarity != CellSimilarity::translation)
798  {
799  const unsigned int n_q_points = data.contravariant.size();
800 
801  Assert(data.n_shape_functions > 0, ExcInternalError());
802 
803  for (unsigned int point = 0; point < n_q_points; ++point)
804  {
805  const Tensor<1, dim> *data_derv =
806  &data.derivative(point + data_set, 0);
807 
808  Tensor<1, dim> result[spacedim];
809 
810  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
811  {
812  unsigned int comp_k =
813  fe.system_to_component_index(k).first;
814  if (fe_mask[comp_k])
815  result[fe_to_real[comp_k]] +=
816  data.local_dof_values[k] * data_derv[k];
817  }
818 
819  // write result into contravariant data
820  for (unsigned int i = 0; i < spacedim; ++i)
821  {
822  data.contravariant[point][i] = result[i];
823  }
824  }
825  }
826  }
827 
828  if (update_flags & update_covariant_transformation)
829  {
830  AssertDimension(data.covariant.size(), data.contravariant.size());
831  if (cell_similarity != CellSimilarity::translation)
832  for (unsigned int point = 0; point < data.contravariant.size();
833  ++point)
834  data.covariant[point] =
835  (data.contravariant[point]).covariant_form();
836  }
837 
838  if (update_flags & update_volume_elements)
839  {
840  AssertDimension(data.contravariant.size(),
841  data.volume_elements.size());
842  if (cell_similarity != CellSimilarity::translation)
843  for (unsigned int point = 0; point < data.contravariant.size();
844  ++point)
845  data.volume_elements[point] =
846  data.contravariant[point].determinant();
847  }
848  }
849 
856  template <int dim,
857  int spacedim,
858  typename VectorType,
859  typename DoFHandlerType>
860  void
861  maybe_update_jacobian_grads(
862  const CellSimilarity::Similarity cell_similarity,
863  const typename ::QProjector<dim>::DataSetDescriptor data_set,
864  const typename ::
866  InternalData & data,
867  const FiniteElement<dim, spacedim> & fe,
868  const ComponentMask & fe_mask,
869  const std::vector<unsigned int> & fe_to_real,
870  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
871  {
872  const UpdateFlags update_flags = data.update_each;
873  if (update_flags & update_jacobian_grads)
874  {
875  const unsigned int n_q_points = jacobian_grads.size();
876 
877  if (cell_similarity != CellSimilarity::translation)
878  {
879  for (unsigned int point = 0; point < n_q_points; ++point)
880  {
881  const Tensor<2, dim> *second =
882  &data.second_derivative(point + data_set, 0);
883 
885 
886  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
887  {
888  unsigned int comp_k =
889  fe.system_to_component_index(k).first;
890  if (fe_mask[comp_k])
891  for (unsigned int j = 0; j < dim; ++j)
892  for (unsigned int l = 0; l < dim; ++l)
893  result[fe_to_real[comp_k]][j][l] +=
894  (second[k][j][l] * data.local_dof_values[k]);
895  }
896 
897  // never touch any data for j=dim in case dim<spacedim, so
898  // it will always be zero as it was initialized
899  for (unsigned int i = 0; i < spacedim; ++i)
900  for (unsigned int j = 0; j < dim; ++j)
901  for (unsigned int l = 0; l < dim; ++l)
902  jacobian_grads[point][i][j][l] = result[i][j][l];
903  }
904  }
905  }
906  }
907 
914  template <int dim,
915  int spacedim,
916  typename VectorType,
917  typename DoFHandlerType>
918  void
919  maybe_update_jacobian_pushed_forward_grads(
920  const CellSimilarity::Similarity cell_similarity,
921  const typename ::QProjector<dim>::DataSetDescriptor data_set,
922  const typename ::
924  InternalData & data,
926  const ComponentMask & fe_mask,
927  const std::vector<unsigned int> & fe_to_real,
928  std::vector<Tensor<3, spacedim>> & jacobian_pushed_forward_grads)
929  {
930  const UpdateFlags update_flags = data.update_each;
931  if (update_flags & update_jacobian_pushed_forward_grads)
932  {
933  const unsigned int n_q_points =
934  jacobian_pushed_forward_grads.size();
935 
936  if (cell_similarity != CellSimilarity::translation)
937  {
938  double tmp[spacedim][spacedim][spacedim];
939  for (unsigned int point = 0; point < n_q_points; ++point)
940  {
941  const Tensor<2, dim> *second =
942  &data.second_derivative(point + data_set, 0);
943 
945 
946  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
947  {
948  unsigned int comp_k =
949  fe.system_to_component_index(k).first;
950  if (fe_mask[comp_k])
951  for (unsigned int j = 0; j < dim; ++j)
952  for (unsigned int l = 0; l < dim; ++l)
953  result[fe_to_real[comp_k]][j][l] +=
954  (second[k][j][l] * data.local_dof_values[k]);
955  }
956 
957  // first push forward the j-components
958  for (unsigned int i = 0; i < spacedim; ++i)
959  for (unsigned int j = 0; j < spacedim; ++j)
960  for (unsigned int l = 0; l < dim; ++l)
961  {
962  tmp[i][j][l] =
963  result[i][0][l] * data.covariant[point][j][0];
964  for (unsigned int jr = 1; jr < dim; ++jr)
965  {
966  tmp[i][j][l] += result[i][jr][l] *
967  data.covariant[point][j][jr];
968  }
969  }
970 
971  // now, pushing forward the l-components
972  for (unsigned int i = 0; i < spacedim; ++i)
973  for (unsigned int j = 0; j < spacedim; ++j)
974  for (unsigned int l = 0; l < spacedim; ++l)
975  {
976  jacobian_pushed_forward_grads[point][i][j][l] =
977  tmp[i][j][0] * data.covariant[point][l][0];
978  for (unsigned int lr = 1; lr < dim; ++lr)
979  {
980  jacobian_pushed_forward_grads[point][i][j][l] +=
981  tmp[i][j][lr] * data.covariant[point][l][lr];
982  }
983  }
984  }
985  }
986  }
987  }
988 
995  template <int dim,
996  int spacedim,
997  typename VectorType,
998  typename DoFHandlerType>
999  void
1000  maybe_update_jacobian_2nd_derivatives(
1001  const CellSimilarity::Similarity cell_similarity,
1002  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1003  const typename ::
1005  InternalData & data,
1006  const FiniteElement<dim, spacedim> & fe,
1007  const ComponentMask & fe_mask,
1008  const std::vector<unsigned int> & fe_to_real,
1009  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1010  {
1011  const UpdateFlags update_flags = data.update_each;
1012  if (update_flags & update_jacobian_2nd_derivatives)
1013  {
1014  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1015 
1016  if (cell_similarity != CellSimilarity::translation)
1017  {
1018  for (unsigned int point = 0; point < n_q_points; ++point)
1019  {
1020  const Tensor<3, dim> *third =
1021  &data.third_derivative(point + data_set, 0);
1022 
1024 
1025  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1026  {
1027  unsigned int comp_k =
1028  fe.system_to_component_index(k).first;
1029  if (fe_mask[comp_k])
1030  for (unsigned int j = 0; j < dim; ++j)
1031  for (unsigned int l = 0; l < dim; ++l)
1032  for (unsigned int m = 0; m < dim; ++m)
1033  result[fe_to_real[comp_k]][j][l][m] +=
1034  (third[k][j][l][m] *
1035  data.local_dof_values[k]);
1036  }
1037 
1038  // never touch any data for j=dim in case dim<spacedim, so
1039  // it will always be zero as it was initialized
1040  for (unsigned int i = 0; i < spacedim; ++i)
1041  for (unsigned int j = 0; j < dim; ++j)
1042  for (unsigned int l = 0; l < dim; ++l)
1043  for (unsigned int m = 0; m < dim; ++m)
1044  jacobian_2nd_derivatives[point][i][j][l][m] =
1045  result[i][j][l][m];
1046  }
1047  }
1048  }
1049  }
1050 
1058  template <int dim,
1059  int spacedim,
1060  typename VectorType,
1061  typename DoFHandlerType>
1062  void
1063  maybe_update_jacobian_pushed_forward_2nd_derivatives(
1064  const CellSimilarity::Similarity cell_similarity,
1065  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1066  const typename ::
1068  InternalData & data,
1069  const FiniteElement<dim, spacedim> &fe,
1070  const ComponentMask & fe_mask,
1071  const std::vector<unsigned int> & fe_to_real,
1072  std::vector<Tensor<4, spacedim>>
1073  &jacobian_pushed_forward_2nd_derivatives)
1074  {
1075  const UpdateFlags update_flags = data.update_each;
1077  {
1078  const unsigned int n_q_points =
1079  jacobian_pushed_forward_2nd_derivatives.size();
1080 
1081  if (cell_similarity != CellSimilarity::translation)
1082  {
1083  double tmp[spacedim][spacedim][spacedim][spacedim];
1084  for (unsigned int point = 0; point < n_q_points; ++point)
1085  {
1086  const Tensor<3, dim> *third =
1087  &data.third_derivative(point + data_set, 0);
1088 
1090 
1091  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1092  {
1093  unsigned int comp_k =
1094  fe.system_to_component_index(k).first;
1095  if (fe_mask[comp_k])
1096  for (unsigned int j = 0; j < dim; ++j)
1097  for (unsigned int l = 0; l < dim; ++l)
1098  for (unsigned int m = 0; m < dim; ++m)
1099  result[fe_to_real[comp_k]][j][l][m] +=
1100  (third[k][j][l][m] *
1101  data.local_dof_values[k]);
1102  }
1103 
1104  // push forward the j-coordinate
1105  for (unsigned int i = 0; i < spacedim; ++i)
1106  for (unsigned int j = 0; j < spacedim; ++j)
1107  for (unsigned int l = 0; l < dim; ++l)
1108  for (unsigned int m = 0; m < dim; ++m)
1109  {
1110  jacobian_pushed_forward_2nd_derivatives
1111  [point][i][j][l][m] =
1112  result[i][0][l][m] *
1113  data.covariant[point][j][0];
1114  for (unsigned int jr = 1; jr < dim; ++jr)
1115  jacobian_pushed_forward_2nd_derivatives[point]
1116  [i][j][l]
1117  [m] +=
1118  result[i][jr][l][m] *
1119  data.covariant[point][j][jr];
1120  }
1121 
1122  // push forward the l-coordinate
1123  for (unsigned int i = 0; i < spacedim; ++i)
1124  for (unsigned int j = 0; j < spacedim; ++j)
1125  for (unsigned int l = 0; l < spacedim; ++l)
1126  for (unsigned int m = 0; m < dim; ++m)
1127  {
1128  tmp[i][j][l][m] =
1129  jacobian_pushed_forward_2nd_derivatives[point]
1130  [i][j][0]
1131  [m] *
1132  data.covariant[point][l][0];
1133  for (unsigned int lr = 1; lr < dim; ++lr)
1134  tmp[i][j][l][m] +=
1135  jacobian_pushed_forward_2nd_derivatives
1136  [point][i][j][lr][m] *
1137  data.covariant[point][l][lr];
1138  }
1139 
1140  // push forward the m-coordinate
1141  for (unsigned int i = 0; i < spacedim; ++i)
1142  for (unsigned int j = 0; j < spacedim; ++j)
1143  for (unsigned int l = 0; l < spacedim; ++l)
1144  for (unsigned int m = 0; m < spacedim; ++m)
1145  {
1146  jacobian_pushed_forward_2nd_derivatives
1147  [point][i][j][l][m] =
1148  tmp[i][j][l][0] * data.covariant[point][m][0];
1149  for (unsigned int mr = 1; mr < dim; ++mr)
1150  jacobian_pushed_forward_2nd_derivatives[point]
1151  [i][j][l]
1152  [m] +=
1153  tmp[i][j][l][mr] *
1154  data.covariant[point][m][mr];
1155  }
1156  }
1157  }
1158  }
1159  }
1160 
1167  template <int dim,
1168  int spacedim,
1169  typename VectorType,
1170  typename DoFHandlerType>
1171  void
1172  maybe_update_jacobian_3rd_derivatives(
1173  const CellSimilarity::Similarity cell_similarity,
1174  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1175  const typename ::
1177  InternalData & data,
1178  const FiniteElement<dim, spacedim> & fe,
1179  const ComponentMask & fe_mask,
1180  const std::vector<unsigned int> & fe_to_real,
1181  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1182  {
1183  const UpdateFlags update_flags = data.update_each;
1184  if (update_flags & update_jacobian_3rd_derivatives)
1185  {
1186  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1187 
1188  if (cell_similarity != CellSimilarity::translation)
1189  {
1190  for (unsigned int point = 0; point < n_q_points; ++point)
1191  {
1192  const Tensor<4, dim> *fourth =
1193  &data.fourth_derivative(point + data_set, 0);
1194 
1196 
1197  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1198  {
1199  unsigned int comp_k =
1200  fe.system_to_component_index(k).first;
1201  if (fe_mask[comp_k])
1202  for (unsigned int j = 0; j < dim; ++j)
1203  for (unsigned int l = 0; l < dim; ++l)
1204  for (unsigned int m = 0; m < dim; ++m)
1205  for (unsigned int n = 0; n < dim; ++n)
1206  result[fe_to_real[comp_k]][j][l][m][n] +=
1207  (fourth[k][j][l][m][n] *
1208  data.local_dof_values[k]);
1209  }
1210 
1211  // never touch any data for j,l,m,n=dim in case
1212  // dim<spacedim, so it will always be zero as it was
1213  // initialized
1214  for (unsigned int i = 0; i < spacedim; ++i)
1215  for (unsigned int j = 0; j < dim; ++j)
1216  for (unsigned int l = 0; l < dim; ++l)
1217  for (unsigned int m = 0; m < dim; ++m)
1218  for (unsigned int n = 0; n < dim; ++n)
1219  jacobian_3rd_derivatives[point][i][j][l][m][n] =
1220  result[i][j][l][m][n];
1221  }
1222  }
1223  }
1224  }
1225 
1233  template <int dim,
1234  int spacedim,
1235  typename VectorType,
1236  typename DoFHandlerType>
1237  void
1238  maybe_update_jacobian_pushed_forward_3rd_derivatives(
1239  const CellSimilarity::Similarity cell_similarity,
1240  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1241  const typename ::
1243  InternalData & data,
1244  const FiniteElement<dim, spacedim> &fe,
1245  const ComponentMask & fe_mask,
1246  const std::vector<unsigned int> & fe_to_real,
1247  std::vector<Tensor<5, spacedim>>
1248  &jacobian_pushed_forward_3rd_derivatives)
1249  {
1250  const UpdateFlags update_flags = data.update_each;
1252  {
1253  const unsigned int n_q_points =
1254  jacobian_pushed_forward_3rd_derivatives.size();
1255 
1256  if (cell_similarity != CellSimilarity::translation)
1257  {
1258  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1259  for (unsigned int point = 0; point < n_q_points; ++point)
1260  {
1261  const Tensor<4, dim> *fourth =
1262  &data.fourth_derivative(point + data_set, 0);
1263 
1265 
1266  for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1267  {
1268  unsigned int comp_k =
1269  fe.system_to_component_index(k).first;
1270  if (fe_mask[comp_k])
1271  for (unsigned int j = 0; j < dim; ++j)
1272  for (unsigned int l = 0; l < dim; ++l)
1273  for (unsigned int m = 0; m < dim; ++m)
1274  for (unsigned int n = 0; n < dim; ++n)
1275  result[fe_to_real[comp_k]][j][l][m][n] +=
1276  (fourth[k][j][l][m][n] *
1277  data.local_dof_values[k]);
1278  }
1279 
1280  // push-forward the j-coordinate
1281  for (unsigned int i = 0; i < spacedim; ++i)
1282  for (unsigned int j = 0; j < spacedim; ++j)
1283  for (unsigned int l = 0; l < dim; ++l)
1284  for (unsigned int m = 0; m < dim; ++m)
1285  for (unsigned int n = 0; n < dim; ++n)
1286  {
1287  tmp[i][j][l][m][n] =
1288  result[i][0][l][m][n] *
1289  data.covariant[point][j][0];
1290  for (unsigned int jr = 1; jr < dim; ++jr)
1291  tmp[i][j][l][m][n] +=
1292  result[i][jr][l][m][n] *
1293  data.covariant[point][j][jr];
1294  }
1295 
1296  // push-forward the l-coordinate
1297  for (unsigned int i = 0; i < spacedim; ++i)
1298  for (unsigned int j = 0; j < spacedim; ++j)
1299  for (unsigned int l = 0; l < spacedim; ++l)
1300  for (unsigned int m = 0; m < dim; ++m)
1301  for (unsigned int n = 0; n < dim; ++n)
1302  {
1303  jacobian_pushed_forward_3rd_derivatives
1304  [point][i][j][l][m][n] =
1305  tmp[i][j][0][m][n] *
1306  data.covariant[point][l][0];
1307  for (unsigned int lr = 1; lr < dim; ++lr)
1308  jacobian_pushed_forward_3rd_derivatives
1309  [point][i][j][l][m][n] +=
1310  tmp[i][j][lr][m][n] *
1311  data.covariant[point][l][lr];
1312  }
1313 
1314  // push-forward the m-coordinate
1315  for (unsigned int i = 0; i < spacedim; ++i)
1316  for (unsigned int j = 0; j < spacedim; ++j)
1317  for (unsigned int l = 0; l < spacedim; ++l)
1318  for (unsigned int m = 0; m < spacedim; ++m)
1319  for (unsigned int n = 0; n < dim; ++n)
1320  {
1321  tmp[i][j][l][m][n] =
1322  jacobian_pushed_forward_3rd_derivatives
1323  [point][i][j][l][0][n] *
1324  data.covariant[point][m][0];
1325  for (unsigned int mr = 1; mr < dim; ++mr)
1326  tmp[i][j][l][m][n] +=
1327  jacobian_pushed_forward_3rd_derivatives
1328  [point][i][j][l][mr][n] *
1329  data.covariant[point][m][mr];
1330  }
1331 
1332  // push-forward the n-coordinate
1333  for (unsigned int i = 0; i < spacedim; ++i)
1334  for (unsigned int j = 0; j < spacedim; ++j)
1335  for (unsigned int l = 0; l < spacedim; ++l)
1336  for (unsigned int m = 0; m < spacedim; ++m)
1337  for (unsigned int n = 0; n < spacedim; ++n)
1338  {
1339  jacobian_pushed_forward_3rd_derivatives
1340  [point][i][j][l][m][n] =
1341  tmp[i][j][l][m][0] *
1342  data.covariant[point][n][0];
1343  for (unsigned int nr = 1; nr < dim; ++nr)
1344  jacobian_pushed_forward_3rd_derivatives
1345  [point][i][j][l][m][n] +=
1346  tmp[i][j][l][m][nr] *
1347  data.covariant[point][n][nr];
1348  }
1349  }
1350  }
1351  }
1352  }
1353 
1354 
1364  template <int dim,
1365  int spacedim,
1366  typename VectorType,
1367  typename DoFHandlerType>
1368  void
1369  maybe_compute_face_data(
1370  const ::Mapping<dim, spacedim> &mapping,
1371  const typename ::Triangulation<dim, spacedim>::cell_iterator
1372  & cell,
1373  const unsigned int face_no,
1374  const unsigned int subface_no,
1375  const std::vector<double> &weights,
1376  const typename ::
1378  InternalData &data,
1380  &output_data)
1381  {
1382  const UpdateFlags update_flags = data.update_each;
1383 
1384  if (update_flags & update_boundary_forms)
1385  {
1386  const unsigned int n_q_points = output_data.boundary_forms.size();
1387  if (update_flags & update_normal_vectors)
1388  AssertDimension(output_data.normal_vectors.size(), n_q_points);
1389  if (update_flags & update_JxW_values)
1390  AssertDimension(output_data.JxW_values.size(), n_q_points);
1391 
1392  // map the unit tangentials to the real cell. checking for d!=dim-1
1393  // eliminates compiler warnings regarding unsigned int expressions <
1394  // 0.
1395  for (unsigned int d = 0; d != dim - 1; ++d)
1396  {
1398  data.unit_tangentials.size(),
1399  ExcInternalError());
1400  Assert(
1401  data.aux[d].size() <=
1402  data
1403  .unit_tangentials[face_no +
1405  .size(),
1406  ExcInternalError());
1407 
1408  mapping.transform(
1409  make_array_view(
1410  data
1411  .unit_tangentials[face_no +
1414  data,
1415  make_array_view(data.aux[d]));
1416  }
1417 
1418  // if dim==spacedim, we can use the unit tangentials to compute the
1419  // boundary form by simply taking the cross product
1420  if (dim == spacedim)
1421  {
1422  for (unsigned int i = 0; i < n_q_points; ++i)
1423  switch (dim)
1424  {
1425  case 1:
1426  // in 1d, we don't have access to any of the data.aux
1427  // fields (because it has only dim-1 components), but we
1428  // can still compute the boundary form by simply looking
1429  // at the number of the face
1430  output_data.boundary_forms[i][0] =
1431  (face_no == 0 ? -1 : +1);
1432  break;
1433  case 2:
1434  output_data.boundary_forms[i] =
1435  cross_product_2d(data.aux[0][i]);
1436  break;
1437  case 3:
1438  output_data.boundary_forms[i] =
1439  cross_product_3d(data.aux[0][i], data.aux[1][i]);
1440  break;
1441  default:
1442  Assert(false, ExcNotImplemented());
1443  }
1444  }
1445  else //(dim < spacedim)
1446  {
1447  // in the codim-one case, the boundary form results from the
1448  // cross product of all the face tangential vectors and the cell
1449  // normal vector
1450  //
1451  // to compute the cell normal, use the same method used in
1452  // fill_fe_values for cells above
1453  AssertDimension(data.contravariant.size(), n_q_points);
1454 
1455  for (unsigned int point = 0; point < n_q_points; ++point)
1456  {
1457  if (dim == 1)
1458  {
1459  // J is a tangent vector
1460  output_data.boundary_forms[point] =
1461  data.contravariant[point].transpose()[0];
1462  output_data.boundary_forms[point] /=
1463  (face_no == 0 ? -1. : +1.) *
1464  output_data.boundary_forms[point].norm();
1465  }
1466 
1467  if (dim == 2)
1468  {
1470  data.contravariant[point].transpose();
1471 
1472  Tensor<1, spacedim> cell_normal =
1473  cross_product_3d(DX_t[0], DX_t[1]);
1474  cell_normal /= cell_normal.norm();
1475 
1476  // then compute the face normal from the face tangent
1477  // and the cell normal:
1478  output_data.boundary_forms[point] =
1479  cross_product_3d(data.aux[0][point], cell_normal);
1480  }
1481  }
1482  }
1483 
1484  if (update_flags & (update_normal_vectors | update_JxW_values))
1485  for (unsigned int i = 0; i < output_data.boundary_forms.size();
1486  ++i)
1487  {
1488  if (update_flags & update_JxW_values)
1489  {
1490  output_data.JxW_values[i] =
1491  output_data.boundary_forms[i].norm() * weights[i];
1492 
1493  if (subface_no != numbers::invalid_unsigned_int)
1494  {
1495  const double area_ratio =
1497  cell->subface_case(face_no), subface_no);
1498  output_data.JxW_values[i] *= area_ratio;
1499  }
1500  }
1501 
1502  if (update_flags & update_normal_vectors)
1503  output_data.normal_vectors[i] =
1504  Point<spacedim>(output_data.boundary_forms[i] /
1505  output_data.boundary_forms[i].norm());
1506  }
1507 
1508  if (update_flags & update_jacobians)
1509  for (unsigned int point = 0; point < n_q_points; ++point)
1510  output_data.jacobians[point] = data.contravariant[point];
1511 
1512  if (update_flags & update_inverse_jacobians)
1513  for (unsigned int point = 0; point < n_q_points; ++point)
1514  output_data.inverse_jacobians[point] =
1515  data.covariant[point].transpose();
1516  }
1517  }
1518 
1525  template <int dim,
1526  int spacedim,
1527  typename VectorType,
1528  typename DoFHandlerType>
1529  void
1530  do_fill_fe_face_values(
1531  const ::Mapping<dim, spacedim> &mapping,
1532  const typename ::Triangulation<dim, spacedim>::cell_iterator
1533  & cell,
1534  const unsigned int face_no,
1535  const unsigned int subface_no,
1536  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1537  const Quadrature<dim - 1> & quadrature,
1538  const typename ::
1540  InternalData & data,
1541  const FiniteElement<dim, spacedim> &fe,
1542  const ComponentMask & fe_mask,
1543  const std::vector<unsigned int> & fe_to_real,
1545  &output_data)
1546  {
1547  maybe_compute_q_points<dim, spacedim, VectorType, DoFHandlerType>(
1548  data_set,
1549  data,
1550  fe,
1551  fe_mask,
1552  fe_to_real,
1553  output_data.quadrature_points);
1554 
1555  maybe_update_Jacobians<dim, spacedim, VectorType, DoFHandlerType>(
1556  CellSimilarity::none, data_set, data, fe, fe_mask, fe_to_real);
1557 
1558  maybe_update_jacobian_grads<dim, spacedim, VectorType, DoFHandlerType>(
1560  data_set,
1561  data,
1562  fe,
1563  fe_mask,
1564  fe_to_real,
1565  output_data.jacobian_grads);
1566 
1567  maybe_update_jacobian_pushed_forward_grads<dim,
1568  spacedim,
1569  VectorType,
1570  DoFHandlerType>(
1572  data_set,
1573  data,
1574  fe,
1575  fe_mask,
1576  fe_to_real,
1577  output_data.jacobian_pushed_forward_grads);
1578 
1579  maybe_update_jacobian_2nd_derivatives<dim,
1580  spacedim,
1581  VectorType,
1582  DoFHandlerType>(
1584  data_set,
1585  data,
1586  fe,
1587  fe_mask,
1588  fe_to_real,
1589  output_data.jacobian_2nd_derivatives);
1590 
1591  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1592  spacedim,
1593  VectorType,
1594  DoFHandlerType>(
1596  data_set,
1597  data,
1598  fe,
1599  fe_mask,
1600  fe_to_real,
1602 
1603  maybe_update_jacobian_3rd_derivatives<dim,
1604  spacedim,
1605  VectorType,
1606  DoFHandlerType>(
1608  data_set,
1609  data,
1610  fe,
1611  fe_mask,
1612  fe_to_real,
1613  output_data.jacobian_3rd_derivatives);
1614 
1615  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1616  spacedim,
1617  VectorType,
1618  DoFHandlerType>(
1620  data_set,
1621  data,
1622  fe,
1623  fe_mask,
1624  fe_to_real,
1626 
1627  maybe_compute_face_data<dim, spacedim, VectorType, DoFHandlerType>(
1628  mapping,
1629  cell,
1630  face_no,
1631  subface_no,
1632  quadrature.get_weights(),
1633  data,
1634  output_data);
1635  }
1636  } // namespace
1637  } // namespace MappingFEFieldImplementation
1638 } // namespace internal
1639 
1640 
1641 // Note that the CellSimilarity flag is modifiable, since MappingFEField can
1642 // need to recalculate data even when cells are similar.
1643 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1646  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1647  const CellSimilarity::Similarity cell_similarity,
1648  const Quadrature<dim> & quadrature,
1649  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1651  &output_data) const
1652 {
1653  // convert data object to internal data for this class. fails with an
1654  // exception if that is not possible
1655  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1656  ExcInternalError());
1657  const InternalData &data = static_cast<const InternalData &>(internal_data);
1658 
1659  const unsigned int n_q_points = quadrature.size();
1660  const CellSimilarity::Similarity updated_cell_similarity =
1661  (get_degree() == 1 ? cell_similarity : CellSimilarity::invalid_next_cell);
1662 
1663  update_internal_dofs(cell, data);
1664 
1665  internal::MappingFEFieldImplementation::
1666  maybe_compute_q_points<dim, spacedim, VectorType, DoFHandlerType>(
1668  data,
1669  euler_dof_handler->get_fe(),
1670  fe_mask,
1671  fe_to_real,
1672  output_data.quadrature_points);
1673 
1674  internal::MappingFEFieldImplementation::
1675  maybe_update_Jacobians<dim, spacedim, VectorType, DoFHandlerType>(
1676  cell_similarity,
1678  data,
1679  euler_dof_handler->get_fe(),
1680  fe_mask,
1681  fe_to_real);
1682 
1683  const UpdateFlags update_flags = data.update_each;
1684  const std::vector<double> &weights = quadrature.get_weights();
1685 
1686  // Multiply quadrature weights by absolute value of Jacobian determinants or
1687  // the area element g=sqrt(DX^t DX) in case of codim > 0
1688 
1689  if (update_flags & (update_normal_vectors | update_JxW_values))
1690  {
1691  AssertDimension(output_data.JxW_values.size(), n_q_points);
1692 
1693  Assert(!(update_flags & update_normal_vectors) ||
1694  (output_data.normal_vectors.size() == n_q_points),
1695  ExcDimensionMismatch(output_data.normal_vectors.size(),
1696  n_q_points));
1697 
1698 
1699  if (cell_similarity != CellSimilarity::translation)
1700  for (unsigned int point = 0; point < n_q_points; ++point)
1701  {
1702  if (dim == spacedim)
1703  {
1704  const double det = data.contravariant[point].determinant();
1705 
1706  // check for distorted cells.
1707 
1708  // TODO: this allows for anisotropies of up to 1e6 in 3D and
1709  // 1e12 in 2D. might want to find a finer
1710  // (dimension-independent) criterion
1711  Assert(det >
1712  1e-12 * Utilities::fixed_power<dim>(
1713  cell->diameter() / std::sqrt(double(dim))),
1715  cell->center(), det, point)));
1716  output_data.JxW_values[point] = weights[point] * det;
1717  }
1718  // if dim==spacedim, then there is no cell normal to
1719  // compute. since this is for FEValues (and not FEFaceValues),
1720  // there are also no face normals to compute
1721  else // codim>0 case
1722  {
1723  Tensor<1, spacedim> DX_t[dim];
1724  for (unsigned int i = 0; i < spacedim; ++i)
1725  for (unsigned int j = 0; j < dim; ++j)
1726  DX_t[j][i] = data.contravariant[point][i][j];
1727 
1728  Tensor<2, dim> G; // First fundamental form
1729  for (unsigned int i = 0; i < dim; ++i)
1730  for (unsigned int j = 0; j < dim; ++j)
1731  G[i][j] = DX_t[i] * DX_t[j];
1732 
1733  output_data.JxW_values[point] =
1734  std::sqrt(determinant(G)) * weights[point];
1735 
1736  if (cell_similarity == CellSimilarity::inverted_translation)
1737  {
1738  // we only need to flip the normal
1739  if (update_flags & update_normal_vectors)
1740  output_data.normal_vectors[point] *= -1.;
1741  }
1742  else
1743  {
1744  if (update_flags & update_normal_vectors)
1745  {
1746  Assert(spacedim - dim == 1,
1747  ExcMessage(
1748  "There is no cell normal in codim 2."));
1749 
1750  if (dim == 1)
1751  output_data.normal_vectors[point] =
1752  cross_product_2d(-DX_t[0]);
1753  else // dim == 2
1754  output_data.normal_vectors[point] =
1755  cross_product_3d(DX_t[0], DX_t[1]);
1756 
1757  output_data.normal_vectors[point] /=
1758  output_data.normal_vectors[point].norm();
1759 
1760  if (cell->direction_flag() == false)
1761  output_data.normal_vectors[point] *= -1.;
1762  }
1763  }
1764  } // codim>0 case
1765  }
1766  }
1767 
1768  // copy values from InternalData to vector given by reference
1769  if (update_flags & update_jacobians)
1770  {
1771  AssertDimension(output_data.jacobians.size(), n_q_points);
1772  if (cell_similarity != CellSimilarity::translation)
1773  for (unsigned int point = 0; point < n_q_points; ++point)
1774  output_data.jacobians[point] = data.contravariant[point];
1775  }
1776 
1777  // copy values from InternalData to vector given by reference
1778  if (update_flags & update_inverse_jacobians)
1779  {
1780  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1781  if (cell_similarity != CellSimilarity::translation)
1782  for (unsigned int point = 0; point < n_q_points; ++point)
1783  output_data.inverse_jacobians[point] =
1784  data.covariant[point].transpose();
1785  }
1786 
1787  // calculate derivatives of the Jacobians
1788  internal::MappingFEFieldImplementation::
1789  maybe_update_jacobian_grads<dim, spacedim, VectorType, DoFHandlerType>(
1790  cell_similarity,
1792  data,
1793  euler_dof_handler->get_fe(),
1794  fe_mask,
1795  fe_to_real,
1796  output_data.jacobian_grads);
1797 
1798  // calculate derivatives of the Jacobians pushed forward to real cell
1799  // coordinates
1800  internal::MappingFEFieldImplementation::
1801  maybe_update_jacobian_pushed_forward_grads<dim,
1802  spacedim,
1803  VectorType,
1804  DoFHandlerType>(
1805  cell_similarity,
1807  data,
1808  euler_dof_handler->get_fe(),
1809  fe_mask,
1810  fe_to_real,
1811  output_data.jacobian_pushed_forward_grads);
1812 
1813  // calculate hessians of the Jacobians
1814  internal::MappingFEFieldImplementation::maybe_update_jacobian_2nd_derivatives<
1815  dim,
1816  spacedim,
1817  VectorType,
1818  DoFHandlerType>(cell_similarity,
1820  data,
1821  euler_dof_handler->get_fe(),
1822  fe_mask,
1823  fe_to_real,
1824  output_data.jacobian_2nd_derivatives);
1825 
1826  // calculate hessians of the Jacobians pushed forward to real cell coordinates
1827  internal::MappingFEFieldImplementation::
1828  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1829  spacedim,
1830  VectorType,
1831  DoFHandlerType>(
1832  cell_similarity,
1834  data,
1835  euler_dof_handler->get_fe(),
1836  fe_mask,
1837  fe_to_real,
1839 
1840  // calculate gradients of the hessians of the Jacobians
1841  internal::MappingFEFieldImplementation::maybe_update_jacobian_3rd_derivatives<
1842  dim,
1843  spacedim,
1844  VectorType,
1845  DoFHandlerType>(cell_similarity,
1847  data,
1848  euler_dof_handler->get_fe(),
1849  fe_mask,
1850  fe_to_real,
1851  output_data.jacobian_3rd_derivatives);
1852 
1853  // calculate gradients of the hessians of the Jacobians pushed forward to real
1854  // cell coordinates
1855  internal::MappingFEFieldImplementation::
1856  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1857  spacedim,
1858  VectorType,
1859  DoFHandlerType>(
1860  cell_similarity,
1862  data,
1863  euler_dof_handler->get_fe(),
1864  fe_mask,
1865  fe_to_real,
1867 
1868  return updated_cell_similarity;
1869 }
1870 
1871 
1872 
1873 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1874 void
1876  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1877  const unsigned int face_no,
1878  const Quadrature<dim - 1> & quadrature,
1879  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1881  &output_data) const
1882 {
1883  // convert data object to internal data for this class. fails with an
1884  // exception if that is not possible
1885  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1886  ExcInternalError());
1887  const InternalData &data = static_cast<const InternalData &>(internal_data);
1888 
1889  update_internal_dofs(cell, data);
1890 
1891  internal::MappingFEFieldImplementation::
1892  do_fill_fe_face_values<dim, spacedim, VectorType, DoFHandlerType>(
1893  *this,
1894  cell,
1895  face_no,
1898  cell->face_orientation(face_no),
1899  cell->face_flip(face_no),
1900  cell->face_rotation(face_no),
1901  quadrature.size()),
1902  quadrature,
1903  data,
1904  euler_dof_handler->get_fe(),
1905  fe_mask,
1906  fe_to_real,
1907  output_data);
1908 }
1909 
1910 
1911 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
1912 void
1915  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1916  const unsigned int face_no,
1917  const unsigned int subface_no,
1918  const Quadrature<dim - 1> & quadrature,
1919  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1921  &output_data) const
1922 {
1923  // convert data object to internal data for this class. fails with an
1924  // exception if that is not possible
1925  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1926  ExcInternalError());
1927  const InternalData &data = static_cast<const InternalData &>(internal_data);
1928 
1929  update_internal_dofs(cell, data);
1930 
1931  internal::MappingFEFieldImplementation::
1932  do_fill_fe_face_values<dim, spacedim, VectorType, DoFHandlerType>(
1933  *this,
1934  cell,
1935  face_no,
1938  subface_no,
1939  cell->face_orientation(
1940  face_no),
1941  cell->face_flip(face_no),
1942  cell->face_rotation(face_no),
1943  quadrature.size(),
1944  cell->subface_case(face_no)),
1945  quadrature,
1946  data,
1947  euler_dof_handler->get_fe(),
1948  fe_mask,
1949  fe_to_real,
1950  output_data);
1951 }
1952 
1953 
1954 namespace internal
1955 {
1956  namespace MappingFEFieldImplementation
1957  {
1958  namespace
1959  {
1960  template <int dim,
1961  int spacedim,
1962  int rank,
1963  typename VectorType,
1964  typename DoFHandlerType>
1965  void
1966  transform_fields(
1967  const ArrayView<const Tensor<rank, dim>> & input,
1968  const MappingKind mapping_kind,
1969  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1970  const ArrayView<Tensor<rank, spacedim>> & output)
1971  {
1972  AssertDimension(input.size(), output.size());
1973  Assert((dynamic_cast<
1974  const typename ::
1975  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1976  InternalData *>(&mapping_data) != nullptr),
1977  ExcInternalError());
1978  const typename ::MappingFEField<dim,
1979  spacedim,
1980  VectorType,
1981  DoFHandlerType>::InternalData
1982  &data = static_cast<
1983  const typename ::
1984  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
1985  InternalData &>(mapping_data);
1986 
1987  switch (mapping_kind)
1988  {
1989  case mapping_contravariant:
1990  {
1991  Assert(
1992  data.update_each & update_contravariant_transformation,
1994  "update_contravariant_transformation"));
1995 
1996  for (unsigned int i = 0; i < output.size(); ++i)
1997  output[i] =
1998  apply_transformation(data.contravariant[i], input[i]);
1999 
2000  return;
2001  }
2002 
2003  case mapping_piola:
2004  {
2005  Assert(
2006  data.update_each & update_contravariant_transformation,
2008  "update_contravariant_transformation"));
2009  Assert(
2010  data.update_each & update_volume_elements,
2012  "update_volume_elements"));
2013  Assert(rank == 1, ExcMessage("Only for rank 1"));
2014  for (unsigned int i = 0; i < output.size(); ++i)
2015  {
2016  output[i] =
2017  apply_transformation(data.contravariant[i], input[i]);
2018  output[i] /= data.volume_elements[i];
2019  }
2020  return;
2021  }
2022 
2023 
2024  // We still allow this operation as in the
2025  // reference cell Derivatives are Tensor
2026  // rather than DerivativeForm
2027  case mapping_covariant:
2028  {
2029  Assert(
2030  data.update_each & update_contravariant_transformation,
2032  "update_contravariant_transformation"));
2033 
2034  for (unsigned int i = 0; i < output.size(); ++i)
2035  output[i] = apply_transformation(data.covariant[i], input[i]);
2036 
2037  return;
2038  }
2039 
2040  default:
2041  Assert(false, ExcNotImplemented());
2042  }
2043  }
2044 
2045 
2046  template <int dim,
2047  int spacedim,
2048  int rank,
2049  typename VectorType,
2050  typename DoFHandlerType>
2051  void
2052  transform_differential_forms(
2053  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
2054  const MappingKind mapping_kind,
2055  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2056  const ArrayView<Tensor<rank + 1, spacedim>> & output)
2057  {
2058  AssertDimension(input.size(), output.size());
2059  Assert((dynamic_cast<
2060  const typename ::
2061  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2062  InternalData *>(&mapping_data) != nullptr),
2063  ExcInternalError());
2064  const typename ::MappingFEField<dim,
2065  spacedim,
2066  VectorType,
2067  DoFHandlerType>::InternalData
2068  &data = static_cast<
2069  const typename ::
2070  MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::
2071  InternalData &>(mapping_data);
2072 
2073  switch (mapping_kind)
2074  {
2075  case mapping_covariant:
2076  {
2077  Assert(
2078  data.update_each & update_contravariant_transformation,
2080  "update_contravariant_transformation"));
2081 
2082  for (unsigned int i = 0; i < output.size(); ++i)
2083  output[i] = apply_transformation(data.covariant[i], input[i]);
2084 
2085  return;
2086  }
2087  default:
2088  Assert(false, ExcNotImplemented());
2089  }
2090  }
2091  } // namespace
2092  } // namespace MappingFEFieldImplementation
2093 } // namespace internal
2094 
2095 
2096 
2097 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2098 void
2100  const ArrayView<const Tensor<1, dim>> & input,
2101  const MappingKind mapping_kind,
2102  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2103  const ArrayView<Tensor<1, spacedim>> & output) const
2104 {
2105  AssertDimension(input.size(), output.size());
2106 
2107  internal::MappingFEFieldImplementation::
2108  transform_fields<dim, spacedim, 1, VectorType, DoFHandlerType>(input,
2109  mapping_kind,
2110  mapping_data,
2111  output);
2112 }
2113 
2114 
2115 
2116 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2117 void
2119  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2120  const MappingKind mapping_kind,
2121  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2122  const ArrayView<Tensor<2, spacedim>> & output) const
2123 {
2124  AssertDimension(input.size(), output.size());
2125 
2126  internal::MappingFEFieldImplementation::
2127  transform_differential_forms<dim, spacedim, 1, VectorType, DoFHandlerType>(
2128  input, mapping_kind, mapping_data, output);
2129 }
2130 
2131 
2132 
2133 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2134 void
2136  const ArrayView<const Tensor<2, dim>> &input,
2137  const MappingKind,
2138  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2139  const ArrayView<Tensor<2, spacedim>> & output) const
2140 {
2141  (void)input;
2142  (void)output;
2143  (void)mapping_data;
2144  AssertDimension(input.size(), output.size());
2145 
2146  AssertThrow(false, ExcNotImplemented());
2147 }
2148 
2149 
2150 
2151 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2152 void
2154  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2155  const MappingKind mapping_kind,
2156  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2157  const ArrayView<Tensor<3, spacedim>> & output) const
2158 {
2159  AssertDimension(input.size(), output.size());
2160  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2161  ExcInternalError());
2162  const InternalData &data = static_cast<const InternalData &>(mapping_data);
2163 
2164  switch (mapping_kind)
2165  {
2167  {
2170  "update_covariant_transformation"));
2171 
2172  for (unsigned int q = 0; q < output.size(); ++q)
2173  for (unsigned int i = 0; i < spacedim; ++i)
2174  for (unsigned int j = 0; j < spacedim; ++j)
2175  for (unsigned int k = 0; k < spacedim; ++k)
2176  {
2177  output[q][i][j][k] = data.covariant[q][j][0] *
2178  data.covariant[q][k][0] *
2179  input[q][i][0][0];
2180  for (unsigned int J = 0; J < dim; ++J)
2181  {
2182  const unsigned int K0 = (0 == J) ? 1 : 0;
2183  for (unsigned int K = K0; K < dim; ++K)
2184  output[q][i][j][k] += data.covariant[q][j][J] *
2185  data.covariant[q][k][K] *
2186  input[q][i][J][K];
2187  }
2188  }
2189  return;
2190  }
2191 
2192  default:
2193  Assert(false, ExcNotImplemented());
2194  }
2195 }
2196 
2197 
2198 
2199 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2200 void
2202  const ArrayView<const Tensor<3, dim>> &input,
2203  const MappingKind /*mapping_kind*/,
2204  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2205  const ArrayView<Tensor<3, spacedim>> & output) const
2206 {
2207  (void)input;
2208  (void)output;
2209  (void)mapping_data;
2210  AssertDimension(input.size(), output.size());
2211 
2212  AssertThrow(false, ExcNotImplemented());
2213 }
2214 
2215 
2216 
2217 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2221  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2222  const Point<dim> & p) const
2223 {
2224  // Use the get_data function to create an InternalData with data vectors of
2225  // the right size and transformation shape values already computed at point
2226  // p.
2227  const Quadrature<dim> point_quadrature(p);
2228  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2229  get_data(update_quadrature_points | update_jacobians, point_quadrature));
2230  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2231  ExcInternalError());
2232 
2233  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2234 
2235  return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
2236 }
2237 
2238 
2239 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2243 {
2244  Point<spacedim> p_real;
2245 
2246  for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2247  {
2248  unsigned int comp_i =
2249  euler_dof_handler->get_fe().system_to_component_index(i).first;
2250  if (fe_mask[comp_i])
2251  p_real[fe_to_real[comp_i]] +=
2252  data.local_dof_values[i] * data.shape(0, i);
2253  }
2254 
2255  return p_real;
2256 }
2257 
2258 
2259 
2260 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2261 Point<dim>
2264  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2265  const Point<spacedim> & p) const
2266 {
2267  // first a Newton iteration based on the real mapping. It uses the center
2268  // point of the cell as a starting point
2269  Point<dim> initial_p_unit;
2270  try
2271  {
2272  initial_p_unit =
2273  StaticMappingQ1<dim, spacedim>::mapping.transform_real_to_unit_cell(
2274  cell, p);
2275  }
2276  catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2277  {
2278  // mirror the conditions of the code below to determine if we need to
2279  // use an arbitrary starting point or if we just need to rethrow the
2280  // exception
2281  for (unsigned int d = 0; d < dim; ++d)
2282  initial_p_unit[d] = 0.5;
2283  }
2284 
2285  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2286 
2287  // for (unsigned int d=0; d<dim; ++d)
2288  // initial_p_unit[d] = 0.;
2289 
2290  const Quadrature<dim> point_quadrature(initial_p_unit);
2291 
2293  if (spacedim > dim)
2294  update_flags |= update_jacobian_grads;
2295  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2296  get_data(update_flags, point_quadrature));
2297  Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2298  ExcInternalError());
2299 
2300  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2301 
2302  return do_transform_real_to_unit_cell(cell,
2303  p,
2304  initial_p_unit,
2305  dynamic_cast<InternalData &>(*mdata));
2306 }
2307 
2308 
2309 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2310 Point<dim>
2313  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2314  const Point<spacedim> & p,
2315  const Point<dim> & initial_p_unit,
2316  InternalData & mdata) const
2317 {
2318  const unsigned int n_shapes = mdata.shape_values.size();
2319  (void)n_shapes;
2320  Assert(n_shapes != 0, ExcInternalError());
2321  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2322 
2323 
2324  // Newton iteration to solve
2325  // f(x)=p(x)-p=0
2326  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2327  // The start value was set to be the
2328  // linear approximation to the cell
2329  // The shape values and derivatives
2330  // of the mapping at this point are
2331  // previously computed.
2332  // f(x)
2333  Point<dim> p_unit = initial_p_unit;
2334  Point<dim> f;
2335  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
2337  Tensor<1, spacedim> p_minus_F = p - p_real;
2338  const double eps = 1.e-12 * cell->diameter();
2339  const unsigned int newton_iteration_limit = 20;
2340  unsigned int newton_iteration = 0;
2341  while (p_minus_F.norm_square() > eps * eps)
2342  {
2343  // f'(x)
2344  Point<spacedim> DF[dim];
2345  Tensor<2, dim> df;
2346  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2347  {
2348  const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2349  unsigned int comp_k =
2350  euler_dof_handler->get_fe().system_to_component_index(k).first;
2351  if (fe_mask[comp_k])
2352  for (unsigned int j = 0; j < dim; ++j)
2353  DF[j][fe_to_real[comp_k]] +=
2354  mdata.local_dof_values[k] * grad_k[j];
2355  }
2356  for (unsigned int j = 0; j < dim; ++j)
2357  {
2358  f[j] = DF[j] * p_minus_F;
2359  for (unsigned int l = 0; l < dim; ++l)
2360  df[j][l] = -DF[j] * DF[l];
2361  }
2362  // Solve [f'(x)]d=f(x)
2363  const Tensor<1, dim> delta =
2364  invert(df) * static_cast<const Tensor<1, dim> &>(f);
2365  // do a line search
2366  double step_length = 1;
2367  do
2368  {
2369  // update of p_unit. The
2370  // spacedimth component of
2371  // transformed point is simply
2372  // ignored in codimension one
2373  // case. When this component is
2374  // not zero, then we are
2375  // projecting the point to the
2376  // surface or curve identified
2377  // by the cell.
2378  Point<dim> p_unit_trial = p_unit;
2379  for (unsigned int i = 0; i < dim; ++i)
2380  p_unit_trial[i] -= step_length * delta[i];
2381  // shape values and derivatives
2382  // at new p_unit point
2383  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
2384  mdata);
2385  // f(x)
2386  Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
2387  const Tensor<1, spacedim> f_trial = p - p_real_trial;
2388  // see if we are making progress with the current step length
2389  // and if not, reduce it by a factor of two and try again
2390  if (f_trial.norm() < p_minus_F.norm())
2391  {
2392  p_real = p_real_trial;
2393  p_unit = p_unit_trial;
2394  p_minus_F = f_trial;
2395  break;
2396  }
2397  else if (step_length > 0.05)
2398  step_length /= 2;
2399  else
2400  goto failure;
2401  }
2402  while (true);
2403  ++newton_iteration;
2404  if (newton_iteration > newton_iteration_limit)
2405  goto failure;
2406  }
2407  return p_unit;
2408  // if we get to the following label, then we have either run out
2409  // of Newton iterations, or the line search has not converged.
2410  // in either case, we need to give up, so throw an exception that
2411  // can then be caught
2412 failure:
2413  AssertThrow(false,
2415  // ...the compiler wants us to return something, though we can
2416  // of course never get here...
2417  return Point<dim>();
2418 }
2419 
2420 
2421 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2422 unsigned int
2424 {
2425  return euler_dof_handler->get_fe().degree;
2426 }
2427 
2428 
2429 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2432  const
2433 {
2434  return this->fe_mask;
2435 }
2436 
2437 
2438 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2439 std::unique_ptr<Mapping<dim, spacedim>>
2441 {
2442  return std_cxx14::make_unique<
2444 }
2445 
2446 
2447 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
2448 void
2450  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2452  InternalData &data) const
2453 {
2454  Assert(euler_dof_handler != nullptr,
2455  ExcMessage("euler_dof_handler is empty"));
2456 
2457  typename DoFHandlerType::cell_iterator dof_cell(*cell, euler_dof_handler);
2458  Assert(uses_level_dofs || dof_cell->active() == true, ExcInactiveCell());
2459  if (uses_level_dofs)
2460  {
2461  AssertIndexRange(cell->level(), euler_vector.size());
2462  AssertDimension(euler_vector[cell->level()]->size(),
2463  euler_dof_handler->n_dofs(cell->level()));
2464  }
2465  else
2466  AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2467 
2468  if (uses_level_dofs)
2469  dof_cell->get_mg_dof_indices(data.local_dof_indices);
2470  else
2471  dof_cell->get_dof_indices(data.local_dof_indices);
2472 
2473  const VectorType &vector =
2474  uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2475 
2476  for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2477  data.local_dof_values[i] =
2478  internal::ElementAccess<VectorType>::get(vector,
2479  data.local_dof_indices[i]);
2480 }
2481 
2482 // explicit instantiations
2483 #define SPLIT_INSTANTIATIONS_COUNT 2
2484 #ifndef SPLIT_INSTANTIATIONS_INDEX
2485 # define SPLIT_INSTANTIATIONS_INDEX 0
2486 #endif
2487 #include "mapping_fe_field.inst"
2488 
2489 
2490 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
unsigned int max_level() const
Shape function values.
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
static const unsigned int invalid_unsigned_int
Definition: types.h:187
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Contravariant transformation.
const std::vector< Point< dim > > & get_points() const
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
const std::vector< double > & get_weights() const
virtual void compute_shapes_virtual(const std::vector< Point< dim >> &unit_points, typename MappingFEField< dim, spacedim, VectorType, DoFHandlerType >::InternalData &data) const
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
std::vector< unsigned int > fe_to_real
Determinant of the Jacobian.
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1416
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:344
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Point< spacedim > do_transform_unit_to_real_cell(const InternalData &mdata) const
unsigned int min_level() const
unsigned int get_degree() const
FEValues< dim, spacedim > fe_values
std::vector< Tensor< 3, dim > > shape_third_derivatives
std::vector< Tensor< 1, dim > > shape_derivatives
static ::ExceptionBase & ExcMessage(std::string arg1)
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< Tensor< 2, dim > > shape_second_derivatives
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
Threads::Mutex fe_values_mutex
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:302
std::vector< double > volume_elements
unsigned int size() const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
Gradient of volume element.
std::vector< double > local_dof_values
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
unsigned int size() const
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< double > shape_values
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< std::vector< Tensor< 1, spacedim > > > aux
std::vector< Point< spacedim > > quadrature_points
static Point< dim > project_to_unit_cell(const Point< dim > &p)
DEAL_II_CONSTEXPR numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1425
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3110
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
const bool uses_level_dofs
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &)
Normal vectors.
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
SmartPointer< const DoFHandlerType, MappingFEField< dim, spacedim, VectorType, DoFHandlerType > > euler_dof_handler
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
std::vector< SmartPointer< const VectorType, MappingFEField< dim, spacedim, VectorType, DoFHandlerType > > > euler_vector
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
static ::ExceptionBase & ExcNotImplemented()
ComponentMask get_component_mask() const
void update_internal_dofs(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename MappingFEField< dim, spacedim, VectorType, DoFHandlerType >::InternalData &data) const
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:356
virtual bool preserves_vertex_locations() const override
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
static ::ExceptionBase & ExcInactiveCell()
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
Point< dim > do_transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_p_unit, InternalData &mdata) const
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Covariant transformation.
std::vector< Tensor< 1, spacedim > > normal_vectors
virtual std::size_t memory_consumption() const override
InternalData(const FiniteElement< dim, spacedim > &fe, const ComponentMask &mask)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
UpdateFlags update_each
Definition: mapping.h:629
static ::ExceptionBase & ExcInternalError()
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: quadrature.cc:1139