Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20: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 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)
126 {
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 }
133 }
134 }
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 }
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)
163{
164 reinit(update_flags, quadrature);
165
166 quadrature_points = quadrature.get_points();
167
168 if (dim > 1 && tensor_product_quadrature)
169 {
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 &
185 {
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)
199 {
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(),
209 }
210 }
211 }
212}
213
214
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(
232 this->polynomial_degree,
233 dim))
235 internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
236 this->polynomial_degree))
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>
246MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
247 : MappingQ<dim, spacedim>(p)
248{}
250
251
252template <int dim, int spacedim>
254 : polynomial_degree(mapping.polynomial_degree)
255 , line_support_points(mapping.line_support_points)
256 , polynomials_1d(mapping.polynomials_1d)
257 , renumber_lexicographic_to_hierarchic(
258 mapping.renumber_lexicographic_to_hierarchic)
259 , unit_cell_support_points(mapping.unit_cell_support_points)
260 , support_point_weights_perimeter_to_interior(
261 mapping.support_point_weights_perimeter_to_interior)
262 , support_point_weights_cell(mapping.support_point_weights_cell)
263{}
264
265
266
267template <int dim, int spacedim>
268std::unique_ptr<Mapping<dim, spacedim>>
270{
271 return std::make_unique<MappingQ<dim, spacedim>>(*this);
272}
273
274
275
276template <int dim, int spacedim>
277unsigned int
280 return polynomial_degree;
281}
282
283
284
285template <int dim, int spacedim>
289 const Point<dim> &p) const
290{
291 if (polynomial_degree == 1)
292 {
293 const auto vertices = this->get_vertices(cell);
294 return Point<spacedim>(
296 }
297 else
299 polynomials_1d,
300 make_const_array_view(this->compute_mapping_support_points(cell)),
301 p,
302 polynomials_1d.size() == 2,
303 renumber_lexicographic_to_hierarchic));
304}
305
306
307// In the code below, GCC tries to instantiate MappingQ<3,4> when
308// seeing which of the overloaded versions of
309// do_transform_real_to_unit_cell_internal() to call. This leads to bad
310// error messages and, generally, nothing very good. Avoid this by ensuring
311// that this class exists, but does not have an inner InternalData
312// type, thereby ruling out the codim-1 version of the function
313// below when doing overload resolution.
314template <>
315class MappingQ<3, 4>
316{};
317
318
319
320// visual studio freaks out when trying to determine if
321// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
322// candidate. So instead of letting the compiler pick the correct overload, we
323// use template specialization to make sure we pick up the right function to
324// call:
325
326template <int dim, int spacedim>
330 const Point<spacedim> &,
331 const Point<dim> &) const
332{
333 // default implementation (should never be called)
335 return {};
336}
337
338
339
340template <>
344 const Point<1> &p,
345 const Point<1> &initial_p_unit) const
346{
347 // dispatch to the various specializations for spacedim=dim,
348 // spacedim=dim+1, etc
349 if (polynomial_degree == 1)
350 {
351 const auto vertices = this->get_vertices(cell);
352 return internal::MappingQImplementation::
353 do_transform_real_to_unit_cell_internal<1>(
354 p,
355 initial_p_unit,
356 ArrayView<const Point<1>>(vertices.data(), vertices.size()),
357 polynomials_1d,
358 renumber_lexicographic_to_hierarchic);
359 }
360 else
361 return internal::MappingQImplementation::
362 do_transform_real_to_unit_cell_internal<1>(
363 p,
364 initial_p_unit,
365 make_const_array_view(this->compute_mapping_support_points(cell)),
366 polynomials_1d,
367 renumber_lexicographic_to_hierarchic);
368}
369
370
371
372template <>
376 const Point<2> &p,
377 const Point<2> &initial_p_unit) const
378{
379 if (polynomial_degree == 1)
380 {
381 const auto vertices = this->get_vertices(cell);
382 return internal::MappingQImplementation::
383 do_transform_real_to_unit_cell_internal<2>(
384 p,
385 initial_p_unit,
386 ArrayView<const Point<2>>(vertices.data(), vertices.size()),
387 polynomials_1d,
388 renumber_lexicographic_to_hierarchic);
389 }
390 else
391 return internal::MappingQImplementation::
392 do_transform_real_to_unit_cell_internal<2>(
393 p,
394 initial_p_unit,
395 make_const_array_view(this->compute_mapping_support_points(cell)),
396 polynomials_1d,
397 renumber_lexicographic_to_hierarchic);
398}
399
400
401
402template <>
406 const Point<3> &p,
407 const Point<3> &initial_p_unit) const
408{
409 if (polynomial_degree == 1)
410 {
411 const auto vertices = this->get_vertices(cell);
412 return internal::MappingQImplementation::
413 do_transform_real_to_unit_cell_internal<3>(
414 p,
415 initial_p_unit,
416 ArrayView<const Point<3>>(vertices.data(), vertices.size()),
417 polynomials_1d,
418 renumber_lexicographic_to_hierarchic);
419 }
420 else
421 return internal::MappingQImplementation::
422 do_transform_real_to_unit_cell_internal<3>(
423 p,
424 initial_p_unit,
425 make_const_array_view(this->compute_mapping_support_points(cell)),
426 polynomials_1d,
427 renumber_lexicographic_to_hierarchic);
428}
429
430
431
432template <>
436 const Point<2> &p,
437 const Point<1> &initial_p_unit) const
438{
439 const int dim = 1;
440 const int spacedim = 2;
441
442 const Quadrature<dim> point_quadrature(initial_p_unit);
443
445 if (spacedim > dim)
446 update_flags |= update_jacobian_grads;
447 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
448 get_data(update_flags, point_quadrature));
449
450 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
451
452 // dispatch to the various specializations for spacedim=dim,
453 // spacedim=dim+1, etc
454 return internal::MappingQImplementation::
455 do_transform_real_to_unit_cell_internal_codim1<1>(
456 p,
457 initial_p_unit,
458 make_const_array_view(mdata->mapping_support_points),
459 polynomials_1d,
460 renumber_lexicographic_to_hierarchic);
461}
462
463
464
465template <>
469 const Point<3> &p,
470 const Point<2> &initial_p_unit) const
472 const int dim = 2;
473 const int spacedim = 3;
474
475 const Quadrature<dim> point_quadrature(initial_p_unit);
478 if (spacedim > dim)
479 update_flags |= update_jacobian_grads;
480 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
481 get_data(update_flags, point_quadrature));
482
483 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
484
485 // dispatch to the various specializations for spacedim=dim,
486 // spacedim=dim+1, etc
487 return internal::MappingQImplementation::
488 do_transform_real_to_unit_cell_internal_codim1<2>(
489 p,
490 initial_p_unit,
491 make_const_array_view(mdata->mapping_support_points),
492 polynomials_1d,
493 renumber_lexicographic_to_hierarchic);
494}
495
496
497
498template <>
508
510
511template <int dim, int spacedim>
515 const Point<spacedim> &p) const
516{
517 // Use an exact formula if one is available. this is only the case
518 // for Q1 mappings in 1d, and in 2d if dim==spacedim
519 if ((polynomial_degree == 1) &&
520 ((dim == 1) || ((dim == 2) && (dim == spacedim))))
521 {
522 // The dimension-dependent algorithms are much faster (about 25-45x in
523 // 2d) but fail most of the time when the given point (p) is not in the
524 // cell. The dimension-independent Newton algorithm given below is
525 // slower, but more robust (though it still sometimes fails). Therefore
526 // this function implements the following strategy based on the
527 // p's dimension:
528 //
529 // * In 1d this mapping is linear, so the mapping is always invertible
530 // (and the exact formula is known) as long as the cell has non-zero
531 // length.
532 // * In 2d the exact (quadratic) formula is called first. If either the
533 // exact formula does not succeed (negative discriminant in the
534 // quadratic formula) or succeeds but finds a solution outside of the
535 // unit cell, then the Newton solver is called. The rationale for the
536 // second choice is that the exact formula may provide two different
537 // answers when mapping a point outside of the real cell, but the
538 // Newton solver (if it converges) will only return one answer.
539 // Otherwise the exact formula successfully found a point in the unit
540 // cell and that value is returned.
541 // * In 3d there is no (known to the authors) exact formula, so the Newton
542 // algorithm is used.
543 const auto vertices_ = this->get_vertices(cell);
544
545 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
546 vertices;
547 for (unsigned int i = 0; i < vertices.size(); ++i)
548 vertices[i] = vertices_[i];
549
550 try
551 {
552 switch (dim)
553 {
554 case 1:
555 {
556 // formula not subject to any issues in 1d
557 if (spacedim == 1)
559 vertices, p);
560 else
561 break;
562 }
563
564 case 2:
565 {
566 const Point<dim> point =
568 p);
569
570 // formula not guaranteed to work for points outside of
571 // the cell. only take the computed point if it lies
572 // inside the reference cell
573 const double eps = 1e-15;
574 if (-eps <= point[1] && point[1] <= 1 + eps &&
575 -eps <= point[0] && point[0] <= 1 + eps)
576 {
577 return point;
578 }
579 else
580 break;
581 }
582
583 default:
584 {
585 // we should get here, based on the if-condition at the top
587 }
588 }
589 }
590 catch (
592 {
593 // simply fall through and continue on to the standard Newton code
594 }
595 }
596 else
597 {
598 // we can't use an explicit formula,
599 }
600
601
602 // Find the initial value for the Newton iteration by a normal
603 // projection to the least square plane determined by the vertices
604 // of the cell
605 Point<dim> initial_p_unit;
606 if (this->preserves_vertex_locations())
607 {
608 initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
609 // in 1d with spacedim > 1 the affine approximation is exact
610 if (dim == 1 && polynomial_degree == 1)
611 return initial_p_unit;
612 }
613 else
614 {
615 // else, we simply use the mid point
616 for (unsigned int d = 0; d < dim; ++d)
617 initial_p_unit[d] = 0.5;
618 }
619
620 // perform the Newton iteration and return the result. note that this
621 // statement may throw an exception, which we simply pass up to the caller
622 const Point<dim> p_unit =
623 this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
626 return p_unit;
627}
629
630
631template <int dim, int spacedim>
632void
635 const ArrayView<const Point<spacedim>> &real_points,
636 const ArrayView<Point<dim>> &unit_points) const
637{
638 // Go to base class functions for dim < spacedim because it is not yet
639 // implemented with optimized code.
640 if (dim < spacedim)
641 {
643 real_points,
644 unit_points);
645 return;
646 }
647
648 AssertDimension(real_points.size(), unit_points.size());
649 std::vector<Point<spacedim>> support_points_higher_order;
650 boost::container::small_vector<Point<spacedim>,
652 vertices;
653 if (polynomial_degree == 1)
654 vertices = this->get_vertices(cell);
655 else
656 support_points_higher_order = this->compute_mapping_support_points(cell);
657 const ArrayView<const Point<spacedim>> support_points(
658 polynomial_degree == 1 ? vertices.data() :
659 support_points_higher_order.data(),
660 Utilities::pow(polynomial_degree + 1, dim));
661
662 // From the given (high-order) support points, now only pick the first
663 // 2^dim points and construct an affine approximation from those.
665 inverse_approximation(support_points, unit_cell_support_points);
666
667 const unsigned int n_points = real_points.size();
668 const unsigned int n_lanes = VectorizedArray<double>::size();
669
670 // Use the more heavy VectorizedArray code path if there is more than
671 // one point left to compute
672 for (unsigned int i = 0; i < n_points; i += n_lanes)
673 if (n_points - i > 1)
676 for (unsigned int j = 0; j < n_lanes; ++j)
677 if (i + j < n_points)
678 for (unsigned int d = 0; d < spacedim; ++d)
679 p_vec[d][j] = real_points[i + j][d];
680 else
681 for (unsigned int d = 0; d < spacedim; ++d)
682 p_vec[d][j] = real_points[i][d];
683
685 internal::MappingQImplementation::
686 do_transform_real_to_unit_cell_internal<dim, spacedim>(
687 p_vec,
688 inverse_approximation.compute(p_vec),
689 support_points,
690 polynomials_1d,
691 renumber_lexicographic_to_hierarchic);
692
693 // If the vectorized computation failed, it could be that only some of
694 // the lanes failed but others would have succeeded if we had let them
695 // compute alone without interference (like negative Jacobian
696 // determinants) from other SIMD lanes. Repeat the computation in this
697 // unlikely case with scalar arguments.
698 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
699 if (numbers::is_finite(unit_point[0][j]))
700 for (unsigned int d = 0; d < dim; ++d)
701 unit_points[i + j][d] = unit_point[d][j];
702 else
703 unit_points[i + j] = internal::MappingQImplementation::
704 do_transform_real_to_unit_cell_internal<dim, spacedim>(
705 real_points[i + j],
706 inverse_approximation.compute(real_points[i + j]),
707 support_points,
708 polynomials_1d,
709 renumber_lexicographic_to_hierarchic);
710 }
711 else
712 unit_points[i] = internal::MappingQImplementation::
713 do_transform_real_to_unit_cell_internal<dim, spacedim>(
714 real_points[i],
715 inverse_approximation.compute(real_points[i]),
716 support_points,
717 polynomials_1d,
718 renumber_lexicographic_to_hierarchic);
719}
720
721
722
723template <int dim, int spacedim>
726{
727 // add flags if the respective quantities are necessary to compute
728 // what we need. note that some flags appear in both the conditions
729 // and in subsequent set operations. this leads to some circular
730 // logic. the only way to treat this is to iterate. since there are
731 // 5 if-clauses in the loop, it will take at most 5 iterations to
732 // converge. do them:
733 UpdateFlags out = in;
734 for (unsigned int i = 0; i < 5; ++i)
735 {
736 // The following is a little incorrect:
737 // If not applied on a face,
738 // update_boundary_forms does not
739 // make sense. On the other hand,
740 // it is necessary on a
741 // face. Currently,
742 // update_boundary_forms is simply
743 // ignored for the interior of a
744 // cell.
747
748 if (out &
752
753 if (out &
758
759 // The contravariant transformation is used in the Piola
760 // transformation, which requires the determinant of the Jacobi
761 // matrix of the transformation. Because we have no way of
762 // knowing here whether the finite element wants to use the
763 // contravariant or the Piola transforms, we add the JxW values
764 // to the list of flags to be updated for each cell.
767
768 // the same is true when computing normal vectors: they require
769 // the determinant of the Jacobian
770 if (out & update_normal_vectors)
772 }
773
774 return out;
775}
776
777
778
779template <int dim, int spacedim>
780std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
782 const Quadrature<dim> &q) const
783{
784 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
785 std::make_unique<InternalData>(polynomial_degree);
786 data_ptr->reinit(update_flags, q);
787 return data_ptr;
788}
789
790
791
792template <int dim, int spacedim>
793std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
795 const UpdateFlags update_flags,
796 const hp::QCollection<dim - 1> &quadrature) const
797{
798 AssertDimension(quadrature.size(), 1);
799
800 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
801 std::make_unique<InternalData>(polynomial_degree);
802 auto &data = dynamic_cast<InternalData &>(*data_ptr);
803 data.initialize_face(this->requires_update_flags(update_flags),
805 ReferenceCells::get_hypercube<dim>(), quadrature[0]),
806 quadrature[0].size());
807
808 return data_ptr;
809}
810
811
812
813template <int dim, int spacedim>
814std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
816 const UpdateFlags update_flags,
817 const Quadrature<dim - 1> &quadrature) const
818{
819 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
820 std::make_unique<InternalData>(polynomial_degree);
821 auto &data = dynamic_cast<InternalData &>(*data_ptr);
822 data.initialize_face(this->requires_update_flags(update_flags),
824 ReferenceCells::get_hypercube<dim>(), quadrature),
825 quadrature.size());
826
827 return data_ptr;
828}
829
830
831
832template <int dim, int spacedim>
836 const CellSimilarity::Similarity cell_similarity,
837 const Quadrature<dim> &quadrature,
838 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
840 &output_data) const
841{
842 // ensure that the following static_cast is really correct:
843 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
845 const InternalData &data = static_cast<const InternalData &>(internal_data);
846 data.output_data = &output_data;
847
848 const unsigned int n_q_points = quadrature.size();
849
850 // recompute the support points of the transformation of this
851 // cell. we tried to be clever here in an earlier version of the
852 // library by checking whether the cell is the same as the one we
853 // had visited last, but it turns out to be difficult to determine
854 // that because a cell for the purposes of a mapping is
855 // characterized not just by its (triangulation, level, index)
856 // triple, but also by the locations of its vertices, the manifold
857 // object attached to the cell and all of its bounding faces/edges,
858 // etc. to reliably test that the "cell" we are on is, therefore,
859 // not easily done
860 if (polynomial_degree == 1)
861 {
863 const auto vertices = this->get_vertices(cell);
864 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
865 data.mapping_support_points[i] = vertices[i];
866 }
867 else
868 data.mapping_support_points = this->compute_mapping_support_points(cell);
869
871
872 // if the order of the mapping is greater than 1, then do not reuse any cell
873 // similarity information. This is necessary because the cell similarity
874 // value is computed with just cell vertices and does not take into account
875 // cell curvature.
876 const CellSimilarity::Similarity computed_cell_similarity =
877 (polynomial_degree == 1 && this->preserves_vertex_locations() ?
878 cell_similarity :
880
881 if (dim > 1 && data.tensor_product_quadrature)
882 {
883 internal::MappingQImplementation::
884 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
885 computed_cell_similarity,
886 data,
887 output_data.quadrature_points,
888 output_data.jacobians,
889 output_data.inverse_jacobians,
890 output_data.jacobian_grads);
891 }
892 else
893 {
895 computed_cell_similarity,
896 data,
897 make_array_view(quadrature.get_points()),
898 polynomials_1d,
899 renumber_lexicographic_to_hierarchic,
900 output_data.quadrature_points,
901 output_data.jacobians,
902 output_data.inverse_jacobians);
903
905 spacedim>(
906 computed_cell_similarity,
907 data,
908 make_array_view(quadrature.get_points()),
909 polynomials_1d,
910 renumber_lexicographic_to_hierarchic,
911 output_data.jacobian_grads);
912 }
913
915 dim,
916 spacedim>(computed_cell_similarity,
917 data,
918 make_array_view(quadrature.get_points()),
919 polynomials_1d,
920 renumber_lexicographic_to_hierarchic,
922
924 dim,
925 spacedim>(computed_cell_similarity,
926 data,
927 make_array_view(quadrature.get_points()),
928 polynomials_1d,
929 renumber_lexicographic_to_hierarchic,
930 output_data.jacobian_2nd_derivatives);
931
932 internal::MappingQImplementation::
933 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
934 computed_cell_similarity,
935 data,
936 make_array_view(quadrature.get_points()),
937 polynomials_1d,
938 renumber_lexicographic_to_hierarchic,
940
942 dim,
943 spacedim>(computed_cell_similarity,
944 data,
945 make_array_view(quadrature.get_points()),
946 polynomials_1d,
947 renumber_lexicographic_to_hierarchic,
948 output_data.jacobian_3rd_derivatives);
949
950 internal::MappingQImplementation::
951 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
952 computed_cell_similarity,
953 data,
954 make_array_view(quadrature.get_points()),
955 polynomials_1d,
956 renumber_lexicographic_to_hierarchic,
958
959 const UpdateFlags update_flags = data.update_each;
960 const std::vector<double> &weights = quadrature.get_weights();
961
962 // Multiply quadrature weights by absolute value of Jacobian determinants or
963 // the area element g=sqrt(DX^t DX) in case of codim > 0
964
965 if (update_flags & (update_normal_vectors | update_JxW_values))
966 {
967 Assert(!(update_flags & update_JxW_values) ||
968 (output_data.JxW_values.size() == n_q_points),
969 ExcDimensionMismatch(output_data.JxW_values.size(), n_q_points));
970
971 Assert(!(update_flags & update_normal_vectors) ||
972 (output_data.normal_vectors.size() == n_q_points),
973 ExcDimensionMismatch(output_data.normal_vectors.size(),
974 n_q_points));
975
976
977 if (computed_cell_similarity != CellSimilarity::translation)
978 for (unsigned int point = 0; point < n_q_points; ++point)
979 {
980 if (dim == spacedim)
981 {
982 const double det = data.volume_elements[point];
983
984 // check for distorted cells.
985
986 // TODO: this allows for anisotropies of up to 1e6 in 3d and
987 // 1e12 in 2d. might want to find a finer
988 // (dimension-independent) criterion
989 Assert(det >
990 1e-12 * Utilities::fixed_power<dim>(
991 cell->diameter() / std::sqrt(double(dim))),
993 cell->center(), det, point)));
994
995 output_data.JxW_values[point] = weights[point] * det;
996 }
997 // if dim==spacedim, then there is no cell normal to
998 // compute. since this is for FEValues (and not FEFaceValues),
999 // there are also no face normals to compute
1000 else // codim>0 case
1001 {
1002 Tensor<1, spacedim> DX_t[dim];
1003 for (unsigned int i = 0; i < spacedim; ++i)
1004 for (unsigned int j = 0; j < dim; ++j)
1005 DX_t[j][i] = output_data.jacobians[point][i][j];
1006
1007 Tensor<2, dim> G; // First fundamental form
1008 for (unsigned int i = 0; i < dim; ++i)
1009 for (unsigned int j = 0; j < dim; ++j)
1010 G[i][j] = DX_t[i] * DX_t[j];
1011
1012 if (update_flags & update_JxW_values)
1013 output_data.JxW_values[point] =
1014 std::sqrt(determinant(G)) * weights[point];
1015
1016 if (computed_cell_similarity ==
1018 {
1019 // we only need to flip the normal
1020 if (update_flags & update_normal_vectors)
1021 output_data.normal_vectors[point] *= -1.;
1022 }
1023 else
1024 {
1025 if (update_flags & update_normal_vectors)
1026 {
1027 Assert(spacedim == dim + 1,
1028 ExcMessage(
1029 "There is no (unique) cell normal for " +
1031 "-dimensional cells in " +
1032 Utilities::int_to_string(spacedim) +
1033 "-dimensional space. This only works if the "
1034 "space dimension is one greater than the "
1035 "dimensionality of the mesh cells."));
1036
1037 if (dim == 1)
1038 output_data.normal_vectors[point] =
1039 cross_product_2d(-DX_t[0]);
1040 else // dim == 2
1041 output_data.normal_vectors[point] =
1042 cross_product_3d(DX_t[0], DX_t[1]);
1043
1044 output_data.normal_vectors[point] /=
1045 output_data.normal_vectors[point].norm();
1046
1047 if (cell->direction_flag() == false)
1048 output_data.normal_vectors[point] *= -1.;
1049 }
1050 }
1051 } // codim>0 case
1052 }
1053 }
1054
1055 return computed_cell_similarity;
1056}
1057
1058
1059
1060template <int dim, int spacedim>
1061void
1064 const unsigned int face_no,
1065 const hp::QCollection<dim - 1> &quadrature,
1066 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1068 &output_data) const
1069{
1070 AssertDimension(quadrature.size(), 1);
1071
1072 // ensure that the following cast is really correct:
1073 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1075 const InternalData &data = static_cast<const InternalData &>(internal_data);
1076 data.output_data = &output_data;
1077
1078 // if necessary, recompute the support points of the transformation of this
1079 // cell (note that we need to first check the triangulation pointer, since
1080 // otherwise the second test might trigger an exception if the
1081 // triangulations are not the same)
1082 if ((data.mapping_support_points.empty()) ||
1083 (&cell->get_triangulation() !=
1085 (cell != data.cell_of_current_support_points))
1086 {
1087 if (polynomial_degree == 1)
1088 {
1089 data.mapping_support_points.resize(
1091 const auto vertices = this->get_vertices(cell);
1092 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
1093 ++i)
1094 data.mapping_support_points[i] = vertices[i];
1095 }
1096 else
1098 this->compute_mapping_support_points(cell);
1100 }
1101
1103 *this,
1104 cell,
1105 face_no,
1108 ReferenceCells::get_hypercube<dim>(),
1109 face_no,
1110 cell->face_orientation(face_no),
1111 cell->face_flip(face_no),
1112 cell->face_rotation(face_no),
1113 quadrature[0].size()),
1114 quadrature[0],
1115 data,
1116 polynomials_1d,
1117 renumber_lexicographic_to_hierarchic,
1118 output_data);
1119}
1120
1121
1122
1123template <int dim, int spacedim>
1124void
1127 const unsigned int face_no,
1128 const unsigned int subface_no,
1129 const Quadrature<dim - 1> &quadrature,
1130 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1132 &output_data) const
1133{
1134 // ensure that the following cast is really correct:
1135 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1137 const InternalData &data = static_cast<const InternalData &>(internal_data);
1138 data.output_data = &output_data;
1139
1140 // if necessary, recompute the support points of the transformation of this
1141 // cell (note that we need to first check the triangulation pointer, since
1142 // otherwise the second test might trigger an exception if the
1143 // triangulations are not the same)
1144 if ((data.mapping_support_points.empty()) ||
1145 (&cell->get_triangulation() !=
1147 (cell != data.cell_of_current_support_points))
1148 {
1149 if (polynomial_degree == 1)
1150 {
1151 data.mapping_support_points.resize(
1153 const auto vertices = this->get_vertices(cell);
1154 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
1155 ++i)
1156 data.mapping_support_points[i] = vertices[i];
1157 }
1158 else
1160 this->compute_mapping_support_points(cell);
1162 }
1163
1165 *this,
1166 cell,
1167 face_no,
1168 subface_no,
1170 ReferenceCells::get_hypercube<dim>(),
1171 face_no,
1172 subface_no,
1173 cell->combined_face_orientation(face_no),
1174 quadrature.size(),
1175 cell->subface_case(face_no)),
1176 quadrature,
1177 data,
1178 polynomials_1d,
1179 renumber_lexicographic_to_hierarchic,
1180 output_data);
1181}
1182
1183
1184
1185template <int dim, int spacedim>
1186void
1190 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1192 &output_data) const
1193{
1194 Assert(dim == spacedim, ExcNotImplemented());
1195
1196 // ensure that the following static_cast is really correct:
1197 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1199 const InternalData &data = static_cast<const InternalData &>(internal_data);
1200 data.output_data = &output_data;
1201
1202 const unsigned int n_q_points = quadrature.size();
1203
1204 if (polynomial_degree == 1)
1205 {
1207 const auto vertices = this->get_vertices(cell);
1208 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1209 data.mapping_support_points[i] = vertices[i];
1210 }
1211 else
1212 data.mapping_support_points = this->compute_mapping_support_points(cell);
1214
1217 data,
1218 make_array_view(quadrature.get_points()),
1219 polynomials_1d,
1220 renumber_lexicographic_to_hierarchic,
1221 output_data.quadrature_points,
1222 output_data.jacobians,
1223 output_data.inverse_jacobians);
1224
1225 internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1227 data,
1228 make_array_view(quadrature.get_points()),
1229 polynomials_1d,
1230 renumber_lexicographic_to_hierarchic,
1231 output_data.jacobian_grads);
1232
1234 dim,
1235 spacedim>(CellSimilarity::none,
1236 data,
1237 make_array_view(quadrature.get_points()),
1238 polynomials_1d,
1239 renumber_lexicographic_to_hierarchic,
1240 output_data.jacobian_pushed_forward_grads);
1241
1243 dim,
1244 spacedim>(CellSimilarity::none,
1245 data,
1246 make_array_view(quadrature.get_points()),
1247 polynomials_1d,
1248 renumber_lexicographic_to_hierarchic,
1249 output_data.jacobian_2nd_derivatives);
1250
1251 internal::MappingQImplementation::
1252 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1254 data,
1255 make_array_view(quadrature.get_points()),
1256 polynomials_1d,
1257 renumber_lexicographic_to_hierarchic,
1259
1261 dim,
1262 spacedim>(CellSimilarity::none,
1263 data,
1264 make_array_view(quadrature.get_points()),
1265 polynomials_1d,
1266 renumber_lexicographic_to_hierarchic,
1267 output_data.jacobian_3rd_derivatives);
1268
1269 internal::MappingQImplementation::
1270 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1272 data,
1273 make_array_view(quadrature.get_points()),
1274 polynomials_1d,
1275 renumber_lexicographic_to_hierarchic,
1277
1278 const UpdateFlags update_flags = data.update_each;
1279 const std::vector<double> &weights = quadrature.get_weights();
1280
1281 if ((update_flags & (update_normal_vectors | update_JxW_values)) != 0u)
1282 {
1283 AssertDimension(output_data.JxW_values.size(), n_q_points);
1284
1285 Assert(!(update_flags & update_normal_vectors) ||
1286 (output_data.normal_vectors.size() == n_q_points),
1287 ExcDimensionMismatch(output_data.normal_vectors.size(),
1288 n_q_points));
1289
1290
1291 for (unsigned int point = 0; point < n_q_points; ++point)
1292 {
1293 const double det = data.volume_elements[point];
1294
1295 // check for distorted cells.
1296
1297 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1298 // 1e12 in 2d. might want to find a finer
1299 // (dimension-independent) criterion
1300 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1301 cell->diameter() / std::sqrt(double(dim))),
1303 cell->center(), det, point)));
1304
1305 // The normals are n = J^{-T} * \hat{n} before normalizing.
1306 Tensor<1, spacedim> normal;
1307 for (unsigned int d = 0; d < spacedim; d++)
1308 normal[d] = output_data.inverse_jacobians[point].transpose()[d] *
1309 quadrature.normal_vector(point);
1310
1311 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1312
1313 if ((update_flags & update_normal_vectors) != 0u)
1314 {
1315 normal /= normal.norm();
1316 output_data.normal_vectors[point] = normal;
1317 }
1318 }
1319 }
1320}
1321
1322
1323
1324template <int dim, int spacedim>
1325void
1328 const ArrayView<const Point<dim>> &unit_points,
1329 const UpdateFlags update_flags,
1331 &output_data) const
1332{
1333 if (update_flags == update_default)
1334 return;
1335
1336 Assert(update_flags & update_inverse_jacobians ||
1337 update_flags & update_jacobians ||
1338 update_flags & update_quadrature_points,
1340
1341 output_data.initialize(unit_points.size(), update_flags);
1342
1343 auto internal_data =
1344 this->get_data(update_flags,
1345 Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
1346 unit_points.end())));
1347 const InternalData &data = static_cast<const InternalData &>(*internal_data);
1348 data.output_data = &output_data;
1349 if (polynomial_degree == 1)
1350 {
1352 const auto vertices = this->get_vertices(cell);
1353 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1354 data.mapping_support_points[i] = vertices[i];
1355 }
1356 else
1357 data.mapping_support_points = this->compute_mapping_support_points(cell);
1358
1361 data,
1362 unit_points,
1363 polynomials_1d,
1364 renumber_lexicographic_to_hierarchic,
1365 output_data.quadrature_points,
1366 output_data.jacobians,
1367 output_data.inverse_jacobians);
1368}
1369
1370
1371
1372template <int dim, int spacedim>
1373void
1376 const unsigned int face_no,
1377 const Quadrature<dim - 1> &face_quadrature,
1378 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1380 &output_data) const
1381{
1382 if (face_quadrature.get_points().empty())
1383 return;
1384
1385 // ensure that the following static_cast is really correct:
1386 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1388 const InternalData &data = static_cast<const InternalData &>(internal_data);
1389
1390 if (polynomial_degree == 1)
1391 {
1393 const auto vertices = this->get_vertices(cell);
1394 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1395 data.mapping_support_points[i] = vertices[i];
1396 }
1397 else
1398 data.mapping_support_points = this->compute_mapping_support_points(cell);
1399 data.output_data = &output_data;
1400
1402 *this,
1403 cell,
1404 face_no,
1407 face_quadrature,
1408 data,
1409 polynomials_1d,
1410 renumber_lexicographic_to_hierarchic,
1411 output_data);
1412}
1413
1414
1415
1416template <int dim, int spacedim>
1417void
1419 const ArrayView<const Tensor<1, dim>> &input,
1420 const MappingKind mapping_kind,
1421 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1422 const ArrayView<Tensor<1, spacedim>> &output) const
1423{
1425 mapping_kind,
1426 mapping_data,
1427 output);
1428}
1429
1430
1431
1432template <int dim, int spacedim>
1433void
1435 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1436 const MappingKind mapping_kind,
1437 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1438 const ArrayView<Tensor<2, spacedim>> &output) const
1439{
1441 mapping_kind,
1442 mapping_data,
1443 output);
1444}
1445
1446
1447
1448template <int dim, int spacedim>
1449void
1451 const ArrayView<const Tensor<2, dim>> &input,
1452 const MappingKind mapping_kind,
1453 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1454 const ArrayView<Tensor<2, spacedim>> &output) const
1455{
1456 switch (mapping_kind)
1457 {
1460 mapping_kind,
1461 mapping_data,
1462 output);
1463 return;
1464
1469 mapping_kind,
1470 mapping_data,
1471 output);
1472 return;
1473 default:
1475 }
1476}
1477
1478
1479
1480template <int dim, int spacedim>
1481void
1483 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1484 const MappingKind mapping_kind,
1485 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1486 const ArrayView<Tensor<3, spacedim>> &output) const
1487{
1488 AssertDimension(input.size(), output.size());
1489 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1492 &data = *static_cast<const InternalData &>(mapping_data).output_data;
1493
1494 switch (mapping_kind)
1495 {
1497 {
1498 Assert(!data.inverse_jacobians.empty(),
1500 "update_covariant_transformation"));
1501
1502 for (unsigned int q = 0; q < output.size(); ++q)
1503 for (unsigned int i = 0; i < spacedim; ++i)
1504 for (unsigned int j = 0; j < spacedim; ++j)
1505 {
1506 double tmp[dim];
1507 const DerivativeForm<1, dim, spacedim> covariant =
1508 data.inverse_jacobians[q].transpose();
1509 for (unsigned int K = 0; K < dim; ++K)
1510 {
1511 tmp[K] = covariant[j][0] * input[q][i][0][K];
1512 for (unsigned int J = 1; J < dim; ++J)
1513 tmp[K] += covariant[j][J] * input[q][i][J][K];
1514 }
1515 for (unsigned int k = 0; k < spacedim; ++k)
1516 {
1517 output[q][i][j][k] = covariant[k][0] * tmp[0];
1518 for (unsigned int K = 1; K < dim; ++K)
1519 output[q][i][j][k] += covariant[k][K] * tmp[K];
1520 }
1521 }
1522 return;
1523 }
1524
1525 default:
1527 }
1528}
1529
1530
1531
1532template <int dim, int spacedim>
1533void
1535 const ArrayView<const Tensor<3, dim>> &input,
1536 const MappingKind mapping_kind,
1537 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1538 const ArrayView<Tensor<3, spacedim>> &output) const
1539{
1540 switch (mapping_kind)
1541 {
1546 mapping_kind,
1547 mapping_data,
1548 output);
1549 return;
1550 default:
1552 }
1553}
1554
1555
1556
1557template <int dim, int spacedim>
1558void
1561 std::vector<Point<spacedim>> &a) const
1562{
1563 // if we only need the midpoint, then ask for it.
1564 if (this->polynomial_degree == 2)
1565 {
1566 for (unsigned int line_no = 0;
1567 line_no < GeometryInfo<dim>::lines_per_cell;
1568 ++line_no)
1569 {
1571 (dim == 1 ?
1572 static_cast<
1574 cell->line(line_no));
1575
1576 const Manifold<dim, spacedim> &manifold =
1577 ((line->manifold_id() == numbers::flat_manifold_id) &&
1578 (dim < spacedim) ?
1579 cell->get_manifold() :
1580 line->get_manifold());
1581 a.push_back(manifold.get_new_point_on_line(line));
1582 }
1583 }
1584 else
1585 // otherwise call the more complicated functions and ask for inner points
1586 // from the manifold description
1587 {
1588 std::vector<Point<spacedim>> tmp_points;
1589 for (unsigned int line_no = 0;
1590 line_no < GeometryInfo<dim>::lines_per_cell;
1591 ++line_no)
1592 {
1594 (dim == 1 ?
1595 static_cast<
1597 cell->line(line_no));
1598
1599 const Manifold<dim, spacedim> &manifold =
1600 ((line->manifold_id() == numbers::flat_manifold_id) &&
1601 (dim < spacedim) ?
1602 cell->get_manifold() :
1603 line->get_manifold());
1604
1605 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1606 const std::array<Point<spacedim>, 2> vertices{
1607 {cell->vertex(reference_cell.line_to_cell_vertices(line_no, 0)),
1608 cell->vertex(reference_cell.line_to_cell_vertices(line_no, 1))}};
1609
1610 const std::size_t n_rows =
1611 support_point_weights_perimeter_to_interior[0].size(0);
1612 a.resize(a.size() + n_rows);
1613 auto a_view = make_array_view(a.end() - n_rows, a.end());
1614 manifold.get_new_points(
1615 make_array_view(vertices.begin(), vertices.end()),
1616 support_point_weights_perimeter_to_interior[0],
1617 a_view);
1618 }
1619 }
1620}
1621
1622
1623
1624template <>
1625void
1628 std::vector<Point<3>> &a) const
1629{
1630 const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1631
1632 // used if face quad at boundary or entirely in the interior of the domain
1633 std::vector<Point<3>> tmp_points;
1634
1635 // loop over all faces and collect points on them
1636 for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1637 {
1638 const Triangulation<3>::face_iterator face = cell->face(face_no);
1639
1640#ifdef DEBUG
1641 const bool face_orientation = cell->face_orientation(face_no),
1642 face_flip = cell->face_flip(face_no),
1643 face_rotation = cell->face_rotation(face_no);
1644 const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1645 lines_per_face = GeometryInfo<3>::lines_per_face;
1646
1647 // some sanity checks up front
1648 for (unsigned int i = 0; i < vertices_per_face; ++i)
1649 Assert(face->vertex_index(i) ==
1650 cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1651 face_no, i, face_orientation, face_flip, face_rotation)),
1653
1654 // indices of the lines that bound a face are given by GeometryInfo<3>::
1655 // face_to_cell_lines
1656 for (unsigned int i = 0; i < lines_per_face; ++i)
1657 Assert(face->line(i) ==
1659 face_no, i, face_orientation, face_flip, face_rotation)),
1661#endif
1662 // extract the points surrounding a quad from the points
1663 // already computed. First get the 4 vertices and then the points on
1664 // the four lines
1665 boost::container::small_vector<Point<3>, 200> tmp_points(
1667 GeometryInfo<2>::lines_per_cell * (polynomial_degree - 1));
1668 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1669 tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1670 if (polynomial_degree > 1)
1671 for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1672 ++line)
1673 for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1674 tmp_points[4 + line * (polynomial_degree - 1) + i] =
1676 (polynomial_degree - 1) *
1678 i];
1679
1680 const std::size_t n_rows =
1681 support_point_weights_perimeter_to_interior[1].size(0);
1682 a.resize(a.size() + n_rows);
1683 auto a_view = make_array_view(a.end() - n_rows, a.end());
1684 face->get_manifold().get_new_points(
1685 make_array_view(tmp_points.begin(), tmp_points.end()),
1686 support_point_weights_perimeter_to_interior[1],
1687 a_view);
1688 }
1689}
1690
1691
1692
1693template <>
1694void
1697 std::vector<Point<3>> &a) const
1698{
1699 std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1700 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1701 vertices[i] = cell->vertex(i);
1702
1703 Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1705 for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1706 for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1707 {
1708 Point<2> point(line_support_points[q1 + 1][0],
1709 line_support_points[q2 + 1][0]);
1710 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1711 weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
1712 }
1713
1714 const std::size_t n_rows = weights.size(0);
1715 a.resize(a.size() + n_rows);
1716 auto a_view = make_array_view(a.end() - n_rows, a.end());
1717 cell->get_manifold().get_new_points(
1718 make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1719}
1720
1721
1722
1723template <int dim, int spacedim>
1724void
1731
1732
1733
1734template <int dim, int spacedim>
1735std::vector<Point<spacedim>>
1737 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1738{
1739 // get the vertices first
1740 std::vector<Point<spacedim>> a;
1741 a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1742 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1743 a.push_back(cell->vertex(i));
1744
1745 if (this->polynomial_degree > 1)
1746 {
1747 // check if all entities have the same manifold id which is when we can
1748 // simply ask the manifold for all points. the transfinite manifold can
1749 // do the interpolation better than this class, so if we detect that we
1750 // do not have to change anything here
1751 Assert(dim <= 3, ExcImpossibleInDim(dim));
1752 bool all_manifold_ids_are_equal = (dim == spacedim);
1753 if (all_manifold_ids_are_equal &&
1755 &cell->get_manifold()) == nullptr)
1756 {
1757 for (auto f : GeometryInfo<dim>::face_indices())
1758 if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1759 all_manifold_ids_are_equal = false;
1760
1761 if (dim == 3)
1762 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1763 if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1764 all_manifold_ids_are_equal = false;
1765 }
1766
1767 if (all_manifold_ids_are_equal)
1768 {
1769 const std::size_t n_rows = support_point_weights_cell.size(0);
1770 a.resize(a.size() + n_rows);
1771 auto a_view = make_array_view(a.end() - n_rows, a.end());
1772 cell->get_manifold().get_new_points(make_array_view(a.begin(),
1773 a.end() - n_rows),
1774 support_point_weights_cell,
1775 a_view);
1776 }
1777 else
1778 switch (dim)
1779 {
1780 case 1:
1781 add_line_support_points(cell, a);
1782 break;
1783 case 2:
1784 // in 2d, add the points on the four bounding lines to the
1785 // exterior (outer) points
1786 add_line_support_points(cell, a);
1787
1788 // then get the interior support points
1789 if (dim != spacedim)
1790 add_quad_support_points(cell, a);
1791 else
1792 {
1793 const std::size_t n_rows =
1794 support_point_weights_perimeter_to_interior[1].size(0);
1795 a.resize(a.size() + n_rows);
1796 auto a_view = make_array_view(a.end() - n_rows, a.end());
1797 cell->get_manifold().get_new_points(
1798 make_array_view(a.begin(), a.end() - n_rows),
1799 support_point_weights_perimeter_to_interior[1],
1800 a_view);
1801 }
1802 break;
1803
1804 case 3:
1805 // in 3d also add the points located on the boundary faces
1806 add_line_support_points(cell, a);
1807 add_quad_support_points(cell, a);
1808
1809 // then compute the interior points
1810 {
1811 const std::size_t n_rows =
1812 support_point_weights_perimeter_to_interior[2].size(0);
1813 a.resize(a.size() + n_rows);
1814 auto a_view = make_array_view(a.end() - n_rows, a.end());
1815 cell->get_manifold().get_new_points(
1816 make_array_view(a.begin(), a.end() - n_rows),
1817 support_point_weights_perimeter_to_interior[2],
1818 a_view);
1819 }
1820 break;
1821
1822 default:
1824 break;
1825 }
1826 }
1827
1828 return a;
1829}
1830
1831
1832
1833template <int dim, int spacedim>
1836 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1837{
1838 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
1839}
1840
1841
1842
1843template <int dim, int spacedim>
1844bool
1846 const ReferenceCell &reference_cell) const
1847{
1848 Assert(dim == reference_cell.get_dimension(),
1849 ExcMessage("The dimension of your mapping (" +
1851 ") and the reference cell cell_type (" +
1852 Utilities::to_string(reference_cell.get_dimension()) +
1853 " ) do not agree."));
1854
1855 return reference_cell.is_hyper_cube();
1856}
1857
1858
1859
1860//--------------------------- Explicit instantiations -----------------------
1861#include "mapping_q.inst"
1862
1863
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
auto make_const_array_view(const Container &container) -> decltype(make_array_view(container))
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:436
virtual std::size_t memory_consumption() const override
Definition mapping_q.cc:63
std::vector< Point< spacedim > > mapping_support_points
Definition mapping_q.h:430
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:450
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
AlignedVector< double > volume_elements
Definition mapping_q.h:442
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition mapping_q.h:548
const Table< 2, double > support_point_weights_cell
Definition mapping_q.h:596
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:781
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:815
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:834
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:524
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:633
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:794
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:328
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition mapping_q.h:582
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:534
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:287
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition mapping_q.h:541
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition mapping_q.cc:269
const std::vector< Point< dim > > unit_cell_support_points
Definition mapping_q.h:560
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:513
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:725
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:278
Abstract base class for mapping classes.
Definition mapping.h:318
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
#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)
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1637
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:79
@ mapping_covariant_gradient
Definition mapping.h:100
@ mapping_contravariant
Definition mapping.h:94
@ mapping_contravariant_hessian
Definition mapping.h:156
@ mapping_covariant_hessian
Definition mapping.h:150
@ mapping_contravariant_gradient
Definition mapping.h:106
@ mapping_piola_gradient
Definition mapping.h:120
@ mapping_piola_hessian
Definition mapping.h:162
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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: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:966
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={})
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 > &)