Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3049-gbbcd74f990 2025-04-08 08:40:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
mapping_cartesian.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
20#include <deal.II/base/tensor.h>
21
23
26
27#include <deal.II/grid/tria.h>
29
31
32#include <algorithm>
33#include <cmath>
34#include <memory>
35
36
38
41 "You are using MappingCartesian, but the incoming cell is not Cartesian.");
42
43
44
50template <typename CellType>
51bool
52is_cartesian(const CellType &cell)
53{
54 if (!cell->reference_cell().is_hyper_cube())
55 return false;
56
57 // The tolerances here are somewhat larger than the square of the machine
58 // epsilon, because we are going to compare the square of distances (to
59 // avoid computing square roots).
60 const double abs_tol = 1e-30;
61 const double rel_tol = 1e-28;
62 const auto bounding_box = cell->bounding_box();
63 const auto &bounding_vertices = bounding_box.get_boundary_points();
64 const auto bb_diagonal_length_squared =
65 bounding_vertices.first.distance_square(bounding_vertices.second);
66
67 for (const unsigned int v : cell->vertex_indices())
68 {
69 // Choose a tolerance that takes into account both that vertices far
70 // away from the origin have only a finite number of digits
71 // that are considered correct (an "absolute tolerance"), as well as that
72 // vertices are supposed to be close to the corresponding vertices of the
73 // bounding box (a tolerance that is "relative" to the size of the cell).
74 //
75 // We need to do it this way because when a vertex is far away from
76 // the origin, computing the difference between two vertices is subject
77 // to cancellation.
78 const double tolerance = std::max(abs_tol * cell->vertex(v).norm_square(),
79 rel_tol * bb_diagonal_length_squared);
80
81 if (cell->vertex(v).distance_square(bounding_box.vertex(v)) > tolerance)
82 return false;
83 }
84
85 return true;
86}
87
88
89
90template <int dim, int spacedim>
92 const Quadrature<dim> &q)
93 : cell_extents(numbers::signaling_nan<Tensor<1, dim>>())
94 , inverse_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
101template <int dim, int spacedim>
102void
104 const UpdateFlags update_flags,
105 const Quadrature<dim> &)
106{
107 // store the flags in the internal data object so we can access them
108 // in fill_fe_*_values(). use the transitive hull of the required
109 // flags
110 this->update_each = update_flags;
111}
112
113
114
115template <int dim, int spacedim>
116std::size_t
124
125
126
127template <int dim, int spacedim>
128bool
133
134
135
136template <int dim, int spacedim>
137bool
139 const ReferenceCell &reference_cell) const
140{
141 Assert(dim == reference_cell.get_dimension(),
142 ExcMessage("The dimension of your mapping (" +
144 ") and the reference cell cell_type (" +
145 Utilities::to_string(reference_cell.get_dimension()) +
146 " ) do not agree."));
147
148 return reference_cell.is_hyper_cube();
149}
150
151
152
153template <int dim, int spacedim>
156 const UpdateFlags in) const
157{
158 // this mapping is pretty simple in that it can basically compute
159 // every piece of information wanted by FEValues without requiring
160 // computing any other quantities. boundary forms are one exception
161 // since they can be computed from the normal vectors without much
162 // further ado
163 UpdateFlags out = in;
164 if (out & update_boundary_forms)
166
167 return out;
168}
169
170
171
172template <int dim, int spacedim>
173std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
175 const Quadrature<dim> &q) const
176{
177 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
178 std::make_unique<InternalData>();
179 data_ptr->reinit(requires_update_flags(update_flags), q);
180
181 return data_ptr;
182}
183
184
185
186template <int dim, int spacedim>
187std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
189 const UpdateFlags update_flags,
190 const hp::QCollection<dim - 1> &quadrature) const
191{
192 AssertDimension(quadrature.size(), 1);
193
194 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
195 std::make_unique<InternalData>(QProjector<dim>::project_to_all_faces(
196 ReferenceCells::get_hypercube<dim>(), quadrature[0]));
197 auto &data = dynamic_cast<InternalData &>(*data_ptr);
198
199 // verify that we have computed the transitive hull of the required
200 // flags and that FEValues has faithfully passed them on to us
201 Assert(update_flags == requires_update_flags(update_flags),
203
204 // store the flags in the internal data object so we can access them
205 // in fill_fe_*_values()
206 data.update_each = update_flags;
207
208 return data_ptr;
209}
210
211
212
213template <int dim, int spacedim>
214std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
216 const UpdateFlags update_flags,
217 const Quadrature<dim - 1> &quadrature) const
218{
219 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
220 std::make_unique<InternalData>(QProjector<dim>::project_to_all_subfaces(
221 ReferenceCells::get_hypercube<dim>(), quadrature));
222 auto &data = dynamic_cast<InternalData &>(*data_ptr);
223
224 // verify that we have computed the transitive hull of the required
225 // flags and that FEValues has faithfully passed them on to us
226 Assert(update_flags == requires_update_flags(update_flags),
228
229 // store the flags in the internal data object so we can access them
230 // in fill_fe_*_values()
231 data.update_each = update_flags;
232
233 return data_ptr;
234}
235
236
237
238template <int dim, int spacedim>
239void
242 const CellSimilarity::Similarity cell_similarity,
243 const InternalData &data) const
244{
245 // Compute start point and sizes along axes. The vertices to be looked at
246 // are 1, 2, 4 compared to the base vertex 0.
247 if (cell_similarity != CellSimilarity::translation)
248 {
249 const Point<dim> start = cell->vertex(0);
250 for (unsigned int d = 0; d < dim; ++d)
251 {
252 const double cell_extent_d = cell->vertex(1 << d)[d] - start[d];
253 data.cell_extents[d] = cell_extent_d;
254 Assert(cell_extent_d != 0.,
255 ExcMessage("Cell does not appear to be Cartesian!"));
256 data.inverse_cell_extents[d] = 1. / cell_extent_d;
257 }
258 }
259}
260
261
262
263namespace
264{
265 template <int dim>
266 void
267 transform_quadrature_points(
268 const Tensor<1, dim> first_vertex,
269 const Tensor<1, dim> cell_extents,
270 const ArrayView<const Point<dim>> &unit_quadrature_points,
271 const typename QProjector<dim>::DataSetDescriptor &offset,
272 std::vector<Point<dim>> &quadrature_points)
273 {
274 for (unsigned int i = 0; i < quadrature_points.size(); ++i)
275 {
276 quadrature_points[i] = first_vertex;
277 for (unsigned int d = 0; d < dim; ++d)
278 quadrature_points[i][d] +=
279 cell_extents[d] * unit_quadrature_points[i + offset][d];
280 }
281 }
282} // namespace
283
284
285
286template <int dim, int spacedim>
287void
290 const InternalData &data,
291 const ArrayView<const Point<dim>> &unit_quadrature_points,
292 std::vector<Point<dim>> &quadrature_points) const
293{
294 if (data.update_each & update_quadrature_points)
295 {
296 const auto offset = QProjector<dim>::DataSetDescriptor::cell();
297
298 transform_quadrature_points(cell->vertex(0),
299 data.cell_extents,
300 unit_quadrature_points,
301 offset,
302 quadrature_points);
303 }
304}
305
306
307
308template <int dim, int spacedim>
309void
312 const unsigned int face_no,
313 const InternalData &data,
314 std::vector<Point<dim>> &quadrature_points) const
315{
317
318 if (data.update_each & update_quadrature_points)
319 {
321 ReferenceCells::get_hypercube<dim>(),
322 face_no,
323 cell->combined_face_orientation(face_no),
324 quadrature_points.size());
325
326
327 transform_quadrature_points(cell->vertex(0),
328 data.cell_extents,
329 make_array_view(data.quadrature_points),
330 offset,
331 quadrature_points);
332 }
333}
334
335
336
337template <int dim, int spacedim>
338void
341 const unsigned int face_no,
342 const unsigned int sub_no,
343 const InternalData &data,
344 std::vector<Point<dim>> &quadrature_points) const
345{
348 if (cell->face(face_no)->has_children())
349 {
350 AssertIndexRange(sub_no, cell->face(face_no)->n_children());
351 }
352
353 if (data.update_each & update_quadrature_points)
354 {
356 ReferenceCells::get_hypercube<dim>(),
357 face_no,
358 sub_no,
359 cell->combined_face_orientation(face_no),
360 quadrature_points.size(),
361 cell->subface_case(face_no));
362
363 transform_quadrature_points(cell->vertex(0),
364 data.cell_extents,
365 make_array_view(data.quadrature_points),
366 offset,
367 quadrature_points);
368 }
369}
370
371
372
373template <int dim, int spacedim>
374void
376 const unsigned int face_no,
377 const InternalData &data,
378 std::vector<Tensor<1, dim>> &normal_vectors) const
379{
380 // compute normal vectors. All normals on a face have the same value.
381 if (data.update_each & update_normal_vectors)
382 {
384 std::fill(normal_vectors.begin(),
385 normal_vectors.end(),
386 ReferenceCells::get_hypercube<dim>()
387 .template face_normal_vector<dim>(face_no));
388 }
389}
390
391
392
393template <int dim, int spacedim>
394void
396 const InternalData &data,
397 const CellSimilarity::Similarity cell_similarity,
399 &output_data) const
400{
401 if (cell_similarity != CellSimilarity::translation)
402 {
403 if (data.update_each & update_jacobian_grads)
404 for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
406
408 for (unsigned int i = 0;
409 i < output_data.jacobian_pushed_forward_grads.size();
410 ++i)
412
413 if (data.update_each & update_jacobian_2nd_derivatives)
414 for (unsigned int i = 0;
415 i < output_data.jacobian_2nd_derivatives.size();
416 ++i)
417 output_data.jacobian_2nd_derivatives[i] =
419
421 for (unsigned int i = 0;
422 i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
423 ++i)
426
427 if (data.update_each & update_jacobian_3rd_derivatives)
428 for (unsigned int i = 0;
429 i < output_data.jacobian_3rd_derivatives.size();
430 ++i)
431 output_data.jacobian_3rd_derivatives[i] =
433
435 for (unsigned int i = 0;
436 i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
437 ++i)
440 }
441}
442
443
444
445template <int dim, int spacedim>
446void
448 const InternalData &data) const
449{
450 if (data.update_each & update_volume_elements)
451 {
452 double volume = data.cell_extents[0];
453 for (unsigned int d = 1; d < dim; ++d)
454 volume *= data.cell_extents[d];
455 data.volume_element = volume;
456 }
457}
458
459
460
461template <int dim, int spacedim>
462void
464 const InternalData &data,
465 const CellSimilarity::Similarity cell_similarity,
467 &output_data) const
468{
469 // "compute" Jacobian at the quadrature points, which are all the
470 // same
471 if (data.update_each & update_jacobians)
472 if (cell_similarity != CellSimilarity::translation)
473 for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
474 {
476 for (unsigned int j = 0; j < dim; ++j)
477 output_data.jacobians[i][j][j] = data.cell_extents[j];
478 }
479}
480
481
482
483template <int dim, int spacedim>
484void
486 const InternalData &data,
487 const CellSimilarity::Similarity cell_similarity,
489 &output_data) const
490{
491 // "compute" inverse Jacobian at the quadrature points, which are
492 // all the same
493 if (data.update_each & update_inverse_jacobians)
494 if (cell_similarity != CellSimilarity::translation)
495 for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
496 {
497 output_data.inverse_jacobians[i] = Tensor<2, dim>();
498 for (unsigned int j = 0; j < dim; ++j)
499 output_data.inverse_jacobians[i][j][j] =
500 data.inverse_cell_extents[j];
501 }
502}
503
504
505
506template <int dim, int spacedim>
510 const CellSimilarity::Similarity cell_similarity,
511 const Quadrature<dim> &quadrature,
512 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
514 &output_data) const
515{
517
518 // convert data object to internal data for this class. fails with
519 // an exception if that is not possible
520 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
522 const InternalData &data = static_cast<const InternalData &>(internal_data);
523
524
525 update_cell_extents(cell, cell_similarity, data);
526
528 data,
529 quadrature.get_points(),
530 output_data.quadrature_points);
531
532 // compute Jacobian determinant. all values are equal and are the
533 // product of the local lengths in each coordinate direction
534 if (data.update_each & (update_JxW_values | update_volume_elements))
535 if (cell_similarity != CellSimilarity::translation)
536 {
537 double J = data.cell_extents[0];
538 for (unsigned int d = 1; d < dim; ++d)
539 J *= data.cell_extents[d];
540 data.volume_element = J;
541 if (data.update_each & update_JxW_values)
542 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
543 output_data.JxW_values[i] = J * quadrature.weight(i);
544 }
545
546
547 maybe_update_jacobians(data, cell_similarity, output_data);
548 maybe_update_jacobian_derivatives(data, cell_similarity, output_data);
549 maybe_update_inverse_jacobians(data, cell_similarity, output_data);
550
551 return cell_similarity;
552}
553
554
555
556template <int dim, int spacedim>
557void
560 const ArrayView<const Point<dim>> &unit_points,
561 const UpdateFlags update_flags,
563 &output_data) const
564{
565 if (update_flags == update_default)
566 return;
567
569
570 Assert(update_flags & update_inverse_jacobians ||
571 update_flags & update_jacobians ||
572 update_flags & update_quadrature_points,
574
575 output_data.initialize(unit_points.size(), update_flags);
576
578 data.update_each = update_flags;
579
581
583 data,
584 unit_points,
585 output_data.quadrature_points);
586
589}
590
591
592
593template <int dim, int spacedim>
594void
597 const unsigned int face_no,
598 const hp::QCollection<dim - 1> &quadrature,
599 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
601 &output_data) const
602{
604 AssertDimension(quadrature.size(), 1);
605
606 // convert data object to internal
607 // data for this class. fails with
608 // an exception if that is not
609 // possible
610 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
612 const InternalData &data = static_cast<const InternalData &>(internal_data);
613
615
617 face_no,
618 data,
619 output_data.quadrature_points);
620
621 maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
622
623 // first compute Jacobian determinant, which is simply the product
624 // of the local lengths since the jacobian is diagonal
625 double J = 1.;
626 for (unsigned int d = 0; d < dim; ++d)
628 J *= data.cell_extents[d];
629
630 if (data.update_each & update_JxW_values)
631 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
632 output_data.JxW_values[i] = J * quadrature[0].weight(i);
633
634 if (data.update_each & update_boundary_forms)
635 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
636 output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
637
642}
643
644
645
646template <int dim, int spacedim>
647void
650 const unsigned int face_no,
651 const unsigned int subface_no,
652 const Quadrature<dim - 1> &quadrature,
653 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
655 &output_data) const
656{
658
659 // convert data object to internal data for this class. fails with
660 // an exception if that is not possible
661 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
663 const InternalData &data = static_cast<const InternalData &>(internal_data);
664
666
668 cell, face_no, subface_no, data, output_data.quadrature_points);
669
670 maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
671
672 // first compute Jacobian determinant, which is simply the product
673 // of the local lengths since the jacobian is diagonal
674 double J = 1.;
675 for (unsigned int d = 0; d < dim; ++d)
677 J *= data.cell_extents[d];
678
679 if (data.update_each & update_JxW_values)
680 {
681 // Here, cell->face(face_no)->n_children() would be the right
682 // choice, but unfortunately the current function is also called
683 // for faces without children (see tests/fe/mapping.cc). Add
684 // following switch to avoid diffs in tests/fe/mapping.OK
685 const unsigned int n_subfaces =
686 cell->face(face_no)->has_children() ?
687 cell->face(face_no)->n_children() :
689 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
690 output_data.JxW_values[i] = J * quadrature.weight(i) / n_subfaces;
691 }
692
693 if (data.update_each & update_boundary_forms)
694 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
695 output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
696
701}
702
703
704
705template <int dim, int spacedim>
706void
710 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
712 &output_data) const
713{
714 AssertDimension(dim, spacedim);
716
717 // Convert data object to internal data for this class. Fails with an
718 // exception if that is not possible.
719 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
721 const InternalData &data = static_cast<const InternalData &>(internal_data);
722
723
725
727 data,
728 quadrature.get_points(),
729 output_data.quadrature_points);
730
731 if (data.update_each & update_normal_vectors)
732 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
733 {
734 // The normals are n = J^{-T} * \hat{n} before normalizing.
735 Tensor<1, dim> normal;
736 const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
737 for (unsigned int d = 0; d < dim; ++d)
738 {
739 normal[d] = ref_space_normal[d] * data.inverse_cell_extents[d];
740 }
741 normal /= normal.norm();
742 output_data.normal_vectors[i] = normal;
743 }
744
745 if (data.update_each & update_JxW_values)
746 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
747 {
748 const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
749
750 // J^{-T} \times \hat{n}
751 Tensor<1, dim> invJTxNormal;
752 double det_jacobian = 1.;
753 for (unsigned int d = 0; d < dim; ++d)
754 {
755 det_jacobian *= data.cell_extents[d];
756 invJTxNormal[d] =
757 ref_space_normal[d] * data.inverse_cell_extents[d];
758 }
759 output_data.JxW_values[i] =
760 det_jacobian * invJTxNormal.norm() * quadrature.weight(i);
761 }
762
767}
768
769
770
771template <int dim, int spacedim>
772void
774 const ArrayView<const Tensor<1, dim>> &input,
775 const MappingKind mapping_kind,
776 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
777 const ArrayView<Tensor<1, spacedim>> &output) const
778{
779 AssertDimension(input.size(), output.size());
780 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
782 const InternalData &data = static_cast<const InternalData &>(mapping_data);
783
784 switch (mapping_kind)
785 {
787 {
790 "update_covariant_transformation"));
791
792 for (unsigned int i = 0; i < output.size(); ++i)
793 for (unsigned int d = 0; d < dim; ++d)
794 output[i][d] = input[i][d] * data.inverse_cell_extents[d];
795 return;
796 }
797
799 {
802 "update_contravariant_transformation"));
803
804 for (unsigned int i = 0; i < output.size(); ++i)
805 for (unsigned int d = 0; d < dim; ++d)
806 output[i][d] = input[i][d] * data.cell_extents[d];
807 return;
808 }
809 case mapping_piola:
810 {
813 "update_contravariant_transformation"));
814 Assert(data.update_each & update_volume_elements,
816 "update_volume_elements"));
817
818 for (unsigned int i = 0; i < output.size(); ++i)
819 for (unsigned int d = 0; d < dim; ++d)
820 output[i][d] =
821 input[i][d] * data.cell_extents[d] / data.volume_element;
822 return;
823 }
824 default:
826 }
827}
828
829
830
831template <int dim, int spacedim>
832void
835 const MappingKind mapping_kind,
836 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
837 const ArrayView<Tensor<2, spacedim>> &output) const
838{
839 AssertDimension(input.size(), output.size());
840 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
842 const InternalData &data = static_cast<const InternalData &>(mapping_data);
843
844 switch (mapping_kind)
845 {
847 {
850 "update_covariant_transformation"));
851
852 for (unsigned int i = 0; i < output.size(); ++i)
853 for (unsigned int d1 = 0; d1 < dim; ++d1)
854 for (unsigned int d2 = 0; d2 < dim; ++d2)
855 output[i][d1][d2] =
856 input[i][d1][d2] * data.inverse_cell_extents[d2];
857 return;
858 }
859
861 {
864 "update_contravariant_transformation"));
865
866 for (unsigned int i = 0; i < output.size(); ++i)
867 for (unsigned int d1 = 0; d1 < dim; ++d1)
868 for (unsigned int d2 = 0; d2 < dim; ++d2)
869 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
870 return;
871 }
872
874 {
877 "update_covariant_transformation"));
878
879 for (unsigned int i = 0; i < output.size(); ++i)
880 for (unsigned int d1 = 0; d1 < dim; ++d1)
881 for (unsigned int d2 = 0; d2 < dim; ++d2)
882 output[i][d1][d2] = input[i][d1][d2] *
883 data.inverse_cell_extents[d2] *
884 data.inverse_cell_extents[d1];
885 return;
886 }
887
889 {
892 "update_contravariant_transformation"));
893
894 for (unsigned int i = 0; i < output.size(); ++i)
895 for (unsigned int d1 = 0; d1 < dim; ++d1)
896 for (unsigned int d2 = 0; d2 < dim; ++d2)
897 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
898 data.inverse_cell_extents[d1];
899 return;
900 }
901
902 case mapping_piola:
903 {
906 "update_contravariant_transformation"));
907 Assert(data.update_each & update_volume_elements,
909 "update_volume_elements"));
910
911 for (unsigned int i = 0; i < output.size(); ++i)
912 for (unsigned int d1 = 0; d1 < dim; ++d1)
913 for (unsigned int d2 = 0; d2 < dim; ++d2)
914 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
915 data.volume_element;
916 return;
917 }
918
920 {
923 "update_contravariant_transformation"));
924 Assert(data.update_each & update_volume_elements,
926 "update_volume_elements"));
927
928 for (unsigned int i = 0; i < output.size(); ++i)
929 for (unsigned int d1 = 0; d1 < dim; ++d1)
930 for (unsigned int d2 = 0; d2 < dim; ++d2)
931 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
932 data.inverse_cell_extents[d1] /
933 data.volume_element;
934 return;
935 }
936
937 default:
939 }
940}
941
942
943
944template <int dim, int spacedim>
945void
947 const ArrayView<const Tensor<2, dim>> &input,
948 const MappingKind mapping_kind,
949 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
950 const ArrayView<Tensor<2, spacedim>> &output) const
951{
952 AssertDimension(input.size(), output.size());
953 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
955 const InternalData &data = static_cast<const InternalData &>(mapping_data);
956
957 switch (mapping_kind)
958 {
960 {
963 "update_covariant_transformation"));
964
965 for (unsigned int i = 0; i < output.size(); ++i)
966 for (unsigned int d1 = 0; d1 < dim; ++d1)
967 for (unsigned int d2 = 0; d2 < dim; ++d2)
968 output[i][d1][d2] =
969 input[i][d1][d2] * data.inverse_cell_extents[d2];
970 return;
971 }
972
974 {
977 "update_contravariant_transformation"));
978
979 for (unsigned int i = 0; i < output.size(); ++i)
980 for (unsigned int d1 = 0; d1 < dim; ++d1)
981 for (unsigned int d2 = 0; d2 < dim; ++d2)
982 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
983 return;
984 }
985
987 {
990 "update_covariant_transformation"));
991
992 for (unsigned int i = 0; i < output.size(); ++i)
993 for (unsigned int d1 = 0; d1 < dim; ++d1)
994 for (unsigned int d2 = 0; d2 < dim; ++d2)
995 output[i][d1][d2] = input[i][d1][d2] *
996 data.inverse_cell_extents[d2] *
997 data.inverse_cell_extents[d1];
998 return;
999 }
1000
1002 {
1005 "update_contravariant_transformation"));
1006
1007 for (unsigned int i = 0; i < output.size(); ++i)
1008 for (unsigned int d1 = 0; d1 < dim; ++d1)
1009 for (unsigned int d2 = 0; d2 < dim; ++d2)
1010 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
1011 data.inverse_cell_extents[d1];
1012 return;
1013 }
1014
1015 case mapping_piola:
1016 {
1019 "update_contravariant_transformation"));
1020 Assert(data.update_each & update_volume_elements,
1022 "update_volume_elements"));
1023
1024 for (unsigned int i = 0; i < output.size(); ++i)
1025 for (unsigned int d1 = 0; d1 < dim; ++d1)
1026 for (unsigned int d2 = 0; d2 < dim; ++d2)
1027 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
1028 data.volume_element;
1029 return;
1030 }
1031
1033 {
1036 "update_contravariant_transformation"));
1037 Assert(data.update_each & update_volume_elements,
1039 "update_volume_elements"));
1040
1041 for (unsigned int i = 0; i < output.size(); ++i)
1042 for (unsigned int d1 = 0; d1 < dim; ++d1)
1043 for (unsigned int d2 = 0; d2 < dim; ++d2)
1044 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
1045 data.inverse_cell_extents[d1] /
1046 data.volume_element;
1047 return;
1048 }
1049
1050 default:
1052 }
1053}
1054
1055
1056
1057template <int dim, int spacedim>
1058void
1060 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1061 const MappingKind mapping_kind,
1062 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1063 const ArrayView<Tensor<3, spacedim>> &output) const
1064{
1065 AssertDimension(input.size(), output.size());
1066 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1068 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1069
1070 switch (mapping_kind)
1071 {
1073 {
1076 "update_covariant_transformation"));
1077
1078 for (unsigned int q = 0; q < output.size(); ++q)
1079 for (unsigned int i = 0; i < spacedim; ++i)
1080 for (unsigned int j = 0; j < spacedim; ++j)
1081 for (unsigned int k = 0; k < spacedim; ++k)
1082 {
1083 output[q][i][j][k] = input[q][i][j][k] *
1084 data.inverse_cell_extents[j] *
1085 data.inverse_cell_extents[k];
1086 }
1087 return;
1088 }
1089 default:
1091 }
1092}
1093
1094
1095
1096template <int dim, int spacedim>
1097void
1099 const ArrayView<const Tensor<3, dim>> &input,
1100 const MappingKind mapping_kind,
1101 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1102 const ArrayView<Tensor<3, spacedim>> &output) const
1103{
1104 AssertDimension(input.size(), output.size());
1105 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1107 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1108
1109 switch (mapping_kind)
1110 {
1112 {
1115 "update_covariant_transformation"));
1118 "update_contravariant_transformation"));
1119
1120 for (unsigned int q = 0; q < output.size(); ++q)
1121 for (unsigned int i = 0; i < spacedim; ++i)
1122 for (unsigned int j = 0; j < spacedim; ++j)
1123 for (unsigned int k = 0; k < spacedim; ++k)
1124 {
1125 output[q][i][j][k] = input[q][i][j][k] *
1126 data.cell_extents[i] *
1127 data.inverse_cell_extents[j] *
1128 data.inverse_cell_extents[k];
1129 }
1130 return;
1131 }
1132
1134 {
1137 "update_covariant_transformation"));
1138
1139 for (unsigned int q = 0; q < output.size(); ++q)
1140 for (unsigned int i = 0; i < spacedim; ++i)
1141 for (unsigned int j = 0; j < spacedim; ++j)
1142 for (unsigned int k = 0; k < spacedim; ++k)
1143 {
1144 output[q][i][j][k] = input[q][i][j][k] *
1145 (data.inverse_cell_extents[i] *
1146 data.inverse_cell_extents[j]) *
1147 data.inverse_cell_extents[k];
1148 }
1149
1150 return;
1151 }
1152
1154 {
1157 "update_covariant_transformation"));
1160 "update_contravariant_transformation"));
1161 Assert(data.update_each & update_volume_elements,
1163 "update_volume_elements"));
1164
1165 for (unsigned int q = 0; q < output.size(); ++q)
1166 for (unsigned int i = 0; i < spacedim; ++i)
1167 for (unsigned int j = 0; j < spacedim; ++j)
1168 for (unsigned int k = 0; k < spacedim; ++k)
1169 {
1170 output[q][i][j][k] =
1171 input[q][i][j][k] *
1172 (data.cell_extents[i] / data.volume_element *
1173 data.inverse_cell_extents[j]) *
1174 data.inverse_cell_extents[k];
1175 }
1176
1177 return;
1178 }
1179
1180 default:
1182 }
1183}
1184
1185
1186
1187template <int dim, int spacedim>
1191 const Point<dim> &p) const
1192{
1194 Assert(dim == spacedim, ExcNotImplemented());
1195
1196 Point<dim> unit = cell->vertex(0);
1197
1198 // Go through vertices with numbers 1, 2, 4
1199 for (unsigned int d = 0; d < dim; ++d)
1200 unit[d] += (cell->vertex(1 << d)[d] - unit[d]) * p[d];
1201
1202 return unit;
1203}
1204
1205
1206
1207template <int dim, int spacedim>
1211 const Point<spacedim> &p) const
1212{
1214 Assert(dim == spacedim, ExcNotImplemented());
1215
1216 const Point<dim> start = cell->vertex(0);
1217 Point<dim> real = p;
1218
1219 // Go through vertices with numbers 1, 2, 4
1220 for (unsigned int d = 0; d < dim; ++d)
1221 real[d] = (real[d] - start[d]) / (cell->vertex(1 << d)[d] - start[d]);
1222
1223 return real;
1224}
1225
1226
1227
1228template <int dim, int spacedim>
1229void
1232 const ArrayView<const Point<spacedim>> &real_points,
1233 const ArrayView<Point<dim>> &unit_points) const
1234{
1236 AssertDimension(real_points.size(), unit_points.size());
1237
1238 if (dim != spacedim)
1240
1241 const Point<dim> start = cell->vertex(0);
1242
1243 // Go through vertices with numbers 1, 2, 4
1244 std::array<double, dim> inverse_lengths;
1245 for (unsigned int d = 0; d < dim; ++d)
1246 inverse_lengths[d] = 1. / (cell->vertex(1 << d)[d] - start[d]);
1247
1248 for (unsigned int i = 0; i < real_points.size(); ++i)
1249 for (unsigned int d = 0; d < dim; ++d)
1250 unit_points[i][d] = (real_points[i][d] - start[d]) * inverse_lengths[d];
1251}
1252
1253
1254
1255template <int dim, int spacedim>
1256std::unique_ptr<Mapping<dim, spacedim>>
1258{
1259 return std::make_unique<MappingCartesian<dim, spacedim>>(*this);
1260}
1261
1262
1263//---------------------------------------------------------------------------
1264// explicit instantiations
1265#include "fe/mapping_cartesian.inst"
1266
1267
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
virtual std::size_t memory_consumption() const override
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_cell_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const InternalData &data, const ArrayView< const Point< dim > > &unit_quadrature_points, std::vector< Point< dim > > &quadrature_points) const
void maybe_update_volume_elements(const InternalData &data) 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
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 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 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 is_compatible_with(const ReferenceCell &reference_cell) const override
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
virtual void transform_points_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< spacedim > > &real_points, const ArrayView< Point< dim > > &unit_points) const override
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
void maybe_update_jacobian_derivatives(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
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 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:320
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition point.h:113
Class storing the offset index into a Quadrature rule created by project_to_all_faces() or project_to...
Definition qprojector.h:340
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 cell()
Definition qprojector.h:522
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)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:69
double weight(const unsigned int i) const
const std::vector< Point< dim > > & get_points() const
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition collection.h:316
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
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< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcCellNotCartesian()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInternalError()
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:81
@ mapping_piola
Definition mapping.h:116
@ mapping_covariant_gradient
Definition mapping.h:102
@ mapping_covariant
Definition mapping.h:91
@ mapping_contravariant
Definition mapping.h:96
@ mapping_contravariant_hessian
Definition mapping.h:158
@ mapping_covariant_hessian
Definition mapping.h:152
@ mapping_contravariant_gradient
Definition mapping.h:108
@ mapping_piola_gradient
Definition mapping.h:122
@ mapping_piola_hessian
Definition mapping.h:164
bool is_cartesian(const CellType &cell)
std::vector< index_type > data
Definition mpi.cc:746
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:475
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)