Reference documentation for deal.II version GIT relicensing-31-gf283a0fc39 2024-02-29 06:50:02+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_q.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
22#include <deal.II/base/table.h>
24
25#include <deal.II/fe/fe_dgq.h>
26#include <deal.II/fe/fe_tools.h>
31
33#include <deal.II/grid/tria.h>
35
36#include <boost/container/small_vector.hpp>
37
38#include <algorithm>
39#include <array>
40#include <cmath>
41#include <limits>
42#include <memory>
43#include <numeric>
44
45
47
48
49template <int dim, int spacedim>
51 const unsigned int polynomial_degree)
52 : polynomial_degree(polynomial_degree)
53 , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
54 , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
55 , tensor_product_quadrature(false)
56 , output_data(nullptr)
57{}
58
59
60
61template <int dim, int spacedim>
62std::size_t
76
77
78
79template <int dim, int spacedim>
80void
82 const UpdateFlags update_flags,
83 const Quadrature<dim> &quadrature,
84 const unsigned int n_original_q_points)
85{
86 // store the flags in the internal data object so we can access them
87 // in fill_fe_*_values()
88 this->update_each = update_flags;
89
90 const unsigned int n_q_points = quadrature.size();
91
92 if (this->update_each & update_volume_elements)
93 volume_elements.resize(n_original_q_points);
94
95 tensor_product_quadrature = quadrature.is_tensor_product();
96
97 // use of MatrixFree only for higher order elements and with more than one
98 // point where tensor products do not make sense
99 if (polynomial_degree < 2 || n_q_points == 1)
100 tensor_product_quadrature = false;
101
102 if (dim > 1)
103 {
104 // find out if the one-dimensional formula is the same
105 // in all directions
106 if (tensor_product_quadrature)
107 {
108 const std::array<Quadrature<1>, dim> &quad_array =
109 quadrature.get_tensor_basis();
110 for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
111 {
112 if (quad_array[i - 1].size() != quad_array[i].size())
113 {
114 tensor_product_quadrature = false;
115 break;
116 }
117 else
118 {
119 const std::vector<Point<1>> &points_1 =
120 quad_array[i - 1].get_points();
121 const std::vector<Point<1>> &points_2 =
122 quad_array[i].get_points();
123 const std::vector<double> &weights_1 =
124 quad_array[i - 1].get_weights();
125 const std::vector<double> &weights_2 =
126 quad_array[i].get_weights();
127 for (unsigned int j = 0; j < quad_array[i].size(); ++j)
128 {
129 if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
130 std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
132 tensor_product_quadrature = false;
133 break;
134 }
136 }
137 }
138
139 if (tensor_product_quadrature)
140 {
141 // use a 1d FE_DGQ and adjust the hierarchic -> lexicographic
142 // numbering manually (building an FE_Q<dim> is relatively
143 // expensive due to constraints)
145 shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
146 shape_info.lexicographic_numbering =
147 FETools::lexicographic_to_hierarchic_numbering<dim>(
149 shape_info.n_q_points = n_q_points;
150 shape_info.dofs_per_component_on_cell =
152 }
154 }
155}
156
158
159template <int dim, int spacedim>
160void
162 const UpdateFlags update_flags,
163 const Quadrature<dim> &quadrature,
164 const unsigned int n_original_q_points)
165{
166 initialize(update_flags, quadrature, n_original_q_points);
167
168 quadrature_points = quadrature.get_points();
169
170 if (dim > 1 && tensor_product_quadrature)
171 {
172 constexpr unsigned int facedim = dim - 1;
174 shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
175 shape_info.lexicographic_numbering =
176 FETools::lexicographic_to_hierarchic_numbering<facedim>(
178 shape_info.n_q_points = n_original_q_points;
179 shape_info.dofs_per_component_on_cell =
181 }
182
183 if (dim > 1)
184 {
185 if (this->update_each &
187 {
188 aux.resize(dim - 1);
189 aux[0].resize(n_original_q_points);
190 if (dim > 2)
191 aux[1].resize(n_original_q_points);
192
193 // Compute tangentials to the unit cell.
194 for (const unsigned int i : GeometryInfo<dim>::face_indices())
195 {
196 unit_tangentials[i].resize(n_original_q_points);
197 std::fill(unit_tangentials[i].begin(),
198 unit_tangentials[i].end(),
200 if (dim > 2)
202 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
203 .resize(n_original_q_points);
204 std::fill(
205 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
206 .begin(),
207 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
208 .end(),
210 }
211 }
212 }
213 }
214}
216
217
218template <int dim, int spacedim>
222 QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
224 Polynomials::generate_complete_Lagrange_basis(line_support_points))
226 FETools::lexicographic_to_hierarchic_numbering<dim>(p))
228 internal::MappingQImplementation::unit_support_points<dim>(
232 internal::MappingQImplementation::
233 compute_support_point_weights_perimeter_to_interior(
234 this->polynomial_degree,
235 dim))
237 internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
238 this->polynomial_degree))
239{
240 Assert(p >= 1,
241 ExcMessage("It only makes sense to create polynomial mappings "
242 "with a polynomial degree greater or equal to one."));
243}
244
245
246
247template <int dim, int spacedim>
248MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
249 : MappingQ<dim, spacedim>(p)
250{}
251
252
253
254template <int dim, int spacedim>
256 : polynomial_degree(mapping.polynomial_degree)
257 , line_support_points(mapping.line_support_points)
258 , polynomials_1d(mapping.polynomials_1d)
259 , renumber_lexicographic_to_hierarchic(
260 mapping.renumber_lexicographic_to_hierarchic)
261 , unit_cell_support_points(mapping.unit_cell_support_points)
262 , support_point_weights_perimeter_to_interior(
263 mapping.support_point_weights_perimeter_to_interior)
264 , support_point_weights_cell(mapping.support_point_weights_cell)
265{}
266
267
268
269template <int dim, int spacedim>
270std::unique_ptr<Mapping<dim, spacedim>>
272{
273 return std::make_unique<MappingQ<dim, spacedim>>(*this);
274}
275
276
277
278template <int dim, int spacedim>
279unsigned int
281{
282 return polynomial_degree;
283}
284
285
286
287template <int dim, int spacedim>
291 const Point<dim> &p) const
292{
294 polynomials_1d,
295 this->compute_mapping_support_points(cell),
296 p,
297 polynomials_1d.size() == 2,
298 renumber_lexicographic_to_hierarchic));
299}
300
301
302// In the code below, GCC tries to instantiate MappingQ<3,4> when
303// seeing which of the overloaded versions of
304// do_transform_real_to_unit_cell_internal() to call. This leads to bad
305// error messages and, generally, nothing very good. Avoid this by ensuring
306// that this class exists, but does not have an inner InternalData
307// type, thereby ruling out the codim-1 version of the function
308// below when doing overload resolution.
309template <>
310class MappingQ<3, 4>
311{};
312
313
314
315// visual studio freaks out when trying to determine if
316// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
317// candidate. So instead of letting the compiler pick the correct overload, we
318// use template specialization to make sure we pick up the right function to
319// call:
320
321template <int dim, int spacedim>
325 const Point<spacedim> &,
326 const Point<dim> &) const
327{
328 // default implementation (should never be called)
330 return {};
332
333
334
335template <>
339 const Point<1> &p,
340 const Point<1> &initial_p_unit) const
341{
342 // dispatch to the various specializations for spacedim=dim,
343 // spacedim=dim+1, etc
344 return internal::MappingQImplementation::
345 do_transform_real_to_unit_cell_internal<1>(
346 p,
347 initial_p_unit,
348 this->compute_mapping_support_points(cell),
349 polynomials_1d,
350 renumber_lexicographic_to_hierarchic);
351}
352
353
354
355template <>
359 const Point<2> &p,
360 const Point<2> &initial_p_unit) const
361{
362 return internal::MappingQImplementation::
363 do_transform_real_to_unit_cell_internal<2>(
364 p,
365 initial_p_unit,
366 this->compute_mapping_support_points(cell),
367 polynomials_1d,
368 renumber_lexicographic_to_hierarchic);
369}
370
371
372
373template <>
377 const Point<3> &p,
378 const Point<3> &initial_p_unit) const
379{
380 return internal::MappingQImplementation::
381 do_transform_real_to_unit_cell_internal<3>(
382 p,
383 initial_p_unit,
384 this->compute_mapping_support_points(cell),
385 polynomials_1d,
386 renumber_lexicographic_to_hierarchic);
387}
388
389
390
391template <>
395 const Point<2> &p,
396 const Point<1> &initial_p_unit) const
397{
398 const int dim = 1;
399 const int spacedim = 2;
400
401 const Quadrature<dim> point_quadrature(initial_p_unit);
402
404 if (spacedim > dim)
405 update_flags |= update_jacobian_grads;
406 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
407 get_data(update_flags, point_quadrature));
408
409 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
410
411 // dispatch to the various specializations for spacedim=dim,
412 // spacedim=dim+1, etc
413 return internal::MappingQImplementation::
414 do_transform_real_to_unit_cell_internal_codim1<1>(
415 p,
416 initial_p_unit,
417 mdata->mapping_support_points,
418 polynomials_1d,
419 renumber_lexicographic_to_hierarchic);
420}
421
422
423
424template <>
428 const Point<3> &p,
429 const Point<2> &initial_p_unit) const
430{
431 const int dim = 2;
432 const int spacedim = 3;
433
434 const Quadrature<dim> point_quadrature(initial_p_unit);
435
437 if (spacedim > dim)
438 update_flags |= update_jacobian_grads;
439 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
440 get_data(update_flags, point_quadrature));
441
442 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
443
444 // dispatch to the various specializations for spacedim=dim,
445 // spacedim=dim+1, etc
446 return internal::MappingQImplementation::
447 do_transform_real_to_unit_cell_internal_codim1<2>(
448 p,
449 initial_p_unit,
450 mdata->mapping_support_points,
451 polynomials_1d,
452 renumber_lexicographic_to_hierarchic);
453}
454
455
456
457template <>
461 const Point<3> &,
462 const Point<1> &) const
463{
465 return {};
466}
467
469
470template <int dim, int spacedim>
474 const Point<spacedim> &p) const
475{
476 // Use an exact formula if one is available. this is only the case
477 // for Q1 mappings in 1d, and in 2d if dim==spacedim
478 if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
479 ((dim == 1) || ((dim == 2) && (dim == spacedim))))
480 {
481 // The dimension-dependent algorithms are much faster (about 25-45x in
482 // 2d) but fail most of the time when the given point (p) is not in the
483 // cell. The dimension-independent Newton algorithm given below is
484 // slower, but more robust (though it still sometimes fails). Therefore
485 // this function implements the following strategy based on the
486 // p's dimension:
487 //
488 // * In 1d this mapping is linear, so the mapping is always invertible
489 // (and the exact formula is known) as long as the cell has non-zero
490 // length.
491 // * In 2d the exact (quadratic) formula is called first. If either the
492 // exact formula does not succeed (negative discriminant in the
493 // quadratic formula) or succeeds but finds a solution outside of the
494 // unit cell, then the Newton solver is called. The rationale for the
495 // second choice is that the exact formula may provide two different
496 // answers when mapping a point outside of the real cell, but the
497 // Newton solver (if it converges) will only return one answer.
498 // Otherwise the exact formula successfully found a point in the unit
499 // cell and that value is returned.
500 // * In 3d there is no (known to the authors) exact formula, so the Newton
501 // algorithm is used.
502 const auto vertices_ = this->get_vertices(cell);
503
504 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
505 vertices;
506 for (unsigned int i = 0; i < vertices.size(); ++i)
507 vertices[i] = vertices_[i];
508
509 try
510 {
511 switch (dim)
512 {
513 case 1:
514 {
515 // formula not subject to any issues in 1d
516 if (spacedim == 1)
518 vertices, p);
519 else
520 break;
521 }
522
523 case 2:
524 {
525 const Point<dim> point =
527 p);
528
529 // formula not guaranteed to work for points outside of
530 // the cell. only take the computed point if it lies
531 // inside the reference cell
532 const double eps = 1e-15;
533 if (-eps <= point[1] && point[1] <= 1 + eps &&
534 -eps <= point[0] && point[0] <= 1 + eps)
535 {
536 return point;
537 }
538 else
539 break;
540 }
541
542 default:
543 {
544 // we should get here, based on the if-condition at the top
546 }
547 }
548 }
549 catch (
551 {
552 // simply fall through and continue on to the standard Newton code
553 }
554 }
555 else
556 {
557 // we can't use an explicit formula,
558 }
559
560
561 // Find the initial value for the Newton iteration by a normal
562 // projection to the least square plane determined by the vertices
563 // of the cell
564 Point<dim> initial_p_unit;
565 if (this->preserves_vertex_locations())
566 {
567 initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
568 // in 1d with spacedim > 1 the affine approximation is exact
569 if (dim == 1 && polynomial_degree == 1)
570 return initial_p_unit;
571 }
572 else
573 {
574 // else, we simply use the mid point
575 for (unsigned int d = 0; d < dim; ++d)
576 initial_p_unit[d] = 0.5;
577 }
578
579 // perform the Newton iteration and return the result. note that this
580 // statement may throw an exception, which we simply pass up to the caller
581 const Point<dim> p_unit =
582 this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
585 return p_unit;
586}
587
588
589
590template <int dim, int spacedim>
591void
594 const ArrayView<const Point<spacedim>> &real_points,
595 const ArrayView<Point<dim>> &unit_points) const
596{
597 // Go to base class functions for dim < spacedim because it is not yet
598 // implemented with optimized code.
599 if (dim < spacedim)
600 {
602 real_points,
603 unit_points);
604 return;
605 }
606
607 AssertDimension(real_points.size(), unit_points.size());
608 const std::vector<Point<spacedim>> support_points =
609 this->compute_mapping_support_points(cell);
610
611 // From the given (high-order) support points, now only pick the first
612 // 2^dim points and construct an affine approximation from those.
614 inverse_approximation(support_points, unit_cell_support_points);
615
616 const unsigned int n_points = real_points.size();
617 const unsigned int n_lanes = VectorizedArray<double>::size();
618
619 // Use the more heavy VectorizedArray code path if there is more than
620 // one point left to compute
621 for (unsigned int i = 0; i < n_points; i += n_lanes)
622 if (n_points - i > 1)
623 {
625 for (unsigned int j = 0; j < n_lanes; ++j)
626 if (i + j < n_points)
627 for (unsigned int d = 0; d < spacedim; ++d)
628 p_vec[d][j] = real_points[i + j][d];
629 else
630 for (unsigned int d = 0; d < spacedim; ++d)
631 p_vec[d][j] = real_points[i][d];
632
634 internal::MappingQImplementation::
635 do_transform_real_to_unit_cell_internal<dim, spacedim>(
636 p_vec,
637 inverse_approximation.compute(p_vec),
638 support_points,
639 polynomials_1d,
640 renumber_lexicographic_to_hierarchic);
641
642 // If the vectorized computation failed, it could be that only some of
643 // the lanes failed but others would have succeeded if we had let them
644 // compute alone without interference (like negative Jacobian
645 // determinants) from other SIMD lanes. Repeat the computation in this
646 // unlikely case with scalar arguments.
647 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
648 if (numbers::is_finite(unit_point[0][j]))
649 for (unsigned int d = 0; d < dim; ++d)
650 unit_points[i + j][d] = unit_point[d][j];
651 else
652 unit_points[i + j] = internal::MappingQImplementation::
653 do_transform_real_to_unit_cell_internal<dim, spacedim>(
654 real_points[i + j],
655 inverse_approximation.compute(real_points[i + j]),
656 support_points,
657 polynomials_1d,
658 renumber_lexicographic_to_hierarchic);
659 }
660 else
661 unit_points[i] = internal::MappingQImplementation::
662 do_transform_real_to_unit_cell_internal<dim, spacedim>(
663 real_points[i],
664 inverse_approximation.compute(real_points[i]),
665 support_points,
666 polynomials_1d,
667 renumber_lexicographic_to_hierarchic);
668}
669
670
671
672template <int dim, int spacedim>
675{
676 // add flags if the respective quantities are necessary to compute
677 // what we need. note that some flags appear in both the conditions
678 // and in subsequent set operations. this leads to some circular
679 // logic. the only way to treat this is to iterate. since there are
680 // 5 if-clauses in the loop, it will take at most 5 iterations to
681 // converge. do them:
682 UpdateFlags out = in;
683 for (unsigned int i = 0; i < 5; ++i)
684 {
685 // The following is a little incorrect:
686 // If not applied on a face,
687 // update_boundary_forms does not
688 // make sense. On the other hand,
689 // it is necessary on a
690 // face. Currently,
691 // update_boundary_forms is simply
692 // ignored for the interior of a
693 // cell.
696
697 if (out &
701
702 if (out &
707
708 // The contravariant transformation is used in the Piola
709 // transformation, which requires the determinant of the Jacobi
710 // matrix of the transformation. Because we have no way of
711 // knowing here whether the finite element wants to use the
712 // contravariant or the Piola transforms, we add the JxW values
713 // to the list of flags to be updated for each cell.
716
717 // the same is true when computing normal vectors: they require
718 // the determinant of the Jacobian
719 if (out & update_normal_vectors)
721 }
722
723 return out;
724}
725
726
727
728template <int dim, int spacedim>
729std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
731 const Quadrature<dim> &q) const
732{
733 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
734 std::make_unique<InternalData>(polynomial_degree);
735 auto &data = dynamic_cast<InternalData &>(*data_ptr);
736 data.initialize(this->requires_update_flags(update_flags), q, q.size());
737
738 return data_ptr;
739}
740
741
742
743template <int dim, int spacedim>
744std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
746 const UpdateFlags update_flags,
747 const hp::QCollection<dim - 1> &quadrature) const
748{
749 AssertDimension(quadrature.size(), 1);
750
751 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
752 std::make_unique<InternalData>(polynomial_degree);
753 auto &data = dynamic_cast<InternalData &>(*data_ptr);
754 data.initialize_face(this->requires_update_flags(update_flags),
756 ReferenceCells::get_hypercube<dim>(), quadrature[0]),
757 quadrature[0].size());
758
759 return data_ptr;
760}
761
762
763
764template <int dim, int spacedim>
765std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
767 const UpdateFlags update_flags,
768 const Quadrature<dim - 1> &quadrature) const
769{
770 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
771 std::make_unique<InternalData>(polynomial_degree);
772 auto &data = dynamic_cast<InternalData &>(*data_ptr);
773 data.initialize_face(this->requires_update_flags(update_flags),
775 ReferenceCells::get_hypercube<dim>(), quadrature),
776 quadrature.size());
777
778 return data_ptr;
779}
780
781
782
783template <int dim, int spacedim>
787 const CellSimilarity::Similarity cell_similarity,
788 const Quadrature<dim> &quadrature,
789 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
791 &output_data) const
792{
793 // ensure that the following static_cast is really correct:
794 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
796 const InternalData &data = static_cast<const InternalData &>(internal_data);
797 data.output_data = &output_data;
798
799 const unsigned int n_q_points = quadrature.size();
800
801 // recompute the support points of the transformation of this
802 // cell. we tried to be clever here in an earlier version of the
803 // library by checking whether the cell is the same as the one we
804 // had visited last, but it turns out to be difficult to determine
805 // that because a cell for the purposes of a mapping is
806 // characterized not just by its (triangulation, level, index)
807 // triple, but also by the locations of its vertices, the manifold
808 // object attached to the cell and all of its bounding faces/edges,
809 // etc. to reliably test that the "cell" we are on is, therefore,
810 // not easily done
811 data.mapping_support_points = this->compute_mapping_support_points(cell);
813
814 // if the order of the mapping is greater than 1, then do not reuse any cell
815 // similarity information. This is necessary because the cell similarity
816 // value is computed with just cell vertices and does not take into account
817 // cell curvature.
818 const CellSimilarity::Similarity computed_cell_similarity =
819 (polynomial_degree == 1 && this->preserves_vertex_locations() ?
820 cell_similarity :
822
823 if (dim > 1 && data.tensor_product_quadrature)
824 {
825 internal::MappingQImplementation::
826 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
827 computed_cell_similarity,
828 data,
829 output_data.quadrature_points,
830 output_data.jacobians,
831 output_data.inverse_jacobians,
832 output_data.jacobian_grads);
833 }
834 else
835 {
837 computed_cell_similarity,
838 data,
839 make_array_view(quadrature.get_points()),
840 polynomials_1d,
841 renumber_lexicographic_to_hierarchic,
842 output_data.quadrature_points,
843 output_data.jacobians,
844 output_data.inverse_jacobians);
845
847 spacedim>(
848 computed_cell_similarity,
849 data,
850 make_array_view(quadrature.get_points()),
851 polynomials_1d,
852 renumber_lexicographic_to_hierarchic,
853 output_data.jacobian_grads);
854 }
855
857 dim,
858 spacedim>(computed_cell_similarity,
859 data,
860 make_array_view(quadrature.get_points()),
861 polynomials_1d,
862 renumber_lexicographic_to_hierarchic,
864
866 dim,
867 spacedim>(computed_cell_similarity,
868 data,
869 make_array_view(quadrature.get_points()),
870 polynomials_1d,
871 renumber_lexicographic_to_hierarchic,
872 output_data.jacobian_2nd_derivatives);
873
874 internal::MappingQImplementation::
875 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
876 computed_cell_similarity,
877 data,
878 make_array_view(quadrature.get_points()),
879 polynomials_1d,
880 renumber_lexicographic_to_hierarchic,
882
884 dim,
885 spacedim>(computed_cell_similarity,
886 data,
887 make_array_view(quadrature.get_points()),
888 polynomials_1d,
889 renumber_lexicographic_to_hierarchic,
890 output_data.jacobian_3rd_derivatives);
891
892 internal::MappingQImplementation::
893 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
894 computed_cell_similarity,
895 data,
896 make_array_view(quadrature.get_points()),
897 polynomials_1d,
898 renumber_lexicographic_to_hierarchic,
900
901 const UpdateFlags update_flags = data.update_each;
902 const std::vector<double> &weights = quadrature.get_weights();
903
904 // Multiply quadrature weights by absolute value of Jacobian determinants or
905 // the area element g=sqrt(DX^t DX) in case of codim > 0
906
907 if (update_flags & (update_normal_vectors | update_JxW_values))
908 {
909 Assert(!(update_flags & update_JxW_values) ||
910 (output_data.JxW_values.size() == n_q_points),
911 ExcDimensionMismatch(output_data.JxW_values.size(), n_q_points));
912
913 Assert(!(update_flags & update_normal_vectors) ||
914 (output_data.normal_vectors.size() == n_q_points),
915 ExcDimensionMismatch(output_data.normal_vectors.size(),
916 n_q_points));
917
918
919 if (computed_cell_similarity != CellSimilarity::translation)
920 for (unsigned int point = 0; point < n_q_points; ++point)
921 {
922 if (dim == spacedim)
923 {
924 const double det = data.volume_elements[point];
925
926 // check for distorted cells.
927
928 // TODO: this allows for anisotropies of up to 1e6 in 3d and
929 // 1e12 in 2d. might want to find a finer
930 // (dimension-independent) criterion
931 Assert(det >
932 1e-12 * Utilities::fixed_power<dim>(
933 cell->diameter() / std::sqrt(double(dim))),
935 cell->center(), det, point)));
936
937 output_data.JxW_values[point] = weights[point] * det;
938 }
939 // if dim==spacedim, then there is no cell normal to
940 // compute. since this is for FEValues (and not FEFaceValues),
941 // there are also no face normals to compute
942 else // codim>0 case
943 {
944 Tensor<1, spacedim> DX_t[dim];
945 for (unsigned int i = 0; i < spacedim; ++i)
946 for (unsigned int j = 0; j < dim; ++j)
947 DX_t[j][i] = output_data.jacobians[point][i][j];
948
949 Tensor<2, dim> G; // First fundamental form
950 for (unsigned int i = 0; i < dim; ++i)
951 for (unsigned int j = 0; j < dim; ++j)
952 G[i][j] = DX_t[i] * DX_t[j];
953
954 if (update_flags & update_JxW_values)
955 output_data.JxW_values[point] =
956 std::sqrt(determinant(G)) * weights[point];
957
958 if (computed_cell_similarity ==
960 {
961 // we only need to flip the normal
962 if (update_flags & update_normal_vectors)
963 output_data.normal_vectors[point] *= -1.;
964 }
965 else
966 {
967 if (update_flags & update_normal_vectors)
968 {
969 Assert(spacedim == dim + 1,
971 "There is no (unique) cell normal for " +
973 "-dimensional cells in " +
974 Utilities::int_to_string(spacedim) +
975 "-dimensional space. This only works if the "
976 "space dimension is one greater than the "
977 "dimensionality of the mesh cells."));
978
979 if (dim == 1)
980 output_data.normal_vectors[point] =
981 cross_product_2d(-DX_t[0]);
982 else // dim == 2
983 output_data.normal_vectors[point] =
984 cross_product_3d(DX_t[0], DX_t[1]);
985
986 output_data.normal_vectors[point] /=
987 output_data.normal_vectors[point].norm();
988
989 if (cell->direction_flag() == false)
990 output_data.normal_vectors[point] *= -1.;
991 }
992 }
993 } // codim>0 case
994 }
995 }
996
997 return computed_cell_similarity;
998}
999
1000
1001
1002template <int dim, int spacedim>
1003void
1006 const unsigned int face_no,
1007 const hp::QCollection<dim - 1> &quadrature,
1008 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1010 &output_data) const
1011{
1012 AssertDimension(quadrature.size(), 1);
1013
1014 // ensure that the following cast is really correct:
1015 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1017 const InternalData &data = static_cast<const InternalData &>(internal_data);
1018 data.output_data = &output_data;
1019
1020 // if necessary, recompute the support points of the transformation of this
1021 // cell (note that we need to first check the triangulation pointer, since
1022 // otherwise the second test might trigger an exception if the
1023 // triangulations are not the same)
1024 if ((data.mapping_support_points.empty()) ||
1025 (&cell->get_triangulation() !=
1027 (cell != data.cell_of_current_support_points))
1028 {
1029 data.mapping_support_points = this->compute_mapping_support_points(cell);
1031 }
1032
1034 *this,
1035 cell,
1036 face_no,
1039 ReferenceCells::get_hypercube<dim>(),
1040 face_no,
1041 cell->face_orientation(face_no),
1042 cell->face_flip(face_no),
1043 cell->face_rotation(face_no),
1044 quadrature[0].size()),
1045 quadrature[0],
1046 data,
1047 polynomials_1d,
1048 renumber_lexicographic_to_hierarchic,
1049 output_data);
1050}
1051
1052
1053
1054template <int dim, int spacedim>
1055void
1058 const unsigned int face_no,
1059 const unsigned int subface_no,
1060 const Quadrature<dim - 1> &quadrature,
1061 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1063 &output_data) const
1064{
1065 // ensure that the following cast is really correct:
1066 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1068 const InternalData &data = static_cast<const InternalData &>(internal_data);
1069 data.output_data = &output_data;
1070
1071 // if necessary, recompute the support points of the transformation of this
1072 // cell (note that we need to first check the triangulation pointer, since
1073 // otherwise the second test might trigger an exception if the
1074 // triangulations are not the same)
1075 if ((data.mapping_support_points.empty()) ||
1076 (&cell->get_triangulation() !=
1078 (cell != data.cell_of_current_support_points))
1079 {
1080 data.mapping_support_points = this->compute_mapping_support_points(cell);
1082 }
1083
1085 *this,
1086 cell,
1087 face_no,
1088 subface_no,
1090 ReferenceCells::get_hypercube<dim>(),
1091 face_no,
1092 subface_no,
1093 cell->face_orientation(face_no),
1094 cell->face_flip(face_no),
1095 cell->face_rotation(face_no),
1096 quadrature.size(),
1097 cell->subface_case(face_no)),
1098 quadrature,
1099 data,
1100 polynomials_1d,
1101 renumber_lexicographic_to_hierarchic,
1102 output_data);
1103}
1104
1105
1106
1107template <int dim, int spacedim>
1108void
1112 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1114 &output_data) const
1115{
1116 Assert(dim == spacedim, ExcNotImplemented());
1117
1118 // ensure that the following static_cast is really correct:
1119 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1121 const InternalData &data = static_cast<const InternalData &>(internal_data);
1122 data.output_data = &output_data;
1123
1124 const unsigned int n_q_points = quadrature.size();
1125
1126 data.mapping_support_points = this->compute_mapping_support_points(cell);
1128
1131 data,
1132 make_array_view(quadrature.get_points()),
1133 polynomials_1d,
1134 renumber_lexicographic_to_hierarchic,
1135 output_data.quadrature_points,
1136 output_data.jacobians,
1137 output_data.inverse_jacobians);
1138
1139 internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1141 data,
1142 make_array_view(quadrature.get_points()),
1143 polynomials_1d,
1144 renumber_lexicographic_to_hierarchic,
1145 output_data.jacobian_grads);
1146
1148 dim,
1149 spacedim>(CellSimilarity::none,
1150 data,
1151 make_array_view(quadrature.get_points()),
1152 polynomials_1d,
1153 renumber_lexicographic_to_hierarchic,
1154 output_data.jacobian_pushed_forward_grads);
1155
1157 dim,
1158 spacedim>(CellSimilarity::none,
1159 data,
1160 make_array_view(quadrature.get_points()),
1161 polynomials_1d,
1162 renumber_lexicographic_to_hierarchic,
1163 output_data.jacobian_2nd_derivatives);
1164
1165 internal::MappingQImplementation::
1166 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1168 data,
1169 make_array_view(quadrature.get_points()),
1170 polynomials_1d,
1171 renumber_lexicographic_to_hierarchic,
1173
1175 dim,
1176 spacedim>(CellSimilarity::none,
1177 data,
1178 make_array_view(quadrature.get_points()),
1179 polynomials_1d,
1180 renumber_lexicographic_to_hierarchic,
1181 output_data.jacobian_3rd_derivatives);
1182
1183 internal::MappingQImplementation::
1184 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1186 data,
1187 make_array_view(quadrature.get_points()),
1188 polynomials_1d,
1189 renumber_lexicographic_to_hierarchic,
1191
1192 const UpdateFlags update_flags = data.update_each;
1193 const std::vector<double> &weights = quadrature.get_weights();
1194
1195 if ((update_flags & (update_normal_vectors | update_JxW_values)) != 0u)
1196 {
1197 AssertDimension(output_data.JxW_values.size(), n_q_points);
1198
1199 Assert(!(update_flags & update_normal_vectors) ||
1200 (output_data.normal_vectors.size() == n_q_points),
1201 ExcDimensionMismatch(output_data.normal_vectors.size(),
1202 n_q_points));
1203
1204
1205 for (unsigned int point = 0; point < n_q_points; ++point)
1206 {
1207 const double det = data.volume_elements[point];
1208
1209 // check for distorted cells.
1210
1211 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1212 // 1e12 in 2d. might want to find a finer
1213 // (dimension-independent) criterion
1214 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1215 cell->diameter() / std::sqrt(double(dim))),
1217 cell->center(), det, point)));
1218
1219 // The normals are n = J^{-T} * \hat{n} before normalizing.
1220 Tensor<1, spacedim> normal;
1221 for (unsigned int d = 0; d < spacedim; d++)
1222 normal[d] = output_data.inverse_jacobians[point].transpose()[d] *
1223 quadrature.normal_vector(point);
1224
1225 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1226
1227 if ((update_flags & update_normal_vectors) != 0u)
1228 {
1229 normal /= normal.norm();
1230 output_data.normal_vectors[point] = normal;
1231 }
1232 }
1233 }
1234}
1235
1236
1237
1238template <int dim, int spacedim>
1239void
1242 const ArrayView<const Point<dim>> &unit_points,
1243 const UpdateFlags update_flags,
1245 &output_data) const
1246{
1247 if (update_flags == update_default)
1248 return;
1249
1250 Assert(update_flags & update_inverse_jacobians ||
1251 update_flags & update_jacobians ||
1252 update_flags & update_quadrature_points,
1254
1255 output_data.initialize(unit_points.size(), update_flags);
1256
1257 auto internal_data =
1258 this->get_data(update_flags,
1259 Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
1260 unit_points.end())));
1261 const InternalData &data = static_cast<const InternalData &>(*internal_data);
1262 data.output_data = &output_data;
1263 data.mapping_support_points = this->compute_mapping_support_points(cell);
1264
1267 data,
1268 unit_points,
1269 polynomials_1d,
1270 renumber_lexicographic_to_hierarchic,
1271 output_data.quadrature_points,
1272 output_data.jacobians,
1273 output_data.inverse_jacobians);
1274}
1275
1276
1277
1278template <int dim, int spacedim>
1279void
1282 const unsigned int face_no,
1283 const Quadrature<dim - 1> &face_quadrature,
1284 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1286 &output_data) const
1287{
1288 if (face_quadrature.get_points().empty())
1289 return;
1290
1291 // ensure that the following static_cast is really correct:
1292 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1294 const InternalData &data = static_cast<const InternalData &>(internal_data);
1295
1296 data.mapping_support_points = this->compute_mapping_support_points(cell);
1297 data.output_data = &output_data;
1298
1300 *this,
1301 cell,
1302 face_no,
1305 face_quadrature,
1306 data,
1307 polynomials_1d,
1308 renumber_lexicographic_to_hierarchic,
1309 output_data);
1310}
1311
1312
1313
1314template <int dim, int spacedim>
1315void
1317 const ArrayView<const Tensor<1, dim>> &input,
1318 const MappingKind mapping_kind,
1319 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1320 const ArrayView<Tensor<1, spacedim>> &output) const
1321{
1323 mapping_kind,
1324 mapping_data,
1325 output);
1326}
1327
1328
1329
1330template <int dim, int spacedim>
1331void
1333 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1334 const MappingKind mapping_kind,
1335 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1336 const ArrayView<Tensor<2, spacedim>> &output) const
1337{
1339 mapping_kind,
1340 mapping_data,
1341 output);
1342}
1343
1344
1345
1346template <int dim, int spacedim>
1347void
1349 const ArrayView<const Tensor<2, dim>> &input,
1350 const MappingKind mapping_kind,
1351 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1352 const ArrayView<Tensor<2, spacedim>> &output) const
1353{
1354 switch (mapping_kind)
1355 {
1358 mapping_kind,
1359 mapping_data,
1360 output);
1361 return;
1362
1367 mapping_kind,
1368 mapping_data,
1369 output);
1370 return;
1371 default:
1373 }
1374}
1375
1376
1377
1378template <int dim, int spacedim>
1379void
1381 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1382 const MappingKind mapping_kind,
1383 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1384 const ArrayView<Tensor<3, spacedim>> &output) const
1385{
1386 AssertDimension(input.size(), output.size());
1387 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1390 &data = *static_cast<const InternalData &>(mapping_data).output_data;
1391
1392 switch (mapping_kind)
1393 {
1395 {
1396 Assert(!data.inverse_jacobians.empty(),
1398 "update_covariant_transformation"));
1399
1400 for (unsigned int q = 0; q < output.size(); ++q)
1401 for (unsigned int i = 0; i < spacedim; ++i)
1402 for (unsigned int j = 0; j < spacedim; ++j)
1403 {
1404 double tmp[dim];
1405 const DerivativeForm<1, dim, spacedim> covariant =
1406 data.inverse_jacobians[q].transpose();
1407 for (unsigned int K = 0; K < dim; ++K)
1408 {
1409 tmp[K] = covariant[j][0] * input[q][i][0][K];
1410 for (unsigned int J = 1; J < dim; ++J)
1411 tmp[K] += covariant[j][J] * input[q][i][J][K];
1412 }
1413 for (unsigned int k = 0; k < spacedim; ++k)
1414 {
1415 output[q][i][j][k] = covariant[k][0] * tmp[0];
1416 for (unsigned int K = 1; K < dim; ++K)
1417 output[q][i][j][k] += covariant[k][K] * tmp[K];
1418 }
1419 }
1420 return;
1421 }
1422
1423 default:
1425 }
1426}
1427
1428
1429
1430template <int dim, int spacedim>
1431void
1433 const ArrayView<const Tensor<3, dim>> &input,
1434 const MappingKind mapping_kind,
1435 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1436 const ArrayView<Tensor<3, spacedim>> &output) const
1437{
1438 switch (mapping_kind)
1439 {
1444 mapping_kind,
1445 mapping_data,
1446 output);
1447 return;
1448 default:
1450 }
1451}
1452
1453
1454
1455template <int dim, int spacedim>
1456void
1459 std::vector<Point<spacedim>> &a) const
1460{
1461 // if we only need the midpoint, then ask for it.
1462 if (this->polynomial_degree == 2)
1463 {
1464 for (unsigned int line_no = 0;
1465 line_no < GeometryInfo<dim>::lines_per_cell;
1466 ++line_no)
1467 {
1469 (dim == 1 ?
1470 static_cast<
1472 cell->line(line_no));
1473
1474 const Manifold<dim, spacedim> &manifold =
1475 ((line->manifold_id() == numbers::flat_manifold_id) &&
1476 (dim < spacedim) ?
1477 cell->get_manifold() :
1478 line->get_manifold());
1479 a.push_back(manifold.get_new_point_on_line(line));
1480 }
1481 }
1482 else
1483 // otherwise call the more complicated functions and ask for inner points
1484 // from the manifold description
1485 {
1486 std::vector<Point<spacedim>> tmp_points;
1487 for (unsigned int line_no = 0;
1488 line_no < GeometryInfo<dim>::lines_per_cell;
1489 ++line_no)
1490 {
1492 (dim == 1 ?
1493 static_cast<
1495 cell->line(line_no));
1496
1497 const Manifold<dim, spacedim> &manifold =
1498 ((line->manifold_id() == numbers::flat_manifold_id) &&
1499 (dim < spacedim) ?
1500 cell->get_manifold() :
1501 line->get_manifold());
1502
1503 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1504 const std::array<Point<spacedim>, 2> vertices{
1505 {cell->vertex(reference_cell.line_to_cell_vertices(line_no, 0)),
1506 cell->vertex(reference_cell.line_to_cell_vertices(line_no, 1))}};
1507
1508 const std::size_t n_rows =
1509 support_point_weights_perimeter_to_interior[0].size(0);
1510 a.resize(a.size() + n_rows);
1511 auto a_view = make_array_view(a.end() - n_rows, a.end());
1512 manifold.get_new_points(
1513 make_array_view(vertices.begin(), vertices.end()),
1514 support_point_weights_perimeter_to_interior[0],
1515 a_view);
1516 }
1517 }
1518}
1519
1520
1521
1522template <>
1523void
1526 std::vector<Point<3>> &a) const
1527{
1528 const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1529
1530 // used if face quad at boundary or entirely in the interior of the domain
1531 std::vector<Point<3>> tmp_points;
1532
1533 // loop over all faces and collect points on them
1534 for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1535 {
1536 const Triangulation<3>::face_iterator face = cell->face(face_no);
1537
1538#ifdef DEBUG
1539 const bool face_orientation = cell->face_orientation(face_no),
1540 face_flip = cell->face_flip(face_no),
1541 face_rotation = cell->face_rotation(face_no);
1542 const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1543 lines_per_face = GeometryInfo<3>::lines_per_face;
1544
1545 // some sanity checks up front
1546 for (unsigned int i = 0; i < vertices_per_face; ++i)
1547 Assert(face->vertex_index(i) ==
1548 cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1549 face_no, i, face_orientation, face_flip, face_rotation)),
1551
1552 // indices of the lines that bound a face are given by GeometryInfo<3>::
1553 // face_to_cell_lines
1554 for (unsigned int i = 0; i < lines_per_face; ++i)
1555 Assert(face->line(i) ==
1557 face_no, i, face_orientation, face_flip, face_rotation)),
1559#endif
1560 // extract the points surrounding a quad from the points
1561 // already computed. First get the 4 vertices and then the points on
1562 // the four lines
1563 boost::container::small_vector<Point<3>, 200> tmp_points(
1565 GeometryInfo<2>::lines_per_cell * (polynomial_degree - 1));
1566 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1567 tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1568 if (polynomial_degree > 1)
1569 for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1570 ++line)
1571 for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1572 tmp_points[4 + line * (polynomial_degree - 1) + i] =
1574 (polynomial_degree - 1) *
1576 i];
1577
1578 const std::size_t n_rows =
1579 support_point_weights_perimeter_to_interior[1].size(0);
1580 a.resize(a.size() + n_rows);
1581 auto a_view = make_array_view(a.end() - n_rows, a.end());
1583 make_array_view(tmp_points.begin(), tmp_points.end()),
1584 support_point_weights_perimeter_to_interior[1],
1585 a_view);
1586 }
1587}
1588
1589
1590
1591template <>
1592void
1595 std::vector<Point<3>> &a) const
1596{
1597 std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1598 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1599 vertices[i] = cell->vertex(i);
1600
1601 Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1603 for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1604 for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1605 {
1606 Point<2> point(line_support_points[q1 + 1][0],
1607 line_support_points[q2 + 1][0]);
1608 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1609 weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1610 }
1611
1612 const std::size_t n_rows = weights.size(0);
1613 a.resize(a.size() + n_rows);
1614 auto a_view = make_array_view(a.end() - n_rows, a.end());
1615 cell->get_manifold().get_new_points(
1616 make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1617}
1618
1619
1620
1621template <int dim, int spacedim>
1622void
1629
1630
1631
1632template <int dim, int spacedim>
1633std::vector<Point<spacedim>>
1635 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1636{
1637 // get the vertices first
1638 std::vector<Point<spacedim>> a;
1639 a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1640 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1641 a.push_back(cell->vertex(i));
1642
1643 if (this->polynomial_degree > 1)
1644 {
1645 // check if all entities have the same manifold id which is when we can
1646 // simply ask the manifold for all points. the transfinite manifold can
1647 // do the interpolation better than this class, so if we detect that we
1648 // do not have to change anything here
1649 Assert(dim <= 3, ExcImpossibleInDim(dim));
1650 bool all_manifold_ids_are_equal = (dim == spacedim);
1651 if (all_manifold_ids_are_equal &&
1653 &cell->get_manifold()) == nullptr)
1654 {
1655 for (auto f : GeometryInfo<dim>::face_indices())
1656 if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1657 all_manifold_ids_are_equal = false;
1658
1659 if (dim == 3)
1660 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1661 if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1662 all_manifold_ids_are_equal = false;
1663 }
1664
1665 if (all_manifold_ids_are_equal)
1666 {
1667 const std::size_t n_rows = support_point_weights_cell.size(0);
1668 a.resize(a.size() + n_rows);
1669 auto a_view = make_array_view(a.end() - n_rows, a.end());
1670 cell->get_manifold().get_new_points(make_array_view(a.begin(),
1671 a.end() - n_rows),
1672 support_point_weights_cell,
1673 a_view);
1674 }
1675 else
1676 switch (dim)
1677 {
1678 case 1:
1679 add_line_support_points(cell, a);
1680 break;
1681 case 2:
1682 // in 2d, add the points on the four bounding lines to the
1683 // exterior (outer) points
1684 add_line_support_points(cell, a);
1685
1686 // then get the interior support points
1687 if (dim != spacedim)
1688 add_quad_support_points(cell, a);
1689 else
1690 {
1691 const std::size_t n_rows =
1692 support_point_weights_perimeter_to_interior[1].size(0);
1693 a.resize(a.size() + n_rows);
1694 auto a_view = make_array_view(a.end() - n_rows, a.end());
1695 cell->get_manifold().get_new_points(
1696 make_array_view(a.begin(), a.end() - n_rows),
1697 support_point_weights_perimeter_to_interior[1],
1698 a_view);
1699 }
1700 break;
1701
1702 case 3:
1703 // in 3d also add the points located on the boundary faces
1704 add_line_support_points(cell, a);
1705 add_quad_support_points(cell, a);
1706
1707 // then compute the interior points
1708 {
1709 const std::size_t n_rows =
1710 support_point_weights_perimeter_to_interior[2].size(0);
1711 a.resize(a.size() + n_rows);
1712 auto a_view = make_array_view(a.end() - n_rows, a.end());
1713 cell->get_manifold().get_new_points(
1714 make_array_view(a.begin(), a.end() - n_rows),
1715 support_point_weights_perimeter_to_interior[2],
1716 a_view);
1717 }
1718 break;
1719
1720 default:
1722 break;
1723 }
1724 }
1725
1726 return a;
1727}
1728
1729
1730
1731template <int dim, int spacedim>
1734 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1735{
1736 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
1737}
1738
1739
1740
1741template <int dim, int spacedim>
1742bool
1744 const ReferenceCell &reference_cell) const
1745{
1746 Assert(dim == reference_cell.get_dimension(),
1747 ExcMessage("The dimension of your mapping (" +
1749 ") and the reference cell cell_type (" +
1750 Utilities::to_string(reference_cell.get_dimension()) +
1751 " ) do not agree."));
1752
1753 return reference_cell.is_hyper_cube();
1754}
1755
1756
1757
1758//--------------------------- Explicit instantiations -----------------------
1759#include "mapping_q.inst"
1760
1761
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:949
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition mapping_q.h:444
virtual std::size_t memory_consumption() const override
Definition mapping_q.cc:63
std::vector< Point< spacedim > > mapping_support_points
Definition mapping_q.h:438
internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > * output_data
Definition mapping_q.h:458
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_q.cc:161
InternalData(const unsigned int polynomial_degree)
Definition mapping_q.cc:50
AlignedVector< double > volume_elements
Definition mapping_q.h:450
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_q.cc:81
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition mapping_q.h:556
const Table< 2, double > support_point_weights_cell
Definition mapping_q.h:604
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition mapping_q.cc:730
void fill_mapping_data_for_face_quadrature(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_number, const Quadrature< dim - 1 > &face_quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition mapping_q.cc:766
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 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 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
Definition mapping_q.cc:785
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
const unsigned int polynomial_degree
Definition mapping_q.h:532
virtual void transform_points_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< spacedim > > &real_points, const ArrayView< Point< dim > > &unit_points) const override
Definition mapping_q.cc:592
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
Definition mapping_q.cc:745
Point< dim > transform_real_to_unit_cell_internal(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_p_unit) const
Definition mapping_q.cc:323
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition mapping_q.h:590
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
const std::vector< Point< 1 > > line_support_points
Definition mapping_q.h:542
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
Definition mapping_q.cc:289
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition mapping_q.h:549
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition mapping_q.cc:271
const std::vector< Point< dim > > unit_cell_support_points
Definition mapping_q.h:568
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
Definition mapping_q.cc:472
MappingQ(const unsigned int polynomial_degree)
Definition mapping_q.cc:219
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim > > &a) const
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim > > &a) const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition mapping_q.cc:674
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 bool is_compatible_with(const ReferenceCell &reference_cell) const override
unsigned int get_degree() const
Definition mapping_q.cc:280
Abstract base class for mapping classes.
Definition mapping.h:316
virtual void transform_points_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const ArrayView< const Point< spacedim > > &real_points, const ArrayView< Point< dim > > &unit_points) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition point.h:111
bool is_tensor_product() const
const std::vector< double > & get_weights() const
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
numbers::NumberTraits< Number >::real_type norm() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int size() const
Definition collection.h:264
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
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
Point< dim, Number > compute(const Point< spacedim, Number > &p) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
TriaIterator< TriaAccessor< dim - 1, dim, spacedim > > face_iterator
Definition tria.h:1609
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1657
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_default
No update.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
MappingKind
Definition mapping.h:77
@ mapping_covariant_gradient
Definition mapping.h:98
@ mapping_contravariant
Definition mapping.h:92
@ mapping_contravariant_hessian
Definition mapping.h:154
@ mapping_covariant_hessian
Definition mapping.h:148
@ mapping_contravariant_gradient
Definition mapping.h:104
@ mapping_piola_gradient
Definition mapping.h:118
@ mapping_piola_hessian
Definition mapping.h:160
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:479
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
constexpr T pow(const T base, const int iexp)
Definition utilities.h:963
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
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 maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 2, dim, spacedim > > &jacobian_grads)
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 maybe_update_jacobian_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 4, dim, spacedim > > &jacobian_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Tensor< 3, spacedim > > &jacobian_pushed_forward_grads)
void maybe_update_q_points_Jacobians_generic(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< Point< spacedim > > &quadrature_points, std::vector< DerivativeForm< 1, dim, spacedim > > &jacobians, std::vector< DerivativeForm< 1, spacedim, dim > > &inverse_jacobians)
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_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename ::MappingQ< dim, spacedim >::InternalData &data, const ArrayView< const Point< dim > > &unit_points, const std::vector< Polynomials::Polynomial< double > > &polynomials_1d, const std::vector< unsigned int > &renumber_lexicographic_to_hierarchic, std::vector< DerivativeForm< 3, dim, spacedim > > &jacobian_2nd_derivatives)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::manifold_id flat_manifold_id
Definition types.h:325
bool is_finite(const double x)
Definition numbers.h:538
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
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 > &)