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