Reference documentation for deal.II version 9.2.0
\(\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_manifold.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 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 
16 
25 
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
32 #include <deal.II/fe/mapping_q1.h>
33 
34 #include <deal.II/grid/manifold.h>
35 #include <deal.II/grid/tria.h>
37 
39 
40 #include <algorithm>
41 #include <array>
42 #include <cmath>
43 #include <memory>
44 #include <numeric>
45 
46 
48 
49 template <int dim, int spacedim>
50 std::size_t
52 {
53  return (
66 }
67 
68 
69 template <int dim, int spacedim>
70 void
72  const UpdateFlags update_flags,
73  const Quadrature<dim> &q,
74  const unsigned int n_original_q_points)
75 {
76  // store the flags in the internal data object so we can access them
77  // in fill_fe_*_values()
78  this->update_each = update_flags;
79 
80  // Store the quadrature
81  this->quad = q;
82 
83  // Resize the weights
85 
86  // see if we need the (transformation) shape function values
87  // and/or gradients and resize the necessary arrays
88  if (this->update_each &
91 
93  covariant.resize(n_original_q_points);
94 
96  contravariant.resize(n_original_q_points);
97 
99  volume_elements.resize(n_original_q_points);
100 }
101 
102 
103 template <int dim, int spacedim>
104 void
106  const UpdateFlags update_flags,
107  const Quadrature<dim> &q,
108  const unsigned int n_original_q_points)
109 {
110  initialize(update_flags, q, n_original_q_points);
111 
112  if (dim > 1)
113  {
115  {
116  aux.resize(dim - 1,
117  std::vector<Tensor<1, spacedim>>(n_original_q_points));
118 
119  // Compute tangentials to the unit cell.
120  for (const unsigned int i : GeometryInfo<dim>::face_indices())
121  {
122  unit_tangentials[i].resize(n_original_q_points);
123  std::fill(unit_tangentials[i].begin(),
124  unit_tangentials[i].end(),
126  if (dim > 2)
127  {
129  .resize(n_original_q_points);
130  std::fill(
132  .begin(),
134  .end(),
136  }
137  }
138  }
139  }
140 }
141 
142 
143 
144 template <int dim, int spacedim>
147 {}
148 
149 
150 
151 template <int dim, int spacedim>
152 std::unique_ptr<Mapping<dim, spacedim>>
154 {
155  return std_cxx14::make_unique<MappingManifold<dim, spacedim>>(*this);
156 }
157 
158 
159 
160 template <int dim, int spacedim>
164  const Point<spacedim> &) const
165 {
166  Assert(false, ExcNotImplemented());
167  return Point<dim>();
168 }
169 
170 
171 
172 template <int dim, int spacedim>
176  const Point<dim> & p) const
177 {
178  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
179  std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;
180 
181  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
182  {
183  vertices[v] = cell->vertex(v);
185  }
186  return cell->get_manifold().get_new_point(
187  make_array_view(vertices.begin(), vertices.end()),
188  make_array_view(weights.begin(), weights.end()));
189 }
190 
191 
192 
193 // In the code below, GCC tries to instantiate MappingManifold<3,4> when
194 // seeing which of the overloaded versions of
195 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
196 // error messages and, generally, nothing very good. Avoid this by ensuring
197 // that this class exists, but does not have an inner InternalData
198 // type, thereby ruling out the codim-1 version of the function
199 // below when doing overload resolution.
200 template <>
201 class MappingManifold<3, 4>
202 {};
203 
204 
205 
206 template <int dim, int spacedim>
209  const UpdateFlags in) const
210 {
211  // add flags if the respective quantities are necessary to compute
212  // what we need. note that some flags appear in both the conditions
213  // and in subsequent set operations. this leads to some circular
214  // logic. the only way to treat this is to iterate. since there are
215  // 5 if-clauses in the loop, it will take at most 5 iterations to
216  // converge. do them:
217  UpdateFlags out = in;
218  for (unsigned int i = 0; i < 5; ++i)
219  {
220  // The following is a little incorrect:
221  // If not applied on a face,
222  // update_boundary_forms does not
223  // make sense. On the other hand,
224  // it is necessary on a
225  // face. Currently,
226  // update_boundary_forms is simply
227  // ignored for the interior of a
228  // cell.
230  out |= update_boundary_forms;
231 
236 
237  if (out &
242 
243  // The contravariant transformation used in the Piola
244  // transformation, which requires the determinant of the Jacobi
245  // matrix of the transformation. Because we have no way of
246  // knowing here whether the finite elements wants to use the
247  // contravariant of the Piola transforms, we add the JxW values
248  // to the list of flags to be updated for each cell.
250  out |= update_JxW_values;
251 
252  if (out & update_normal_vectors)
253  out |= update_JxW_values;
254  }
255 
256  // Now throw an exception if we stumble upon something that was not
257  // implemented yet
262 
263  return out;
264 }
265 
266 
267 
268 template <int dim, int spacedim>
269 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
271  const Quadrature<dim> &q) const
272 {
273  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
274  std_cxx14::make_unique<InternalData>();
275  auto &data = dynamic_cast<InternalData &>(*data_ptr);
276  data.initialize(this->requires_update_flags(update_flags), q, q.size());
277 
278  return data_ptr;
279 }
280 
281 
282 
283 template <int dim, int spacedim>
284 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
286  const UpdateFlags update_flags,
287  const Quadrature<dim - 1> &quadrature) const
288 {
289  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
290  std_cxx14::make_unique<InternalData>();
291  auto &data = dynamic_cast<InternalData &>(*data_ptr);
292  data.initialize_face(this->requires_update_flags(update_flags),
294  quadrature.size());
295 
296  return data_ptr;
297 }
298 
299 
300 
301 template <int dim, int spacedim>
302 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
304  const UpdateFlags update_flags,
305  const Quadrature<dim - 1> &quadrature) const
306 {
307  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
308  std_cxx14::make_unique<InternalData>();
309  auto &data = dynamic_cast<InternalData &>(*data_ptr);
310  data.initialize_face(this->requires_update_flags(update_flags),
312  quadrature.size());
313 
314  return data_ptr;
315 }
316 
317 
318 
319 namespace internal
320 {
321  namespace MappingManifoldImplementation
322  {
323  namespace
324  {
331  template <int dim, int spacedim>
332  void
333  maybe_compute_q_points(
334  const typename QProjector<dim>::DataSetDescriptor data_set,
335  const typename ::MappingManifold<dim, spacedim>::InternalData
336  & data,
337  std::vector<Point<spacedim>> &quadrature_points)
338  {
339  const UpdateFlags update_flags = data.update_each;
340 
341  AssertDimension(data.vertices.size(),
343 
344  if (update_flags & update_quadrature_points)
345  {
346  for (unsigned int point = 0; point < quadrature_points.size();
347  ++point)
348  {
349  quadrature_points[point] = data.manifold->get_new_point(
350  make_array_view(data.vertices),
352  data.cell_manifold_quadrature_weights[point + data_set]));
353  }
354  }
355  }
356 
357 
358 
365  template <int dim, int spacedim>
366  void
367  maybe_update_Jacobians(
368  const typename ::QProjector<dim>::DataSetDescriptor data_set,
369  const typename ::MappingManifold<dim, spacedim>::InternalData
370  &data)
371  {
372  const UpdateFlags update_flags = data.update_each;
373 
374  if (update_flags & update_contravariant_transformation)
375  {
376  const unsigned int n_q_points = data.contravariant.size();
377 
378  std::fill(data.contravariant.begin(),
379  data.contravariant.end(),
381 
383  data.vertices.size());
384  for (unsigned int point = 0; point < n_q_points; ++point)
385  {
386  // Start by figuring out how to compute the direction in
387  // the reference space:
388  const Point<dim> &p = data.quad.point(point + data_set);
389 
390  // And get its image on the manifold:
391  const Point<spacedim> P = data.manifold->get_new_point(
392  make_array_view(data.vertices),
394  data.cell_manifold_quadrature_weights[point + data_set]));
395 
396  // To compute the Jacobian, we choose dim points aligned
397  // with the dim reference axes, which are still in the
398  // given cell, and ask for the tangent vector in these
399  // directions. Choosing the points is somewhat arbitrary,
400  // so we try to be smart and we pick points which are
401  // on the opposite quadrant w.r.t. the evaluation
402  // point.
403  for (unsigned int i = 0; i < dim; ++i)
404  {
405  const Point<dim> ei = Point<dim>::unit_vector(i);
406  const double pi = p[i];
407  Assert(pi >= 0 && pi <= 1.0,
409  "Was expecting a quadrature point "
410  "inside the unit reference element."));
411 
412  // In the length L, we store also the direction sign,
413  // which is positive, if the coordinate is < .5,
414  const double L = pi > .5 ? -pi : 1 - pi;
415 
416  const Point<dim> np(p + L * ei);
417 
418  // Get the weights to compute the np point in real space
419  for (const unsigned int j :
421  data.vertex_weights[j] =
423 
424  const Point<spacedim> NP = data.manifold->get_new_point(
425  make_array_view(data.vertices),
426  make_array_view(data.vertex_weights));
427 
428  const Tensor<1, spacedim> T =
429  data.manifold->get_tangent_vector(P, NP);
430 
431  for (unsigned int d = 0; d < spacedim; ++d)
432  data.contravariant[point][d][i] = T[d] / L;
433  }
434  }
435 
436  if (update_flags & update_covariant_transformation)
437  {
438  const unsigned int n_q_points = data.contravariant.size();
439  for (unsigned int point = 0; point < n_q_points; ++point)
440  {
441  data.covariant[point] =
442  (data.contravariant[point]).covariant_form();
443  }
444  }
445 
446  if (update_flags & update_volume_elements)
447  {
448  const unsigned int n_q_points = data.contravariant.size();
449  for (unsigned int point = 0; point < n_q_points; ++point)
450  data.volume_elements[point] =
451  data.contravariant[point].determinant();
452  }
453  }
454  }
455  } // namespace
456  } // namespace MappingManifoldImplementation
457 } // namespace internal
458 
459 
460 
461 template <int dim, int spacedim>
466  const Quadrature<dim> & quadrature,
467  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
469  &output_data) const
470 {
471  // ensure that the following static_cast is really correct:
472  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
473  ExcInternalError());
474  const InternalData &data = static_cast<const InternalData &>(internal_data);
475 
476  const unsigned int n_q_points = quadrature.size();
477 
478  data.store_vertices(cell);
479  data.manifold = &(cell->get_manifold());
480 
481  internal::MappingManifoldImplementation::maybe_compute_q_points<dim,
482  spacedim>(
484  data,
485  output_data.quadrature_points);
486 
487  internal::MappingManifoldImplementation::maybe_update_Jacobians<dim,
488  spacedim>(
490 
491  const UpdateFlags update_flags = data.update_each;
492  const std::vector<double> &weights = quadrature.get_weights();
493 
494  // Multiply quadrature weights by absolute value of Jacobian determinants or
495  // the area element g=sqrt(DX^t DX) in case of codim > 0
496 
497  if (update_flags & (update_normal_vectors | update_JxW_values))
498  {
499  AssertDimension(output_data.JxW_values.size(), n_q_points);
500 
501  Assert(!(update_flags & update_normal_vectors) ||
502  (output_data.normal_vectors.size() == n_q_points),
503  ExcDimensionMismatch(output_data.normal_vectors.size(),
504  n_q_points));
505 
506 
507  for (unsigned int point = 0; point < n_q_points; ++point)
508  {
509  if (dim == spacedim)
510  {
511  const double det = data.contravariant[point].determinant();
512 
513  // check for distorted cells.
514 
515  // TODO: this allows for anisotropies of up to 1e6 in 3D and
516  // 1e12 in 2D. might want to find a finer
517  // (dimension-independent) criterion
518  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
519  cell->diameter() / std::sqrt(double(dim))),
521  cell->center(), det, point)));
522 
523  output_data.JxW_values[point] = weights[point] * det;
524  }
525  // if dim==spacedim, then there is no cell normal to
526  // compute. since this is for FEValues (and not FEFaceValues),
527  // there are also no face normals to compute
528  else // codim>0 case
529  {
530  Tensor<1, spacedim> DX_t[dim];
531  for (unsigned int i = 0; i < spacedim; ++i)
532  for (unsigned int j = 0; j < dim; ++j)
533  DX_t[j][i] = data.contravariant[point][i][j];
534 
535  Tensor<2, dim> G; // First fundamental form
536  for (unsigned int i = 0; i < dim; ++i)
537  for (unsigned int j = 0; j < dim; ++j)
538  G[i][j] = DX_t[i] * DX_t[j];
539 
540  output_data.JxW_values[point] =
541  std::sqrt(determinant(G)) * weights[point];
542 
543  if (update_flags & update_normal_vectors)
544  {
545  Assert(spacedim == dim + 1,
546  ExcMessage(
547  "There is no (unique) cell normal for " +
549  "-dimensional cells in " +
550  Utilities::int_to_string(spacedim) +
551  "-dimensional space. This only works if the "
552  "space dimension is one greater than the "
553  "dimensionality of the mesh cells."));
554 
555  if (dim == 1)
556  output_data.normal_vectors[point] =
557  cross_product_2d(-DX_t[0]);
558  else // dim == 2
559  output_data.normal_vectors[point] =
560  cross_product_3d(DX_t[0], DX_t[1]);
561 
562  output_data.normal_vectors[point] /=
563  output_data.normal_vectors[point].norm();
564 
565  if (cell->direction_flag() == false)
566  output_data.normal_vectors[point] *= -1.;
567  }
568  } // codim>0 case
569  }
570  }
571 
572 
573 
574  // copy values from InternalData to vector given by reference
575  if (update_flags & update_jacobians)
576  {
577  AssertDimension(output_data.jacobians.size(), n_q_points);
578  for (unsigned int point = 0; point < n_q_points; ++point)
579  output_data.jacobians[point] = data.contravariant[point];
580  }
581 
582  // copy values from InternalData to vector given by reference
583  if (update_flags & update_inverse_jacobians)
584  {
585  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
586  for (unsigned int point = 0; point < n_q_points; ++point)
587  output_data.inverse_jacobians[point] =
588  data.covariant[point].transpose();
589  }
590 
592 }
593 
594 
595 
596 namespace internal
597 {
598  namespace MappingManifoldImplementation
599  {
600  namespace
601  {
611  template <int dim, int spacedim>
612  void
613  maybe_compute_face_data(
614  const ::MappingManifold<dim, spacedim> &mapping,
615  const typename ::Triangulation<dim, spacedim>::cell_iterator
616  & cell,
617  const unsigned int face_no,
618  const unsigned int subface_no,
619  const unsigned int n_q_points,
620  const std::vector<double> &weights,
621  const typename ::MappingManifold<dim, spacedim>::InternalData
622  &data,
624  &output_data)
625  {
626  const UpdateFlags update_flags = data.update_each;
627 
628  if (update_flags & update_boundary_forms)
629  {
630  AssertDimension(output_data.boundary_forms.size(), n_q_points);
631  if (update_flags & update_normal_vectors)
632  AssertDimension(output_data.normal_vectors.size(), n_q_points);
633  if (update_flags & update_JxW_values)
634  AssertDimension(output_data.JxW_values.size(), n_q_points);
635 
636  // map the unit tangentials to the real cell. checking for d!=dim-1
637  // eliminates compiler warnings regarding unsigned int expressions <
638  // 0.
639  for (unsigned int d = 0; d != dim - 1; ++d)
640  {
642  data.unit_tangentials.size(),
643  ExcInternalError());
644  Assert(
645  data.aux[d].size() <=
646  data
647  .unit_tangentials[face_no +
649  .size(),
650  ExcInternalError());
651 
652  mapping.transform(
654  data
655  .unit_tangentials[face_no +
658  data,
659  make_array_view(data.aux[d]));
660  }
661 
662  // if dim==spacedim, we can use the unit tangentials to compute the
663  // boundary form by simply taking the cross product
664  if (dim == spacedim)
665  {
666  for (unsigned int i = 0; i < n_q_points; ++i)
667  switch (dim)
668  {
669  case 1:
670  // in 1d, we don't have access to any of the data.aux
671  // fields (because it has only dim-1 components), but we
672  // can still compute the boundary form by simply
673  // looking at the number of the face
674  output_data.boundary_forms[i][0] =
675  (face_no == 0 ? -1 : +1);
676  break;
677  case 2:
678  output_data.boundary_forms[i] =
679  cross_product_2d(data.aux[0][i]);
680  break;
681  case 3:
682  output_data.boundary_forms[i] =
683  cross_product_3d(data.aux[0][i], data.aux[1][i]);
684  break;
685  default:
686  Assert(false, ExcNotImplemented());
687  }
688  }
689  else //(dim < spacedim)
690  {
691  // in the codim-one case, the boundary form results from the
692  // cross product of all the face tangential vectors and the cell
693  // normal vector
694  //
695  // to compute the cell normal, use the same method used in
696  // fill_fe_values for cells above
697  AssertDimension(data.contravariant.size(), n_q_points);
698 
699  for (unsigned int point = 0; point < n_q_points; ++point)
700  {
701  switch (dim)
702  {
703  case 1:
704  {
705  // J is a tangent vector
706  output_data.boundary_forms[point] =
707  data.contravariant[point].transpose()[0];
708  output_data.boundary_forms[point] /=
709  (face_no == 0 ? -1. : +1.) *
710  output_data.boundary_forms[point].norm();
711 
712  break;
713  }
714 
715  case 2:
716  {
718  data.contravariant[point].transpose();
719 
720  Tensor<1, spacedim> cell_normal =
721  cross_product_3d(DX_t[0], DX_t[1]);
722  cell_normal /= cell_normal.norm();
723 
724  // then compute the face normal from the face
725  // tangent and the cell normal:
726  output_data.boundary_forms[point] =
727  cross_product_3d(data.aux[0][point], cell_normal);
728 
729  break;
730  }
731 
732  default:
733  Assert(false, ExcNotImplemented());
734  }
735  }
736  }
737 
738  if (update_flags & (update_normal_vectors | update_JxW_values))
739  for (unsigned int i = 0; i < output_data.boundary_forms.size();
740  ++i)
741  {
742  if (update_flags & update_JxW_values)
743  {
744  output_data.JxW_values[i] =
745  output_data.boundary_forms[i].norm() * weights[i];
746 
747  if (subface_no != numbers::invalid_unsigned_int)
748  {
749  const double area_ratio =
751  cell->subface_case(face_no), subface_no);
752  output_data.JxW_values[i] *= area_ratio;
753  }
754  }
755 
756  if (update_flags & update_normal_vectors)
757  output_data.normal_vectors[i] =
758  Point<spacedim>(output_data.boundary_forms[i] /
759  output_data.boundary_forms[i].norm());
760  }
761 
762  if (update_flags & update_jacobians)
763  for (unsigned int point = 0; point < n_q_points; ++point)
764  output_data.jacobians[point] = data.contravariant[point];
765 
766  if (update_flags & update_inverse_jacobians)
767  for (unsigned int point = 0; point < n_q_points; ++point)
768  output_data.inverse_jacobians[point] =
769  data.covariant[point].transpose();
770  }
771  }
772 
773 
780  template <int dim, int spacedim>
781  void
782  do_fill_fe_face_values(
783  const ::MappingManifold<dim, spacedim> &mapping,
784  const typename ::Triangulation<dim, spacedim>::cell_iterator
785  & cell,
786  const unsigned int face_no,
787  const unsigned int subface_no,
788  const typename QProjector<dim>::DataSetDescriptor data_set,
789  const Quadrature<dim - 1> & quadrature,
790  const typename ::MappingManifold<dim, spacedim>::InternalData
791  &data,
793  &output_data)
794  {
795  data.store_vertices(cell);
796 
797  data.manifold = &cell->face(face_no)->get_manifold();
798 
799  maybe_compute_q_points<dim, spacedim>(data_set,
800  data,
801  output_data.quadrature_points);
802  maybe_update_Jacobians<dim, spacedim>(data_set, data);
803 
804  maybe_compute_face_data(mapping,
805  cell,
806  face_no,
807  subface_no,
808  quadrature.size(),
809  quadrature.get_weights(),
810  data,
811  output_data);
812  }
813 
814  template <int dim, int spacedim, int rank>
815  void
816  transform_fields(
817  const ArrayView<const Tensor<rank, dim>> & input,
818  const MappingKind mapping_kind,
819  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
820  const ArrayView<Tensor<rank, spacedim>> & output)
821  {
822  AssertDimension(input.size(), output.size());
823  Assert((dynamic_cast<const typename ::
824  MappingManifold<dim, spacedim>::InternalData *>(
825  &mapping_data) != nullptr),
826  ExcInternalError());
827  const typename ::MappingManifold<dim, spacedim>::InternalData
828  &data =
829  static_cast<const typename ::MappingManifold<dim, spacedim>::
830  InternalData &>(mapping_data);
831 
832  switch (mapping_kind)
833  {
835  {
836  Assert(
837  data.update_each & update_contravariant_transformation,
839  "update_contravariant_transformation"));
840 
841  for (unsigned int i = 0; i < output.size(); ++i)
842  output[i] =
843  apply_transformation(data.contravariant[i], input[i]);
844 
845  return;
846  }
847 
848  case mapping_piola:
849  {
850  Assert(
851  data.update_each & update_contravariant_transformation,
853  "update_contravariant_transformation"));
854  Assert(
855  data.update_each & update_volume_elements,
857  "update_volume_elements"));
858  Assert(rank == 1, ExcMessage("Only for rank 1"));
859  if (rank != 1)
860  return;
861 
862  for (unsigned int i = 0; i < output.size(); ++i)
863  {
864  output[i] =
865  apply_transformation(data.contravariant[i], input[i]);
866  output[i] /= data.volume_elements[i];
867  }
868  return;
869  }
870  // We still allow this operation as in the
871  // reference cell Derivatives are Tensor
872  // rather than DerivativeForm
873  case mapping_covariant:
874  {
875  Assert(
876  data.update_each & update_contravariant_transformation,
878  "update_covariant_transformation"));
879 
880  for (unsigned int i = 0; i < output.size(); ++i)
881  output[i] = apply_transformation(data.covariant[i], input[i]);
882 
883  return;
884  }
885 
886  default:
887  Assert(false, ExcNotImplemented());
888  }
889  }
890 
891 
892  template <int dim, int spacedim, int rank>
893  void
894  transform_gradients(
895  const ArrayView<const Tensor<rank, dim>> & input,
896  const MappingKind mapping_kind,
897  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
898  const ArrayView<Tensor<rank, spacedim>> & output)
899  {
900  AssertDimension(input.size(), output.size());
901  Assert((dynamic_cast<const typename ::
902  MappingManifold<dim, spacedim>::InternalData *>(
903  &mapping_data) != nullptr),
904  ExcInternalError());
905  const typename ::MappingManifold<dim, spacedim>::InternalData
906  &data =
907  static_cast<const typename ::MappingManifold<dim, spacedim>::
908  InternalData &>(mapping_data);
909 
910  switch (mapping_kind)
911  {
913  {
914  Assert(
915  data.update_each & update_covariant_transformation,
917  "update_covariant_transformation"));
918  Assert(
919  data.update_each & update_contravariant_transformation,
921  "update_contravariant_transformation"));
922  Assert(rank == 2, ExcMessage("Only for rank 2"));
923 
924  for (unsigned int i = 0; i < output.size(); ++i)
925  {
927  apply_transformation(data.contravariant[i],
928  transpose(input[i]));
929  output[i] =
930  apply_transformation(data.covariant[i], A.transpose());
931  }
932 
933  return;
934  }
935 
937  {
938  Assert(
939  data.update_each & update_covariant_transformation,
941  "update_covariant_transformation"));
942  Assert(rank == 2, ExcMessage("Only for rank 2"));
943 
944  for (unsigned int i = 0; i < output.size(); ++i)
945  {
947  apply_transformation(data.covariant[i],
948  transpose(input[i]));
949  output[i] =
950  apply_transformation(data.covariant[i], A.transpose());
951  }
952 
953  return;
954  }
955 
957  {
958  Assert(
959  data.update_each & update_covariant_transformation,
961  "update_covariant_transformation"));
962  Assert(
963  data.update_each & update_contravariant_transformation,
965  "update_contravariant_transformation"));
966  Assert(
967  data.update_each & update_volume_elements,
969  "update_volume_elements"));
970  Assert(rank == 2, ExcMessage("Only for rank 2"));
971 
972  for (unsigned int i = 0; i < output.size(); ++i)
973  {
975  apply_transformation(data.covariant[i], input[i]);
977  apply_transformation(data.contravariant[i],
978  A.transpose());
979 
980  output[i] = transpose(T);
981  output[i] /= data.volume_elements[i];
982  }
983 
984  return;
985  }
986 
987  default:
988  Assert(false, ExcNotImplemented());
989  }
990  }
991 
992 
993 
994  template <int dim, int spacedim>
995  void
996  transform_hessians(
997  const ArrayView<const Tensor<3, dim>> & input,
998  const MappingKind mapping_kind,
999  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1000  const ArrayView<Tensor<3, spacedim>> & output)
1001  {
1002  AssertDimension(input.size(), output.size());
1003  Assert((dynamic_cast<const typename ::
1004  MappingManifold<dim, spacedim>::InternalData *>(
1005  &mapping_data) != nullptr),
1006  ExcInternalError());
1007  const typename ::MappingManifold<dim, spacedim>::InternalData
1008  &data =
1009  static_cast<const typename ::MappingManifold<dim, spacedim>::
1010  InternalData &>(mapping_data);
1011 
1012  switch (mapping_kind)
1013  {
1015  {
1016  Assert(
1017  data.update_each & update_covariant_transformation,
1019  "update_covariant_transformation"));
1020  Assert(
1021  data.update_each & update_contravariant_transformation,
1023  "update_contravariant_transformation"));
1024 
1025  for (unsigned int q = 0; q < output.size(); ++q)
1026  for (unsigned int i = 0; i < spacedim; ++i)
1027  {
1028  double tmp1[dim][dim];
1029  for (unsigned int J = 0; J < dim; ++J)
1030  for (unsigned int K = 0; K < dim; ++K)
1031  {
1032  tmp1[J][K] =
1033  data.contravariant[q][i][0] * input[q][0][J][K];
1034  for (unsigned int I = 1; I < dim; ++I)
1035  tmp1[J][K] +=
1036  data.contravariant[q][i][I] * input[q][I][J][K];
1037  }
1038  for (unsigned int j = 0; j < spacedim; ++j)
1039  {
1040  double tmp2[dim];
1041  for (unsigned int K = 0; K < dim; ++K)
1042  {
1043  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1044  for (unsigned int J = 1; J < dim; ++J)
1045  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1046  }
1047  for (unsigned int k = 0; k < spacedim; ++k)
1048  {
1049  output[q][i][j][k] =
1050  data.covariant[q][k][0] * tmp2[0];
1051  for (unsigned int K = 1; K < dim; ++K)
1052  output[q][i][j][k] +=
1053  data.covariant[q][k][K] * tmp2[K];
1054  }
1055  }
1056  }
1057  return;
1058  }
1059 
1061  {
1062  Assert(
1063  data.update_each & update_covariant_transformation,
1065  "update_covariant_transformation"));
1066 
1067  for (unsigned int q = 0; q < output.size(); ++q)
1068  for (unsigned int i = 0; i < spacedim; ++i)
1069  {
1070  double tmp1[dim][dim];
1071  for (unsigned int J = 0; J < dim; ++J)
1072  for (unsigned int K = 0; K < dim; ++K)
1073  {
1074  tmp1[J][K] =
1075  data.covariant[q][i][0] * input[q][0][J][K];
1076  for (unsigned int I = 1; I < dim; ++I)
1077  tmp1[J][K] +=
1078  data.covariant[q][i][I] * input[q][I][J][K];
1079  }
1080  for (unsigned int j = 0; j < spacedim; ++j)
1081  {
1082  double tmp2[dim];
1083  for (unsigned int K = 0; K < dim; ++K)
1084  {
1085  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1086  for (unsigned int J = 1; J < dim; ++J)
1087  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1088  }
1089  for (unsigned int k = 0; k < spacedim; ++k)
1090  {
1091  output[q][i][j][k] =
1092  data.covariant[q][k][0] * tmp2[0];
1093  for (unsigned int K = 1; K < dim; ++K)
1094  output[q][i][j][k] +=
1095  data.covariant[q][k][K] * tmp2[K];
1096  }
1097  }
1098  }
1099 
1100  return;
1101  }
1102 
1103  case mapping_piola_hessian:
1104  {
1105  Assert(
1106  data.update_each & update_covariant_transformation,
1108  "update_covariant_transformation"));
1109  Assert(
1110  data.update_each & update_contravariant_transformation,
1112  "update_contravariant_transformation"));
1113  Assert(
1114  data.update_each & update_volume_elements,
1116  "update_volume_elements"));
1117 
1118  for (unsigned int q = 0; q < output.size(); ++q)
1119  for (unsigned int i = 0; i < spacedim; ++i)
1120  {
1121  double factor[dim];
1122  for (unsigned int I = 0; I < dim; ++I)
1123  factor[I] =
1124  data.contravariant[q][i][I] / data.volume_elements[q];
1125  double tmp1[dim][dim];
1126  for (unsigned int J = 0; J < dim; ++J)
1127  for (unsigned int K = 0; K < dim; ++K)
1128  {
1129  tmp1[J][K] = factor[0] * input[q][0][J][K];
1130  for (unsigned int I = 1; I < dim; ++I)
1131  tmp1[J][K] += factor[I] * input[q][I][J][K];
1132  }
1133  for (unsigned int j = 0; j < spacedim; ++j)
1134  {
1135  double tmp2[dim];
1136  for (unsigned int K = 0; K < dim; ++K)
1137  {
1138  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1139  for (unsigned int J = 1; J < dim; ++J)
1140  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1141  }
1142  for (unsigned int k = 0; k < spacedim; ++k)
1143  {
1144  output[q][i][j][k] =
1145  data.covariant[q][k][0] * tmp2[0];
1146  for (unsigned int K = 1; K < dim; ++K)
1147  output[q][i][j][k] +=
1148  data.covariant[q][k][K] * tmp2[K];
1149  }
1150  }
1151  }
1152 
1153  return;
1154  }
1155 
1156  default:
1157  Assert(false, ExcNotImplemented());
1158  }
1159  }
1160 
1161 
1162 
1163  template <int dim, int spacedim, int rank>
1164  void
1165  transform_differential_forms(
1166  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
1167  const MappingKind mapping_kind,
1168  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1169  const ArrayView<Tensor<rank + 1, spacedim>> & output)
1170  {
1171  AssertDimension(input.size(), output.size());
1172  Assert((dynamic_cast<const typename ::
1173  MappingManifold<dim, spacedim>::InternalData *>(
1174  &mapping_data) != nullptr),
1175  ExcInternalError());
1176  const typename ::MappingManifold<dim, spacedim>::InternalData
1177  &data =
1178  static_cast<const typename ::MappingManifold<dim, spacedim>::
1179  InternalData &>(mapping_data);
1180 
1181  switch (mapping_kind)
1182  {
1183  case mapping_covariant:
1184  {
1185  Assert(
1186  data.update_each & update_contravariant_transformation,
1188  "update_covariant_transformation"));
1189 
1190  for (unsigned int i = 0; i < output.size(); ++i)
1191  output[i] = apply_transformation(data.covariant[i], input[i]);
1192 
1193  return;
1194  }
1195  default:
1196  Assert(false, ExcNotImplemented());
1197  }
1198  }
1199  } // namespace
1200  } // namespace MappingManifoldImplementation
1201 } // namespace internal
1202 
1203 
1204 
1205 template <int dim, int spacedim>
1206 void
1209  const unsigned int face_no,
1210  const Quadrature<dim - 1> & quadrature,
1211  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1213  &output_data) const
1214 {
1215  // ensure that the following cast is really correct:
1216  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1217  ExcInternalError());
1218  const InternalData &data = static_cast<const InternalData &>(internal_data);
1219 
1220  internal::MappingManifoldImplementation::do_fill_fe_face_values(
1221  *this,
1222  cell,
1223  face_no,
1226  cell->face_orientation(face_no),
1227  cell->face_flip(face_no),
1228  cell->face_rotation(face_no),
1229  quadrature.size()),
1230  quadrature,
1231  data,
1232  output_data);
1233 }
1234 
1235 
1236 
1237 template <int dim, int spacedim>
1238 void
1241  const unsigned int face_no,
1242  const unsigned int subface_no,
1243  const Quadrature<dim - 1> & quadrature,
1244  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1246  &output_data) const
1247 {
1248  // ensure that the following cast is really correct:
1249  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1250  ExcInternalError());
1251  const InternalData &data = static_cast<const InternalData &>(internal_data);
1252 
1253  internal::MappingManifoldImplementation::do_fill_fe_face_values(
1254  *this,
1255  cell,
1256  face_no,
1257  subface_no,
1259  subface_no,
1260  cell->face_orientation(face_no),
1261  cell->face_flip(face_no),
1262  cell->face_rotation(face_no),
1263  quadrature.size(),
1264  cell->subface_case(face_no)),
1265  quadrature,
1266  data,
1267  output_data);
1268 }
1269 
1270 
1271 
1272 template <int dim, int spacedim>
1273 void
1275  const ArrayView<const Tensor<1, dim>> & input,
1276  const MappingKind mapping_kind,
1277  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1278  const ArrayView<Tensor<1, spacedim>> & output) const
1279 {
1280  internal::MappingManifoldImplementation::transform_fields(input,
1281  mapping_kind,
1282  mapping_data,
1283  output);
1284 }
1285 
1286 
1287 
1288 template <int dim, int spacedim>
1289 void
1291  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1292  const MappingKind mapping_kind,
1293  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1294  const ArrayView<Tensor<2, spacedim>> & output) const
1295 {
1296  internal::MappingManifoldImplementation::transform_differential_forms(
1297  input, mapping_kind, mapping_data, output);
1298 }
1299 
1300 
1301 
1302 template <int dim, int spacedim>
1303 void
1305  const ArrayView<const Tensor<2, dim>> & input,
1306  const MappingKind mapping_kind,
1307  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1308  const ArrayView<Tensor<2, spacedim>> & output) const
1309 {
1310  switch (mapping_kind)
1311  {
1312  case mapping_contravariant:
1313  internal::MappingManifoldImplementation::transform_fields(input,
1314  mapping_kind,
1315  mapping_data,
1316  output);
1317  return;
1318 
1322  internal::MappingManifoldImplementation::transform_gradients(
1323  input, mapping_kind, mapping_data, output);
1324  return;
1325  default:
1326  Assert(false, ExcNotImplemented());
1327  }
1328 }
1329 
1330 
1331 
1332 template <int dim, int spacedim>
1333 void
1335  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1336  const MappingKind mapping_kind,
1337  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1338  const ArrayView<Tensor<3, spacedim>> & output) const
1339 {
1340  AssertDimension(input.size(), output.size());
1341  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1342  ExcInternalError());
1343  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1344 
1345  switch (mapping_kind)
1346  {
1348  {
1351  "update_covariant_transformation"));
1352 
1353  for (unsigned int q = 0; q < output.size(); ++q)
1354  for (unsigned int i = 0; i < spacedim; ++i)
1355  for (unsigned int j = 0; j < spacedim; ++j)
1356  {
1357  double tmp[dim];
1358  for (unsigned int K = 0; K < dim; ++K)
1359  {
1360  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1361  for (unsigned int J = 1; J < dim; ++J)
1362  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1363  }
1364  for (unsigned int k = 0; k < spacedim; ++k)
1365  {
1366  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1367  for (unsigned int K = 1; K < dim; ++K)
1368  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1369  }
1370  }
1371  return;
1372  }
1373 
1374  default:
1375  Assert(false, ExcNotImplemented());
1376  }
1377 }
1378 
1379 
1380 
1381 template <int dim, int spacedim>
1382 void
1384  const ArrayView<const Tensor<3, dim>> & input,
1385  const MappingKind mapping_kind,
1386  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1387  const ArrayView<Tensor<3, spacedim>> & output) const
1388 {
1389  switch (mapping_kind)
1390  {
1391  case mapping_piola_hessian:
1394  internal::MappingManifoldImplementation::transform_hessians(
1395  input, mapping_kind, mapping_data, output);
1396  return;
1397  default:
1398  Assert(false, ExcNotImplemented());
1399  }
1400 }
1401 
1402 //--------------------------- Explicit instantiations -----------------------
1403 #include "mapping_manifold.inst"
1404 
1405 
Transformed quadrature weights.
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Tensor< 1, spacedim, Number > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number > &grad_F, const Tensor< 1, dim, Number > &d_x)
static const char L
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
SmartPointer< const Manifold< dim, spacedim > > manifold
Contravariant transformation.
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
const std::vector< double > & get_weights() const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
static Point< dim, Number > unit_vector(const unsigned int i)
Outer normal vector, not normalized.
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:442
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Determinant of the Jacobian.
Transformed quadrature points.
void store_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual std::size_t memory_consumption() const override
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:346
Triangulation< dim, spacedim >::cell_iterator cell
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:2407
std::vector< Point< spacedim > > vertices
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
std::vector< std::vector< Tensor< 1, spacedim > > > aux
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char T
#define Assert(cond, exc)
Definition: exceptions.h:1419
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:302
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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:607
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
DerivativeForm< 1, spacedim, dim, Number > transpose() const
Gradient of volume element.
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
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)
std::vector< double > volume_elements
unsigned int size() const
std::vector< double > vertex_weights
void compute_manifold_quadrature_weights(const Quadrature< dim > &quadrature)
std::vector< Point< spacedim > > quadrature_points
static const char A
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
Normal vectors.
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
static ::ExceptionBase & ExcNotImplemented()
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
numbers::NumberTraits< Number >::real_type norm() const
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
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:2381
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
Covariant transformation.
std::vector< std::vector< double > > cell_manifold_quadrature_weights
std::vector< Tensor< 1, spacedim > > normal_vectors
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
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
MappingManifold()=default
UpdateFlags update_each
Definition: mapping.h:630
static ::ExceptionBase & ExcInternalError()