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_manifold.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2023 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
16
24
26
27#include <deal.II/fe/fe.h>
28#include <deal.II/fe/fe_tools.h>
31
33#include <deal.II/grid/tria.h>
35
37
38#include <algorithm>
39#include <array>
40#include <cmath>
41#include <memory>
42#include <numeric>
43
44
46
47template <int dim, int spacedim>
48std::size_t
50{
51 return (
64}
65
66
67template <int dim, int spacedim>
68void
70 const UpdateFlags update_flags,
71 const Quadrature<dim> &q,
72 const unsigned int n_original_q_points)
73{
74 // store the flags in the internal data object so we can access them
75 // in fill_fe_*_values()
76 this->update_each = update_flags;
77
78 // Store the quadrature
79 this->quad = q;
80
81 // Resize the weights
82 this->vertex_weights.resize(GeometryInfo<dim>::vertices_per_cell);
83
84 // see if we need the (transformation) shape function values
85 // and/or gradients and resize the necessary arrays
86 if (this->update_each &
88 compute_manifold_quadrature_weights(q);
89
90 if (this->update_each & update_covariant_transformation)
91 covariant.resize(n_original_q_points);
92
93 if (this->update_each & update_contravariant_transformation)
94 contravariant.resize(n_original_q_points);
95
96 if (this->update_each & update_volume_elements)
97 volume_elements.resize(n_original_q_points);
98}
99
100
101template <int dim, int spacedim>
102void
104 const UpdateFlags update_flags,
105 const Quadrature<dim> &q,
106 const unsigned int n_original_q_points)
107{
108 initialize(update_flags, q, n_original_q_points);
109
110 if (dim > 1)
111 {
112 if (this->update_each & update_boundary_forms)
113 {
114 aux.resize(dim - 1,
115 std::vector<Tensor<1, spacedim>>(n_original_q_points));
116
117 // Compute tangentials to the unit cell.
118 for (const unsigned int i : GeometryInfo<dim>::face_indices())
119 {
120 unit_tangentials[i].resize(n_original_q_points);
121 std::fill(unit_tangentials[i].begin(),
122 unit_tangentials[i].end(),
124 if (dim > 2)
125 {
126 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
127 .resize(n_original_q_points);
128 std::fill(
129 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
130 .begin(),
131 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
132 .end(),
134 }
135 }
136 }
137 }
138}
139
140
141
142template <int dim, int spacedim>
145{}
146
147
148
149template <int dim, int spacedim>
150std::unique_ptr<Mapping<dim, spacedim>>
152{
153 return std::make_unique<MappingManifold<dim, spacedim>>(*this);
154}
155
156
157
158template <int dim, int spacedim>
162 const Point<spacedim> &) const
163{
164 Assert(false, ExcNotImplemented());
165 return {};
166}
167
168
169
170template <int dim, int spacedim>
174 const Point<dim> & p) const
175{
176 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
177 std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;
178
179 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
180 {
181 vertices[v] = cell->vertex(v);
183 }
184 return cell->get_manifold().get_new_point(
185 make_array_view(vertices.begin(), vertices.end()),
186 make_array_view(weights.begin(), weights.end()));
187}
188
189
190
191// In the code below, GCC tries to instantiate MappingManifold<3,4> when
192// seeing which of the overloaded versions of
193// do_transform_real_to_unit_cell_internal() to call. This leads to bad
194// error messages and, generally, nothing very good. Avoid this by ensuring
195// that this class exists, but does not have an inner InternalData
196// type, thereby ruling out the codim-1 version of the function
197// below when doing overload resolution.
198template <>
200{};
201
202
203
204template <int dim, int spacedim>
207 const UpdateFlags in) const
208{
209 // add flags if the respective quantities are necessary to compute
210 // what we need. note that some flags appear in both the conditions
211 // and in subsequent set operations. this leads to some circular
212 // logic. the only way to treat this is to iterate. since there are
213 // 5 if-clauses in the loop, it will take at most 5 iterations to
214 // converge. do them:
215 UpdateFlags out = in;
216 for (unsigned int i = 0; i < 5; ++i)
217 {
218 // The following is a little incorrect:
219 // If not applied on a face,
220 // update_boundary_forms does not
221 // make sense. On the other hand,
222 // it is necessary on a
223 // face. Currently,
224 // update_boundary_forms is simply
225 // ignored for the interior of a
226 // cell.
229
234
235 if (out &
240
241 // The contravariant transformation used in the Piola
242 // transformation, which requires the determinant of the Jacobi
243 // matrix of the transformation. Because we have no way of
244 // knowing here whether the finite elements wants to use the
245 // contravariant of the Piola transforms, we add the JxW values
246 // to the list of flags to be updated for each cell.
248 out |= update_JxW_values;
249
250 if (out & update_normal_vectors)
251 out |= update_JxW_values;
252 }
253
254 // Now throw an exception if we stumble upon something that was not
255 // implemented yet
260
261 return out;
262}
263
264
265
266template <int dim, int spacedim>
267std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
269 const Quadrature<dim> &q) const
270{
271 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
272 std::make_unique<InternalData>();
273 auto &data = dynamic_cast<InternalData &>(*data_ptr);
274 data.initialize(this->requires_update_flags(update_flags), q, q.size());
275
276 return data_ptr;
277}
278
279
280
281template <int dim, int spacedim>
282std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
284 const UpdateFlags update_flags,
285 const hp::QCollection<dim - 1> &quadrature) const
286{
287 AssertDimension(quadrature.size(), 1);
288
289 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
290 std::make_unique<InternalData>();
291 auto &data = dynamic_cast<InternalData &>(*data_ptr);
292 data.initialize_face(this->requires_update_flags(update_flags),
294 ReferenceCells::get_hypercube<dim>(), quadrature[0]),
295 quadrature[0].size());
296
297 return data_ptr;
298}
299
300
301
302template <int dim, int spacedim>
303std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
305 const UpdateFlags update_flags,
306 const Quadrature<dim - 1> &quadrature) const
307{
308 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
309 std::make_unique<InternalData>();
310 auto &data = dynamic_cast<InternalData &>(*data_ptr);
311 data.initialize_face(this->requires_update_flags(update_flags),
313 ReferenceCells::get_hypercube<dim>(), quadrature),
314 quadrature.size());
315
316 return data_ptr;
317}
318
319
320
321namespace internal
322{
323 namespace MappingManifoldImplementation
324 {
325 namespace
326 {
333 template <int dim, int spacedim>
334 void
335 maybe_compute_q_points(
336 const typename QProjector<dim>::DataSetDescriptor data_set,
337 const typename ::MappingManifold<dim, spacedim>::InternalData
338 & data,
339 std::vector<Point<spacedim>> &quadrature_points)
340 {
341 const UpdateFlags update_flags = data.update_each;
342
343 AssertDimension(data.vertices.size(),
345
346 if (update_flags & update_quadrature_points)
347 {
348 for (unsigned int point = 0; point < quadrature_points.size();
349 ++point)
350 {
351 quadrature_points[point] = data.manifold->get_new_point(
352 make_array_view(data.vertices),
354 data.cell_manifold_quadrature_weights[point + data_set]));
355 }
356 }
357 }
358
359
360
367 template <int dim, int spacedim>
368 void
369 maybe_update_Jacobians(
370 const typename ::QProjector<dim>::DataSetDescriptor data_set,
371 const typename ::MappingManifold<dim, spacedim>::InternalData
372 &data)
373 {
374 const UpdateFlags update_flags = data.update_each;
375
376 if (update_flags & update_contravariant_transformation)
377 {
378 const unsigned int n_q_points = data.contravariant.size();
379
380 std::fill(data.contravariant.begin(),
381 data.contravariant.end(),
383
385 data.vertices.size());
386 for (unsigned int point = 0; point < n_q_points; ++point)
387 {
388 // Start by figuring out how to compute the direction in
389 // the reference space:
390 const Point<dim> &p = data.quad.point(point + data_set);
391
392 // And get its image on the manifold:
393 const Point<spacedim> P = data.manifold->get_new_point(
394 make_array_view(data.vertices),
396 data.cell_manifold_quadrature_weights[point + data_set]));
397
398 // To compute the Jacobian, we choose dim points aligned
399 // with the dim reference axes, which are still in the
400 // given cell, and ask for the tangent vector in these
401 // directions. Choosing the points is somewhat arbitrary,
402 // so we try to be smart and we pick points which are
403 // on the opposite quadrant w.r.t. the evaluation
404 // point.
405 for (unsigned int i = 0; i < dim; ++i)
406 {
408 const double pi = p[i];
409 Assert(pi >= 0 && pi <= 1.0,
411 "Was expecting a quadrature point "
412 "inside the unit reference element."));
413
414 // In the length L, we store also the direction sign,
415 // which is positive, if the coordinate is < .5,
416 const double L = pi > .5 ? -pi : 1 - pi;
417
418 const Point<dim> np(p + L * ei);
419
420 // Get the weights to compute the np point in real space
421 for (const unsigned int j :
423 data.vertex_weights[j] =
425
426 const Point<spacedim> NP = data.manifold->get_new_point(
427 make_array_view(data.vertices),
428 make_array_view(data.vertex_weights));
429
430 const Tensor<1, spacedim> T =
431 data.manifold->get_tangent_vector(P, NP);
432
433 for (unsigned int d = 0; d < spacedim; ++d)
434 data.contravariant[point][d][i] = T[d] / L;
435 }
436 }
437
438 if (update_flags & update_covariant_transformation)
439 {
440 const unsigned int n_q_points = data.contravariant.size();
441 for (unsigned int point = 0; point < n_q_points; ++point)
442 {
443 data.covariant[point] =
444 (data.contravariant[point]).covariant_form();
445 }
446 }
447
448 if (update_flags & update_volume_elements)
449 {
450 const unsigned int n_q_points = data.contravariant.size();
451 for (unsigned int point = 0; point < n_q_points; ++point)
452 data.volume_elements[point] =
453 data.contravariant[point].determinant();
454 }
455 }
456 }
457 } // namespace
458 } // namespace MappingManifoldImplementation
459} // namespace internal
460
461
462
463template <int dim, int spacedim>
468 const Quadrature<dim> & quadrature,
469 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
471 &output_data) const
472{
473 // ensure that the following static_cast is really correct:
474 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
476 const InternalData &data = static_cast<const InternalData &>(internal_data);
477
478 const unsigned int n_q_points = quadrature.size();
479
480 data.store_vertices(cell);
481 data.manifold = &(cell->get_manifold());
482
483 internal::MappingManifoldImplementation::maybe_compute_q_points<dim,
484 spacedim>(
486 data,
487 output_data.quadrature_points);
488
489 internal::MappingManifoldImplementation::maybe_update_Jacobians<dim,
490 spacedim>(
492
493 const UpdateFlags update_flags = data.update_each;
494 const std::vector<double> &weights = quadrature.get_weights();
495
496 // Multiply quadrature weights by absolute value of Jacobian determinants or
497 // the area element g=sqrt(DX^t DX) in case of codim > 0
498
499 if (update_flags & (update_normal_vectors | update_JxW_values))
500 {
501 AssertDimension(output_data.JxW_values.size(), n_q_points);
502
503 Assert(!(update_flags & update_normal_vectors) ||
504 (output_data.normal_vectors.size() == n_q_points),
505 ExcDimensionMismatch(output_data.normal_vectors.size(),
506 n_q_points));
507
508
509 for (unsigned int point = 0; point < n_q_points; ++point)
510 {
511 if (dim == spacedim)
512 {
513 const double det = data.contravariant[point].determinant();
514
515 // check for distorted cells.
516
517 // TODO: this allows for anisotropies of up to 1e6 in 3d and
518 // 1e12 in 2d. might want to find a finer
519 // (dimension-independent) criterion
520 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
521 cell->diameter() / std::sqrt(double(dim))),
523 cell->center(), det, point)));
524
525 output_data.JxW_values[point] = weights[point] * det;
526 }
527 // if dim==spacedim, then there is no cell normal to
528 // compute. since this is for FEValues (and not FEFaceValues),
529 // there are also no face normals to compute
530 else // codim>0 case
531 {
532 Tensor<1, spacedim> DX_t[dim];
533 for (unsigned int i = 0; i < spacedim; ++i)
534 for (unsigned int j = 0; j < dim; ++j)
535 DX_t[j][i] = data.contravariant[point][i][j];
536
537 Tensor<2, dim> G; // First fundamental form
538 for (unsigned int i = 0; i < dim; ++i)
539 for (unsigned int j = 0; j < dim; ++j)
540 G[i][j] = DX_t[i] * DX_t[j];
541
542 output_data.JxW_values[point] =
543 std::sqrt(determinant(G)) * weights[point];
544
545 if (update_flags & update_normal_vectors)
546 {
547 Assert(spacedim == dim + 1,
549 "There is no (unique) cell normal for " +
551 "-dimensional cells in " +
552 Utilities::int_to_string(spacedim) +
553 "-dimensional space. This only works if the "
554 "space dimension is one greater than the "
555 "dimensionality of the mesh cells."));
556
557 if (dim == 1)
558 output_data.normal_vectors[point] =
559 cross_product_2d(-DX_t[0]);
560 else // dim == 2
561 output_data.normal_vectors[point] =
562 cross_product_3d(DX_t[0], DX_t[1]);
563
564 output_data.normal_vectors[point] /=
565 output_data.normal_vectors[point].norm();
566
567 if (cell->direction_flag() == false)
568 output_data.normal_vectors[point] *= -1.;
569 }
570 } // codim>0 case
571 }
572 }
573
574
575
576 // copy values from InternalData to vector given by reference
577 if (update_flags & update_jacobians)
578 {
579 AssertDimension(output_data.jacobians.size(), n_q_points);
580 for (unsigned int point = 0; point < n_q_points; ++point)
581 output_data.jacobians[point] = data.contravariant[point];
582 }
583
584 // copy values from InternalData to vector given by reference
585 if (update_flags & update_inverse_jacobians)
586 {
587 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
588 for (unsigned int point = 0; point < n_q_points; ++point)
589 output_data.inverse_jacobians[point] =
590 data.covariant[point].transpose();
591 }
592
594}
595
596
597
598namespace internal
599{
600 namespace MappingManifoldImplementation
601 {
602 namespace
603 {
613 template <int dim, int spacedim>
614 void
615 maybe_compute_face_data(
616 const ::MappingManifold<dim, spacedim> &mapping,
617 const typename ::Triangulation<dim, spacedim>::cell_iterator
618 & cell,
619 const unsigned int face_no,
620 const unsigned int subface_no,
621 const unsigned int n_q_points,
622 const std::vector<double> &weights,
623 const typename ::MappingManifold<dim, spacedim>::InternalData
624 &data,
626 &output_data)
627 {
628 const UpdateFlags update_flags = data.update_each;
629
630 if (update_flags & update_boundary_forms)
631 {
632 AssertDimension(output_data.boundary_forms.size(), n_q_points);
633 if (update_flags & update_normal_vectors)
634 AssertDimension(output_data.normal_vectors.size(), n_q_points);
635 if (update_flags & update_JxW_values)
636 AssertDimension(output_data.JxW_values.size(), n_q_points);
637
638 // map the unit tangentials to the real cell. checking for d!=dim-1
639 // eliminates compiler warnings regarding unsigned int expressions <
640 // 0.
641 for (unsigned int d = 0; d != dim - 1; ++d)
642 {
644 data.unit_tangentials.size(),
646 Assert(
647 data.aux[d].size() <=
648 data
649 .unit_tangentials[face_no +
651 .size(),
653
654 mapping.transform(
656 data
657 .unit_tangentials[face_no +
660 data,
661 make_array_view(data.aux[d]));
662 }
663
664 // if dim==spacedim, we can use the unit tangentials to compute the
665 // boundary form by simply taking the cross product
666 if (dim == spacedim)
667 {
668 for (unsigned int i = 0; i < n_q_points; ++i)
669 switch (dim)
670 {
671 case 1:
672 // in 1d, we don't have access to any of the data.aux
673 // fields (because it has only dim-1 components), but we
674 // can still compute the boundary form by simply
675 // looking at the number of the face
676 output_data.boundary_forms[i][0] =
677 (face_no == 0 ? -1 : +1);
678 break;
679 case 2:
680 output_data.boundary_forms[i] =
681 cross_product_2d(data.aux[0][i]);
682 break;
683 case 3:
684 output_data.boundary_forms[i] =
685 cross_product_3d(data.aux[0][i], data.aux[1][i]);
686 break;
687 default:
688 Assert(false, ExcNotImplemented());
689 }
690 }
691 else //(dim < spacedim)
692 {
693 // in the codim-one case, the boundary form results from the
694 // cross product of all the face tangential vectors and the cell
695 // normal vector
696 //
697 // to compute the cell normal, use the same method used in
698 // fill_fe_values for cells above
699 AssertDimension(data.contravariant.size(), n_q_points);
700
701 for (unsigned int point = 0; point < n_q_points; ++point)
702 {
703 switch (dim)
704 {
705 case 1:
706 {
707 // J is a tangent vector
708 output_data.boundary_forms[point] =
709 data.contravariant[point].transpose()[0];
710 output_data.boundary_forms[point] /=
711 (face_no == 0 ? -1. : +1.) *
712 output_data.boundary_forms[point].norm();
713
714 break;
715 }
716
717 case 2:
718 {
720 data.contravariant[point].transpose();
721
722 Tensor<1, spacedim> cell_normal =
723 cross_product_3d(DX_t[0], DX_t[1]);
724 cell_normal /= cell_normal.norm();
725
726 // then compute the face normal from the face
727 // tangent and the cell normal:
728 output_data.boundary_forms[point] =
729 cross_product_3d(data.aux[0][point], cell_normal);
730
731 break;
732 }
733
734 default:
735 Assert(false, ExcNotImplemented());
736 }
737 }
738 }
739
740 if (update_flags & (update_normal_vectors | update_JxW_values))
741 for (unsigned int i = 0; i < output_data.boundary_forms.size();
742 ++i)
743 {
744 if (update_flags & update_JxW_values)
745 {
746 output_data.JxW_values[i] =
747 output_data.boundary_forms[i].norm() * weights[i];
748
749 if (subface_no != numbers::invalid_unsigned_int)
750 {
751 const double area_ratio =
753 cell->subface_case(face_no), subface_no);
754 output_data.JxW_values[i] *= area_ratio;
755 }
756 }
757
758 if (update_flags & update_normal_vectors)
759 output_data.normal_vectors[i] =
760 Point<spacedim>(output_data.boundary_forms[i] /
761 output_data.boundary_forms[i].norm());
762 }
763
764 if (update_flags & update_jacobians)
765 for (unsigned int point = 0; point < n_q_points; ++point)
766 output_data.jacobians[point] = data.contravariant[point];
767
768 if (update_flags & update_inverse_jacobians)
769 for (unsigned int point = 0; point < n_q_points; ++point)
770 output_data.inverse_jacobians[point] =
771 data.covariant[point].transpose();
772 }
773 }
774
775
782 template <int dim, int spacedim>
783 void
785 const ::MappingManifold<dim, spacedim> &mapping,
786 const typename ::Triangulation<dim, spacedim>::cell_iterator
787 & cell,
788 const unsigned int face_no,
789 const unsigned int subface_no,
790 const typename QProjector<dim>::DataSetDescriptor data_set,
791 const Quadrature<dim - 1> & quadrature,
792 const typename ::MappingManifold<dim, spacedim>::InternalData
793 &data,
795 &output_data)
796 {
797 data.store_vertices(cell);
798
799 data.manifold = &cell->face(face_no)->get_manifold();
800
801 maybe_compute_q_points<dim, spacedim>(data_set,
802 data,
803 output_data.quadrature_points);
804 maybe_update_Jacobians<dim, spacedim>(data_set, data);
805
807 cell,
808 face_no,
809 subface_no,
810 quadrature.size(),
811 quadrature.get_weights(),
812 data,
813 output_data);
814 }
815
816 template <int dim, int spacedim, int rank>
817 void
819 const ArrayView<const Tensor<rank, dim>> & input,
820 const MappingKind mapping_kind,
821 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
822 const ArrayView<Tensor<rank, spacedim>> & output)
823 {
824 AssertDimension(input.size(), output.size());
825 Assert((dynamic_cast<const typename ::
826 MappingManifold<dim, spacedim>::InternalData *>(
827 &mapping_data) != nullptr),
829 const typename ::MappingManifold<dim, spacedim>::InternalData
830 &data =
831 static_cast<const typename ::MappingManifold<dim, spacedim>::
832 InternalData &>(mapping_data);
833
834 switch (mapping_kind)
835 {
837 {
838 Assert(
839 data.update_each & update_contravariant_transformation,
841 "update_contravariant_transformation"));
842
843 for (unsigned int i = 0; i < output.size(); ++i)
844 output[i] =
845 apply_transformation(data.contravariant[i], input[i]);
846
847 return;
848 }
849
850 case mapping_piola:
851 {
852 Assert(
853 data.update_each & update_contravariant_transformation,
855 "update_contravariant_transformation"));
856 Assert(
857 data.update_each & update_volume_elements,
859 "update_volume_elements"));
860 Assert(rank == 1, ExcMessage("Only for rank 1"));
861 if (rank != 1)
862 return;
863
864 for (unsigned int i = 0; i < output.size(); ++i)
865 {
866 output[i] =
867 apply_transformation(data.contravariant[i], input[i]);
868 output[i] /= data.volume_elements[i];
869 }
870 return;
871 }
872 // We still allow this operation as in the
873 // reference cell Derivatives are Tensor
874 // rather than DerivativeForm
876 {
877 Assert(
878 data.update_each & update_contravariant_transformation,
880 "update_covariant_transformation"));
881
882 for (unsigned int i = 0; i < output.size(); ++i)
883 output[i] = apply_transformation(data.covariant[i], input[i]);
884
885 return;
886 }
887
888 default:
889 Assert(false, ExcNotImplemented());
890 }
891 }
892
893
894 template <int dim, int spacedim, int rank>
895 void
897 const ArrayView<const Tensor<rank, dim>> & input,
898 const MappingKind mapping_kind,
899 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
900 const ArrayView<Tensor<rank, spacedim>> & output)
901 {
902 AssertDimension(input.size(), output.size());
903 Assert((dynamic_cast<const typename ::
904 MappingManifold<dim, spacedim>::InternalData *>(
905 &mapping_data) != nullptr),
907 const typename ::MappingManifold<dim, spacedim>::InternalData
908 &data =
909 static_cast<const typename ::MappingManifold<dim, spacedim>::
910 InternalData &>(mapping_data);
911
912 switch (mapping_kind)
913 {
915 {
916 Assert(
917 data.update_each & update_covariant_transformation,
919 "update_covariant_transformation"));
920 Assert(
921 data.update_each & update_contravariant_transformation,
923 "update_contravariant_transformation"));
924 Assert(rank == 2, ExcMessage("Only for rank 2"));
925
926 for (unsigned int i = 0; i < output.size(); ++i)
927 {
929 apply_transformation(data.contravariant[i],
930 transpose(input[i]));
931 output[i] =
932 apply_transformation(data.covariant[i], A.transpose());
933 }
934
935 return;
936 }
937
939 {
940 Assert(
941 data.update_each & update_covariant_transformation,
943 "update_covariant_transformation"));
944 Assert(rank == 2, ExcMessage("Only for rank 2"));
945
946 for (unsigned int i = 0; i < output.size(); ++i)
947 {
949 apply_transformation(data.covariant[i],
950 transpose(input[i]));
951 output[i] =
952 apply_transformation(data.covariant[i], A.transpose());
953 }
954
955 return;
956 }
957
959 {
960 Assert(
961 data.update_each & update_covariant_transformation,
963 "update_covariant_transformation"));
964 Assert(
965 data.update_each & update_contravariant_transformation,
967 "update_contravariant_transformation"));
968 Assert(
969 data.update_each & update_volume_elements,
971 "update_volume_elements"));
972 Assert(rank == 2, ExcMessage("Only for rank 2"));
973
974 for (unsigned int i = 0; i < output.size(); ++i)
975 {
977 apply_transformation(data.covariant[i], input[i]);
979 apply_transformation(data.contravariant[i],
980 A.transpose());
981
982 output[i] = transpose(T);
983 output[i] /= data.volume_elements[i];
984 }
985
986 return;
987 }
988
989 default:
990 Assert(false, ExcNotImplemented());
991 }
992 }
993
994
995
996 template <int dim, int spacedim>
997 void
999 const ArrayView<const Tensor<3, dim>> & input,
1000 const MappingKind mapping_kind,
1001 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1002 const ArrayView<Tensor<3, spacedim>> & output)
1003 {
1004 AssertDimension(input.size(), output.size());
1005 Assert((dynamic_cast<const typename ::
1006 MappingManifold<dim, spacedim>::InternalData *>(
1007 &mapping_data) != nullptr),
1009 const typename ::MappingManifold<dim, spacedim>::InternalData
1010 &data =
1011 static_cast<const typename ::MappingManifold<dim, spacedim>::
1012 InternalData &>(mapping_data);
1013
1014 switch (mapping_kind)
1015 {
1017 {
1018 Assert(
1019 data.update_each & update_covariant_transformation,
1021 "update_covariant_transformation"));
1022 Assert(
1023 data.update_each & update_contravariant_transformation,
1025 "update_contravariant_transformation"));
1026
1027 for (unsigned int q = 0; q < output.size(); ++q)
1028 for (unsigned int i = 0; i < spacedim; ++i)
1029 {
1030 double tmp1[dim][dim];
1031 for (unsigned int J = 0; J < dim; ++J)
1032 for (unsigned int K = 0; K < dim; ++K)
1033 {
1034 tmp1[J][K] =
1035 data.contravariant[q][i][0] * input[q][0][J][K];
1036 for (unsigned int I = 1; I < dim; ++I)
1037 tmp1[J][K] +=
1038 data.contravariant[q][i][I] * input[q][I][J][K];
1039 }
1040 for (unsigned int j = 0; j < spacedim; ++j)
1041 {
1042 double tmp2[dim];
1043 for (unsigned int K = 0; K < dim; ++K)
1044 {
1045 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1046 for (unsigned int J = 1; J < dim; ++J)
1047 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1048 }
1049 for (unsigned int k = 0; k < spacedim; ++k)
1050 {
1051 output[q][i][j][k] =
1052 data.covariant[q][k][0] * tmp2[0];
1053 for (unsigned int K = 1; K < dim; ++K)
1054 output[q][i][j][k] +=
1055 data.covariant[q][k][K] * tmp2[K];
1056 }
1057 }
1058 }
1059 return;
1060 }
1061
1063 {
1064 Assert(
1065 data.update_each & update_covariant_transformation,
1067 "update_covariant_transformation"));
1068
1069 for (unsigned int q = 0; q < output.size(); ++q)
1070 for (unsigned int i = 0; i < spacedim; ++i)
1071 {
1072 double tmp1[dim][dim];
1073 for (unsigned int J = 0; J < dim; ++J)
1074 for (unsigned int K = 0; K < dim; ++K)
1075 {
1076 tmp1[J][K] =
1077 data.covariant[q][i][0] * input[q][0][J][K];
1078 for (unsigned int I = 1; I < dim; ++I)
1079 tmp1[J][K] +=
1080 data.covariant[q][i][I] * input[q][I][J][K];
1081 }
1082 for (unsigned int j = 0; j < spacedim; ++j)
1083 {
1084 double tmp2[dim];
1085 for (unsigned int K = 0; K < dim; ++K)
1086 {
1087 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1088 for (unsigned int J = 1; J < dim; ++J)
1089 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1090 }
1091 for (unsigned int k = 0; k < spacedim; ++k)
1092 {
1093 output[q][i][j][k] =
1094 data.covariant[q][k][0] * tmp2[0];
1095 for (unsigned int K = 1; K < dim; ++K)
1096 output[q][i][j][k] +=
1097 data.covariant[q][k][K] * tmp2[K];
1098 }
1099 }
1100 }
1101
1102 return;
1103 }
1104
1106 {
1107 Assert(
1108 data.update_each & update_covariant_transformation,
1110 "update_covariant_transformation"));
1111 Assert(
1112 data.update_each & update_contravariant_transformation,
1114 "update_contravariant_transformation"));
1115 Assert(
1116 data.update_each & update_volume_elements,
1118 "update_volume_elements"));
1119
1120 for (unsigned int q = 0; q < output.size(); ++q)
1121 for (unsigned int i = 0; i < spacedim; ++i)
1122 {
1123 double factor[dim];
1124 for (unsigned int I = 0; I < dim; ++I)
1125 factor[I] =
1126 data.contravariant[q][i][I] / data.volume_elements[q];
1127 double tmp1[dim][dim];
1128 for (unsigned int J = 0; J < dim; ++J)
1129 for (unsigned int K = 0; K < dim; ++K)
1130 {
1131 tmp1[J][K] = factor[0] * input[q][0][J][K];
1132 for (unsigned int I = 1; I < dim; ++I)
1133 tmp1[J][K] += factor[I] * input[q][I][J][K];
1134 }
1135 for (unsigned int j = 0; j < spacedim; ++j)
1136 {
1137 double tmp2[dim];
1138 for (unsigned int K = 0; K < dim; ++K)
1139 {
1140 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1141 for (unsigned int J = 1; J < dim; ++J)
1142 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1143 }
1144 for (unsigned int k = 0; k < spacedim; ++k)
1145 {
1146 output[q][i][j][k] =
1147 data.covariant[q][k][0] * tmp2[0];
1148 for (unsigned int K = 1; K < dim; ++K)
1149 output[q][i][j][k] +=
1150 data.covariant[q][k][K] * tmp2[K];
1151 }
1152 }
1153 }
1154
1155 return;
1156 }
1157
1158 default:
1159 Assert(false, ExcNotImplemented());
1160 }
1161 }
1162
1163
1164
1165 template <int dim, int spacedim, int rank>
1166 void
1169 const MappingKind mapping_kind,
1170 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1171 const ArrayView<Tensor<rank + 1, spacedim>> & output)
1172 {
1173 AssertDimension(input.size(), output.size());
1174 Assert((dynamic_cast<const typename ::
1175 MappingManifold<dim, spacedim>::InternalData *>(
1176 &mapping_data) != nullptr),
1178 const typename ::MappingManifold<dim, spacedim>::InternalData
1179 &data =
1180 static_cast<const typename ::MappingManifold<dim, spacedim>::
1181 InternalData &>(mapping_data);
1182
1183 switch (mapping_kind)
1184 {
1185 case mapping_covariant:
1186 {
1187 Assert(
1188 data.update_each & update_contravariant_transformation,
1190 "update_covariant_transformation"));
1191
1192 for (unsigned int i = 0; i < output.size(); ++i)
1193 output[i] = apply_transformation(data.covariant[i], input[i]);
1194
1195 return;
1196 }
1197 default:
1198 Assert(false, ExcNotImplemented());
1199 }
1200 }
1201 } // namespace
1202 } // namespace MappingManifoldImplementation
1203} // namespace internal
1204
1205
1206
1207template <int dim, int spacedim>
1208void
1211 const unsigned int face_no,
1212 const hp::QCollection<dim - 1> & quadrature,
1213 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1215 &output_data) const
1216{
1217 AssertDimension(quadrature.size(), 1);
1218
1219 // ensure that the following cast is really correct:
1220 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1222 const InternalData &data = static_cast<const InternalData &>(internal_data);
1223
1224 internal::MappingManifoldImplementation::do_fill_fe_face_values(
1225 *this,
1226 cell,
1227 face_no,
1230 ReferenceCells::get_hypercube<dim>(),
1231 face_no,
1232 cell->face_orientation(face_no),
1233 cell->face_flip(face_no),
1234 cell->face_rotation(face_no),
1235 quadrature[0].size()),
1236 quadrature[0],
1237 data,
1238 output_data);
1239}
1240
1241
1242
1243template <int dim, int spacedim>
1244void
1247 const unsigned int face_no,
1248 const unsigned int subface_no,
1249 const Quadrature<dim - 1> & quadrature,
1250 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1252 &output_data) const
1253{
1254 // ensure that the following cast is really correct:
1255 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1257 const InternalData &data = static_cast<const InternalData &>(internal_data);
1258
1259 internal::MappingManifoldImplementation::do_fill_fe_face_values(
1260 *this,
1261 cell,
1262 face_no,
1263 subface_no,
1265 ReferenceCells::get_hypercube<dim>(),
1266 face_no,
1267 subface_no,
1268 cell->face_orientation(face_no),
1269 cell->face_flip(face_no),
1270 cell->face_rotation(face_no),
1271 quadrature.size(),
1272 cell->subface_case(face_no)),
1273 quadrature,
1274 data,
1275 output_data);
1276}
1277
1278
1279
1280template <int dim, int spacedim>
1281void
1283 const ArrayView<const Tensor<1, dim>> & input,
1284 const MappingKind mapping_kind,
1285 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1286 const ArrayView<Tensor<1, spacedim>> & output) const
1287{
1288 internal::MappingManifoldImplementation::transform_fields(input,
1289 mapping_kind,
1290 mapping_data,
1291 output);
1292}
1293
1294
1295
1296template <int dim, int spacedim>
1297void
1299 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1300 const MappingKind mapping_kind,
1301 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1302 const ArrayView<Tensor<2, spacedim>> & output) const
1303{
1304 internal::MappingManifoldImplementation::transform_differential_forms(
1305 input, mapping_kind, mapping_data, output);
1306}
1307
1308
1309
1310template <int dim, int spacedim>
1311void
1313 const ArrayView<const Tensor<2, dim>> & input,
1314 const MappingKind mapping_kind,
1315 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1316 const ArrayView<Tensor<2, spacedim>> & output) const
1317{
1318 switch (mapping_kind)
1319 {
1321 internal::MappingManifoldImplementation::transform_fields(input,
1322 mapping_kind,
1323 mapping_data,
1324 output);
1325 return;
1326
1330 internal::MappingManifoldImplementation::transform_gradients(
1331 input, mapping_kind, mapping_data, output);
1332 return;
1333 default:
1334 Assert(false, ExcNotImplemented());
1335 }
1336}
1337
1338
1339
1340template <int dim, int spacedim>
1341void
1343 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1344 const MappingKind mapping_kind,
1345 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1346 const ArrayView<Tensor<3, spacedim>> & output) const
1347{
1348 AssertDimension(input.size(), output.size());
1349 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1351 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1352
1353 switch (mapping_kind)
1354 {
1356 {
1359 "update_covariant_transformation"));
1360
1361 for (unsigned int q = 0; q < output.size(); ++q)
1362 for (unsigned int i = 0; i < spacedim; ++i)
1363 for (unsigned int j = 0; j < spacedim; ++j)
1364 {
1365 double tmp[dim];
1366 for (unsigned int K = 0; K < dim; ++K)
1367 {
1368 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1369 for (unsigned int J = 1; J < dim; ++J)
1370 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1371 }
1372 for (unsigned int k = 0; k < spacedim; ++k)
1373 {
1374 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1375 for (unsigned int K = 1; K < dim; ++K)
1376 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1377 }
1378 }
1379 return;
1380 }
1381
1382 default:
1383 Assert(false, ExcNotImplemented());
1384 }
1385}
1386
1387
1388
1389template <int dim, int spacedim>
1390void
1392 const ArrayView<const Tensor<3, dim>> & input,
1393 const MappingKind mapping_kind,
1394 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1395 const ArrayView<Tensor<3, spacedim>> & output) const
1396{
1397 switch (mapping_kind)
1398 {
1402 internal::MappingManifoldImplementation::transform_hessians(
1403 input, mapping_kind, mapping_data, output);
1404 return;
1405 default:
1406 Assert(false, ExcNotImplemented());
1407 }
1408}
1409
1410//--------------------------- Explicit instantiations -----------------------
1411#include "mapping_manifold.inst"
1412
1413
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:704
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
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
SmartPointer< const Manifold< dim, spacedim > > manifold
std::vector< std::vector< Tensor< 1, spacedim > > > aux
virtual std::size_t memory_consumption() const override
std::vector< std::vector< double > > cell_manifold_quadrature_weights
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< Point< spacedim > > vertices
std::vector< double > vertex_weights
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:317
Definition point.h:112
static Point< dim, Number > unit_vector(const unsigned int i)
static DataSetDescriptor cell()
Definition qprojector.h:394
const std::vector< double > & get_weights() const
unsigned int size() const
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
Definition collection.h:265
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
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
Point< 3 > vertices[4]
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: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
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471
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:213
::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 > &)