Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mapping_cartesian.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2001 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
21#include <deal.II/base/tensor.h>
22
24
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
51template <class CellType>
52bool
53is_cartesian(const CellType &cell)
54{
55 if (!cell->reference_cell().is_hyper_cube())
56 return false;
57
58 const double tolerance = 1e-14 * cell->diameter();
59 const auto bounding_box = cell->bounding_box();
60
61 for (const unsigned int v : cell->vertex_indices())
62 if (cell->vertex(v).distance(bounding_box.vertex(v)) > tolerance)
63 return false;
64
65 return true;
66}
67
68
69
70template <int dim, int spacedim>
72 const Quadrature<dim> &q)
73 : cell_extents(numbers::signaling_nan<Tensor<1, dim>>())
74 , volume_element(numbers::signaling_nan<double>())
75 , quadrature_points(q.get_points())
76{}
77
78
79
80template <int dim, int spacedim>
81std::size_t
83{
88}
89
90
91
92template <int dim, int spacedim>
93bool
95{
96 return true;
97}
98
99
100
101template <int dim, int spacedim>
102bool
104 const ReferenceCell &reference_cell) const
105{
106 Assert(dim == reference_cell.get_dimension(),
107 ExcMessage("The dimension of your mapping (" +
109 ") and the reference cell cell_type (" +
110 Utilities::to_string(reference_cell.get_dimension()) +
111 " ) do not agree."));
112
113 return reference_cell.is_hyper_cube();
114}
115
116
117
118template <int dim, int spacedim>
121 const UpdateFlags in) const
122{
123 // this mapping is pretty simple in that it can basically compute
124 // every piece of information wanted by FEValues without requiring
125 // computing any other quantities. boundary forms are one exception
126 // since they can be computed from the normal vectors without much
127 // further ado
128 UpdateFlags out = in;
129 if (out & update_boundary_forms)
131
132 return out;
133}
134
135
136
137template <int dim, int spacedim>
138std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
140 const Quadrature<dim> &q) const
141{
142 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
143 std::make_unique<InternalData>(q);
144 auto &data = dynamic_cast<InternalData &>(*data_ptr);
145
146 // store the flags in the internal data object so we can access them
147 // in fill_fe_*_values(). use the transitive hull of the required
148 // flags
149 data.update_each = requires_update_flags(update_flags);
150
151 return data_ptr;
152}
153
154
155
156template <int dim, int spacedim>
157std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
159 const UpdateFlags update_flags,
160 const hp::QCollection<dim - 1> &quadrature) const
161{
162 AssertDimension(quadrature.size(), 1);
163
164 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
165 std::make_unique<InternalData>(QProjector<dim>::project_to_all_faces(
166 ReferenceCells::get_hypercube<dim>(), quadrature[0]));
167 auto &data = dynamic_cast<InternalData &>(*data_ptr);
168
169 // verify that we have computed the transitive hull of the required
170 // flags and that FEValues has faithfully passed them on to us
171 Assert(update_flags == requires_update_flags(update_flags),
173
174 // store the flags in the internal data object so we can access them
175 // in fill_fe_*_values()
176 data.update_each = update_flags;
177
178 return data_ptr;
179}
180
181
182
183template <int dim, int spacedim>
184std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
186 const UpdateFlags update_flags,
187 const Quadrature<dim - 1> &quadrature) const
188{
189 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
190 std::make_unique<InternalData>(QProjector<dim>::project_to_all_subfaces(
191 ReferenceCells::get_hypercube<dim>(), quadrature));
192 auto &data = dynamic_cast<InternalData &>(*data_ptr);
193
194 // verify that we have computed the transitive hull of the required
195 // flags and that FEValues has faithfully passed them on to us
196 Assert(update_flags == requires_update_flags(update_flags),
198
199 // store the flags in the internal data object so we can access them
200 // in fill_fe_*_values()
201 data.update_each = update_flags;
202
203 return data_ptr;
204}
205
206
207
208template <int dim, int spacedim>
209void
212 const CellSimilarity::Similarity cell_similarity,
213 const InternalData & data) const
214{
215 // Compute start point and sizes
216 // along axes. Strange vertex
217 // numbering makes this complicated
218 // again.
219 if (cell_similarity != CellSimilarity::translation)
220 {
221 const Point<dim> &start = cell->vertex(0);
222 switch (dim)
223 {
224 case 1:
225 data.cell_extents[0] = cell->vertex(1)(0) - start(0);
226 break;
227 case 2:
228 data.cell_extents[0] = cell->vertex(1)(0) - start(0);
229 data.cell_extents[1] = cell->vertex(2)(1) - start(1);
230 break;
231 case 3:
232 data.cell_extents[0] = cell->vertex(1)(0) - start(0);
233 data.cell_extents[1] = cell->vertex(2)(1) - start(1);
234 data.cell_extents[2] = cell->vertex(4)(2) - start(2);
235 break;
236 default:
237 Assert(false, ExcNotImplemented());
238 }
239 }
240}
241
242
243
244template <int dim, int spacedim>
245void
248 const InternalData & data,
249 std::vector<Point<dim>> &quadrature_points) const
250{
251 if (data.update_each & update_quadrature_points)
252 {
253 const auto offset = QProjector<dim>::DataSetDescriptor::cell();
254
256 }
257}
258
259
260
261template <int dim, int spacedim>
262void
265 const unsigned int face_no,
266 const InternalData & data,
267 std::vector<Point<dim>> &quadrature_points) const
268{
270
271 if (data.update_each & update_quadrature_points)
272 {
274 ReferenceCells::get_hypercube<dim>(),
275 face_no,
276 cell->face_orientation(face_no),
277 cell->face_flip(face_no),
278 cell->face_rotation(face_no),
279 quadrature_points.size());
280
281
283 }
284}
285
286
287
288template <int dim, int spacedim>
289void
292 const unsigned int face_no,
293 const unsigned int sub_no,
294 const InternalData & data,
295 std::vector<Point<dim>> &quadrature_points) const
296{
299 if (cell->face(face_no)->has_children())
300 {
301 AssertIndexRange(sub_no, cell->face(face_no)->n_children());
302 }
303
304 if (data.update_each & update_quadrature_points)
305 {
307 ReferenceCells::get_hypercube<dim>(),
308 face_no,
309 sub_no,
310 cell->face_orientation(face_no),
311 cell->face_flip(face_no),
312 cell->face_rotation(face_no),
313 quadrature_points.size(),
314 cell->subface_case(face_no));
315
317 }
318}
319
320
321
322template <int dim, int spacedim>
323void
326 const InternalData & data,
327 const typename QProjector<dim>::DataSetDescriptor & offset,
328 std::vector<Point<dim>> &quadrature_points) const
329{
330 // let @p{start} be the origin of a local coordinate system. it is chosen
331 // as the (lower) left vertex
332 const Point<dim> &start = cell->vertex(0);
333
334 for (unsigned int i = 0; i < quadrature_points.size(); ++i)
335 {
336 quadrature_points[i] = start;
337 for (unsigned int d = 0; d < dim; ++d)
338 quadrature_points[i](d) +=
339 data.cell_extents[d] * data.quadrature_points[i + offset](d);
340 }
341}
342
343
344
345template <int dim, int spacedim>
346void
348 const unsigned int face_no,
349 const InternalData & data,
350 std::vector<Tensor<1, dim>> &normal_vectors) const
351{
352 // compute normal vectors. All normals on a face have the same value.
353 if (data.update_each & update_normal_vectors)
354 {
356 std::fill(normal_vectors.begin(),
357 normal_vectors.end(),
359 }
360}
361
362
363
364template <int dim, int spacedim>
365void
367 const InternalData & data,
368 const CellSimilarity::Similarity cell_similarity,
370 &output_data) const
371{
372 if (cell_similarity != CellSimilarity::translation)
373 {
374 if (data.update_each & update_jacobian_grads)
375 for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
377
378 if (data.update_each & update_jacobian_pushed_forward_grads)
379 for (unsigned int i = 0;
380 i < output_data.jacobian_pushed_forward_grads.size();
381 ++i)
383
384 if (data.update_each & update_jacobian_2nd_derivatives)
385 for (unsigned int i = 0;
386 i < output_data.jacobian_2nd_derivatives.size();
387 ++i)
388 output_data.jacobian_2nd_derivatives[i] =
390
392 for (unsigned int i = 0;
393 i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
394 ++i)
397
398 if (data.update_each & update_jacobian_3rd_derivatives)
399 for (unsigned int i = 0;
400 i < output_data.jacobian_3rd_derivatives.size();
401 ++i)
402 output_data.jacobian_3rd_derivatives[i] =
404
406 for (unsigned int i = 0;
407 i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
408 ++i)
411 }
412}
413
414
415
416template <int dim, int spacedim>
417void
419 const InternalData &data) const
420{
421 if (data.update_each & update_volume_elements)
422 {
423 double volume = data.cell_extents[0];
424 for (unsigned int d = 1; d < dim; ++d)
425 volume *= data.cell_extents[d];
426 data.volume_element = volume;
427 }
428}
429
430
431
432template <int dim, int spacedim>
433void
435 const InternalData & data,
436 const CellSimilarity::Similarity cell_similarity,
438 &output_data) const
439{
440 // "compute" Jacobian at the quadrature points, which are all the
441 // same
442 if (data.update_each & update_jacobians)
443 if (cell_similarity != CellSimilarity::translation)
444 for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
445 {
447 for (unsigned int j = 0; j < dim; ++j)
448 output_data.jacobians[i][j][j] = data.cell_extents[j];
449 }
450}
451
452
453
454template <int dim, int spacedim>
455void
457 const InternalData & data,
458 const CellSimilarity::Similarity cell_similarity,
460 &output_data) const
461{
462 // "compute" inverse Jacobian at the quadrature points, which are
463 // all the same
464 if (data.update_each & update_inverse_jacobians)
465 if (cell_similarity != CellSimilarity::translation)
466 for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
467 {
468 output_data.inverse_jacobians[i] = Tensor<2, dim>();
469 for (unsigned int j = 0; j < dim; ++j)
470 output_data.inverse_jacobians[i][j][j] = 1. / data.cell_extents[j];
471 }
472}
473
474
475
476template <int dim, int spacedim>
480 const CellSimilarity::Similarity cell_similarity,
481 const Quadrature<dim> & quadrature,
482 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
484 &output_data) const
485{
487
488 // convert data object to internal data for this class. fails with
489 // an exception if that is not possible
490 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
492 const InternalData &data = static_cast<const InternalData &>(internal_data);
493
494
495 update_cell_extents(cell, cell_similarity, data);
496
498 data,
499 output_data.quadrature_points);
500
501 // compute Jacobian determinant. all values are equal and are the
502 // product of the local lengths in each coordinate direction
503 if (data.update_each & (update_JxW_values | update_volume_elements))
504 if (cell_similarity != CellSimilarity::translation)
505 {
506 double J = data.cell_extents[0];
507 for (unsigned int d = 1; d < dim; ++d)
508 J *= data.cell_extents[d];
509 data.volume_element = J;
510 if (data.update_each & update_JxW_values)
511 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
512 output_data.JxW_values[i] = J * quadrature.weight(i);
513 }
514
515
516 maybe_update_jacobians(data, cell_similarity, output_data);
517 maybe_update_jacobian_derivatives(data, cell_similarity, output_data);
518 maybe_update_inverse_jacobians(data, cell_similarity, output_data);
519
520 return cell_similarity;
521}
522
523
524
525template <int dim, int spacedim>
526void
529 const ArrayView<const Point<dim>> & unit_points,
530 const UpdateFlags update_flags,
532 &output_data) const
533{
534 if (update_flags == update_default)
535 return;
536
538
539 Assert(update_flags & update_inverse_jacobians ||
540 update_flags & update_jacobians ||
541 update_flags & update_quadrature_points,
543
544 output_data.initialize(unit_points.size(), update_flags);
545
546 InternalData data;
547 data.update_each = update_flags;
548 data.quadrature_points =
549 std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
550
552
554 data,
555 output_data.quadrature_points);
556
559}
560
561
562
563template <int dim, int spacedim>
564void
567 const unsigned int face_no,
568 const hp::QCollection<dim - 1> & quadrature,
569 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
571 &output_data) const
572{
574 AssertDimension(quadrature.size(), 1);
575
576 // convert data object to internal
577 // data for this class. fails with
578 // an exception if that is not
579 // possible
580 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
582 const InternalData &data = static_cast<const InternalData &>(internal_data);
583
585
587 face_no,
588 data,
589 output_data.quadrature_points);
590
591 maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
592
593 // first compute Jacobian determinant, which is simply the product
594 // of the local lengths since the jacobian is diagonal
595 double J = 1.;
596 for (unsigned int d = 0; d < dim; ++d)
598 J *= data.cell_extents[d];
599
600 if (data.update_each & update_JxW_values)
601 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
602 output_data.JxW_values[i] = J * quadrature[0].weight(i);
603
604 if (data.update_each & update_boundary_forms)
605 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
606 output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
607
612}
613
614
615
616template <int dim, int spacedim>
617void
620 const unsigned int face_no,
621 const unsigned int subface_no,
622 const Quadrature<dim - 1> & quadrature,
623 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
625 &output_data) const
626{
628
629 // convert data object to internal data for this class. fails with
630 // an exception if that is not possible
631 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
633 const InternalData &data = static_cast<const InternalData &>(internal_data);
634
636
638 cell, face_no, subface_no, data, output_data.quadrature_points);
639
640 maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
641
642 // first compute Jacobian determinant, which is simply the product
643 // of the local lengths since the jacobian is diagonal
644 double J = 1.;
645 for (unsigned int d = 0; d < dim; ++d)
647 J *= data.cell_extents[d];
648
649 if (data.update_each & update_JxW_values)
650 {
651 // Here, cell->face(face_no)->n_children() would be the right
652 // choice, but unfortunately the current function is also called
653 // for faces without children (see tests/fe/mapping.cc). Add
654 // following switch to avoid diffs in tests/fe/mapping.OK
655 const unsigned int n_subfaces =
656 cell->face(face_no)->has_children() ?
657 cell->face(face_no)->n_children() :
659 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
660 output_data.JxW_values[i] = J * quadrature.weight(i) / n_subfaces;
661 }
662
663 if (data.update_each & update_boundary_forms)
664 for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
665 output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
666
671}
672
673
674
675template <int dim, int spacedim>
676void
680 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
682 &output_data) const
683{
684 AssertDimension(dim, spacedim);
686
687 // Convert data object to internal data for this class. Fails with an
688 // exception if that is not possible.
689 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
691 const InternalData &data = static_cast<const InternalData &>(internal_data);
692
693
695
697 data,
698 output_data.quadrature_points);
699
700 if (data.update_each & update_normal_vectors)
701 for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
702 {
703 // The normals are n = J^{-T} * \hat{n} before normalizing.
704 Tensor<1, dim> normal;
705 const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
706 for (unsigned int d = 0; d < dim; ++d)
707 {
708 normal[d] = ref_space_normal[d] / data.cell_extents[d];
709 }
710 normal /= normal.norm();
711 output_data.normal_vectors[i] = normal;
712 }
713
714 if (data.update_each & update_JxW_values)
715 for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
716 {
717 const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i);
718
719 // J^{-T} \times \hat{n}
720 Tensor<1, dim> invJTxNormal;
721 double det_jacobian = 1.;
722 for (unsigned int d = 0; d < dim; ++d)
723 {
724 det_jacobian *= data.cell_extents[d];
725 invJTxNormal[d] = ref_space_normal[d] / data.cell_extents[d];
726 }
727 output_data.JxW_values[i] =
728 det_jacobian * invJTxNormal.norm() * quadrature.weight(i);
729 }
730
735}
736
737
738
739template <int dim, int spacedim>
740void
742 const ArrayView<const Tensor<1, dim>> & input,
743 const MappingKind mapping_kind,
744 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
745 const ArrayView<Tensor<1, spacedim>> & output) const
746{
747 AssertDimension(input.size(), output.size());
748 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
750 const InternalData &data = static_cast<const InternalData &>(mapping_data);
751
752 switch (mapping_kind)
753 {
755 {
756 Assert(data.update_each & update_covariant_transformation,
758 "update_covariant_transformation"));
759
760 for (unsigned int i = 0; i < output.size(); ++i)
761 for (unsigned int d = 0; d < dim; ++d)
762 output[i][d] = input[i][d] / data.cell_extents[d];
763 return;
764 }
765
767 {
770 "update_contravariant_transformation"));
771
772 for (unsigned int i = 0; i < output.size(); ++i)
773 for (unsigned int d = 0; d < dim; ++d)
774 output[i][d] = input[i][d] * data.cell_extents[d];
775 return;
776 }
777 case mapping_piola:
778 {
781 "update_contravariant_transformation"));
782 Assert(data.update_each & update_volume_elements,
784 "update_volume_elements"));
785
786 for (unsigned int i = 0; i < output.size(); ++i)
787 for (unsigned int d = 0; d < dim; ++d)
788 output[i][d] =
789 input[i][d] * data.cell_extents[d] / data.volume_element;
790 return;
791 }
792 default:
793 Assert(false, ExcNotImplemented());
794 }
795}
796
797
798
799template <int dim, int spacedim>
800void
803 const MappingKind mapping_kind,
804 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
805 const ArrayView<Tensor<2, spacedim>> & output) const
806{
807 AssertDimension(input.size(), output.size());
808 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
810 const InternalData &data = static_cast<const InternalData &>(mapping_data);
811
812 switch (mapping_kind)
813 {
815 {
816 Assert(data.update_each & update_covariant_transformation,
818 "update_covariant_transformation"));
819
820 for (unsigned int i = 0; i < output.size(); ++i)
821 for (unsigned int d1 = 0; d1 < dim; ++d1)
822 for (unsigned int d2 = 0; d2 < dim; ++d2)
823 output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
824 return;
825 }
826
828 {
831 "update_contravariant_transformation"));
832
833 for (unsigned int i = 0; i < output.size(); ++i)
834 for (unsigned int d1 = 0; d1 < dim; ++d1)
835 for (unsigned int d2 = 0; d2 < dim; ++d2)
836 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
837 return;
838 }
839
841 {
842 Assert(data.update_each & update_covariant_transformation,
844 "update_covariant_transformation"));
845
846 for (unsigned int i = 0; i < output.size(); ++i)
847 for (unsigned int d1 = 0; d1 < dim; ++d1)
848 for (unsigned int d2 = 0; d2 < dim; ++d2)
849 output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
850 data.cell_extents[d1];
851 return;
852 }
853
855 {
858 "update_contravariant_transformation"));
859
860 for (unsigned int i = 0; i < output.size(); ++i)
861 for (unsigned int d1 = 0; d1 < dim; ++d1)
862 for (unsigned int d2 = 0; d2 < dim; ++d2)
863 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
864 data.cell_extents[d1];
865 return;
866 }
867
868 case mapping_piola:
869 {
872 "update_contravariant_transformation"));
873 Assert(data.update_each & update_volume_elements,
875 "update_volume_elements"));
876
877 for (unsigned int i = 0; i < output.size(); ++i)
878 for (unsigned int d1 = 0; d1 < dim; ++d1)
879 for (unsigned int d2 = 0; d2 < dim; ++d2)
880 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
881 data.volume_element;
882 return;
883 }
884
886 {
889 "update_contravariant_transformation"));
890 Assert(data.update_each & update_volume_elements,
892 "update_volume_elements"));
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.cell_extents[d1] / data.volume_element;
899 return;
900 }
901
902 default:
903 Assert(false, ExcNotImplemented());
904 }
905}
906
907
908
909template <int dim, int spacedim>
910void
912 const ArrayView<const Tensor<2, dim>> & input,
913 const MappingKind mapping_kind,
914 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
915 const ArrayView<Tensor<2, spacedim>> & output) const
916{
917 AssertDimension(input.size(), output.size());
918 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
920 const InternalData &data = static_cast<const InternalData &>(mapping_data);
921
922 switch (mapping_kind)
923 {
925 {
926 Assert(data.update_each & update_covariant_transformation,
928 "update_covariant_transformation"));
929
930 for (unsigned int i = 0; i < output.size(); ++i)
931 for (unsigned int d1 = 0; d1 < dim; ++d1)
932 for (unsigned int d2 = 0; d2 < dim; ++d2)
933 output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
934 return;
935 }
936
938 {
941 "update_contravariant_transformation"));
942
943 for (unsigned int i = 0; i < output.size(); ++i)
944 for (unsigned int d1 = 0; d1 < dim; ++d1)
945 for (unsigned int d2 = 0; d2 < dim; ++d2)
946 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
947 return;
948 }
949
951 {
952 Assert(data.update_each & update_covariant_transformation,
954 "update_covariant_transformation"));
955
956 for (unsigned int i = 0; i < output.size(); ++i)
957 for (unsigned int d1 = 0; d1 < dim; ++d1)
958 for (unsigned int d2 = 0; d2 < dim; ++d2)
959 output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
960 data.cell_extents[d1];
961 return;
962 }
963
965 {
968 "update_contravariant_transformation"));
969
970 for (unsigned int i = 0; i < output.size(); ++i)
971 for (unsigned int d1 = 0; d1 < dim; ++d1)
972 for (unsigned int d2 = 0; d2 < dim; ++d2)
973 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
974 data.cell_extents[d1];
975 return;
976 }
977
978 case mapping_piola:
979 {
982 "update_contravariant_transformation"));
983 Assert(data.update_each & update_volume_elements,
985 "update_volume_elements"));
986
987 for (unsigned int i = 0; i < output.size(); ++i)
988 for (unsigned int d1 = 0; d1 < dim; ++d1)
989 for (unsigned int d2 = 0; d2 < dim; ++d2)
990 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
991 data.volume_element;
992 return;
993 }
994
996 {
999 "update_contravariant_transformation"));
1000 Assert(data.update_each & update_volume_elements,
1002 "update_volume_elements"));
1003
1004 for (unsigned int i = 0; i < output.size(); ++i)
1005 for (unsigned int d1 = 0; d1 < dim; ++d1)
1006 for (unsigned int d2 = 0; d2 < dim; ++d2)
1007 output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
1008 data.cell_extents[d1] / data.volume_element;
1009 return;
1010 }
1011
1012 default:
1013 Assert(false, ExcNotImplemented());
1014 }
1015}
1016
1017
1018
1019template <int dim, int spacedim>
1020void
1022 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1023 const MappingKind mapping_kind,
1024 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1025 const ArrayView<Tensor<3, spacedim>> & output) const
1026{
1027 AssertDimension(input.size(), output.size());
1028 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1030 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1031
1032 switch (mapping_kind)
1033 {
1035 {
1038 "update_covariant_transformation"));
1039
1040 for (unsigned int q = 0; q < output.size(); ++q)
1041 for (unsigned int i = 0; i < spacedim; ++i)
1042 for (unsigned int j = 0; j < spacedim; ++j)
1043 for (unsigned int k = 0; k < spacedim; ++k)
1044 {
1045 output[q][i][j][k] = input[q][i][j][k] /
1046 data.cell_extents[j] /
1047 data.cell_extents[k];
1048 }
1049 return;
1050 }
1051 default:
1052 Assert(false, ExcNotImplemented());
1053 }
1054}
1055
1056
1057
1058template <int dim, int spacedim>
1059void
1061 const ArrayView<const Tensor<3, dim>> & input,
1062 const MappingKind mapping_kind,
1063 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1064 const ArrayView<Tensor<3, spacedim>> & output) const
1065{
1066 AssertDimension(input.size(), output.size());
1067 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1069 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1070
1071 switch (mapping_kind)
1072 {
1074 {
1075 Assert(data.update_each & update_covariant_transformation,
1077 "update_covariant_transformation"));
1080 "update_contravariant_transformation"));
1081
1082 for (unsigned int q = 0; q < output.size(); ++q)
1083 for (unsigned int i = 0; i < spacedim; ++i)
1084 for (unsigned int j = 0; j < spacedim; ++j)
1085 for (unsigned int k = 0; k < spacedim; ++k)
1086 {
1087 output[q][i][j][k] =
1088 input[q][i][j][k] * data.cell_extents[i] /
1089 data.cell_extents[j] / data.cell_extents[k];
1090 }
1091 return;
1092 }
1093
1095 {
1096 Assert(data.update_each & update_covariant_transformation,
1098 "update_covariant_transformation"));
1099
1100 for (unsigned int q = 0; q < output.size(); ++q)
1101 for (unsigned int i = 0; i < spacedim; ++i)
1102 for (unsigned int j = 0; j < spacedim; ++j)
1103 for (unsigned int k = 0; k < spacedim; ++k)
1104 {
1105 output[q][i][j][k] =
1106 input[q][i][j][k] / data.cell_extents[i] /
1107 data.cell_extents[j] / data.cell_extents[k];
1108 }
1109
1110 return;
1111 }
1112
1114 {
1115 Assert(data.update_each & update_covariant_transformation,
1117 "update_covariant_transformation"));
1120 "update_contravariant_transformation"));
1121 Assert(data.update_each & update_volume_elements,
1123 "update_volume_elements"));
1124
1125 for (unsigned int q = 0; q < output.size(); ++q)
1126 for (unsigned int i = 0; i < spacedim; ++i)
1127 for (unsigned int j = 0; j < spacedim; ++j)
1128 for (unsigned int k = 0; k < spacedim; ++k)
1129 {
1130 output[q][i][j][k] =
1131 input[q][i][j][k] * data.cell_extents[i] /
1132 data.volume_element / data.cell_extents[j] /
1133 data.cell_extents[k];
1134 }
1135
1136 return;
1137 }
1138
1139 default:
1140 Assert(false, ExcNotImplemented());
1141 }
1142}
1143
1144
1145
1146template <int dim, int spacedim>
1150 const Point<dim> & p) const
1151{
1153
1154 Tensor<1, dim> length;
1155 const Point<dim> start = cell->vertex(0);
1156 switch (dim)
1157 {
1158 case 1:
1159 length[0] = cell->vertex(1)(0) - start(0);
1160 break;
1161 case 2:
1162 length[0] = cell->vertex(1)(0) - start(0);
1163 length[1] = cell->vertex(2)(1) - start(1);
1164 break;
1165 case 3:
1166 length[0] = cell->vertex(1)(0) - start(0);
1167 length[1] = cell->vertex(2)(1) - start(1);
1168 length[2] = cell->vertex(4)(2) - start(2);
1169 break;
1170 default:
1171 Assert(false, ExcNotImplemented());
1172 }
1173
1174 Point<dim> p_real = cell->vertex(0);
1175 for (unsigned int d = 0; d < dim; ++d)
1176 p_real(d) += length[d] * p(d);
1177
1178 return p_real;
1179}
1180
1181
1182
1183template <int dim, int spacedim>
1187 const Point<spacedim> & p) const
1188{
1190
1191 if (dim != spacedim)
1192 Assert(false, ExcNotImplemented());
1193 const Point<dim> &start = cell->vertex(0);
1194 Point<dim> real = p;
1195 real -= start;
1196
1197 switch (dim)
1198 {
1199 case 1:
1200 real(0) /= cell->vertex(1)(0) - start(0);
1201 break;
1202 case 2:
1203 real(0) /= cell->vertex(1)(0) - start(0);
1204 real(1) /= cell->vertex(2)(1) - start(1);
1205 break;
1206 case 3:
1207 real(0) /= cell->vertex(1)(0) - start(0);
1208 real(1) /= cell->vertex(2)(1) - start(1);
1209 real(2) /= cell->vertex(4)(2) - start(2);
1210 break;
1211 default:
1212 Assert(false, ExcNotImplemented());
1213 }
1214 return real;
1215}
1216
1217
1218
1219template <int dim, int spacedim>
1220std::unique_ptr<Mapping<dim, spacedim>>
1222{
1223 return std::make_unique<MappingCartesian<dim, spacedim>>(*this);
1224}
1225
1226
1227//---------------------------------------------------------------------------
1228// explicit instantiations
1229#include "mapping_cartesian.inst"
1230
1231
std::vector< Point< dim > > quadrature_points
virtual std::size_t memory_consumption() 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
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
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
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
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_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
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 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 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 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:311
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition: point.h:111
static DataSetDescriptor cell()
Definition: qprojector.h:563
static DataSetDescriptor subface(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 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:1365
double weight(const unsigned int i) const
Definition: tensor.h:503
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition: collection.h:264
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
Definition: fe_values.cc:2704
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< Point< spacedim > > quadrature_points
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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
static ::ExceptionBase & ExcCellNotCartesian()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
MappingKind
Definition: mapping.h:72
@ mapping_piola
Definition: mapping.h:107
@ mapping_covariant_gradient
Definition: mapping.h:93
@ mapping_covariant
Definition: mapping.h:82
@ mapping_contravariant
Definition: mapping.h:87
@ mapping_contravariant_hessian
Definition: mapping.h:149
@ mapping_covariant_hessian
Definition: mapping.h:143
@ mapping_contravariant_gradient
Definition: mapping.h:99
@ mapping_piola_gradient
Definition: mapping.h:113
@ mapping_piola_hessian
Definition: mapping.h:155
bool is_cartesian(const CellType &cell)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:139
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type 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:493
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
T signaling_nan()