Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40:02+00:00
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
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>
94  : cell_extents(numbers::signaling_nan<Tensor<1, dim>>())
95  , volume_element(numbers::signaling_nan<double>())
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
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>
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 {
184
185  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
186  std::make_unique<InternalData>(QProjector<dim>::project_to_all_faces(
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,
209 {
210  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
211  std::make_unique<InternalData>(QProjector<dim>::project_to_all_subfaces(
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,
271 {
273  {
274  const auto offset = QProjector<dim>::DataSetDescriptor::cell();
275
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,
289 {
291
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),
301
302
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,
317 {
320  if (cell->face(face_no)->has_children())
321  {
322  AssertIndexRange(sub_no, cell->face(face_no)->n_children());
323  }
324
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),
335  cell->subface_case(face_no));
336
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,
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  {
358  for (unsigned int d = 0; d < dim; ++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  {
396  for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
398
400  for (unsigned int i = 0;
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,
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,
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 ||
564
565  output_data.initialize(unit_points.size(), update_flags);
566
567  InternalData data;
568  data.update_each = update_flags;
570  std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
571
573
575  data,
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 {
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,
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,
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,
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_3rd_derivatives
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
@ update_covariant_transformation
Covariant transformation.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ 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
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
Definition: mapping.h:105
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()