Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mapping_cartesian.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
21 #include <deal.II/base/tensor.h>
22 
24 
25 #include <deal.II/fe/fe_values.h>
27 
28 #include <deal.II/grid/tria.h>
30 
32 
33 #include <algorithm>
34 #include <cmath>
35 #include <memory>
36 
37 
39 
42  "You are using MappingCartesian, but the incoming cell is not Cartesian.");
43 
44 
45 
51 template <typename CellType>
52 bool
53 is_cartesian(const CellType &cell)
54 {
55  if (!cell->reference_cell().is_hyper_cube())
56  return false;
57 
58  // The tolerances here are somewhat larger than the square of the machine
59  // epsilon, because we are going to compare the square of distances (to
60  // avoid computing square roots).
61  const double abs_tol = 1e-30;
62  const double rel_tol = 1e-28;
63  const auto bounding_box = cell->bounding_box();
64  const auto &bounding_vertices = bounding_box.get_boundary_points();
65  const auto bb_diagonal_length_squared =
66  bounding_vertices.first.distance_square(bounding_vertices.second);
67 
68  for (const unsigned int v : cell->vertex_indices())
69  {
70  // Choose a tolerance that takes into account both that vertices far
71  // away from the origin have only a finite number of digits
72  // that are considered correct (an "absolute tolerance"), as well as that
73  // vertices are supposed to be close to the corresponding vertices of the
74  // bounding box (a tolerance that is "relative" to the size of the cell).
75  //
76  // We need to do it this way because when a vertex is far away from
77  // the origin, computing the difference between two vertices is subject
78  // to cancellation.
79  const double tolerance = std::max(abs_tol * cell->vertex(v).norm_square(),
80  rel_tol * bb_diagonal_length_squared);
81 
82  if (cell->vertex(v).distance_square(bounding_box.vertex(v)) > tolerance)
83  return false;
84  }
85 
86  return true;
87 }
88 
89 
90 
91 template <int dim, int spacedim>
93  const Quadrature<dim> &q)
94  : cell_extents(numbers::signaling_nan<Tensor<1, dim>>())
95  , volume_element(numbers::signaling_nan<double>())
96  , quadrature_points(q.get_points())
97 {}
98 
99 
100 
101 template <int dim, int spacedim>
102 std::size_t
104 {
107  MemoryConsumption::memory_consumption(volume_element) +
109 }
110 
111 
112 
113 template <int dim, int spacedim>
114 bool
116 {
117  return true;
118 }
119 
120 
121 
122 template <int dim, int spacedim>
123 bool
125  const ReferenceCell &reference_cell) const
126 {
127  Assert(dim == reference_cell.get_dimension(),
128  ExcMessage("The dimension of your mapping (" +
129  Utilities::to_string(dim) +
130  ") and the reference cell cell_type (" +
131  Utilities::to_string(reference_cell.get_dimension()) +
132  " ) do not agree."));
133 
134  return reference_cell.is_hyper_cube();
135 }
136 
137 
138 
139 template <int dim, int spacedim>
142  const UpdateFlags in) const
143 {
144  // this mapping is pretty simple in that it can basically compute
145  // every piece of information wanted by FEValues without requiring
146  // computing any other quantities. boundary forms are one exception
147  // since they can be computed from the normal vectors without much
148  // further ado
149  UpdateFlags out = in;
150  if (out & update_boundary_forms)
151  out |= update_normal_vectors;
152 
153  return out;
154 }
155 
156 
157 
158 template <int dim, int spacedim>
159 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
161  const Quadrature<dim> &q) const
162 {
163  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
164  std::make_unique<InternalData>(q);
165  auto &data = dynamic_cast<InternalData &>(*data_ptr);
166 
167  // store the flags in the internal data object so we can access them
168  // in fill_fe_*_values(). use the transitive hull of the required
169  // flags
170  data.update_each = requires_update_flags(update_flags);
171 
172  return data_ptr;
173 }
174 
175 
176 
177 template <int dim, int spacedim>
178 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
180  const UpdateFlags update_flags,
181  const hp::QCollection<dim - 1> &quadrature) const
182 {
183  AssertDimension(quadrature.size(), 1);
184 
185  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
186  std::make_unique<InternalData>(QProjector<dim>::project_to_all_faces(
187  ReferenceCells::get_hypercube<dim>(), quadrature[0]));
188  auto &data = dynamic_cast<InternalData &>(*data_ptr);
189 
190  // verify that we have computed the transitive hull of the required
191  // flags and that FEValues has faithfully passed them on to us
192  Assert(update_flags == requires_update_flags(update_flags),
193  ExcInternalError());
194 
195  // store the flags in the internal data object so we can access them
196  // in fill_fe_*_values()
197  data.update_each = update_flags;
198 
199  return data_ptr;
200 }
201 
202 
203 
204 template <int dim, int spacedim>
205 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
207  const UpdateFlags update_flags,
208  const Quadrature<dim - 1> &quadrature) const
209 {
210  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
211  std::make_unique<InternalData>(QProjector<dim>::project_to_all_subfaces(
212  ReferenceCells::get_hypercube<dim>(), quadrature));
213  auto &data = dynamic_cast<InternalData &>(*data_ptr);
214 
215  // verify that we have computed the transitive hull of the required
216  // flags and that FEValues has faithfully passed them on to us
217  Assert(update_flags == requires_update_flags(update_flags),
218  ExcInternalError());
219 
220  // store the flags in the internal data object so we can access them
221  // in fill_fe_*_values()
222  data.update_each = update_flags;
223 
224  return data_ptr;
225 }
226 
227 
228 
229 template <int dim, int spacedim>
230 void
232  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
233  const CellSimilarity::Similarity cell_similarity,
234  const InternalData &data) const
235 {
236  // Compute start point and sizes
237  // along axes. Strange vertex
238  // numbering makes this complicated
239  // again.
240  if (cell_similarity != CellSimilarity::translation)
241  {
242  const Point<dim> &start = cell->vertex(0);
243  switch (dim)
244  {
245  case 1:
246  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
247  break;
248  case 2:
249  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
250  data.cell_extents[1] = cell->vertex(2)(1) - start(1);
251  break;
252  case 3:
253  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
254  data.cell_extents[1] = cell->vertex(2)(1) - start(1);
255  data.cell_extents[2] = cell->vertex(4)(2) - start(2);
256  break;
257  default:
258  Assert(false, ExcNotImplemented());
259  }
260  }
261 }
262 
263 
264 
265 template <int dim, int spacedim>
266 void
268  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
269  const InternalData &data,
270  std::vector<Point<dim>> &quadrature_points) const
271 {
272  if (data.update_each & update_quadrature_points)
273  {
274  const auto offset = QProjector<dim>::DataSetDescriptor::cell();
275 
276  transform_quadrature_points(cell, data, offset, quadrature_points);
277  }
278 }
279 
280 
281 
282 template <int dim, int spacedim>
283 void
285  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
286  const unsigned int face_no,
287  const InternalData &data,
288  std::vector<Point<dim>> &quadrature_points) const
289 {
291 
292  if (data.update_each & update_quadrature_points)
293  {
294  const auto offset = QProjector<dim>::DataSetDescriptor::face(
295  ReferenceCells::get_hypercube<dim>(),
296  face_no,
297  cell->face_orientation(face_no),
298  cell->face_flip(face_no),
299  cell->face_rotation(face_no),
300  quadrature_points.size());
301 
302 
303  transform_quadrature_points(cell, data, offset, quadrature_points);
304  }
305 }
306 
307 
308 
309 template <int dim, int spacedim>
310 void
312  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
313  const unsigned int face_no,
314  const unsigned int sub_no,
315  const InternalData &data,
316  std::vector<Point<dim>> &quadrature_points) const
317 {
320  if (cell->face(face_no)->has_children())
321  {
322  AssertIndexRange(sub_no, cell->face(face_no)->n_children());
323  }
324 
325  if (data.update_each & update_quadrature_points)
326  {
328  ReferenceCells::get_hypercube<dim>(),
329  face_no,
330  sub_no,
331  cell->face_orientation(face_no),
332  cell->face_flip(face_no),
333  cell->face_rotation(face_no),
334  quadrature_points.size(),
335  cell->subface_case(face_no));
336 
337  transform_quadrature_points(cell, data, offset, quadrature_points);
338  }
339 }
340 
341 
342 
343 template <int dim, int spacedim>
344 void
346  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
347  const InternalData &data,
348  const typename QProjector<dim>::DataSetDescriptor &offset,
349  std::vector<Point<dim>> &quadrature_points) const
350 {
351  // let @p{start} be the origin of a local coordinate system. it is chosen
352  // as the (lower) left vertex
353  const Point<dim> &start = cell->vertex(0);
354 
355  for (unsigned int i = 0; i < quadrature_points.size(); ++i)
356  {
357  quadrature_points[i] = start;
358  for (unsigned int d = 0; d < dim; ++d)
359  quadrature_points[i](d) +=
360  data.cell_extents[d] * data.quadrature_points[i + offset](d);
361  }
362 }
363 
364 
365 
366 template <int dim, int spacedim>
367 void
369  const unsigned int face_no,
370  const InternalData &data,
371  std::vector<Tensor<1, dim>> &normal_vectors) const
372 {
373  // compute normal vectors. All normals on a face have the same value.
374  if (data.update_each & update_normal_vectors)
375  {
377  std::fill(normal_vectors.begin(),
378  normal_vectors.end(),
380  }
381 }
382 
383 
384 
385 template <int dim, int spacedim>
386 void
388  const InternalData &data,
389  const CellSimilarity::Similarity cell_similarity,
391  &output_data) const
392 {
393  if (cell_similarity != CellSimilarity::translation)
394  {
395  if (data.update_each & update_jacobian_grads)
396  for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
398 
399  if (data.update_each & update_jacobian_pushed_forward_grads)
400  for (unsigned int i = 0;
401  i < output_data.jacobian_pushed_forward_grads.size();
402  ++i)
404 
405  if (data.update_each & update_jacobian_2nd_derivatives)
406  for (unsigned int i = 0;
407  i < output_data.jacobian_2nd_derivatives.size();
408  ++i)
409  output_data.jacobian_2nd_derivatives[i] =
411 
412  if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
413  for (unsigned int i = 0;
414  i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
415  ++i)
418 
419  if (data.update_each & update_jacobian_3rd_derivatives)
420  for (unsigned int i = 0;
421  i < output_data.jacobian_3rd_derivatives.size();
422  ++i)
423  output_data.jacobian_3rd_derivatives[i] =
425 
426  if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
427  for (unsigned int i = 0;
428  i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
429  ++i)
432  }
433 }
434 
435 
436 
437 template <int dim, int spacedim>
438 void
440  const InternalData &data) const
441 {
442  if (data.update_each & update_volume_elements)
443  {
444  double volume = data.cell_extents[0];
445  for (unsigned int d = 1; d < dim; ++d)
446  volume *= data.cell_extents[d];
447  data.volume_element = volume;
448  }
449 }
450 
451 
452 
453 template <int dim, int spacedim>
454 void
456  const InternalData &data,
457  const CellSimilarity::Similarity cell_similarity,
459  &output_data) const
460 {
461  // "compute" Jacobian at the quadrature points, which are all the
462  // same
463  if (data.update_each & update_jacobians)
464  if (cell_similarity != CellSimilarity::translation)
465  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
466  {
467  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
468  for (unsigned int j = 0; j < dim; ++j)
469  output_data.jacobians[i][j][j] = data.cell_extents[j];
470  }
471 }
472 
473 
474 
475 template <int dim, int spacedim>
476 void
478  const InternalData &data,
479  const CellSimilarity::Similarity cell_similarity,
481  &output_data) const
482 {
483  // "compute" inverse Jacobian at the quadrature points, which are
484  // all the same
485  if (data.update_each & update_inverse_jacobians)
486  if (cell_similarity != CellSimilarity::translation)
487  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
488  {
489  output_data.inverse_jacobians[i] = Tensor<2, dim>();
490  for (unsigned int j = 0; j < dim; ++j)
491  output_data.inverse_jacobians[i][j][j] = 1. / data.cell_extents[j];
492  }
493 }
494 
495 
496 
497 template <int dim, int spacedim>
500  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
501  const CellSimilarity::Similarity cell_similarity,
502  const Quadrature<dim> &quadrature,
503  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
505  &output_data) const
506 {
508 
509  // convert data object to internal data for this class. fails with
510  // an exception if that is not possible
511  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
512  ExcInternalError());
513  const InternalData &data = static_cast<const InternalData &>(internal_data);
514 
515 
516  update_cell_extents(cell, cell_similarity, data);
517 
519  data,
520  output_data.quadrature_points);
521 
522  // compute Jacobian determinant. all values are equal and are the
523  // product of the local lengths in each coordinate direction
524  if (data.update_each & (update_JxW_values | update_volume_elements))
525  if (cell_similarity != CellSimilarity::translation)
526  {
527  double J = data.cell_extents[0];
528  for (unsigned int d = 1; d < dim; ++d)
529  J *= data.cell_extents[d];
530  data.volume_element = J;
531  if (data.update_each & update_JxW_values)
532  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
533  output_data.JxW_values[i] = J * quadrature.weight(i);
534  }
535 
536 
537  maybe_update_jacobians(data, cell_similarity, output_data);
538  maybe_update_jacobian_derivatives(data, cell_similarity, output_data);
539  maybe_update_inverse_jacobians(data, cell_similarity, output_data);
540 
541  return cell_similarity;
542 }
543 
544 
545 
546 template <int dim, int spacedim>
547 void
549  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
550  const ArrayView<const Point<dim>> &unit_points,
551  const UpdateFlags update_flags,
553  &output_data) const
554 {
555  if (update_flags == update_default)
556  return;
557 
559 
560  Assert(update_flags & update_inverse_jacobians ||
561  update_flags & update_jacobians ||
562  update_flags & update_quadrature_points,
564 
565  output_data.initialize(unit_points.size(), update_flags);
566 
567  InternalData data;
568  data.update_each = update_flags;
569  data.quadrature_points =
570  std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
571 
573 
575  data,
576  output_data.quadrature_points);
577 
578  maybe_update_jacobians(data, CellSimilarity::none, output_data);
580 }
581 
582 
583 
584 template <int dim, int spacedim>
585 void
587  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
588  const unsigned int face_no,
589  const hp::QCollection<dim - 1> &quadrature,
590  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
592  &output_data) const
593 {
595  AssertDimension(quadrature.size(), 1);
596 
597  // convert data object to internal
598  // data for this class. fails with
599  // an exception if that is not
600  // possible
601  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
602  ExcInternalError());
603  const InternalData &data = static_cast<const InternalData &>(internal_data);
604 
606 
608  face_no,
609  data,
610  output_data.quadrature_points);
611 
612  maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
613 
614  // first compute Jacobian determinant, which is simply the product
615  // of the local lengths since the jacobian is diagonal
616  double J = 1.;
617  for (unsigned int d = 0; d < dim; ++d)
619  J *= data.cell_extents[d];
620 
621  if (data.update_each & update_JxW_values)
622  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
623  output_data.JxW_values[i] = J * quadrature[0].weight(i);
624 
625  if (data.update_each & update_boundary_forms)
626  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
627  output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
628 
630  maybe_update_jacobians(data, CellSimilarity::none, output_data);
633 }
634 
635 
636 
637 template <int dim, int spacedim>
638 void
640  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
641  const unsigned int face_no,
642  const unsigned int subface_no,
643  const Quadrature<dim - 1> &quadrature,
644  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
646  &output_data) const
647 {
649 
650  // convert data object to internal data for this class. fails with
651  // an exception if that is not possible
652  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
653  ExcInternalError());
654  const InternalData &data = static_cast<const InternalData &>(internal_data);
655 
657 
659  cell, face_no, subface_no, data, output_data.quadrature_points);
660 
661  maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
662 
663  // first compute Jacobian determinant, which is simply the product
664  // of the local lengths since the jacobian is diagonal
665  double J = 1.;
666  for (unsigned int d = 0; d < dim; ++d)
668  J *= data.cell_extents[d];
669 
670  if (data.update_each & update_JxW_values)
671  {
672  // Here, cell->face(face_no)->n_children() would be the right
673  // choice, but unfortunately the current function is also called
674  // for faces without children (see tests/fe/mapping.cc). Add
675  // following switch to avoid diffs in tests/fe/mapping.OK
676  const unsigned int n_subfaces =
677  cell->face(face_no)->has_children() ?
678  cell->face(face_no)->n_children() :
680  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
681  output_data.JxW_values[i] = J * quadrature.weight(i) / n_subfaces;
682  }
683 
684  if (data.update_each & update_boundary_forms)
685  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
686  output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
687 
689  maybe_update_jacobians(data, CellSimilarity::none, output_data);
692 }
693 
694 
695 
696 template <int dim, int spacedim>
697 void
699  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
701  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
703  &output_data) const
704 {
705  AssertDimension(dim, spacedim);
707 
708  // Convert data object to internal data for this class. Fails with an
709  // exception if that is not possible.
710  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
711  ExcInternalError());
712  const InternalData &data = static_cast<const InternalData &>(internal_data);
713 
714 
716 
718  data,
719  output_data.quadrature_points);
720 
721  if (data.update_each & update_normal_vectors)
722  for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
723  {
724  // The normals are n = J^{-T} * \hat{n} before normalizing.
725  Tensor<1, dim> normal;
726  const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
727  for (unsigned int d = 0; d < dim; ++d)
728  {
729  normal[d] = ref_space_normal[d] / data.cell_extents[d];
730  }
731  normal /= normal.norm();
732  output_data.normal_vectors[i] = normal;
733  }
734 
735  if (data.update_each & update_JxW_values)
736  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
737  {
738  const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
739 
740  // J^{-T} \times \hat{n}
741  Tensor<1, dim> invJTxNormal;
742  double det_jacobian = 1.;
743  for (unsigned int d = 0; d < dim; ++d)
744  {
745  det_jacobian *= data.cell_extents[d];
746  invJTxNormal[d] = ref_space_normal[d] / data.cell_extents[d];
747  }
748  output_data.JxW_values[i] =
749  det_jacobian * invJTxNormal.norm() * quadrature.weight(i);
750  }
751 
753  maybe_update_jacobians(data, CellSimilarity::none, output_data);
756 }
757 
758 
759 
760 template <int dim, int spacedim>
761 void
763  const ArrayView<const Tensor<1, dim>> &input,
764  const MappingKind mapping_kind,
765  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
766  const ArrayView<Tensor<1, spacedim>> &output) const
767 {
768  AssertDimension(input.size(), output.size());
769  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
770  ExcInternalError());
771  const InternalData &data = static_cast<const InternalData &>(mapping_data);
772 
773  switch (mapping_kind)
774  {
775  case mapping_covariant:
776  {
777  Assert(data.update_each & update_covariant_transformation,
779  "update_covariant_transformation"));
780 
781  for (unsigned int i = 0; i < output.size(); ++i)
782  for (unsigned int d = 0; d < dim; ++d)
783  output[i][d] = input[i][d] / data.cell_extents[d];
784  return;
785  }
786 
788  {
789  Assert(data.update_each & update_contravariant_transformation,
791  "update_contravariant_transformation"));
792 
793  for (unsigned int i = 0; i < output.size(); ++i)
794  for (unsigned int d = 0; d < dim; ++d)
795  output[i][d] = input[i][d] * data.cell_extents[d];
796  return;
797  }
798  case mapping_piola:
799  {
800  Assert(data.update_each & update_contravariant_transformation,
802  "update_contravariant_transformation"));
803  Assert(data.update_each & update_volume_elements,
805  "update_volume_elements"));
806 
807  for (unsigned int i = 0; i < output.size(); ++i)
808  for (unsigned int d = 0; d < dim; ++d)
809  output[i][d] =
810  input[i][d] * data.cell_extents[d] / data.volume_element;
811  return;
812  }
813  default:
814  Assert(false, ExcNotImplemented());
815  }
816 }
817 
818 
819 
820 template <int dim, int spacedim>
821 void
823  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
824  const MappingKind mapping_kind,
825  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
826  const ArrayView<Tensor<2, spacedim>> &output) const
827 {
828  AssertDimension(input.size(), output.size());
829  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
830  ExcInternalError());
831  const InternalData &data = static_cast<const InternalData &>(mapping_data);
832 
833  switch (mapping_kind)
834  {
835  case mapping_covariant:
836  {
837  Assert(data.update_each & update_covariant_transformation,
839  "update_covariant_transformation"));
840 
841  for (unsigned int i = 0; i < output.size(); ++i)
842  for (unsigned int d1 = 0; d1 < dim; ++d1)
843  for (unsigned int d2 = 0; d2 < dim; ++d2)
844  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
845  return;
846  }
847 
849  {
850  Assert(data.update_each & update_contravariant_transformation,
852  "update_contravariant_transformation"));
853 
854  for (unsigned int i = 0; i < output.size(); ++i)
855  for (unsigned int d1 = 0; d1 < dim; ++d1)
856  for (unsigned int d2 = 0; d2 < dim; ++d2)
857  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
858  return;
859  }
860 
862  {
863  Assert(data.update_each & update_covariant_transformation,
865  "update_covariant_transformation"));
866 
867  for (unsigned int i = 0; i < output.size(); ++i)
868  for (unsigned int d1 = 0; d1 < dim; ++d1)
869  for (unsigned int d2 = 0; d2 < dim; ++d2)
870  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
871  data.cell_extents[d1];
872  return;
873  }
874 
876  {
877  Assert(data.update_each & update_contravariant_transformation,
879  "update_contravariant_transformation"));
880 
881  for (unsigned int i = 0; i < output.size(); ++i)
882  for (unsigned int d1 = 0; d1 < dim; ++d1)
883  for (unsigned int d2 = 0; d2 < dim; ++d2)
884  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
885  data.cell_extents[d1];
886  return;
887  }
888 
889  case mapping_piola:
890  {
891  Assert(data.update_each & update_contravariant_transformation,
893  "update_contravariant_transformation"));
894  Assert(data.update_each & update_volume_elements,
896  "update_volume_elements"));
897 
898  for (unsigned int i = 0; i < output.size(); ++i)
899  for (unsigned int d1 = 0; d1 < dim; ++d1)
900  for (unsigned int d2 = 0; d2 < dim; ++d2)
901  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
902  data.volume_element;
903  return;
904  }
905 
907  {
908  Assert(data.update_each & update_contravariant_transformation,
910  "update_contravariant_transformation"));
911  Assert(data.update_each & update_volume_elements,
913  "update_volume_elements"));
914 
915  for (unsigned int i = 0; i < output.size(); ++i)
916  for (unsigned int d1 = 0; d1 < dim; ++d1)
917  for (unsigned int d2 = 0; d2 < dim; ++d2)
918  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
919  data.cell_extents[d1] / data.volume_element;
920  return;
921  }
922 
923  default:
924  Assert(false, ExcNotImplemented());
925  }
926 }
927 
928 
929 
930 template <int dim, int spacedim>
931 void
933  const ArrayView<const Tensor<2, dim>> &input,
934  const MappingKind mapping_kind,
935  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
936  const ArrayView<Tensor<2, spacedim>> &output) const
937 {
938  AssertDimension(input.size(), output.size());
939  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
940  ExcInternalError());
941  const InternalData &data = static_cast<const InternalData &>(mapping_data);
942 
943  switch (mapping_kind)
944  {
945  case mapping_covariant:
946  {
947  Assert(data.update_each & update_covariant_transformation,
949  "update_covariant_transformation"));
950 
951  for (unsigned int i = 0; i < output.size(); ++i)
952  for (unsigned int d1 = 0; d1 < dim; ++d1)
953  for (unsigned int d2 = 0; d2 < dim; ++d2)
954  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
955  return;
956  }
957 
959  {
960  Assert(data.update_each & update_contravariant_transformation,
962  "update_contravariant_transformation"));
963 
964  for (unsigned int i = 0; i < output.size(); ++i)
965  for (unsigned int d1 = 0; d1 < dim; ++d1)
966  for (unsigned int d2 = 0; d2 < dim; ++d2)
967  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
968  return;
969  }
970 
972  {
973  Assert(data.update_each & update_covariant_transformation,
975  "update_covariant_transformation"));
976 
977  for (unsigned int i = 0; i < output.size(); ++i)
978  for (unsigned int d1 = 0; d1 < dim; ++d1)
979  for (unsigned int d2 = 0; d2 < dim; ++d2)
980  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
981  data.cell_extents[d1];
982  return;
983  }
984 
986  {
987  Assert(data.update_each & update_contravariant_transformation,
989  "update_contravariant_transformation"));
990 
991  for (unsigned int i = 0; i < output.size(); ++i)
992  for (unsigned int d1 = 0; d1 < dim; ++d1)
993  for (unsigned int d2 = 0; d2 < dim; ++d2)
994  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
995  data.cell_extents[d1];
996  return;
997  }
998 
999  case mapping_piola:
1000  {
1001  Assert(data.update_each & update_contravariant_transformation,
1003  "update_contravariant_transformation"));
1004  Assert(data.update_each & update_volume_elements,
1006  "update_volume_elements"));
1007 
1008  for (unsigned int i = 0; i < output.size(); ++i)
1009  for (unsigned int d1 = 0; d1 < dim; ++d1)
1010  for (unsigned int d2 = 0; d2 < dim; ++d2)
1011  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
1012  data.volume_element;
1013  return;
1014  }
1015 
1017  {
1018  Assert(data.update_each & update_contravariant_transformation,
1020  "update_contravariant_transformation"));
1021  Assert(data.update_each & update_volume_elements,
1023  "update_volume_elements"));
1024 
1025  for (unsigned int i = 0; i < output.size(); ++i)
1026  for (unsigned int d1 = 0; d1 < dim; ++d1)
1027  for (unsigned int d2 = 0; d2 < dim; ++d2)
1028  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
1029  data.cell_extents[d1] / data.volume_element;
1030  return;
1031  }
1032 
1033  default:
1034  Assert(false, ExcNotImplemented());
1035  }
1036 }
1037 
1038 
1039 
1040 template <int dim, int spacedim>
1041 void
1043  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1044  const MappingKind mapping_kind,
1045  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1046  const ArrayView<Tensor<3, spacedim>> &output) const
1047 {
1048  AssertDimension(input.size(), output.size());
1049  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1050  ExcInternalError());
1051  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1052 
1053  switch (mapping_kind)
1054  {
1056  {
1057  Assert(data.update_each & update_contravariant_transformation,
1059  "update_covariant_transformation"));
1060 
1061  for (unsigned int q = 0; q < output.size(); ++q)
1062  for (unsigned int i = 0; i < spacedim; ++i)
1063  for (unsigned int j = 0; j < spacedim; ++j)
1064  for (unsigned int k = 0; k < spacedim; ++k)
1065  {
1066  output[q][i][j][k] = input[q][i][j][k] /
1067  data.cell_extents[j] /
1068  data.cell_extents[k];
1069  }
1070  return;
1071  }
1072  default:
1073  Assert(false, ExcNotImplemented());
1074  }
1075 }
1076 
1077 
1078 
1079 template <int dim, int spacedim>
1080 void
1082  const ArrayView<const Tensor<3, dim>> &input,
1083  const MappingKind mapping_kind,
1084  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1085  const ArrayView<Tensor<3, spacedim>> &output) const
1086 {
1087  AssertDimension(input.size(), output.size());
1088  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1089  ExcInternalError());
1090  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1091 
1092  switch (mapping_kind)
1093  {
1095  {
1096  Assert(data.update_each & update_covariant_transformation,
1098  "update_covariant_transformation"));
1099  Assert(data.update_each & update_contravariant_transformation,
1101  "update_contravariant_transformation"));
1102 
1103  for (unsigned int q = 0; q < output.size(); ++q)
1104  for (unsigned int i = 0; i < spacedim; ++i)
1105  for (unsigned int j = 0; j < spacedim; ++j)
1106  for (unsigned int k = 0; k < spacedim; ++k)
1107  {
1108  output[q][i][j][k] =
1109  input[q][i][j][k] * data.cell_extents[i] /
1110  data.cell_extents[j] / data.cell_extents[k];
1111  }
1112  return;
1113  }
1114 
1116  {
1117  Assert(data.update_each & update_covariant_transformation,
1119  "update_covariant_transformation"));
1120 
1121  for (unsigned int q = 0; q < output.size(); ++q)
1122  for (unsigned int i = 0; i < spacedim; ++i)
1123  for (unsigned int j = 0; j < spacedim; ++j)
1124  for (unsigned int k = 0; k < spacedim; ++k)
1125  {
1126  output[q][i][j][k] =
1127  input[q][i][j][k] / data.cell_extents[i] /
1128  data.cell_extents[j] / data.cell_extents[k];
1129  }
1130 
1131  return;
1132  }
1133 
1134  case mapping_piola_hessian:
1135  {
1136  Assert(data.update_each & update_covariant_transformation,
1138  "update_covariant_transformation"));
1139  Assert(data.update_each & update_contravariant_transformation,
1141  "update_contravariant_transformation"));
1142  Assert(data.update_each & update_volume_elements,
1144  "update_volume_elements"));
1145 
1146  for (unsigned int q = 0; q < output.size(); ++q)
1147  for (unsigned int i = 0; i < spacedim; ++i)
1148  for (unsigned int j = 0; j < spacedim; ++j)
1149  for (unsigned int k = 0; k < spacedim; ++k)
1150  {
1151  output[q][i][j][k] =
1152  input[q][i][j][k] * data.cell_extents[i] /
1153  data.volume_element / data.cell_extents[j] /
1154  data.cell_extents[k];
1155  }
1156 
1157  return;
1158  }
1159 
1160  default:
1161  Assert(false, ExcNotImplemented());
1162  }
1163 }
1164 
1165 
1166 
1167 template <int dim, int spacedim>
1170  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1171  const Point<dim> &p) const
1172 {
1174 
1175  Tensor<1, dim> length;
1176  const Point<dim> start = cell->vertex(0);
1177  switch (dim)
1178  {
1179  case 1:
1180  length[0] = cell->vertex(1)(0) - start(0);
1181  break;
1182  case 2:
1183  length[0] = cell->vertex(1)(0) - start(0);
1184  length[1] = cell->vertex(2)(1) - start(1);
1185  break;
1186  case 3:
1187  length[0] = cell->vertex(1)(0) - start(0);
1188  length[1] = cell->vertex(2)(1) - start(1);
1189  length[2] = cell->vertex(4)(2) - start(2);
1190  break;
1191  default:
1192  Assert(false, ExcNotImplemented());
1193  }
1194 
1195  Point<dim> p_real = cell->vertex(0);
1196  for (unsigned int d = 0; d < dim; ++d)
1197  p_real(d) += length[d] * p(d);
1198 
1199  return p_real;
1200 }
1201 
1202 
1203 
1204 template <int dim, int spacedim>
1205 Point<dim>
1207  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1208  const Point<spacedim> &p) const
1209 {
1211 
1212  if (dim != spacedim)
1213  Assert(false, ExcNotImplemented());
1214  const Point<dim> &start = cell->vertex(0);
1215  Point<dim> real = p;
1216  real -= start;
1217 
1218  switch (dim)
1219  {
1220  case 1:
1221  real(0) /= cell->vertex(1)(0) - start(0);
1222  break;
1223  case 2:
1224  real(0) /= cell->vertex(1)(0) - start(0);
1225  real(1) /= cell->vertex(2)(1) - start(1);
1226  break;
1227  case 3:
1228  real(0) /= cell->vertex(1)(0) - start(0);
1229  real(1) /= cell->vertex(2)(1) - start(1);
1230  real(2) /= cell->vertex(4)(2) - start(2);
1231  break;
1232  default:
1233  Assert(false, ExcNotImplemented());
1234  }
1235  return real;
1236 }
1237 
1238 
1239 
1240 template <int dim, int spacedim>
1241 std::unique_ptr<Mapping<dim, spacedim>>
1243 {
1244  return std::make_unique<MappingCartesian<dim, spacedim>>(*this);
1245 }
1246 
1247 
1248 //---------------------------------------------------------------------------
1249 // explicit instantiations
1250 #include "mapping_cartesian.inst"
1251 
1252 
std::vector< Point< dim > > quadrature_points
virtual std::size_t memory_consumption() const override
void maybe_update_cell_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const InternalData &data, std::vector< Point< dim >> &quadrature_points) const
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
void maybe_update_volume_elements(const InternalData &data) const
void maybe_update_subface_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const InternalData &data, std::vector< Point< dim >> &quadrature_points) const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
void update_cell_extents(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const InternalData &data) const
virtual void fill_fe_immersed_surface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const NonMatching::ImmersedSurfaceQuadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
void maybe_update_jacobians(const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
void transform_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const InternalData &data, const typename QProjector< dim >::DataSetDescriptor &offset, std::vector< Point< dim >> &quadrature_points) const
void fill_mapping_data_for_generic_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< dim >> &unit_points, const UpdateFlags update_flags, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
void maybe_update_face_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const InternalData &data, std::vector< Point< dim >> &quadrature_points) const
virtual bool preserves_vertex_locations() const override
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
void maybe_update_jacobian_derivatives(const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
void maybe_update_normal_vectors(const unsigned int face_no, const InternalData &data, std::vector< Tensor< 1, dim >> &normal_vectors) const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
void maybe_update_inverse_jacobians(const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Abstract base class for mapping classes.
Definition: mapping.h:317
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition: point.h:112
static DataSetDescriptor cell()
Definition: qprojector.h:394
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: qprojector.cc:1319
double weight(const unsigned int i) const
Definition: tensor.h:516
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition: collection.h:265
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
std::vector< Tensor< 1, spacedim > > normal_vectors
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< Tensor< 1, spacedim > > boundary_forms
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:495
static ::ExceptionBase & ExcCellNotCartesian()
static ::ExceptionBase & ExcMessage(std::string arg1)
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_volume_elements
Determinant of the Jacobian.
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_jacobian_3rd_derivatives
@ update_jacobian_grads
Gradient of volume element.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_quadrature_points
Transformed quadrature points.
@ update_default
No update.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
MappingKind
Definition: mapping.h:78
@ mapping_piola
Definition: mapping.h:113
@ mapping_covariant_gradient
Definition: mapping.h:99
@ mapping_covariant
Definition: mapping.h:88
@ mapping_contravariant
Definition: mapping.h:93
@ mapping_contravariant_hessian
Definition: mapping.h:155
@ mapping_covariant_hessian
Definition: mapping.h:149
@ mapping_contravariant_gradient
Definition: mapping.h:105
@ mapping_piola_gradient
Definition: mapping.h:119
@ mapping_piola_hessian
Definition: mapping.h:161
bool is_cartesian(const CellType &cell)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
double volume(const Triangulation< dim, spacedim > &tria)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:490
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:480
T signaling_nan()