deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+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 Quadrature<dim> &quadrature)
83{
84 // store the flags in the internal data object so we can access them
85 // in fill_fe_*_values()
86 this->update_each = update_flags;
87
88 const unsigned int n_q_points = quadrature.size();
89
90 if (this->update_each & update_volume_elements)
91 volume_elements.resize(n_q_points);
92
93 tensor_product_quadrature = quadrature.is_tensor_product();
94
95 // use of MatrixFree only for higher order elements and with more than one
96 // point where tensor products do not make sense
97 if (polynomial_degree < 2 || n_q_points == 1)
98 tensor_product_quadrature = false;
99
100 if (dim > 1)
101 {
102 // find out if the one-dimensional formula is the same
103 // in all directions
104 if (tensor_product_quadrature)
105 {
106 const std::array<Quadrature<1>, dim> &quad_array =
107 quadrature.get_tensor_basis();
108 for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
109 {
110 if (quad_array[i - 1].size() != quad_array[i].size())
111 {
112 tensor_product_quadrature = false;
113 break;
114 }
115 else
116 {
117 const std::vector<Point<1>> &points_1 =
118 quad_array[i - 1].get_points();
119 const std::vector<Point<1>> &points_2 =
120 quad_array[i].get_points();
121 const std::vector<double> &weights_1 =
122 quad_array[i - 1].get_weights();
123 const std::vector<double> &weights_2 =
124 quad_array[i].get_weights();
125 for (unsigned int j = 0; j < quad_array[i].size(); ++j)
127 if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
128 std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
129 {
130 tensor_product_quadrature = false;
131 break;
132 }
134 }
135 }
136
137 if (tensor_product_quadrature)
138 {
139 // use a 1d FE_DGQ and adjust the hierarchic -> lexicographic
140 // numbering manually (building an FE_Q<dim> is relatively
141 // expensive due to constraints)
143 shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
144 shape_info.lexicographic_numbering =
145 FETools::lexicographic_to_hierarchic_numbering<dim>(
147 shape_info.n_q_points = n_q_points;
148 shape_info.dofs_per_component_on_cell =
150 }
151 }
152 }
153}
154
155
156
157template <int dim, int spacedim>
158void
160 const UpdateFlags update_flags,
161 const Quadrature<dim> &quadrature,
162 const unsigned int n_original_q_points)
164 reinit(update_flags, quadrature);
165
166 quadrature_points = quadrature.get_points();
167
168 if (dim > 1 && tensor_product_quadrature)
170 constexpr unsigned int facedim = dim - 1;
172 shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
173 shape_info.lexicographic_numbering =
174 FETools::lexicographic_to_hierarchic_numbering<facedim>(
176 shape_info.n_q_points = n_original_q_points;
177 shape_info.dofs_per_component_on_cell =
179 }
180
181 if (dim > 1)
182 {
183 if (this->update_each &
186 aux.resize(dim - 1);
187 aux[0].resize(n_original_q_points);
188 if (dim > 2)
189 aux[1].resize(n_original_q_points);
190
191 // Compute tangentials to the unit cell.
192 for (const unsigned int i : GeometryInfo<dim>::face_indices())
193 {
194 unit_tangentials[i].resize(n_original_q_points);
195 std::fill(unit_tangentials[i].begin(),
196 unit_tangentials[i].end(),
198 if (dim > 2)
200 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
201 .resize(n_original_q_points);
202 std::fill(
203 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
204 .begin(),
205 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
206 .end(),
208 }
209 }
210 }
211 }
212}
214
215
216template <int dim, int spacedim>
220 QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
222 Polynomials::generate_complete_Lagrange_basis(line_support_points))
224 FETools::lexicographic_to_hierarchic_numbering<dim>(p))
226 internal::MappingQImplementation::unit_support_points<dim>(
230 internal::MappingQImplementation::
231 compute_support_point_weights_perimeter_to_interior(
233 dim))
235 internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
237{
238 Assert(p >= 1,
239 ExcMessage("It only makes sense to create polynomial mappings "
240 "with a polynomial degree greater or equal to one."));
241}
242
243
244
245template <int dim, int spacedim>
247 : polynomial_degree(mapping.polynomial_degree)
248 , line_support_points(mapping.line_support_points)
249 , polynomials_1d(mapping.polynomials_1d)
250 , renumber_lexicographic_to_hierarchic(
251 mapping.renumber_lexicographic_to_hierarchic)
252 , unit_cell_support_points(mapping.unit_cell_support_points)
253 , support_point_weights_perimeter_to_interior(
254 mapping.support_point_weights_perimeter_to_interior)
255 , support_point_weights_cell(mapping.support_point_weights_cell)
256{}
257
258
259
260template <int dim, int spacedim>
261std::unique_ptr<Mapping<dim, spacedim>>
263{
264 return std::make_unique<MappingQ<dim, spacedim>>(*this);
265}
266
267
268
269template <int dim, int spacedim>
270unsigned int
272{
273 return polynomial_degree;
274}
275
276
277
278template <int dim, int spacedim>
282 const Point<dim> &p) const
283{
284 if (polynomial_degree == 1)
285 {
286 const auto vertices = this->get_vertices(cell);
287 return Point<spacedim>(
289 }
290 else
292 polynomials_1d,
293 make_const_array_view(this->compute_mapping_support_points(cell)),
294 p,
295 polynomials_1d.size() == 2,
296 renumber_lexicographic_to_hierarchic));
297}
298
299
300// In the code below, GCC tries to instantiate MappingQ<3,4> when
301// seeing which of the overloaded versions of
302// do_transform_real_to_unit_cell_internal() to call. This leads to bad
303// error messages and, generally, nothing very good. Avoid this by ensuring
304// that this class exists, but does not have an inner InternalData
305// type, thereby ruling out the codim-1 version of the function
306// below when doing overload resolution.
307template <>
308class MappingQ<3, 4>
309{};
310
311
312
313// visual studio freaks out when trying to determine if
314// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
315// candidate. So instead of letting the compiler pick the correct overload, we
316// use template specialization to make sure we pick up the right function to
317// call:
318
319template <int dim, int spacedim>
323 const Point<spacedim> &,
324 const Point<dim> &) const
325{
326 // default implementation (should never be called)
328 return {};
329}
330
331
332
333template <>
337 const Point<1> &p,
338 const Point<1> &initial_p_unit) const
339{
340 // dispatch to the various specializations for spacedim=dim,
341 // spacedim=dim+1, etc
342 if (polynomial_degree == 1)
343 {
344 const auto vertices = this->get_vertices(cell);
345 return internal::MappingQImplementation::
346 do_transform_real_to_unit_cell_internal<1>(
347 p,
349 ArrayView<const Point<1>>(vertices.data(), vertices.size()),
350 polynomials_1d,
351 renumber_lexicographic_to_hierarchic);
352 }
353 else
354 return internal::MappingQImplementation::
355 do_transform_real_to_unit_cell_internal<1>(
356 p,
358 make_const_array_view(this->compute_mapping_support_points(cell)),
359 polynomials_1d,
360 renumber_lexicographic_to_hierarchic);
361}
362
363
364
365template <>
369 const Point<2> &p,
370 const Point<2> &initial_p_unit) const
371{
372 if (polynomial_degree == 1)
373 {
374 const auto vertices = this->get_vertices(cell);
375 return internal::MappingQImplementation::
376 do_transform_real_to_unit_cell_internal<2>(
377 p,
379 ArrayView<const Point<2>>(vertices.data(), vertices.size()),
380 polynomials_1d,
381 renumber_lexicographic_to_hierarchic);
382 }
383 else
384 return internal::MappingQImplementation::
385 do_transform_real_to_unit_cell_internal<2>(
386 p,
388 make_const_array_view(this->compute_mapping_support_points(cell)),
389 polynomials_1d,
390 renumber_lexicographic_to_hierarchic);
391}
392
393
394
395template <>
399 const Point<3> &p,
400 const Point<3> &initial_p_unit) const
401{
402 if (polynomial_degree == 1)
403 {
404 const auto vertices = this->get_vertices(cell);
405 return internal::MappingQImplementation::
406 do_transform_real_to_unit_cell_internal<3>(
407 p,
409 ArrayView<const Point<3>>(vertices.data(), vertices.size()),
410 polynomials_1d,
411 renumber_lexicographic_to_hierarchic);
412 }
413 else
414 return internal::MappingQImplementation::
415 do_transform_real_to_unit_cell_internal<3>(
416 p,
418 make_const_array_view(this->compute_mapping_support_points(cell)),
419 polynomials_1d,
420 renumber_lexicographic_to_hierarchic);
421}
422
423
424
425template <>
429 const Point<2> &p,
430 const Point<1> &initial_p_unit) const
431{
432 const int dim = 1;
433 const int spacedim = 2;
434
436
438 if (spacedim > dim)
439 update_flags |= update_jacobian_grads;
440 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
441 get_data(update_flags, point_quadrature));
442
443 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
444
445 // dispatch to the various specializations for spacedim=dim,
446 // spacedim=dim+1, etc
447 return internal::MappingQImplementation::
448 do_transform_real_to_unit_cell_internal_codim1<1>(
449 p,
451 make_const_array_view(mdata->mapping_support_points),
452 polynomials_1d,
453 renumber_lexicographic_to_hierarchic);
454}
455
456
458template <>
462 const Point<3> &p,
463 const Point<2> &initial_p_unit) const
464{
465 const int dim = 2;
466 const int spacedim = 3;
469
471 if (spacedim > dim)
472 update_flags |= update_jacobian_grads;
473 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
474 get_data(update_flags, point_quadrature));
475
476 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
477
478 // dispatch to the various specializations for spacedim=dim,
479 // spacedim=dim+1, etc
480 return internal::MappingQImplementation::
481 do_transform_real_to_unit_cell_internal_codim1<2>(
482 p,
484 make_const_array_view(mdata->mapping_support_points),
485 polynomials_1d,
486 renumber_lexicographic_to_hierarchic);
487}
488
490
491template <>
501
502
503
504template <int dim, int spacedim>
508 const Point<spacedim> &p) const
509{
510 // Use an exact formula if one is available. this is only the case
511 // for Q1 mappings in 1d, and in 2d if dim==spacedim
512 if ((polynomial_degree == 1) &&
513 ((dim == 1) || ((dim == 2) && (dim == spacedim))))
514 {
515 // The dimension-dependent algorithms are much faster (about 25-45x in
516 // 2d) but fail most of the time when the given point (p) is not in the
517 // cell. The dimension-independent Newton algorithm given below is
518 // slower, but more robust (though it still sometimes fails). Therefore
519 // this function implements the following strategy based on the
520 // p's dimension:
521 //
522 // * In 1d this mapping is linear, so the mapping is always invertible
523 // (and the exact formula is known) as long as the cell has non-zero
524 // length.
525 // * In 2d the exact (quadratic) formula is called first. If either the
526 // exact formula does not succeed (negative discriminant in the
527 // quadratic formula) or succeeds but finds a solution outside of the
528 // unit cell, then the Newton solver is called. The rationale for the
529 // second choice is that the exact formula may provide two different
530 // answers when mapping a point outside of the real cell, but the
531 // Newton solver (if it converges) will only return one answer.
532 // Otherwise the exact formula successfully found a point in the unit
533 // cell and that value is returned.
534 // * In 3d there is no (known to the authors) exact formula, so the Newton
535 // algorithm is used.
536 const auto vertices_ = this->get_vertices(cell);
537
538 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
539 vertices;
540 for (unsigned int i = 0; i < vertices.size(); ++i)
541 vertices[i] = vertices_[i];
542
543 try
544 {
545 switch (dim)
546 {
547 case 1:
548 {
549 // formula not subject to any issues in 1d
550 if (spacedim == 1)
552 vertices, p);
553 else
554 break;
555 }
556
557 case 2:
558 {
559 const Point<dim> point =
561 p);
562
563 // formula not guaranteed to work for points outside of
564 // the cell. only take the computed point if it lies
565 // inside the reference cell
566 const double eps = 1e-15;
567 if (-eps <= point[1] && point[1] <= 1 + eps &&
568 -eps <= point[0] && point[0] <= 1 + eps)
569 {
570 return point;
571 }
572 else
573 break;
574 }
575
576 default:
577 {
578 // we should not get here, based on the if-condition at the
579 // top
581 }
582 }
583 }
584 catch (
586 {
587 // simply fall through and continue on to the standard Newton code
588 }
589 }
590 else
591 {
592 // we can't use an explicit formula,
593 }
594
595
596 // Find the initial value for the Newton iteration by a normal
597 // projection to the least square plane determined by the vertices
598 // of the cell
600 if (this->preserves_vertex_locations())
601 {
602 initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
603 // in 1d with spacedim > 1 the affine approximation is exact
604 if (dim == 1 && polynomial_degree == 1)
605 return initial_p_unit;
606 }
607 else
608 {
609 // else, we simply use the mid point
610 for (unsigned int d = 0; d < dim; ++d)
611 initial_p_unit[d] = 0.5;
612 }
613
614 // perform the Newton iteration and return the result. note that this
615 // statement may throw an exception, which we simply pass up to the caller
616 const Point<dim> p_unit =
617 this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
618 AssertThrow(p_unit[0] != std::numeric_limits<double>::lowest(),
620 return p_unit;
621}
622
623
624
625template <int dim, int spacedim>
626void
629 const ArrayView<const Point<spacedim>> &real_points,
630 const ArrayView<Point<dim>> &unit_points) const
631{
632 // Go to base class functions for dim < spacedim because it is not yet
633 // implemented with optimized code.
634 if (dim < spacedim)
635 {
637 real_points,
638 unit_points);
639 return;
640 }
641
642 AssertDimension(real_points.size(), unit_points.size());
643 std::vector<Point<spacedim>> support_points_higher_order;
644 boost::container::small_vector<Point<spacedim>,
645#ifndef _MSC_VER
646 ReferenceCells::max_n_vertices<dim>()
647#else
649#endif
650 >
651 vertices;
652 if (polynomial_degree == 1)
653 vertices = this->get_vertices(cell);
654 else
655 support_points_higher_order = this->compute_mapping_support_points(cell);
657 polynomial_degree == 1 ? vertices.data() :
659 Utilities::pow(polynomial_degree + 1, dim));
660
661 // From the given (high-order) support points, now only pick the first
662 // 2^dim points and construct an affine approximation from those.
664 inverse_approximation(support_points, unit_cell_support_points);
666 const unsigned int n_points = real_points.size();
667 const unsigned int n_lanes = VectorizedArray<double>::size();
668
669 // Use the more heavy VectorizedArray code path if there is more than
670 // one point left to compute
671 for (unsigned int i = 0; i < n_points; i += n_lanes)
672 if (n_points - i > 1)
673 {
675 for (unsigned int j = 0; j < n_lanes; ++j)
676 if (i + j < n_points)
677 for (unsigned int d = 0; d < spacedim; ++d)
678 p_vec[d][j] = real_points[i + j][d];
679 else
680 for (unsigned int d = 0; d < spacedim; ++d)
681 p_vec[d][j] = real_points[i][d];
682
684 internal::MappingQImplementation::
685 do_transform_real_to_unit_cell_internal<dim, spacedim>(
686 p_vec,
689 polynomials_1d,
690 renumber_lexicographic_to_hierarchic);
691
692 // If the vectorized computation failed, it could be that only some of
693 // the lanes failed but others would have succeeded if we had let them
694 // compute alone without interference (like negative Jacobian
695 // determinants) from other SIMD lanes. Repeat the computation in this
696 // unlikely case with scalar arguments.
697 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
698 if (unit_point[0][j] != std::numeric_limits<double>::lowest())
699 for (unsigned int d = 0; d < dim; ++d)
700 unit_points[i + j][d] = unit_point[d][j];
701 else
702 unit_points[i + j] = internal::MappingQImplementation::
703 do_transform_real_to_unit_cell_internal<dim, spacedim>(
704 real_points[i + j],
705 inverse_approximation.compute(real_points[i + j]),
707 polynomials_1d,
708 renumber_lexicographic_to_hierarchic);
709 }
710 else
711 unit_points[i] = internal::MappingQImplementation::
712 do_transform_real_to_unit_cell_internal<dim, spacedim>(
713 real_points[i],
714 inverse_approximation.compute(real_points[i]),
716 polynomials_1d,
717 renumber_lexicographic_to_hierarchic);
718}
719
720
721
722template <int dim, int spacedim>
725{
726 // add flags if the respective quantities are necessary to compute
727 // what we need. note that some flags appear in both the conditions
728 // and in subsequent set operations. this leads to some circular
729 // logic. the only way to treat this is to iterate. since there are
730 // 5 if-clauses in the loop, it will take at most 5 iterations to
731 // converge. do them:
732 UpdateFlags out = in;
733 for (unsigned int i = 0; i < 5; ++i)
734 {
735 // The following is a little incorrect:
736 // If not applied on a face,
737 // update_boundary_forms does not
738 // make sense. On the other hand,
739 // it is necessary on a
740 // face. Currently,
741 // update_boundary_forms is simply
742 // ignored for the interior of a
743 // cell.
746
747 if (out &
751
752 if (out &
757
758 // The contravariant transformation is used in the Piola
759 // transformation, which requires the determinant of the Jacobi
760 // matrix of the transformation. Because we have no way of
761 // knowing here whether the finite element wants to use the
762 // contravariant or the Piola transforms, we add the JxW values
763 // to the list of flags to be updated for each cell.
766
767 // the same is true when computing normal vectors: they require
768 // the determinant of the Jacobian
769 if (out & update_normal_vectors)
771 }
772
773 return out;
774}
775
776
777
778template <int dim, int spacedim>
779std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
781 const Quadrature<dim> &q) const
782{
783 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
784 std::make_unique<InternalData>(polynomial_degree);
785 data_ptr->reinit(update_flags, q);
786 return data_ptr;
787}
788
789
790
791template <int dim, int spacedim>
792std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
794 const UpdateFlags update_flags,
795 const hp::QCollection<dim - 1> &quadrature) const
796{
797 AssertDimension(quadrature.size(), 1);
798
799 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
800 std::make_unique<InternalData>(polynomial_degree);
801 auto &data = dynamic_cast<InternalData &>(*data_ptr);
802 data.initialize_face(this->requires_update_flags(update_flags),
804 ReferenceCells::get_hypercube<dim>(), quadrature[0]),
805 quadrature[0].size());
806
807 return data_ptr;
808}
809
810
811
812template <int dim, int spacedim>
813std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
815 const UpdateFlags update_flags,
816 const Quadrature<dim - 1> &quadrature) const
817{
818 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
819 std::make_unique<InternalData>(polynomial_degree);
820 auto &data = dynamic_cast<InternalData &>(*data_ptr);
821 data.initialize_face(this->requires_update_flags(update_flags),
823 ReferenceCells::get_hypercube<dim>(), quadrature),
824 quadrature.size());
825
826 return data_ptr;
827}
828
829
830
831template <int dim, int spacedim>
835 const CellSimilarity::Similarity cell_similarity,
836 const Quadrature<dim> &quadrature,
839 &output_data) const
840{
841 // ensure that the following static_cast is really correct:
842 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
844 const InternalData &data = static_cast<const InternalData &>(internal_data);
845 data.output_data = &output_data;
846
847 const unsigned int n_q_points = quadrature.size();
848
849 // recompute the support points of the transformation of this
850 // cell. we tried to be clever here in an earlier version of the
851 // library by checking whether the cell is the same as the one we
852 // had visited last, but it turns out to be difficult to determine
853 // that because a cell for the purposes of a mapping is
854 // characterized not just by its (triangulation, level, index)
855 // triple, but also by the locations of its vertices, the manifold
856 // object attached to the cell and all of its bounding faces/edges,
857 // etc. to reliably test that the "cell" we are on is, therefore,
858 // not easily done
859 if (polynomial_degree == 1)
860 {
861 data.mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
862 const auto vertices = this->get_vertices(cell);
863 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
864 data.mapping_support_points[i] = vertices[i];
865 }
866 else
867 data.mapping_support_points = this->compute_mapping_support_points(cell);
868
869 data.cell_of_current_support_points = cell;
870
871 // if the order of the mapping is greater than 1, then do not reuse any cell
872 // similarity information. This is necessary because the cell similarity
873 // value is computed with just cell vertices and does not take into account
874 // cell curvature.
876 (polynomial_degree == 1 && this->preserves_vertex_locations() ?
877 cell_similarity :
879
880 if (dim > 1 && data.tensor_product_quadrature)
881 {
882 internal::MappingQImplementation::
883 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
885 data,
886 output_data.quadrature_points,
887 output_data.jacobians,
888 output_data.inverse_jacobians,
889 output_data.jacobian_grads);
890 }
891 else
892 {
895 data,
896 make_array_view(quadrature.get_points()),
897 polynomials_1d,
898 renumber_lexicographic_to_hierarchic,
899 output_data.quadrature_points,
900 output_data.jacobians,
901 output_data.inverse_jacobians);
902
904 spacedim>(
906 data,
907 make_array_view(quadrature.get_points()),
908 polynomials_1d,
909 renumber_lexicographic_to_hierarchic,
910 output_data.jacobian_grads);
911 }
912
914 dim,
915 spacedim>(computed_cell_similarity,
916 data,
917 make_array_view(quadrature.get_points()),
918 polynomials_1d,
919 renumber_lexicographic_to_hierarchic,
921
923 dim,
924 spacedim>(computed_cell_similarity,
925 data,
926 make_array_view(quadrature.get_points()),
927 polynomials_1d,
928 renumber_lexicographic_to_hierarchic,
929 output_data.jacobian_2nd_derivatives);
930
931 internal::MappingQImplementation::
932 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
934 data,
935 make_array_view(quadrature.get_points()),
936 polynomials_1d,
937 renumber_lexicographic_to_hierarchic,
939
941 dim,
942 spacedim>(computed_cell_similarity,
943 data,
944 make_array_view(quadrature.get_points()),
945 polynomials_1d,
946 renumber_lexicographic_to_hierarchic,
947 output_data.jacobian_3rd_derivatives);
948
949 internal::MappingQImplementation::
950 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
952 data,
953 make_array_view(quadrature.get_points()),
954 polynomials_1d,
955 renumber_lexicographic_to_hierarchic,
957
958 const UpdateFlags update_flags = data.update_each;
959 const std::vector<double> &weights = quadrature.get_weights();
960
961 // Multiply quadrature weights by absolute value of Jacobian determinants or
962 // the area element g=sqrt(DX^t DX) in case of codim > 0
963
964 if (update_flags & (update_normal_vectors | update_JxW_values))
965 {
966 Assert(!(update_flags & update_JxW_values) ||
967 (output_data.JxW_values.size() == n_q_points),
968 ExcDimensionMismatch(output_data.JxW_values.size(), n_q_points));
969
970 Assert(!(update_flags & update_normal_vectors) ||
971 (output_data.normal_vectors.size() == n_q_points),
972 ExcDimensionMismatch(output_data.normal_vectors.size(),
973 n_q_points));
974
975
977 for (unsigned int point = 0; point < n_q_points; ++point)
978 {
979 if (dim == spacedim)
980 {
981 const double det = data.volume_elements[point];
982
983 // check for distorted cells.
984
985 // TODO: this allows for anisotropies of up to 1e6 in 3d and
986 // 1e12 in 2d. might want to find a finer
987 // (dimension-independent) criterion
988 Assert(det >
989 1e-12 * Utilities::fixed_power<dim>(
990 cell->diameter() / std::sqrt(double(dim))),
992 cell->center(), det, point)));
993
994 output_data.JxW_values[point] = weights[point] * det;
995 }
996 // if dim==spacedim, then there is no cell normal to
997 // compute. since this is for FEValues (and not FEFaceValues),
998 // there are also no face normals to compute
999 else // codim>0 case
1000 {
1002 for (unsigned int i = 0; i < spacedim; ++i)
1003 for (unsigned int j = 0; j < dim; ++j)
1004 DX_t[j][i] = output_data.jacobians[point][i][j];
1005
1006 Tensor<2, dim> G; // First fundamental form
1007 for (unsigned int i = 0; i < dim; ++i)
1008 for (unsigned int j = 0; j < dim; ++j)
1009 G[i][j] = DX_t[i] * DX_t[j];
1010
1011 if (update_flags & update_JxW_values)
1012 output_data.JxW_values[point] =
1013 std::sqrt(determinant(G)) * weights[point];
1014
1017 {
1018 // we only need to flip the normal
1019 if (update_flags & update_normal_vectors)
1020 output_data.normal_vectors[point] *= -1.;
1021 }
1022 else
1023 {
1024 if (update_flags & update_normal_vectors)
1025 {
1026 Assert(spacedim == dim + 1,
1027 ExcMessage(
1028 "There is no (unique) cell normal for " +
1030 "-dimensional cells in " +
1031 Utilities::int_to_string(spacedim) +
1032 "-dimensional space. This only works if the "
1033 "space dimension is one greater than the "
1034 "dimensionality of the mesh cells."));
1035
1036 if (dim == 1)
1037 output_data.normal_vectors[point] =
1039 else // dim == 2
1040 output_data.normal_vectors[point] =
1041 cross_product_3d(DX_t[0], DX_t[1]);
1042
1043 output_data.normal_vectors[point] /=
1044 output_data.normal_vectors[point].norm();
1045
1046 if (cell->direction_flag() == false)
1047 output_data.normal_vectors[point] *= -1.;
1048 }
1049 }
1050 } // codim>0 case
1051 }
1052 }
1053
1055}
1056
1057
1058
1059template <int dim, int spacedim>
1060void
1063 const unsigned int face_no,
1064 const hp::QCollection<dim - 1> &quadrature,
1067 &output_data) const
1068{
1069 AssertDimension(quadrature.size(), 1);
1070
1071 // ensure that the following cast is really correct:
1072 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1074 const InternalData &data = static_cast<const InternalData &>(internal_data);
1075 data.output_data = &output_data;
1076
1077 // if necessary, recompute the support points of the transformation of this
1078 // cell (note that we need to first check the triangulation pointer, since
1079 // otherwise the second test might trigger an exception if the
1080 // triangulations are not the same)
1081 if ((data.mapping_support_points.empty()) ||
1082 (&cell->get_triangulation() !=
1083 &data.cell_of_current_support_points->get_triangulation()) ||
1084 (cell != data.cell_of_current_support_points))
1085 {
1086 if (polynomial_degree == 1)
1087 {
1088 data.mapping_support_points.resize(
1090 const auto vertices = this->get_vertices(cell);
1091 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
1092 ++i)
1093 data.mapping_support_points[i] = vertices[i];
1094 }
1095 else
1096 data.mapping_support_points =
1097 this->compute_mapping_support_points(cell);
1098 data.cell_of_current_support_points = cell;
1099 }
1100
1102 *this,
1103 cell,
1104 face_no,
1107 ReferenceCells::get_hypercube<dim>(),
1108 face_no,
1109 cell->combined_face_orientation(face_no),
1110 quadrature[0].size()),
1111 quadrature[0],
1112 data,
1113 polynomials_1d,
1114 renumber_lexicographic_to_hierarchic,
1115 output_data);
1116}
1117
1118
1119
1120template <int dim, int spacedim>
1121void
1124 const unsigned int face_no,
1125 const unsigned int subface_no,
1126 const Quadrature<dim - 1> &quadrature,
1129 &output_data) const
1130{
1131 // ensure that the following cast is really correct:
1132 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1134 const InternalData &data = static_cast<const InternalData &>(internal_data);
1135 data.output_data = &output_data;
1136
1137 // if necessary, recompute the support points of the transformation of this
1138 // cell (note that we need to first check the triangulation pointer, since
1139 // otherwise the second test might trigger an exception if the
1140 // triangulations are not the same)
1141 if ((data.mapping_support_points.empty()) ||
1142 (&cell->get_triangulation() !=
1143 &data.cell_of_current_support_points->get_triangulation()) ||
1144 (cell != data.cell_of_current_support_points))
1145 {
1146 if (polynomial_degree == 1)
1147 {
1148 data.mapping_support_points.resize(
1150 const auto vertices = this->get_vertices(cell);
1151 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
1152 ++i)
1153 data.mapping_support_points[i] = vertices[i];
1154 }
1155 else
1156 data.mapping_support_points =
1157 this->compute_mapping_support_points(cell);
1158 data.cell_of_current_support_points = cell;
1159 }
1160
1162 *this,
1163 cell,
1164 face_no,
1165 subface_no,
1167 ReferenceCells::get_hypercube<dim>(),
1168 face_no,
1169 subface_no,
1170 cell->combined_face_orientation(face_no),
1171 quadrature.size(),
1172 cell->subface_case(face_no)),
1173 quadrature,
1174 data,
1175 polynomials_1d,
1176 renumber_lexicographic_to_hierarchic,
1177 output_data);
1178}
1179
1180
1181
1182template <int dim, int spacedim>
1183void
1189 &output_data) const
1190{
1191 Assert(dim == spacedim, ExcNotImplemented());
1192
1193 // ensure that the following static_cast is really correct:
1194 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1196 const InternalData &data = static_cast<const InternalData &>(internal_data);
1197 data.output_data = &output_data;
1198
1199 const unsigned int n_q_points = quadrature.size();
1200
1201 if (polynomial_degree == 1)
1202 {
1203 data.mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
1204 const auto vertices = this->get_vertices(cell);
1205 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1206 data.mapping_support_points[i] = vertices[i];
1207 }
1208 else
1209 data.mapping_support_points = this->compute_mapping_support_points(cell);
1210 data.cell_of_current_support_points = cell;
1211
1214 data,
1215 make_array_view(quadrature.get_points()),
1216 polynomials_1d,
1217 renumber_lexicographic_to_hierarchic,
1218 output_data.quadrature_points,
1219 output_data.jacobians,
1220 output_data.inverse_jacobians);
1221
1222 internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1224 data,
1225 make_array_view(quadrature.get_points()),
1226 polynomials_1d,
1227 renumber_lexicographic_to_hierarchic,
1228 output_data.jacobian_grads);
1229
1231 dim,
1232 spacedim>(CellSimilarity::none,
1233 data,
1234 make_array_view(quadrature.get_points()),
1235 polynomials_1d,
1236 renumber_lexicographic_to_hierarchic,
1237 output_data.jacobian_pushed_forward_grads);
1238
1240 dim,
1241 spacedim>(CellSimilarity::none,
1242 data,
1243 make_array_view(quadrature.get_points()),
1244 polynomials_1d,
1245 renumber_lexicographic_to_hierarchic,
1246 output_data.jacobian_2nd_derivatives);
1247
1248 internal::MappingQImplementation::
1249 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1251 data,
1252 make_array_view(quadrature.get_points()),
1253 polynomials_1d,
1254 renumber_lexicographic_to_hierarchic,
1256
1258 dim,
1259 spacedim>(CellSimilarity::none,
1260 data,
1261 make_array_view(quadrature.get_points()),
1262 polynomials_1d,
1263 renumber_lexicographic_to_hierarchic,
1264 output_data.jacobian_3rd_derivatives);
1265
1266 internal::MappingQImplementation::
1267 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1269 data,
1270 make_array_view(quadrature.get_points()),
1271 polynomials_1d,
1272 renumber_lexicographic_to_hierarchic,
1274
1275 const UpdateFlags update_flags = data.update_each;
1276 const std::vector<double> &weights = quadrature.get_weights();
1277
1278 if ((update_flags & (update_normal_vectors | update_JxW_values)) != 0u)
1279 {
1280 AssertDimension(output_data.JxW_values.size(), n_q_points);
1281
1282 Assert(!(update_flags & update_normal_vectors) ||
1283 (output_data.normal_vectors.size() == n_q_points),
1284 ExcDimensionMismatch(output_data.normal_vectors.size(),
1285 n_q_points));
1286
1287
1288 for (unsigned int point = 0; point < n_q_points; ++point)
1289 {
1290 const double det = data.volume_elements[point];
1291
1292 // check for distorted cells.
1293
1294 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1295 // 1e12 in 2d. might want to find a finer
1296 // (dimension-independent) criterion
1297 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1298 cell->diameter() / std::sqrt(double(dim))),
1300 cell->center(), det, point)));
1301
1302 // The normals are n = J^{-T} * \hat{n} before normalizing.
1303 Tensor<1, spacedim> normal;
1304 for (unsigned int d = 0; d < spacedim; d++)
1305 normal[d] = output_data.inverse_jacobians[point].transpose()[d] *
1306 quadrature.normal_vector(point);
1307
1308 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1309
1310 if ((update_flags & update_normal_vectors) != 0u)
1311 {
1312 normal /= normal.norm();
1313 output_data.normal_vectors[point] = normal;
1314 }
1315 }
1316 }
1317}
1318
1319
1320
1321template <int dim, int spacedim>
1322void
1325 const ArrayView<const Point<dim>> &unit_points,
1326 const UpdateFlags update_flags,
1328 &output_data) const
1329{
1330 if (update_flags == update_default)
1331 return;
1332
1333 Assert(update_flags & update_inverse_jacobians ||
1334 update_flags & update_jacobians ||
1335 update_flags & update_quadrature_points,
1337
1338 output_data.initialize(unit_points.size(), update_flags);
1339
1340 auto internal_data =
1341 this->get_data(update_flags,
1342 Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
1343 unit_points.end())));
1344 const InternalData &data = static_cast<const InternalData &>(*internal_data);
1345 data.output_data = &output_data;
1346 if (polynomial_degree == 1)
1347 {
1348 data.mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
1349 const auto vertices = this->get_vertices(cell);
1350 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1351 data.mapping_support_points[i] = vertices[i];
1352 }
1353 else
1354 data.mapping_support_points = this->compute_mapping_support_points(cell);
1355
1358 data,
1359 unit_points,
1360 polynomials_1d,
1361 renumber_lexicographic_to_hierarchic,
1362 output_data.quadrature_points,
1363 output_data.jacobians,
1364 output_data.inverse_jacobians);
1365}
1366
1367
1368
1369template <int dim, int spacedim>
1370void
1373 const unsigned int face_no,
1374 const Quadrature<dim - 1> &face_quadrature,
1377 &output_data) const
1378{
1379 if (face_quadrature.get_points().empty())
1380 return;
1381
1382 // ensure that the following static_cast is really correct:
1383 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1385 const InternalData &data = static_cast<const InternalData &>(internal_data);
1386
1387 if (polynomial_degree == 1)
1388 {
1390 const auto vertices = this->get_vertices(cell);
1391 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1392 data.mapping_support_points[i] = vertices[i];
1393 }
1394 else
1395 data.mapping_support_points = this->compute_mapping_support_points(cell);
1396 data.output_data = &output_data;
1397
1399 *this,
1400 cell,
1401 face_no,
1404 face_quadrature,
1405 data,
1406 polynomials_1d,
1407 renumber_lexicographic_to_hierarchic,
1408 output_data);
1409}
1410
1411
1412
1413template <int dim, int spacedim>
1414void
1416 const ArrayView<const Tensor<1, dim>> &input,
1417 const MappingKind mapping_kind,
1418 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1419 const ArrayView<Tensor<1, spacedim>> &output) const
1420{
1422 mapping_kind,
1423 mapping_data,
1424 output);
1425}
1426
1427
1428
1429template <int dim, int spacedim>
1430void
1432 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1433 const MappingKind mapping_kind,
1434 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1435 const ArrayView<Tensor<2, spacedim>> &output) const
1436{
1438 mapping_kind,
1439 mapping_data,
1440 output);
1441}
1442
1443
1444
1445template <int dim, int spacedim>
1446void
1448 const ArrayView<const Tensor<2, dim>> &input,
1449 const MappingKind mapping_kind,
1450 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1451 const ArrayView<Tensor<2, spacedim>> &output) const
1452{
1453 switch (mapping_kind)
1454 {
1457 mapping_kind,
1458 mapping_data,
1459 output);
1460 return;
1461
1466 mapping_kind,
1467 mapping_data,
1468 output);
1469 return;
1470 default:
1472 }
1473}
1474
1475
1476
1477template <int dim, int spacedim>
1478void
1480 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1481 const MappingKind mapping_kind,
1482 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1483 const ArrayView<Tensor<3, spacedim>> &output) const
1484{
1485 AssertDimension(input.size(), output.size());
1486 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1489 &data = *static_cast<const InternalData &>(mapping_data).output_data;
1490
1491 switch (mapping_kind)
1492 {
1494 {
1495 Assert(!data.inverse_jacobians.empty(),
1497 "update_covariant_transformation"));
1498
1499 for (unsigned int q = 0; q < output.size(); ++q)
1500 for (unsigned int i = 0; i < spacedim; ++i)
1501 for (unsigned int j = 0; j < spacedim; ++j)
1502 {
1503 double tmp[dim];
1504 const DerivativeForm<1, dim, spacedim> covariant =
1505 data.inverse_jacobians[q].transpose();
1506 for (unsigned int K = 0; K < dim; ++K)
1507 {
1508 tmp[K] = covariant[j][0] * input[q][i][0][K];
1509 for (unsigned int J = 1; J < dim; ++J)
1510 tmp[K] += covariant[j][J] * input[q][i][J][K];
1511 }
1512 for (unsigned int k = 0; k < spacedim; ++k)
1513 {
1514 output[q][i][j][k] = covariant[k][0] * tmp[0];
1515 for (unsigned int K = 1; K < dim; ++K)
1516 output[q][i][j][k] += covariant[k][K] * tmp[K];
1517 }
1518 }
1519 return;
1520 }
1521
1522 default:
1524 }
1525}
1526
1527
1528
1529template <int dim, int spacedim>
1530void
1532 const ArrayView<const Tensor<3, dim>> &input,
1533 const MappingKind mapping_kind,
1534 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1535 const ArrayView<Tensor<3, spacedim>> &output) const
1536{
1537 switch (mapping_kind)
1538 {
1543 mapping_kind,
1544 mapping_data,
1545 output);
1546 return;
1547 default:
1549 }
1550}
1551
1552
1553
1554template <int dim, int spacedim>
1555void
1558 std::vector<Point<spacedim>> &a) const
1559{
1560 // if we only need the midpoint, then ask for it.
1561 if (this->polynomial_degree == 2)
1562 {
1563 for (unsigned int line_no = 0;
1565 ++line_no)
1566 {
1568 (dim == 1 ?
1569 static_cast<
1571 cell->line(line_no));
1572
1573 const Manifold<dim, spacedim> &manifold =
1574 ((line->manifold_id() == numbers::flat_manifold_id) &&
1575 (dim < spacedim) ?
1576 cell->get_manifold() :
1577 line->get_manifold());
1578 a.push_back(manifold.get_new_point_on_line(line));
1579 }
1580 }
1581 else
1582 // otherwise call the more complicated functions and ask for inner points
1583 // from the manifold description
1584 {
1585 std::vector<Point<spacedim>> tmp_points;
1586 for (unsigned int line_no = 0;
1588 ++line_no)
1589 {
1591 (dim == 1 ?
1592 static_cast<
1594 cell->line(line_no));
1595
1596 const Manifold<dim, spacedim> &manifold =
1597 ((line->manifold_id() == numbers::flat_manifold_id) &&
1598 (dim < spacedim) ?
1599 cell->get_manifold() :
1600 line->get_manifold());
1601
1602 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1603 const std::array<Point<spacedim>, 2> vertices{
1604 {cell->vertex(reference_cell.line_to_cell_vertices(line_no, 0)),
1605 cell->vertex(reference_cell.line_to_cell_vertices(line_no, 1))}};
1606
1607 const std::size_t n_rows =
1608 support_point_weights_perimeter_to_interior[0].size(0);
1609 a.resize(a.size() + n_rows);
1610 auto a_view = make_array_view(a.end() - n_rows, a.end());
1611 manifold.get_new_points(
1612 make_array_view(vertices.begin(), vertices.end()),
1613 support_point_weights_perimeter_to_interior[0],
1614 a_view);
1615 }
1616 }
1617}
1618
1619
1620
1621template <>
1622void
1625 std::vector<Point<3>> &a) const
1626{
1627 const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1628
1629 // used if face quad at boundary or entirely in the interior of the domain
1630 std::vector<Point<3>> tmp_points;
1631
1632 // loop over all faces and collect points on them
1633 for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1634 {
1635 const Triangulation<3>::face_iterator face = cell->face(face_no);
1636
1637#ifdef DEBUG
1638 const bool face_orientation = cell->face_orientation(face_no),
1639 face_flip = cell->face_flip(face_no),
1640 face_rotation = cell->face_rotation(face_no);
1641 const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1642 lines_per_face = GeometryInfo<3>::lines_per_face;
1643
1644 // some sanity checks up front
1645 for (unsigned int i = 0; i < vertices_per_face; ++i)
1646 Assert(face->vertex_index(i) ==
1647 cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1648 face_no, i, face_orientation, face_flip, face_rotation)),
1650
1651 // indices of the lines that bound a face are given by GeometryInfo<3>::
1652 // face_to_cell_lines
1653 for (unsigned int i = 0; i < lines_per_face; ++i)
1654 Assert(face->line(i) ==
1656 face_no, i, face_orientation, face_flip, face_rotation)),
1658#endif
1659 // extract the points surrounding a quad from the points
1660 // already computed. First get the 4 vertices and then the points on
1661 // the four lines
1662 boost::container::small_vector<Point<3>, 200> tmp_points(
1664 GeometryInfo<2>::lines_per_cell * (polynomial_degree - 1));
1665 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1667 if (polynomial_degree > 1)
1668 for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1669 ++line)
1670 for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1671 tmp_points[4 + line * (polynomial_degree - 1) + i] =
1673 (polynomial_degree - 1) *
1675 i];
1676
1677 const std::size_t n_rows =
1678 support_point_weights_perimeter_to_interior[1].size(0);
1679 a.resize(a.size() + n_rows);
1680 auto a_view = make_array_view(a.end() - n_rows, a.end());
1681 face->get_manifold().get_new_points(
1682 make_array_view(tmp_points.begin(), tmp_points.end()),
1683 support_point_weights_perimeter_to_interior[1],
1684 a_view);
1685 }
1686}
1687
1688
1689
1690template <>
1691void
1694 std::vector<Point<3>> &a) const
1695{
1696 std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1697 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1698 vertices[i] = cell->vertex(i);
1699
1700 Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1702 for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1703 for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1704 {
1705 Point<2> point(line_support_points[q1 + 1][0],
1706 line_support_points[q2 + 1][0]);
1707 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1708 weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1709 }
1710
1711 const std::size_t n_rows = weights.size(0);
1712 a.resize(a.size() + n_rows);
1713 auto a_view = make_array_view(a.end() - n_rows, a.end());
1714 cell->get_manifold().get_new_points(
1715 make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1716}
1717
1718
1719
1720template <int dim, int spacedim>
1721void
1728
1729
1730
1731template <int dim, int spacedim>
1732std::vector<Point<spacedim>>
1734 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1735{
1736 // get the vertices first
1737 std::vector<Point<spacedim>> a;
1738 a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1739 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1740 a.push_back(cell->vertex(i));
1741
1742 if (this->polynomial_degree > 1)
1743 {
1744 // check if all entities have the same manifold id which is when we can
1745 // simply ask the manifold for all points. the transfinite manifold can
1746 // do the interpolation better than this class, so if we detect that we
1747 // do not have to change anything here
1748 Assert(dim <= 3, ExcImpossibleInDim(dim));
1749 bool all_manifold_ids_are_equal = (dim == spacedim);
1752 &cell->get_manifold()) == nullptr)
1753 {
1754 for (auto f : GeometryInfo<dim>::face_indices())
1755 if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1757
1758 if (dim == 3)
1759 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1760 if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1762 }
1763
1765 {
1766 const std::size_t n_rows = support_point_weights_cell.size(0);
1767 a.resize(a.size() + n_rows);
1768 auto a_view = make_array_view(a.end() - n_rows, a.end());
1769 cell->get_manifold().get_new_points(make_array_view(a.begin(),
1770 a.end() - n_rows),
1771 support_point_weights_cell,
1772 a_view);
1773 }
1774 else
1775 switch (dim)
1776 {
1777 case 1:
1778 add_line_support_points(cell, a);
1779 break;
1780 case 2:
1781 // in 2d, add the points on the four bounding lines to the
1782 // exterior (outer) points
1783 add_line_support_points(cell, a);
1784
1785 // then get the interior support points
1786 if (dim != spacedim)
1787 add_quad_support_points(cell, a);
1788 else
1789 {
1790 const std::size_t n_rows =
1791 support_point_weights_perimeter_to_interior[1].size(0);
1792 a.resize(a.size() + n_rows);
1793 auto a_view = make_array_view(a.end() - n_rows, a.end());
1794 cell->get_manifold().get_new_points(
1795 make_array_view(a.begin(), a.end() - n_rows),
1796 support_point_weights_perimeter_to_interior[1],
1797 a_view);
1798 }
1799 break;
1800
1801 case 3:
1802 // in 3d also add the points located on the boundary faces
1803 add_line_support_points(cell, a);
1804 add_quad_support_points(cell, a);
1805
1806 // then compute the interior points
1807 {
1808 const std::size_t n_rows =
1809 support_point_weights_perimeter_to_interior[2].size(0);
1810 a.resize(a.size() + n_rows);
1811 auto a_view = make_array_view(a.end() - n_rows, a.end());
1812 cell->get_manifold().get_new_points(
1813 make_array_view(a.begin(), a.end() - n_rows),
1814 support_point_weights_perimeter_to_interior[2],
1815 a_view);
1816 }
1817 break;
1818
1819 default:
1821 break;
1822 }
1823 }
1824
1825 return a;
1826}
1827
1828
1829
1830template <int dim, int spacedim>
1833 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1834{
1835 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
1836}
1837
1838
1839
1840template <int dim, int spacedim>
1841bool
1843 const ReferenceCell &reference_cell) const
1844{
1845 Assert(dim == reference_cell.get_dimension(),
1846 ExcMessage("The dimension of your mapping (" +
1848 ") and the reference cell cell_type (" +
1849 Utilities::to_string(reference_cell.get_dimension()) +
1850 " ) do not agree."));
1851
1852 return reference_cell.is_hyper_cube();
1853}
1854
1855
1856
1857//--------------------------- Explicit instantiations -----------------------
1858#include "fe/mapping_q.inst"
1859
1860
auto make_const_array_view(const Container &container) -> decltype(make_array_view(container))
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
virtual std::size_t memory_consumption() const override
Definition mapping_q.cc:63
std::vector< Point< spacedim > > mapping_support_points
Definition mapping_q.h:421
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
Definition mapping_q.cc:81
internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > * output_data
Definition mapping_q.h:441
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition mapping_q.cc:159
InternalData(const unsigned int polynomial_degree)
Definition mapping_q.cc:50
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition mapping_q.h:539
const Table< 2, double > support_point_weights_cell
Definition mapping_q.h:587
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:780
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:814
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:833
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:515
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:627
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:793
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:321
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition mapping_q.h:573
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:525
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:280
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition mapping_q.h:532
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition mapping_q.cc:262
const std::vector< Point< dim > > unit_cell_support_points
Definition mapping_q.h:551
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:506
MappingQ(const unsigned int polynomial_degree)
Definition mapping_q.cc:217
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:724
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:271
Abstract base class for mapping classes.
Definition mapping.h:320
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
numbers::NumberTraits< Number >::real_type norm() const
friend class Tensor
Definition tensor.h:865
unsigned int size() const
Definition collection.h:316
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
void initialize(const unsigned int n_quadrature_points, const UpdateFlags flags)
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
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)
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1644
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:81
@ mapping_covariant_gradient
Definition mapping.h:102
@ mapping_contravariant
Definition mapping.h:96
@ mapping_contravariant_hessian
Definition mapping.h:158
@ mapping_covariant_hessian
Definition mapping.h:152
@ mapping_contravariant_gradient
Definition mapping.h:108
@ mapping_piola_gradient
Definition mapping.h:122
@ mapping_piola_hessian
Definition mapping.h:164
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:483
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
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_linear(const Number *values, const Point< dim, Number2 > &p)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::manifold_id flat_manifold_id
Definition types.h:336
::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 > &)