deal.II version GIT relicensing-1927-g3de9220933 2024-10-03 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\}}\)
Loading...
Searching...
No Matches
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,
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,
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(),
387 }
388}
389
390
391
392template <int dim, int spacedim>
393void
395 const InternalData &data,
396 const CellSimilarity::Similarity cell_similarity,
398 &output_data) const
399{
400 if (cell_similarity != CellSimilarity::translation)
401 {
402 if (data.update_each & update_jacobian_grads)
403 for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
405
406 if (data.update_each & update_jacobian_pushed_forward_grads)
407 for (unsigned int i = 0;
408 i < output_data.jacobian_pushed_forward_grads.size();
409 ++i)
411
412 if (data.update_each & update_jacobian_2nd_derivatives)
413 for (unsigned int i = 0;
414 i < output_data.jacobian_2nd_derivatives.size();
415 ++i)
416 output_data.jacobian_2nd_derivatives[i] =
418
420 for (unsigned int i = 0;
421 i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
422 ++i)
425
426 if (data.update_each & update_jacobian_3rd_derivatives)
427 for (unsigned int i = 0;
428 i < output_data.jacobian_3rd_derivatives.size();
429 ++i)
430 output_data.jacobian_3rd_derivatives[i] =
432
434 for (unsigned int i = 0;
435 i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
436 ++i)
439 }
440}
441
442
443
444template <int dim, int spacedim>
445void
447 const InternalData &data) const
448{
449 if (data.update_each & update_volume_elements)
450 {
451 double volume = data.cell_extents[0];
452 for (unsigned int d = 1; d < dim; ++d)
453 volume *= data.cell_extents[d];
454 data.volume_element = volume;
455 }
456}
457
458
459
460template <int dim, int spacedim>
461void
463 const InternalData &data,
464 const CellSimilarity::Similarity cell_similarity,
466 &output_data) const
467{
468 // "compute" Jacobian at the quadrature points, which are all the
469 // same
470 if (data.update_each & update_jacobians)
471 if (cell_similarity != CellSimilarity::translation)
472 for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
473 {
475 for (unsigned int j = 0; j < dim; ++j)
476 output_data.jacobians[i][j][j] = data.cell_extents[j];
477 }
478}
479
480
481
482template <int dim, int spacedim>
483void
485 const InternalData &data,
486 const CellSimilarity::Similarity cell_similarity,
488 &output_data) const
489{
490 // "compute" inverse Jacobian at the quadrature points, which are
491 // all the same
492 if (data.update_each & update_inverse_jacobians)
493 if (cell_similarity != CellSimilarity::translation)
494 for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
495 {
496 output_data.inverse_jacobians[i] = Tensor<2, dim>();
497 for (unsigned int j = 0; j < dim; ++j)
498 output_data.inverse_jacobians[i][j][j] =
499 data.inverse_cell_extents[j];
500 }
501}
502
503
504
505template <int dim, int spacedim>
509 const CellSimilarity::Similarity cell_similarity,
510 const Quadrature<dim> &quadrature,
511 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
513 &output_data) const
514{
516
517 // convert data object to internal data for this class. fails with
518 // an exception if that is not possible
519 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
521 const InternalData &data = static_cast<const InternalData &>(internal_data);
522
523
524 update_cell_extents(cell, cell_similarity, data);
525
527 data,
528 quadrature.get_points(),
529 output_data.quadrature_points);
530
531 // compute Jacobian determinant. all values are equal and are the
532 // product of the local lengths in each coordinate direction
533 if (data.update_each & (update_JxW_values | update_volume_elements))
534 if (cell_similarity != CellSimilarity::translation)
535 {
536 double J = data.cell_extents[0];
537 for (unsigned int d = 1; d < dim; ++d)
538 J *= data.cell_extents[d];
539 data.volume_element = J;
540 if (data.update_each & update_JxW_values)
541 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
542 output_data.JxW_values[i] = J * quadrature.weight(i);
543 }
544
545
546 maybe_update_jacobians(data, cell_similarity, output_data);
547 maybe_update_jacobian_derivatives(data, cell_similarity, output_data);
548 maybe_update_inverse_jacobians(data, cell_similarity, output_data);
549
550 return cell_similarity;
551}
552
553
554
555template <int dim, int spacedim>
556void
559 const ArrayView<const Point<dim>> &unit_points,
560 const UpdateFlags update_flags,
562 &output_data) const
563{
564 if (update_flags == update_default)
565 return;
566
568
569 Assert(update_flags & update_inverse_jacobians ||
570 update_flags & update_jacobians ||
571 update_flags & update_quadrature_points,
573
574 output_data.initialize(unit_points.size(), update_flags);
575
576 InternalData data;
577 data.update_each = update_flags;
578
580
582 data,
583 unit_points,
584 output_data.quadrature_points);
585
588}
589
590
591
592template <int dim, int spacedim>
593void
596 const unsigned int face_no,
597 const hp::QCollection<dim - 1> &quadrature,
598 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
600 &output_data) const
601{
603 AssertDimension(quadrature.size(), 1);
604
605 // convert data object to internal
606 // data for this class. fails with
607 // an exception if that is not
608 // possible
609 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
611 const InternalData &data = static_cast<const InternalData &>(internal_data);
612
614
616 face_no,
617 data,
618 output_data.quadrature_points);
619
620 maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
621
622 // first compute Jacobian determinant, which is simply the product
623 // of the local lengths since the jacobian is diagonal
624 double J = 1.;
625 for (unsigned int d = 0; d < dim; ++d)
627 J *= data.cell_extents[d];
628
629 if (data.update_each & update_JxW_values)
630 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
631 output_data.JxW_values[i] = J * quadrature[0].weight(i);
632
633 if (data.update_each & update_boundary_forms)
634 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
635 output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
636
641}
642
643
644
645template <int dim, int spacedim>
646void
649 const unsigned int face_no,
650 const unsigned int subface_no,
651 const Quadrature<dim - 1> &quadrature,
652 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
654 &output_data) const
655{
657
658 // convert data object to internal data for this class. fails with
659 // an exception if that is not possible
660 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
662 const InternalData &data = static_cast<const InternalData &>(internal_data);
663
665
667 cell, face_no, subface_no, data, output_data.quadrature_points);
668
669 maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
670
671 // first compute Jacobian determinant, which is simply the product
672 // of the local lengths since the jacobian is diagonal
673 double J = 1.;
674 for (unsigned int d = 0; d < dim; ++d)
676 J *= data.cell_extents[d];
677
678 if (data.update_each & update_JxW_values)
679 {
680 // Here, cell->face(face_no)->n_children() would be the right
681 // choice, but unfortunately the current function is also called
682 // for faces without children (see tests/fe/mapping.cc). Add
683 // following switch to avoid diffs in tests/fe/mapping.OK
684 const unsigned int n_subfaces =
685 cell->face(face_no)->has_children() ?
686 cell->face(face_no)->n_children() :
688 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
689 output_data.JxW_values[i] = J * quadrature.weight(i) / n_subfaces;
690 }
691
692 if (data.update_each & update_boundary_forms)
693 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
694 output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
695
700}
701
702
703
704template <int dim, int spacedim>
705void
709 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
711 &output_data) const
712{
713 AssertDimension(dim, spacedim);
715
716 // Convert data object to internal data for this class. Fails with an
717 // exception if that is not possible.
718 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
720 const InternalData &data = static_cast<const InternalData &>(internal_data);
721
722
724
726 data,
727 quadrature.get_points(),
728 output_data.quadrature_points);
729
730 if (data.update_each & update_normal_vectors)
731 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
732 {
733 // The normals are n = J^{-T} * \hat{n} before normalizing.
734 Tensor<1, dim> normal;
735 const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
736 for (unsigned int d = 0; d < dim; ++d)
737 {
738 normal[d] = ref_space_normal[d] * data.inverse_cell_extents[d];
739 }
740 normal /= normal.norm();
741 output_data.normal_vectors[i] = normal;
742 }
743
744 if (data.update_each & update_JxW_values)
745 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
746 {
747 const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
748
749 // J^{-T} \times \hat{n}
750 Tensor<1, dim> invJTxNormal;
751 double det_jacobian = 1.;
752 for (unsigned int d = 0; d < dim; ++d)
753 {
754 det_jacobian *= data.cell_extents[d];
755 invJTxNormal[d] =
756 ref_space_normal[d] * data.inverse_cell_extents[d];
757 }
758 output_data.JxW_values[i] =
759 det_jacobian * invJTxNormal.norm() * quadrature.weight(i);
760 }
761
766}
767
768
769
770template <int dim, int spacedim>
771void
773 const ArrayView<const Tensor<1, dim>> &input,
774 const MappingKind mapping_kind,
775 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
776 const ArrayView<Tensor<1, spacedim>> &output) const
777{
778 AssertDimension(input.size(), output.size());
779 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
781 const InternalData &data = static_cast<const InternalData &>(mapping_data);
782
783 switch (mapping_kind)
784 {
786 {
787 Assert(data.update_each & update_covariant_transformation,
789 "update_covariant_transformation"));
790
791 for (unsigned int i = 0; i < output.size(); ++i)
792 for (unsigned int d = 0; d < dim; ++d)
793 output[i][d] = input[i][d] * data.inverse_cell_extents[d];
794 return;
795 }
796
798 {
801 "update_contravariant_transformation"));
802
803 for (unsigned int i = 0; i < output.size(); ++i)
804 for (unsigned int d = 0; d < dim; ++d)
805 output[i][d] = input[i][d] * data.cell_extents[d];
806 return;
807 }
808 case mapping_piola:
809 {
812 "update_contravariant_transformation"));
813 Assert(data.update_each & update_volume_elements,
815 "update_volume_elements"));
816
817 for (unsigned int i = 0; i < output.size(); ++i)
818 for (unsigned int d = 0; d < dim; ++d)
819 output[i][d] =
820 input[i][d] * data.cell_extents[d] / data.volume_element;
821 return;
822 }
823 default:
825 }
826}
827
828
829
830template <int dim, int spacedim>
831void
834 const MappingKind mapping_kind,
835 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
836 const ArrayView<Tensor<2, spacedim>> &output) const
837{
838 AssertDimension(input.size(), output.size());
839 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
841 const InternalData &data = static_cast<const InternalData &>(mapping_data);
842
843 switch (mapping_kind)
844 {
846 {
847 Assert(data.update_each & update_covariant_transformation,
849 "update_covariant_transformation"));
850
851 for (unsigned int i = 0; i < output.size(); ++i)
852 for (unsigned int d1 = 0; d1 < dim; ++d1)
853 for (unsigned int d2 = 0; d2 < dim; ++d2)
854 output[i][d1][d2] =
855 input[i][d1][d2] * data.inverse_cell_extents[d2];
856 return;
857 }
858
860 {
863 "update_contravariant_transformation"));
864
865 for (unsigned int i = 0; i < output.size(); ++i)
866 for (unsigned int d1 = 0; d1 < dim; ++d1)
867 for (unsigned int d2 = 0; d2 < dim; ++d2)
868 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
869 return;
870 }
871
873 {
874 Assert(data.update_each & update_covariant_transformation,
876 "update_covariant_transformation"));
877
878 for (unsigned int i = 0; i < output.size(); ++i)
879 for (unsigned int d1 = 0; d1 < dim; ++d1)
880 for (unsigned int d2 = 0; d2 < dim; ++d2)
881 output[i][d1][d2] = input[i][d1][d2] *
882 data.inverse_cell_extents[d2] *
883 data.inverse_cell_extents[d1];
884 return;
885 }
886
888 {
891 "update_contravariant_transformation"));
892
893 for (unsigned int i = 0; i < output.size(); ++i)
894 for (unsigned int d1 = 0; d1 < dim; ++d1)
895 for (unsigned int d2 = 0; d2 < dim; ++d2)
896 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
897 data.inverse_cell_extents[d1];
898 return;
899 }
900
901 case mapping_piola:
902 {
905 "update_contravariant_transformation"));
906 Assert(data.update_each & update_volume_elements,
908 "update_volume_elements"));
909
910 for (unsigned int i = 0; i < output.size(); ++i)
911 for (unsigned int d1 = 0; d1 < dim; ++d1)
912 for (unsigned int d2 = 0; d2 < dim; ++d2)
913 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
914 data.volume_element;
915 return;
916 }
917
919 {
922 "update_contravariant_transformation"));
923 Assert(data.update_each & update_volume_elements,
925 "update_volume_elements"));
926
927 for (unsigned int i = 0; i < output.size(); ++i)
928 for (unsigned int d1 = 0; d1 < dim; ++d1)
929 for (unsigned int d2 = 0; d2 < dim; ++d2)
930 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
931 data.inverse_cell_extents[d1] /
932 data.volume_element;
933 return;
934 }
935
936 default:
938 }
939}
940
941
942
943template <int dim, int spacedim>
944void
946 const ArrayView<const Tensor<2, dim>> &input,
947 const MappingKind mapping_kind,
948 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
949 const ArrayView<Tensor<2, spacedim>> &output) const
950{
951 AssertDimension(input.size(), output.size());
952 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
954 const InternalData &data = static_cast<const InternalData &>(mapping_data);
955
956 switch (mapping_kind)
957 {
959 {
960 Assert(data.update_each & update_covariant_transformation,
962 "update_covariant_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] =
968 input[i][d1][d2] * data.inverse_cell_extents[d2];
969 return;
970 }
971
973 {
976 "update_contravariant_transformation"));
977
978 for (unsigned int i = 0; i < output.size(); ++i)
979 for (unsigned int d1 = 0; d1 < dim; ++d1)
980 for (unsigned int d2 = 0; d2 < dim; ++d2)
981 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
982 return;
983 }
984
986 {
987 Assert(data.update_each & update_covariant_transformation,
989 "update_covariant_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] *
995 data.inverse_cell_extents[d2] *
996 data.inverse_cell_extents[d1];
997 return;
998 }
999
1001 {
1004 "update_contravariant_transformation"));
1005
1006 for (unsigned int i = 0; i < output.size(); ++i)
1007 for (unsigned int d1 = 0; d1 < dim; ++d1)
1008 for (unsigned int d2 = 0; d2 < dim; ++d2)
1009 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
1010 data.inverse_cell_extents[d1];
1011 return;
1012 }
1013
1014 case mapping_piola:
1015 {
1018 "update_contravariant_transformation"));
1019 Assert(data.update_each & update_volume_elements,
1021 "update_volume_elements"));
1022
1023 for (unsigned int i = 0; i < output.size(); ++i)
1024 for (unsigned int d1 = 0; d1 < dim; ++d1)
1025 for (unsigned int d2 = 0; d2 < dim; ++d2)
1026 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
1027 data.volume_element;
1028 return;
1029 }
1030
1032 {
1035 "update_contravariant_transformation"));
1036 Assert(data.update_each & update_volume_elements,
1038 "update_volume_elements"));
1039
1040 for (unsigned int i = 0; i < output.size(); ++i)
1041 for (unsigned int d1 = 0; d1 < dim; ++d1)
1042 for (unsigned int d2 = 0; d2 < dim; ++d2)
1043 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] *
1044 data.inverse_cell_extents[d1] /
1045 data.volume_element;
1046 return;
1047 }
1048
1049 default:
1051 }
1052}
1053
1054
1055
1056template <int dim, int spacedim>
1057void
1059 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1060 const MappingKind mapping_kind,
1061 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1062 const ArrayView<Tensor<3, spacedim>> &output) const
1063{
1064 AssertDimension(input.size(), output.size());
1065 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1067 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1068
1069 switch (mapping_kind)
1070 {
1072 {
1075 "update_covariant_transformation"));
1076
1077 for (unsigned int q = 0; q < output.size(); ++q)
1078 for (unsigned int i = 0; i < spacedim; ++i)
1079 for (unsigned int j = 0; j < spacedim; ++j)
1080 for (unsigned int k = 0; k < spacedim; ++k)
1081 {
1082 output[q][i][j][k] = input[q][i][j][k] *
1083 data.inverse_cell_extents[j] *
1084 data.inverse_cell_extents[k];
1085 }
1086 return;
1087 }
1088 default:
1090 }
1091}
1092
1093
1094
1095template <int dim, int spacedim>
1096void
1098 const ArrayView<const Tensor<3, dim>> &input,
1099 const MappingKind mapping_kind,
1100 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1101 const ArrayView<Tensor<3, spacedim>> &output) const
1102{
1103 AssertDimension(input.size(), output.size());
1104 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1106 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1107
1108 switch (mapping_kind)
1109 {
1111 {
1112 Assert(data.update_each & update_covariant_transformation,
1114 "update_covariant_transformation"));
1117 "update_contravariant_transformation"));
1118
1119 for (unsigned int q = 0; q < output.size(); ++q)
1120 for (unsigned int i = 0; i < spacedim; ++i)
1121 for (unsigned int j = 0; j < spacedim; ++j)
1122 for (unsigned int k = 0; k < spacedim; ++k)
1123 {
1124 output[q][i][j][k] = input[q][i][j][k] *
1125 data.cell_extents[i] *
1126 data.inverse_cell_extents[j] *
1127 data.inverse_cell_extents[k];
1128 }
1129 return;
1130 }
1131
1133 {
1134 Assert(data.update_each & update_covariant_transformation,
1136 "update_covariant_transformation"));
1137
1138 for (unsigned int q = 0; q < output.size(); ++q)
1139 for (unsigned int i = 0; i < spacedim; ++i)
1140 for (unsigned int j = 0; j < spacedim; ++j)
1141 for (unsigned int k = 0; k < spacedim; ++k)
1142 {
1143 output[q][i][j][k] = input[q][i][j][k] *
1144 (data.inverse_cell_extents[i] *
1145 data.inverse_cell_extents[j]) *
1146 data.inverse_cell_extents[k];
1147 }
1148
1149 return;
1150 }
1151
1153 {
1154 Assert(data.update_each & update_covariant_transformation,
1156 "update_covariant_transformation"));
1159 "update_contravariant_transformation"));
1160 Assert(data.update_each & update_volume_elements,
1162 "update_volume_elements"));
1163
1164 for (unsigned int q = 0; q < output.size(); ++q)
1165 for (unsigned int i = 0; i < spacedim; ++i)
1166 for (unsigned int j = 0; j < spacedim; ++j)
1167 for (unsigned int k = 0; k < spacedim; ++k)
1168 {
1169 output[q][i][j][k] =
1170 input[q][i][j][k] *
1171 (data.cell_extents[i] / data.volume_element *
1172 data.inverse_cell_extents[j]) *
1173 data.inverse_cell_extents[k];
1174 }
1175
1176 return;
1177 }
1178
1179 default:
1181 }
1182}
1183
1184
1185
1186template <int dim, int spacedim>
1190 const Point<dim> &p) const
1191{
1193 Assert(dim == spacedim, ExcNotImplemented());
1194
1195 Point<dim> unit = cell->vertex(0);
1196
1197 // Go through vertices with numbers 1, 2, 4
1198 for (unsigned int d = 0; d < dim; ++d)
1199 unit[d] += (cell->vertex(1 << d)[d] - unit[d]) * p[d];
1200
1201 return unit;
1202}
1203
1204
1205
1206template <int dim, int spacedim>
1210 const Point<spacedim> &p) const
1211{
1213 Assert(dim == spacedim, ExcNotImplemented());
1214
1215 const Point<dim> start = cell->vertex(0);
1216 Point<dim> real = p;
1217
1218 // Go through vertices with numbers 1, 2, 4
1219 for (unsigned int d = 0; d < dim; ++d)
1220 real[d] = (real[d] - start[d]) / (cell->vertex(1 << d)[d] - start[d]);
1221
1222 return real;
1223}
1224
1225
1226
1227template <int dim, int spacedim>
1228void
1231 const ArrayView<const Point<spacedim>> &real_points,
1232 const ArrayView<Point<dim>> &unit_points) const
1233{
1235 AssertDimension(real_points.size(), unit_points.size());
1236
1237 if (dim != spacedim)
1239
1240 const Point<dim> start = cell->vertex(0);
1241
1242 // Go through vertices with numbers 1, 2, 4
1243 std::array<double, dim> inverse_lengths;
1244 for (unsigned int d = 0; d < dim; ++d)
1245 inverse_lengths[d] = 1. / (cell->vertex(1 << d)[d] - start[d]);
1246
1247 for (unsigned int i = 0; i < real_points.size(); ++i)
1248 for (unsigned int d = 0; d < dim; ++d)
1249 unit_points[i][d] = (real_points[i][d] - start[d]) * inverse_lengths[d];
1250}
1251
1252
1253
1254template <int dim, int spacedim>
1255std::unique_ptr<Mapping<dim, spacedim>>
1257{
1258 return std::make_unique<MappingCartesian<dim, spacedim>>(*this);
1259}
1260
1261
1262//---------------------------------------------------------------------------
1263// explicit instantiations
1264#include "mapping_cartesian.inst"
1265
1266
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:949
std::vector< Point< dim > > quadrature_points
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:318
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition point.h:111
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:461
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)
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:308
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcCellNotCartesian()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:489
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:79
@ mapping_piola
Definition mapping.h:114
@ mapping_covariant_gradient
Definition mapping.h:100
@ mapping_covariant
Definition mapping.h:89
@ mapping_contravariant
Definition mapping.h:94
@ mapping_contravariant_hessian
Definition mapping.h:156
@ mapping_covariant_hessian
Definition mapping.h:150
@ mapping_contravariant_gradient
Definition mapping.h:106
@ mapping_piola_gradient
Definition mapping.h:120
@ mapping_piola_hessian
Definition mapping.h:162
#define DEAL_II_NOT_IMPLEMENTED()
bool is_cartesian(const CellType &cell)
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:479
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)