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