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