deal.II version GIT relicensing-2330-gf6dfc6c370 2025-01-06 13:10: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_manifold.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 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
15
23
25
26#include <deal.II/fe/fe.h>
27#include <deal.II/fe/fe_tools.h>
30
32#include <deal.II/grid/tria.h>
34
36
37#include <algorithm>
38#include <array>
39#include <cmath>
40#include <memory>
41#include <numeric>
42
43
45
46template <int dim, int spacedim>
47std::size_t
64
65
66
67template <int dim, int spacedim>
68void
70 const UpdateFlags update_flags,
71 const Quadrature<dim> &q)
72{
73 // store the flags in the internal data object so we can access them
74 // in fill_fe_*_values()
75 this->update_each = update_flags;
76
77 const unsigned int n_q_points = q.size();
78
79 // Store the quadrature
80 this->quad.initialize(q.get_points(), q.get_weights());
81
82 // see if we need the (transformation) shape function values
83 // and/or gradients and resize the necessary arrays
84 if (this->update_each &
86 compute_manifold_quadrature_weights(q);
87
88 if (this->update_each & update_covariant_transformation)
89 covariant.resize(n_q_points);
90
91 if (this->update_each & update_contravariant_transformation)
92 contravariant.resize(n_q_points);
93
94 if (this->update_each & update_volume_elements)
95 volume_elements.resize(n_q_points);
96}
97
98
99
100template <int dim, int spacedim>
101void
103 const UpdateFlags update_flags,
104 const Quadrature<dim> &q,
105 const unsigned int n_original_q_points)
106{
107 reinit(update_flags, q);
108
109 // Set to the size of a single quadrature object for faces, as the size set
110 // in in reinit() is for all points
111 if (this->update_each & update_covariant_transformation)
112 covariant.resize(n_original_q_points);
113
114 if (this->update_each & update_contravariant_transformation)
115 contravariant.resize(n_original_q_points);
116
117 if (this->update_each & update_volume_elements)
118 volume_elements.resize(n_original_q_points);
119
120 if (dim > 1)
121 {
122 if (this->update_each & update_boundary_forms)
123 {
124 aux.resize(dim - 1,
125 std::vector<Tensor<1, spacedim>>(n_original_q_points));
126
127 // Compute tangentials to the unit cell.
128 for (const unsigned int i : GeometryInfo<dim>::face_indices())
129 {
130 unit_tangentials[i].resize(n_original_q_points);
131 std::fill(unit_tangentials[i].begin(),
132 unit_tangentials[i].end(),
134 if (dim > 2)
135 {
136 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
137 .resize(n_original_q_points);
138 std::fill(
139 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
140 .begin(),
141 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
142 .end(),
144 }
145 }
146 }
147 }
148}
149
150
151
152template <int dim, int spacedim>
153void
155 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
156{
157 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
158 vertices[i] = cell->vertex(i);
159 this->cell = cell;
160}
161
162
163
164template <int dim, int spacedim>
165void
168{
169 cell_manifold_quadrature_weights.resize(quad.size());
170 for (unsigned int q = 0; q < quad.size(); ++q)
171 {
172 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
173 {
174 cell_manifold_quadrature_weights[q][i] =
176 }
177 }
178}
179
180
181
182template <int dim, int spacedim>
186
187
188
189template <int dim, int spacedim>
190std::unique_ptr<Mapping<dim, spacedim>>
192{
193 return std::make_unique<MappingManifold<dim, spacedim>>(*this);
194}
195
196
197
198template <int dim, int spacedim>
207
208
209
210template <int dim, int spacedim>
214 const Point<dim> &p) const
215{
216 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
217 std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;
218
219 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
220 {
221 vertices[v] = cell->vertex(v);
223 }
224 return cell->get_manifold().get_new_point(
225 make_array_view(vertices.begin(), vertices.end()),
226 make_array_view(weights.begin(), weights.end()));
227}
228
229
230
231// In the code below, GCC tries to instantiate MappingManifold<3,4> when
232// seeing which of the overloaded versions of
233// do_transform_real_to_unit_cell_internal() to call. This leads to bad
234// error messages and, generally, nothing very good. Avoid this by ensuring
235// that this class exists, but does not have an inner InternalData
236// type, thereby ruling out the codim-1 version of the function
237// below when doing overload resolution.
238template <>
240{};
241
242
243
244template <int dim, int spacedim>
247 const UpdateFlags in) const
248{
249 // add flags if the respective quantities are necessary to compute
250 // what we need. note that some flags appear in both the conditions
251 // and in subsequent set operations. this leads to some circular
252 // logic. the only way to treat this is to iterate. since there are
253 // 5 if-clauses in the loop, it will take at most 5 iterations to
254 // converge. do them:
255 UpdateFlags out = in;
256 for (unsigned int i = 0; i < 5; ++i)
257 {
258 // The following is a little incorrect:
259 // If not applied on a face,
260 // update_boundary_forms does not
261 // make sense. On the other hand,
262 // it is necessary on a
263 // face. Currently,
264 // update_boundary_forms is simply
265 // ignored for the interior of a
266 // cell.
269
274
275 if (out &
280
281 // The contravariant transformation used in the Piola
282 // transformation, which requires the determinant of the Jacobi
283 // matrix of the transformation. Because we have no way of
284 // knowing here whether the finite elements wants to use the
285 // contravariant of the Piola transforms, we add the JxW values
286 // to the list of flags to be updated for each cell.
288 out |= update_JxW_values;
289
290 if (out & update_normal_vectors)
291 out |= update_JxW_values;
292 }
293
294 // Now throw an exception if we stumble upon something that was not
295 // implemented yet
300
301 return out;
302}
303
304
305
306template <int dim, int spacedim>
307std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
309 const Quadrature<dim> &q) const
310{
311 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
312 std::make_unique<InternalData>();
313 data_ptr->reinit(this->requires_update_flags(update_flags), q);
314
315 return data_ptr;
316}
317
318
319
320template <int dim, int spacedim>
321std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
323 const UpdateFlags update_flags,
324 const hp::QCollection<dim - 1> &quadrature) const
325{
326 AssertDimension(quadrature.size(), 1);
327
328 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
329 std::make_unique<InternalData>();
330 auto &data = dynamic_cast<InternalData &>(*data_ptr);
331 data.initialize_face(this->requires_update_flags(update_flags),
333 ReferenceCells::get_hypercube<dim>(), quadrature[0]),
334 quadrature[0].size());
335
336 return data_ptr;
337}
338
339
340
341template <int dim, int spacedim>
342std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
344 const UpdateFlags update_flags,
345 const Quadrature<dim - 1> &quadrature) const
346{
347 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
348 std::make_unique<InternalData>();
349 auto &data = dynamic_cast<InternalData &>(*data_ptr);
350 data.initialize_face(this->requires_update_flags(update_flags),
352 ReferenceCells::get_hypercube<dim>(), quadrature),
353 quadrature.size());
354
355 return data_ptr;
356}
357
358
359
360namespace internal
361{
362 namespace MappingManifoldImplementation
363 {
364 namespace
365 {
372 template <int dim, int spacedim>
373 void
374 maybe_compute_q_points(
375 const typename QProjector<dim>::DataSetDescriptor data_set,
376 const typename ::MappingManifold<dim, spacedim>::InternalData
377 &data,
378 std::vector<Point<spacedim>> &quadrature_points)
379 {
380 const UpdateFlags update_flags = data.update_each;
381
382 if (update_flags & update_quadrature_points)
383 {
384 for (unsigned int point = 0; point < quadrature_points.size();
385 ++point)
386 {
387 quadrature_points[point] = data.manifold->get_new_point(
388 make_array_view(data.vertices),
390 data.cell_manifold_quadrature_weights[point + data_set]));
391 }
392 }
393 }
394
395
396
403 template <int dim, int spacedim>
404 void
405 maybe_update_Jacobians(
406 const typename ::QProjector<dim>::DataSetDescriptor data_set,
407 const typename ::MappingManifold<dim, spacedim>::InternalData
408 &data)
409 {
410 const UpdateFlags update_flags = data.update_each;
411
412 if (update_flags & update_contravariant_transformation)
413 {
414 const unsigned int n_q_points = data.contravariant.size();
415
416 std::fill(data.contravariant.begin(),
417 data.contravariant.end(),
419
420 for (unsigned int point = 0; point < n_q_points; ++point)
421 {
422 // Start by figuring out how to compute the direction in
423 // the reference space:
424 const Point<dim> &p = data.quad.point(point + data_set);
425
426 // And get its image on the manifold:
427 const Point<spacedim> P = data.manifold->get_new_point(
428 make_array_view(data.vertices),
430 data.cell_manifold_quadrature_weights[point + data_set]));
431
432 // To compute the Jacobian, we choose dim points aligned
433 // with the dim reference axes, which are still in the
434 // given cell, and ask for the tangent vector in these
435 // directions. Choosing the points is somewhat arbitrary,
436 // so we try to be smart and we pick points which are
437 // on the opposite quadrant w.r.t. the evaluation
438 // point.
439 for (unsigned int i = 0; i < dim; ++i)
440 {
442 const double pi = p[i];
443 Assert(pi >= 0 && pi <= 1.0,
445 "Was expecting a quadrature point "
446 "inside the unit reference element."));
447
448 // In the length L, we store also the direction sign,
449 // which is positive, if the coordinate is < .5,
450 const double L = pi > .5 ? -pi : 1 - pi;
451
452 const Point<dim> np(p + L * ei);
453
454 // Get the weights to compute the np point in real space
455 for (const unsigned int j :
457 data.vertex_weights[j] =
459
460 const Point<spacedim> NP = data.manifold->get_new_point(
461 make_array_view(data.vertices),
462 make_array_view(data.vertex_weights));
463
464 const Tensor<1, spacedim> T =
465 data.manifold->get_tangent_vector(P, NP);
466
467 for (unsigned int d = 0; d < spacedim; ++d)
468 data.contravariant[point][d][i] = T[d] / L;
469 }
470 }
471
472 if (update_flags & update_covariant_transformation)
473 {
474 const unsigned int n_q_points = data.contravariant.size();
475 for (unsigned int point = 0; point < n_q_points; ++point)
476 {
477 data.covariant[point] =
478 (data.contravariant[point]).covariant_form();
479 }
480 }
481
482 if (update_flags & update_volume_elements)
483 {
484 const unsigned int n_q_points = data.contravariant.size();
485 for (unsigned int point = 0; point < n_q_points; ++point)
486 data.volume_elements[point] =
487 data.contravariant[point].determinant();
488 }
489 }
490 }
491 } // namespace
492 } // namespace MappingManifoldImplementation
493} // namespace internal
494
495
496
497template <int dim, int spacedim>
502 const Quadrature<dim> &quadrature,
503 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
505 &output_data) const
506{
507 // ensure that the following static_cast is really correct:
508 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
510 const InternalData &data = static_cast<const InternalData &>(internal_data);
511
512 const unsigned int n_q_points = quadrature.size();
513
514 data.store_vertices(cell);
515 data.manifold = &(cell->get_manifold());
516
517 internal::MappingManifoldImplementation::maybe_compute_q_points<dim,
518 spacedim>(
520 data,
521 output_data.quadrature_points);
522
523 internal::MappingManifoldImplementation::maybe_update_Jacobians<dim,
524 spacedim>(
526
527 const UpdateFlags update_flags = data.update_each;
528 const std::vector<double> &weights = quadrature.get_weights();
529
530 // Multiply quadrature weights by absolute value of Jacobian determinants or
531 // the area element g=sqrt(DX^t DX) in case of codim > 0
532
533 if (update_flags & (update_normal_vectors | update_JxW_values))
534 {
535 AssertDimension(output_data.JxW_values.size(), n_q_points);
536
537 Assert(!(update_flags & update_normal_vectors) ||
538 (output_data.normal_vectors.size() == n_q_points),
539 ExcDimensionMismatch(output_data.normal_vectors.size(),
540 n_q_points));
541
542
543 for (unsigned int point = 0; point < n_q_points; ++point)
544 {
545 if (dim == spacedim)
546 {
547 const double det = data.contravariant[point].determinant();
548
549 // check for distorted cells.
550
551 // TODO: this allows for anisotropies of up to 1e6 in 3d and
552 // 1e12 in 2d. might want to find a finer
553 // (dimension-independent) criterion
554 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
555 cell->diameter() / std::sqrt(double(dim))),
557 cell->center(), det, point)));
558
559 output_data.JxW_values[point] = weights[point] * det;
560 }
561 // if dim==spacedim, then there is no cell normal to
562 // compute. since this is for FEValues (and not FEFaceValues),
563 // there are also no face normals to compute
564 else // codim>0 case
565 {
566 Tensor<1, spacedim> DX_t[dim];
567 for (unsigned int i = 0; i < spacedim; ++i)
568 for (unsigned int j = 0; j < dim; ++j)
569 DX_t[j][i] = data.contravariant[point][i][j];
570
571 Tensor<2, dim> G; // First fundamental form
572 for (unsigned int i = 0; i < dim; ++i)
573 for (unsigned int j = 0; j < dim; ++j)
574 G[i][j] = DX_t[i] * DX_t[j];
575
576 output_data.JxW_values[point] =
577 std::sqrt(determinant(G)) * weights[point];
578
579 if (update_flags & update_normal_vectors)
580 {
581 Assert(spacedim == dim + 1,
583 "There is no (unique) cell normal for " +
585 "-dimensional cells in " +
586 Utilities::int_to_string(spacedim) +
587 "-dimensional space. This only works if the "
588 "space dimension is one greater than the "
589 "dimensionality of the mesh cells."));
590
591 if (dim == 1)
592 output_data.normal_vectors[point] =
593 cross_product_2d(-DX_t[0]);
594 else // dim == 2
595 output_data.normal_vectors[point] =
596 cross_product_3d(DX_t[0], DX_t[1]);
597
598 output_data.normal_vectors[point] /=
599 output_data.normal_vectors[point].norm();
600
601 if (cell->direction_flag() == false)
602 output_data.normal_vectors[point] *= -1.;
603 }
604 } // codim>0 case
605 }
606 }
607
608
609
610 // copy values from InternalData to vector given by reference
611 if (update_flags & update_jacobians)
612 {
613 AssertDimension(output_data.jacobians.size(), n_q_points);
614 for (unsigned int point = 0; point < n_q_points; ++point)
615 output_data.jacobians[point] = data.contravariant[point];
616 }
617
618 // copy values from InternalData to vector given by reference
619 if (update_flags & update_inverse_jacobians)
620 {
621 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
622 for (unsigned int point = 0; point < n_q_points; ++point)
623 output_data.inverse_jacobians[point] =
624 data.covariant[point].transpose();
625 }
626
628}
629
630
631
632namespace internal
633{
634 namespace MappingManifoldImplementation
635 {
636 namespace
637 {
647 template <int dim, int spacedim>
648 void
649 maybe_compute_face_data(
650 const ::MappingManifold<dim, spacedim> &mapping,
651 const typename ::Triangulation<dim, spacedim>::cell_iterator
652 &cell,
653 const unsigned int face_no,
654 const unsigned int subface_no,
655 const unsigned int n_q_points,
656 const std::vector<double> &weights,
657 const typename ::MappingManifold<dim, spacedim>::InternalData
658 &data,
660 &output_data)
661 {
662 const UpdateFlags update_flags = data.update_each;
663
664 if (update_flags & update_boundary_forms)
665 {
666 AssertDimension(output_data.boundary_forms.size(), n_q_points);
667 if (update_flags & update_normal_vectors)
668 AssertDimension(output_data.normal_vectors.size(), n_q_points);
669 if (update_flags & update_JxW_values)
670 AssertDimension(output_data.JxW_values.size(), n_q_points);
671
672 // map the unit tangentials to the real cell. checking for d!=dim-1
673 // eliminates compiler warnings regarding unsigned int expressions <
674 // 0.
675 for (unsigned int d = 0; d != dim - 1; ++d)
676 {
678 data.unit_tangentials.size(),
680 Assert(
681 data.aux[d].size() <=
682 data
683 .unit_tangentials[face_no +
685 .size(),
687
688 mapping.transform(
690 data
691 .unit_tangentials[face_no +
694 data,
695 make_array_view(data.aux[d]));
696 }
697
698 // if dim==spacedim, we can use the unit tangentials to compute the
699 // boundary form by simply taking the cross product
700 if (dim == spacedim)
701 {
702 for (unsigned int i = 0; i < n_q_points; ++i)
703 switch (dim)
704 {
705 case 1:
706 // in 1d, we don't have access to any of the data.aux
707 // fields (because it has only dim-1 components), but we
708 // can still compute the boundary form by simply
709 // looking at the number of the face
710 output_data.boundary_forms[i][0] =
711 (face_no == 0 ? -1 : +1);
712 break;
713 case 2:
714 output_data.boundary_forms[i] =
715 cross_product_2d(data.aux[0][i]);
716 break;
717 case 3:
718 output_data.boundary_forms[i] =
719 cross_product_3d(data.aux[0][i], data.aux[1][i]);
720 break;
721 default:
723 }
724 }
725 else //(dim < spacedim)
726 {
727 // in the codim-one case, the boundary form results from the
728 // cross product of all the face tangential vectors and the cell
729 // normal vector
730 //
731 // to compute the cell normal, use the same method used in
732 // fill_fe_values for cells above
733 AssertDimension(data.contravariant.size(), n_q_points);
734
735 for (unsigned int point = 0; point < n_q_points; ++point)
736 {
737 switch (dim)
738 {
739 case 1:
740 {
741 // J is a tangent vector
742 output_data.boundary_forms[point] =
743 data.contravariant[point].transpose()[0];
744 output_data.boundary_forms[point] /=
745 (face_no == 0 ? -1. : +1.) *
746 output_data.boundary_forms[point].norm();
747
748 break;
749 }
750
751 case 2:
752 {
754 data.contravariant[point].transpose();
755
756 Tensor<1, spacedim> cell_normal =
757 cross_product_3d(DX_t[0], DX_t[1]);
758 cell_normal /= cell_normal.norm();
759
760 // then compute the face normal from the face
761 // tangent and the cell normal:
762 output_data.boundary_forms[point] =
763 cross_product_3d(data.aux[0][point], cell_normal);
764
765 break;
766 }
767
768 default:
770 }
771 }
772 }
773
774 if (update_flags & (update_normal_vectors | update_JxW_values))
775 for (unsigned int i = 0; i < output_data.boundary_forms.size();
776 ++i)
777 {
778 if (update_flags & update_JxW_values)
779 {
780 output_data.JxW_values[i] =
781 output_data.boundary_forms[i].norm() * weights[i];
782
783 if (subface_no != numbers::invalid_unsigned_int)
784 {
785 const double area_ratio =
787 cell->subface_case(face_no), subface_no);
788 output_data.JxW_values[i] *= area_ratio;
789 }
790 }
791
792 if (update_flags & update_normal_vectors)
793 output_data.normal_vectors[i] =
794 Point<spacedim>(output_data.boundary_forms[i] /
795 output_data.boundary_forms[i].norm());
796 }
797
798 if (update_flags & update_jacobians)
799 for (unsigned int point = 0; point < n_q_points; ++point)
800 output_data.jacobians[point] = data.contravariant[point];
801
802 if (update_flags & update_inverse_jacobians)
803 for (unsigned int point = 0; point < n_q_points; ++point)
804 output_data.inverse_jacobians[point] =
805 data.covariant[point].transpose();
806 }
807 }
808
809
816 template <int dim, int spacedim>
817 void
819 const ::MappingManifold<dim, spacedim> &mapping,
820 const typename ::Triangulation<dim, spacedim>::cell_iterator
821 &cell,
822 const unsigned int face_no,
823 const unsigned int subface_no,
824 const typename QProjector<dim>::DataSetDescriptor data_set,
825 const Quadrature<dim - 1> &quadrature,
826 const typename ::MappingManifold<dim, spacedim>::InternalData
827 &data,
829 &output_data)
830 {
831 data.store_vertices(cell);
832
833 data.manifold = &cell->face(face_no)->get_manifold();
834
835 maybe_compute_q_points<dim, spacedim>(data_set,
836 data,
837 output_data.quadrature_points);
838 maybe_update_Jacobians<dim, spacedim>(data_set, data);
839
841 cell,
842 face_no,
843 subface_no,
844 quadrature.size(),
845 quadrature.get_weights(),
846 data,
847 output_data);
848 }
849
850 template <int dim, int spacedim, int rank>
851 void
853 const ArrayView<const Tensor<rank, dim>> &input,
854 const MappingKind mapping_kind,
855 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
856 const ArrayView<Tensor<rank, spacedim>> &output)
857 {
858 AssertDimension(input.size(), output.size());
859 Assert((dynamic_cast<const typename ::
860 MappingManifold<dim, spacedim>::InternalData *>(
861 &mapping_data) != nullptr),
863 const typename ::MappingManifold<dim, spacedim>::InternalData
864 &data =
865 static_cast<const typename ::MappingManifold<dim, spacedim>::
866 InternalData &>(mapping_data);
867
868 switch (mapping_kind)
869 {
871 {
872 Assert(
875 "update_contravariant_transformation"));
876
877 for (unsigned int i = 0; i < output.size(); ++i)
878 output[i] =
879 apply_transformation(data.contravariant[i], input[i]);
880
881 return;
882 }
883
884 case mapping_piola:
885 {
886 Assert(
889 "update_contravariant_transformation"));
890 Assert(
891 data.update_each & update_volume_elements,
893 "update_volume_elements"));
894 Assert(rank == 1, ExcMessage("Only for rank 1"));
895 if (rank != 1)
896 return;
897
898 for (unsigned int i = 0; i < output.size(); ++i)
899 {
900 output[i] =
901 apply_transformation(data.contravariant[i], input[i]);
902 output[i] /= data.volume_elements[i];
903 }
904 return;
905 }
906 // We still allow this operation as in the
907 // reference cell Derivatives are Tensor
908 // rather than DerivativeForm
910 {
911 Assert(
914 "update_covariant_transformation"));
915
916 for (unsigned int i = 0; i < output.size(); ++i)
917 output[i] = apply_transformation(data.covariant[i], input[i]);
918
919 return;
920 }
921
922 default:
924 }
925 }
926
927
928 template <int dim, int spacedim, int rank>
929 void
931 const ArrayView<const Tensor<rank, dim>> &input,
932 const MappingKind mapping_kind,
933 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
934 const ArrayView<Tensor<rank, spacedim>> &output)
935 {
936 AssertDimension(input.size(), output.size());
937 Assert((dynamic_cast<const typename ::
938 MappingManifold<dim, spacedim>::InternalData *>(
939 &mapping_data) != nullptr),
941 const typename ::MappingManifold<dim, spacedim>::InternalData
942 &data =
943 static_cast<const typename ::MappingManifold<dim, spacedim>::
944 InternalData &>(mapping_data);
945
946 switch (mapping_kind)
947 {
949 {
950 Assert(
953 "update_covariant_transformation"));
954 Assert(
957 "update_contravariant_transformation"));
958 Assert(rank == 2, ExcMessage("Only for rank 2"));
959
960 for (unsigned int i = 0; i < output.size(); ++i)
961 {
963 apply_transformation(data.contravariant[i],
964 transpose(input[i]));
965 output[i] =
966 apply_transformation(data.covariant[i], A.transpose());
967 }
968
969 return;
970 }
971
973 {
974 Assert(
977 "update_covariant_transformation"));
978 Assert(rank == 2, ExcMessage("Only for rank 2"));
979
980 for (unsigned int i = 0; i < output.size(); ++i)
981 {
983 apply_transformation(data.covariant[i],
984 transpose(input[i]));
985 output[i] =
986 apply_transformation(data.covariant[i], A.transpose());
987 }
988
989 return;
990 }
991
993 {
994 Assert(
997 "update_covariant_transformation"));
998 Assert(
1001 "update_contravariant_transformation"));
1002 Assert(
1003 data.update_each & update_volume_elements,
1005 "update_volume_elements"));
1006 Assert(rank == 2, ExcMessage("Only for rank 2"));
1007
1008 for (unsigned int i = 0; i < output.size(); ++i)
1009 {
1011 apply_transformation(data.covariant[i], input[i]);
1013 apply_transformation(data.contravariant[i],
1014 A.transpose());
1015
1016 output[i] = transpose(T);
1017 output[i] /= data.volume_elements[i];
1018 }
1019
1020 return;
1021 }
1022
1023 default:
1025 }
1026 }
1027
1028
1029
1030 template <int dim, int spacedim>
1031 void
1033 const ArrayView<const Tensor<3, dim>> &input,
1034 const MappingKind mapping_kind,
1035 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1036 const ArrayView<Tensor<3, spacedim>> &output)
1037 {
1038 AssertDimension(input.size(), output.size());
1039 Assert((dynamic_cast<const typename ::
1040 MappingManifold<dim, spacedim>::InternalData *>(
1041 &mapping_data) != nullptr),
1043 const typename ::MappingManifold<dim, spacedim>::InternalData
1044 &data =
1045 static_cast<const typename ::MappingManifold<dim, spacedim>::
1046 InternalData &>(mapping_data);
1047
1048 switch (mapping_kind)
1049 {
1051 {
1052 Assert(
1055 "update_covariant_transformation"));
1056 Assert(
1059 "update_contravariant_transformation"));
1060
1061 for (unsigned int q = 0; q < output.size(); ++q)
1062 for (unsigned int i = 0; i < spacedim; ++i)
1063 {
1064 double tmp1[dim][dim];
1065 for (unsigned int J = 0; J < dim; ++J)
1066 for (unsigned int K = 0; K < dim; ++K)
1067 {
1068 tmp1[J][K] =
1069 data.contravariant[q][i][0] * input[q][0][J][K];
1070 for (unsigned int I = 1; I < dim; ++I)
1071 tmp1[J][K] +=
1072 data.contravariant[q][i][I] * input[q][I][J][K];
1073 }
1074 for (unsigned int j = 0; j < spacedim; ++j)
1075 {
1076 double tmp2[dim];
1077 for (unsigned int K = 0; K < dim; ++K)
1078 {
1079 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1080 for (unsigned int J = 1; J < dim; ++J)
1081 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1082 }
1083 for (unsigned int k = 0; k < spacedim; ++k)
1084 {
1085 output[q][i][j][k] =
1086 data.covariant[q][k][0] * tmp2[0];
1087 for (unsigned int K = 1; K < dim; ++K)
1088 output[q][i][j][k] +=
1089 data.covariant[q][k][K] * tmp2[K];
1090 }
1091 }
1092 }
1093 return;
1094 }
1095
1097 {
1098 Assert(
1101 "update_covariant_transformation"));
1102
1103 for (unsigned int q = 0; q < output.size(); ++q)
1104 for (unsigned int i = 0; i < spacedim; ++i)
1105 {
1106 double tmp1[dim][dim];
1107 for (unsigned int J = 0; J < dim; ++J)
1108 for (unsigned int K = 0; K < dim; ++K)
1109 {
1110 tmp1[J][K] =
1111 data.covariant[q][i][0] * input[q][0][J][K];
1112 for (unsigned int I = 1; I < dim; ++I)
1113 tmp1[J][K] +=
1114 data.covariant[q][i][I] * input[q][I][J][K];
1115 }
1116 for (unsigned int j = 0; j < spacedim; ++j)
1117 {
1118 double tmp2[dim];
1119 for (unsigned int K = 0; K < dim; ++K)
1120 {
1121 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1122 for (unsigned int J = 1; J < dim; ++J)
1123 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1124 }
1125 for (unsigned int k = 0; k < spacedim; ++k)
1126 {
1127 output[q][i][j][k] =
1128 data.covariant[q][k][0] * tmp2[0];
1129 for (unsigned int K = 1; K < dim; ++K)
1130 output[q][i][j][k] +=
1131 data.covariant[q][k][K] * tmp2[K];
1132 }
1133 }
1134 }
1135
1136 return;
1137 }
1138
1140 {
1141 Assert(
1144 "update_covariant_transformation"));
1145 Assert(
1148 "update_contravariant_transformation"));
1149 Assert(
1150 data.update_each & update_volume_elements,
1152 "update_volume_elements"));
1153
1154 for (unsigned int q = 0; q < output.size(); ++q)
1155 for (unsigned int i = 0; i < spacedim; ++i)
1156 {
1157 double factor[dim];
1158 for (unsigned int I = 0; I < dim; ++I)
1159 factor[I] =
1160 data.contravariant[q][i][I] / data.volume_elements[q];
1161 double tmp1[dim][dim];
1162 for (unsigned int J = 0; J < dim; ++J)
1163 for (unsigned int K = 0; K < dim; ++K)
1164 {
1165 tmp1[J][K] = factor[0] * input[q][0][J][K];
1166 for (unsigned int I = 1; I < dim; ++I)
1167 tmp1[J][K] += factor[I] * input[q][I][J][K];
1168 }
1169 for (unsigned int j = 0; j < spacedim; ++j)
1170 {
1171 double tmp2[dim];
1172 for (unsigned int K = 0; K < dim; ++K)
1173 {
1174 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1175 for (unsigned int J = 1; J < dim; ++J)
1176 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1177 }
1178 for (unsigned int k = 0; k < spacedim; ++k)
1179 {
1180 output[q][i][j][k] =
1181 data.covariant[q][k][0] * tmp2[0];
1182 for (unsigned int K = 1; K < dim; ++K)
1183 output[q][i][j][k] +=
1184 data.covariant[q][k][K] * tmp2[K];
1185 }
1186 }
1187 }
1188
1189 return;
1190 }
1191
1192 default:
1194 }
1195 }
1196
1197
1198
1199 template <int dim, int spacedim, int rank>
1200 void
1203 const MappingKind mapping_kind,
1204 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1206 {
1207 AssertDimension(input.size(), output.size());
1208 Assert((dynamic_cast<const typename ::
1209 MappingManifold<dim, spacedim>::InternalData *>(
1210 &mapping_data) != nullptr),
1212 const typename ::MappingManifold<dim, spacedim>::InternalData
1213 &data =
1214 static_cast<const typename ::MappingManifold<dim, spacedim>::
1215 InternalData &>(mapping_data);
1216
1217 switch (mapping_kind)
1218 {
1219 case mapping_covariant:
1220 {
1221 Assert(
1224 "update_covariant_transformation"));
1225
1226 for (unsigned int i = 0; i < output.size(); ++i)
1227 output[i] = apply_transformation(data.covariant[i], input[i]);
1228
1229 return;
1230 }
1231 default:
1233 }
1234 }
1235 } // namespace
1236 } // namespace MappingManifoldImplementation
1237} // namespace internal
1238
1239
1240
1241template <int dim, int spacedim>
1242void
1245 const unsigned int face_no,
1246 const hp::QCollection<dim - 1> &quadrature,
1247 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1249 &output_data) const
1250{
1251 AssertDimension(quadrature.size(), 1);
1252
1253 // ensure that the following cast is really correct:
1254 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1256 const InternalData &data = static_cast<const InternalData &>(internal_data);
1257
1258 internal::MappingManifoldImplementation::do_fill_fe_face_values(
1259 *this,
1260 cell,
1261 face_no,
1264 ReferenceCells::get_hypercube<dim>(),
1265 face_no,
1266 cell->combined_face_orientation(face_no),
1267 quadrature[0].size()),
1268 quadrature[0],
1269 data,
1270 output_data);
1271}
1272
1273
1274
1275template <int dim, int spacedim>
1276void
1279 const unsigned int face_no,
1280 const unsigned int subface_no,
1281 const Quadrature<dim - 1> &quadrature,
1282 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1284 &output_data) const
1285{
1286 // ensure that the following cast is really correct:
1287 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1289 const InternalData &data = static_cast<const InternalData &>(internal_data);
1290
1291 internal::MappingManifoldImplementation::do_fill_fe_face_values(
1292 *this,
1293 cell,
1294 face_no,
1295 subface_no,
1297 ReferenceCells::get_hypercube<dim>(),
1298 face_no,
1299 subface_no,
1300 cell->combined_face_orientation(face_no),
1301 quadrature.size(),
1302 cell->subface_case(face_no)),
1303 quadrature,
1304 data,
1305 output_data);
1306}
1307
1308
1309
1310template <int dim, int spacedim>
1311void
1313 const ArrayView<const Tensor<1, dim>> &input,
1314 const MappingKind mapping_kind,
1315 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1316 const ArrayView<Tensor<1, spacedim>> &output) const
1317{
1318 internal::MappingManifoldImplementation::transform_fields(input,
1319 mapping_kind,
1320 mapping_data,
1321 output);
1322}
1323
1324
1325
1326template <int dim, int spacedim>
1327void
1329 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1330 const MappingKind mapping_kind,
1331 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1332 const ArrayView<Tensor<2, spacedim>> &output) const
1333{
1334 internal::MappingManifoldImplementation::transform_differential_forms(
1335 input, mapping_kind, mapping_data, output);
1336}
1337
1338
1339
1340template <int dim, int spacedim>
1341void
1343 const ArrayView<const Tensor<2, dim>> &input,
1344 const MappingKind mapping_kind,
1345 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1346 const ArrayView<Tensor<2, spacedim>> &output) const
1347{
1348 switch (mapping_kind)
1349 {
1351 internal::MappingManifoldImplementation::transform_fields(input,
1352 mapping_kind,
1353 mapping_data,
1354 output);
1355 return;
1356
1360 internal::MappingManifoldImplementation::transform_gradients(
1361 input, mapping_kind, mapping_data, output);
1362 return;
1363 default:
1365 }
1366}
1367
1368
1369
1370template <int dim, int spacedim>
1371void
1373 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1374 const MappingKind mapping_kind,
1375 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1376 const ArrayView<Tensor<3, spacedim>> &output) const
1377{
1378 AssertDimension(input.size(), output.size());
1379 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1381 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1382
1383 switch (mapping_kind)
1384 {
1386 {
1389 "update_covariant_transformation"));
1390
1391 for (unsigned int q = 0; q < output.size(); ++q)
1392 for (unsigned int i = 0; i < spacedim; ++i)
1393 for (unsigned int j = 0; j < spacedim; ++j)
1394 {
1395 double tmp[dim];
1396 for (unsigned int K = 0; K < dim; ++K)
1397 {
1398 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1399 for (unsigned int J = 1; J < dim; ++J)
1400 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1401 }
1402 for (unsigned int k = 0; k < spacedim; ++k)
1403 {
1404 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1405 for (unsigned int K = 1; K < dim; ++K)
1406 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1407 }
1408 }
1409 return;
1410 }
1411
1412 default:
1414 }
1415}
1416
1417
1418
1419template <int dim, int spacedim>
1420void
1422 const ArrayView<const Tensor<3, dim>> &input,
1423 const MappingKind mapping_kind,
1424 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1425 const ArrayView<Tensor<3, spacedim>> &output) const
1426{
1427 switch (mapping_kind)
1428 {
1432 internal::MappingManifoldImplementation::transform_hessians(
1433 input, mapping_kind, mapping_data, output);
1434 return;
1435 default:
1437 }
1438}
1439
1440//--------------------------- Explicit instantiations -----------------------
1441#include "mapping_manifold.inst"
1442
1443
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
DerivativeForm< 1, spacedim, dim, Number > transpose() const
void store_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
std::vector< double > volume_elements
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > vertices
std::vector< std::array< double, GeometryInfo< dim >::vertices_per_cell > > cell_manifold_quadrature_weights
std::vector< std::vector< Tensor< 1, spacedim > > > aux
virtual std::size_t memory_consumption() const override
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
ObserverPointer< const Manifold< dim, spacedim > > manifold
std::array< double, GeometryInfo< dim >::vertices_per_cell > vertex_weights
void compute_manifold_quadrature_weights(const Quadrature< dim > &quadrature)
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
Triangulation< dim, spacedim >::cell_iterator cell
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
MappingManifold()=default
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(const ArrayView< const Tensor< 1, dim > > &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim > > &output) const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
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
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
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
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
static constexpr Point< dim, Number > unit_vector(const unsigned int i)
static DataSetDescriptor cell()
Definition qprojector.h:461
const Point< dim > & point(const unsigned int i) const
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition collection.h:315
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
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
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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_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_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
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()
std::vector< index_type > data
Definition mpi.cc:735
std::size_t size
Definition mpi.cc:734
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim > > &output)
void do_fill_fe_face_values(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const typename QProjector< dim >::DataSetDescriptor data_set, const Quadrature< dim - 1 > &quadrature, const typename ::MappingQ< dim, spacedim >::InternalData &data, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void transform_fields(const ArrayView< const Tensor< rank, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim > > &output)
void transform_gradients(const ArrayView< const Tensor< rank, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank, spacedim > > &output)
void transform_hessians(const ArrayView< const Tensor< 3, dim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< 3, spacedim > > &output)
void maybe_compute_face_data(const ::MappingQ< dim, spacedim > &mapping, const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int n_q_points, const std::vector< double > &weights, const typename ::MappingQ< dim, spacedim >::InternalData &data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)