Reference documentation for deal.II version 9.4.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\}}\)
mapping_q.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2022 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
38#include <boost/container/small_vector.hpp>
40
41#include <algorithm>
42#include <array>
43#include <cmath>
44#include <memory>
45#include <numeric>
46
47
49
50
51template <int dim, int spacedim>
53 const unsigned int polynomial_degree)
54 : polynomial_degree(polynomial_degree)
55 , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
56 , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
57 , tensor_product_quadrature(false)
58{}
59
60
61
62template <int dim, int spacedim>
63std::size_t
65{
66 return (
69 MemoryConsumption::memory_consumption(shape_derivatives) +
74 MemoryConsumption::memory_consumption(mapping_support_points) +
75 MemoryConsumption::memory_consumption(cell_of_current_support_points) +
78 MemoryConsumption::memory_consumption(n_shape_functions));
79}
80
81
82
83template <int dim, int spacedim>
84void
86 const UpdateFlags update_flags,
87 const Quadrature<dim> &q,
88 const unsigned int n_original_q_points)
89{
90 // store the flags in the internal data object so we can access them
91 // in fill_fe_*_values()
92 this->update_each = update_flags;
93
94 const unsigned int n_q_points = q.size();
95
96 const bool needs_higher_order_terms =
97 this->update_each &
102
103 if (this->update_each & update_covariant_transformation)
104 covariant.resize(n_original_q_points);
105
106 if (this->update_each & update_contravariant_transformation)
107 contravariant.resize(n_original_q_points);
108
109 if (this->update_each & update_volume_elements)
110 volume_elements.resize(n_original_q_points);
111
112 tensor_product_quadrature = q.is_tensor_product();
113
114 // use of MatrixFree only for higher order elements and with more than one
115 // point where tensor products do not make sense
116 if (polynomial_degree < 2 || n_q_points == 1)
117 tensor_product_quadrature = false;
119 if (dim > 1)
120 {
121 // find out if the one-dimensional formula is the same
122 // in all directions
123 if (tensor_product_quadrature)
124 {
125 const std::array<Quadrature<1>, dim> quad_array =
127 for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
128 {
129 if (quad_array[i - 1].size() != quad_array[i].size())
130 {
131 tensor_product_quadrature = false;
132 break;
133 }
134 else
135 {
136 const std::vector<Point<1>> &points_1 =
137 quad_array[i - 1].get_points();
138 const std::vector<Point<1>> &points_2 =
139 quad_array[i].get_points();
140 const std::vector<double> &weights_1 =
141 quad_array[i - 1].get_weights();
142 const std::vector<double> &weights_2 =
143 quad_array[i].get_weights();
144 for (unsigned int j = 0; j < quad_array[i].size(); ++j)
145 {
146 if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
147 std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
148 {
149 tensor_product_quadrature = false;
150 break;
151 }
152 }
153 }
155
156 if (tensor_product_quadrature)
157 {
158 // use a 1D FE_DGQ and adjust the hierarchic -> lexicographic
159 // numbering manually (building an FE_Q<dim> is relatively
160 // expensive due to constraints)
162 shape_info.reinit(q.get_tensor_basis()[0], fe);
163 shape_info.lexicographic_numbering =
164 FETools::lexicographic_to_hierarchic_numbering<dim>(
166 shape_info.n_q_points = q.size();
167 shape_info.dofs_per_component_on_cell =
169 }
170 }
171 }
172
173 // Only fill the big arrays on demand in case we cannot use the tensor
174 // product quadrature code path
175 if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
176 {
177 // see if we need the (transformation) shape function values
178 // and/or gradients and resize the necessary arrays
179 if (this->update_each & update_quadrature_points)
180 shape_values.resize(n_shape_functions * n_q_points);
181
182 if (this->update_each &
192 shape_derivatives.resize(n_shape_functions * n_q_points);
193
194 if (this->update_each &
196 shape_second_derivatives.resize(n_shape_functions * n_q_points);
197
198 if (this->update_each & (update_jacobian_2nd_derivatives |
200 shape_third_derivatives.resize(n_shape_functions * n_q_points);
201
202 if (this->update_each & (update_jacobian_3rd_derivatives |
204 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
205
206 // now also fill the various fields with their correct values
207 compute_shape_function_values(q.get_points());
208 }
210
211
212
213template <int dim, int spacedim>
214void
216 const UpdateFlags update_flags,
217 const Quadrature<dim> &q,
218 const unsigned int n_original_q_points)
219{
220 initialize(update_flags, q, n_original_q_points);
221
222 if (dim > 1 && tensor_product_quadrature)
224 constexpr unsigned int facedim = dim - 1;
226 shape_info.reinit(q.get_tensor_basis()[0], fe);
227 shape_info.lexicographic_numbering =
228 FETools::lexicographic_to_hierarchic_numbering<facedim>(
230 shape_info.n_q_points = n_original_q_points;
231 shape_info.dofs_per_component_on_cell =
233 }
234
235 if (dim > 1)
236 {
237 if (this->update_each &
240 {
241 aux.resize(dim - 1,
242 AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
243
244 // Compute tangentials to the unit cell.
245 for (const unsigned int i : GeometryInfo<dim>::face_indices())
246 {
247 unit_tangentials[i].resize(n_original_q_points);
248 std::fill(unit_tangentials[i].begin(),
249 unit_tangentials[i].end(),
251 if (dim > 2)
252 {
253 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
254 .resize(n_original_q_points);
255 std::fill(
256 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
257 .begin(),
258 unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
259 .end(),
261 }
262 }
263 }
264 }
265}
266
267
268
269template <int dim, int spacedim>
270void
272 const std::vector<Point<dim>> &unit_points)
273{
274 const unsigned int n_points = unit_points.size();
275
276 // Construct the tensor product polynomials used as shape functions for
277 // the Qp mapping of cells at the boundary.
278 const TensorProductPolynomials<dim> tensor_pols(
280 line_support_points.get_points()));
281 Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
282
283 // then also construct the mapping from lexicographic to the Qp shape
284 // function numbering
285 const std::vector<unsigned int> renumber =
286 FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
287
288 std::vector<double> values;
289 std::vector<Tensor<1, dim>> grads;
290 if (shape_values.size() != 0)
292 Assert(shape_values.size() == n_shape_functions * n_points,
294 values.resize(n_shape_functions);
295 }
296 if (shape_derivatives.size() != 0)
297 {
298 Assert(shape_derivatives.size() == n_shape_functions * n_points,
300 grads.resize(n_shape_functions);
302
303 std::vector<Tensor<2, dim>> grad2;
304 if (shape_second_derivatives.size() != 0)
305 {
306 Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
308 grad2.resize(n_shape_functions);
309 }
310
311 std::vector<Tensor<3, dim>> grad3;
312 if (shape_third_derivatives.size() != 0)
313 {
314 Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
316 grad3.resize(n_shape_functions);
317 }
318
319 std::vector<Tensor<4, dim>> grad4;
320 if (shape_fourth_derivatives.size() != 0)
322 Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
324 grad4.resize(n_shape_functions);
325 }
326
327
328 if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
329 shape_second_derivatives.size() != 0 ||
330 shape_third_derivatives.size() != 0 ||
331 shape_fourth_derivatives.size() != 0)
332 for (unsigned int point = 0; point < n_points; ++point)
333 {
334 tensor_pols.evaluate(
335 unit_points[point], values, grads, grad2, grad3, grad4);
336
337 if (shape_values.size() != 0)
338 for (unsigned int i = 0; i < n_shape_functions; ++i)
339 shape(point, i) = values[renumber[i]];
340
341 if (shape_derivatives.size() != 0)
342 for (unsigned int i = 0; i < n_shape_functions; ++i)
343 derivative(point, i) = grads[renumber[i]];
344
345 if (shape_second_derivatives.size() != 0)
346 for (unsigned int i = 0; i < n_shape_functions; ++i)
347 second_derivative(point, i) = grad2[renumber[i]];
348
349 if (shape_third_derivatives.size() != 0)
350 for (unsigned int i = 0; i < n_shape_functions; ++i)
351 third_derivative(point, i) = grad3[renumber[i]];
352
353 if (shape_fourth_derivatives.size() != 0)
354 for (unsigned int i = 0; i < n_shape_functions; ++i)
355 fourth_derivative(point, i) = grad4[renumber[i]];
356 }
357}
358
359
360
361template <int dim, int spacedim>
365 QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
371 internal::MappingQImplementation::unit_support_points<dim>(
375 internal::MappingQImplementation::
377 this->polynomial_degree,
378 dim))
380 internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
381 this->polynomial_degree))
382{
383 Assert(p >= 1,
384 ExcMessage("It only makes sense to create polynomial mappings "
385 "with a polynomial degree greater or equal to one."));
386}
387
388
389
390template <int dim, int spacedim>
391MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
392 : polynomial_degree(p)
393 , line_support_points(
394 QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
395 , polynomials_1d(
396 Polynomials::generate_complete_Lagrange_basis(line_support_points))
397 , renumber_lexicographic_to_hierarchic(
399 , unit_cell_support_points(
400 internal::MappingQImplementation::unit_support_points<dim>(
401 line_support_points,
402 renumber_lexicographic_to_hierarchic))
403 , support_point_weights_perimeter_to_interior(
404 internal::MappingQImplementation::
406 this->polynomial_degree,
407 dim))
408 , support_point_weights_cell(
409 internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
410 this->polynomial_degree))
411{
412 Assert(p >= 1,
413 ExcMessage("It only makes sense to create polynomial mappings "
414 "with a polynomial degree greater or equal to one."));
415}
416
417
418
419template <int dim, int spacedim>
421 : polynomial_degree(mapping.polynomial_degree)
422 , line_support_points(mapping.line_support_points)
423 , polynomials_1d(mapping.polynomials_1d)
424 , renumber_lexicographic_to_hierarchic(
425 mapping.renumber_lexicographic_to_hierarchic)
426 , support_point_weights_perimeter_to_interior(
427 mapping.support_point_weights_perimeter_to_interior)
428 , support_point_weights_cell(mapping.support_point_weights_cell)
429{}
430
431
432
433template <int dim, int spacedim>
434std::unique_ptr<Mapping<dim, spacedim>>
436{
437 return std::make_unique<MappingQ<dim, spacedim>>(*this);
438}
439
440
441
442template <int dim, int spacedim>
443unsigned int
445{
446 return polynomial_degree;
447}
448
449
450
451template <int dim, int spacedim>
455 const Point<dim> & p) const
456{
458 polynomials_1d,
459 this->compute_mapping_support_points(cell),
460 p,
461 polynomials_1d.size() == 2,
462 renumber_lexicographic_to_hierarchic)
463 .first);
464}
465
466
467// In the code below, GCC tries to instantiate MappingQ<3,4> when
468// seeing which of the overloaded versions of
469// do_transform_real_to_unit_cell_internal() to call. This leads to bad
470// error messages and, generally, nothing very good. Avoid this by ensuring
471// that this class exists, but does not have an inner InternalData
472// type, thereby ruling out the codim-1 version of the function
473// below when doing overload resolution.
474template <>
475class MappingQ<3, 4>
476{};
477
478
479
480// visual studio freaks out when trying to determine if
481// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
482// candidate. So instead of letting the compiler pick the correct overload, we
483// use template specialization to make sure we pick up the right function to
484// call:
485
486template <int dim, int spacedim>
490 const Point<spacedim> &,
491 const Point<dim> &) const
492{
493 // default implementation (should never be called)
494 Assert(false, ExcInternalError());
495 return {};
496}
497
498
499
500template <>
504 const Point<1> & p,
505 const Point<1> & initial_p_unit) const
506{
507 // dispatch to the various specializations for spacedim=dim,
508 // spacedim=dim+1, etc
509 return internal::MappingQImplementation::
510 do_transform_real_to_unit_cell_internal<1>(
511 p,
512 initial_p_unit,
513 this->compute_mapping_support_points(cell),
514 polynomials_1d,
515 renumber_lexicographic_to_hierarchic);
516}
517
518
519
520template <>
524 const Point<2> & p,
525 const Point<2> & initial_p_unit) const
526{
527 return internal::MappingQImplementation::
528 do_transform_real_to_unit_cell_internal<2>(
529 p,
530 initial_p_unit,
531 this->compute_mapping_support_points(cell),
532 polynomials_1d,
533 renumber_lexicographic_to_hierarchic);
534}
535
536
537
538template <>
542 const Point<3> & p,
543 const Point<3> & initial_p_unit) const
544{
545 return internal::MappingQImplementation::
546 do_transform_real_to_unit_cell_internal<3>(
547 p,
548 initial_p_unit,
549 this->compute_mapping_support_points(cell),
550 polynomials_1d,
551 renumber_lexicographic_to_hierarchic);
552}
553
554
555
556template <>
560 const Point<2> & p,
561 const Point<1> & initial_p_unit) const
562{
563 const int dim = 1;
564 const int spacedim = 2;
565
566 const Quadrature<dim> point_quadrature(initial_p_unit);
567
569 if (spacedim > dim)
570 update_flags |= update_jacobian_grads;
571 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
572 get_data(update_flags, point_quadrature));
573
574 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
575
576 // dispatch to the various specializations for spacedim=dim,
577 // spacedim=dim+1, etc
578 return internal::MappingQImplementation::
579 do_transform_real_to_unit_cell_internal_codim1<1>(cell,
580 p,
581 initial_p_unit,
582 *mdata);
583}
584
585
586
587template <>
591 const Point<3> & p,
592 const Point<2> & initial_p_unit) const
593{
594 const int dim = 2;
595 const int spacedim = 3;
596
597 const Quadrature<dim> point_quadrature(initial_p_unit);
598
600 if (spacedim > dim)
601 update_flags |= update_jacobian_grads;
602 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
603 get_data(update_flags, point_quadrature));
604
605 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
606
607 // dispatch to the various specializations for spacedim=dim,
608 // spacedim=dim+1, etc
609 return internal::MappingQImplementation::
610 do_transform_real_to_unit_cell_internal_codim1<2>(cell,
611 p,
612 initial_p_unit,
613 *mdata);
614}
615
616template <>
620 const Point<3> &,
621 const Point<1> &) const
622{
623 Assert(false, ExcNotImplemented());
624 return {};
625}
626
627
628
629template <int dim, int spacedim>
633 const Point<spacedim> & p) const
634{
635 // Use an exact formula if one is available. this is only the case
636 // for Q1 mappings in 1d, and in 2d if dim==spacedim
637 if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
638 ((dim == 1) || ((dim == 2) && (dim == spacedim))))
639 {
640 // The dimension-dependent algorithms are much faster (about 25-45x in
641 // 2D) but fail most of the time when the given point (p) is not in the
642 // cell. The dimension-independent Newton algorithm given below is
643 // slower, but more robust (though it still sometimes fails). Therefore
644 // this function implements the following strategy based on the
645 // p's dimension:
646 //
647 // * In 1D this mapping is linear, so the mapping is always invertible
648 // (and the exact formula is known) as long as the cell has non-zero
649 // length.
650 // * In 2D the exact (quadratic) formula is called first. If either the
651 // exact formula does not succeed (negative discriminant in the
652 // quadratic formula) or succeeds but finds a solution outside of the
653 // unit cell, then the Newton solver is called. The rationale for the
654 // second choice is that the exact formula may provide two different
655 // answers when mapping a point outside of the real cell, but the
656 // Newton solver (if it converges) will only return one answer.
657 // Otherwise the exact formula successfully found a point in the unit
658 // cell and that value is returned.
659 // * In 3D there is no (known to the authors) exact formula, so the Newton
660 // algorithm is used.
661 const auto vertices_ = this->get_vertices(cell);
662
663 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
664 vertices;
665 for (unsigned int i = 0; i < vertices.size(); ++i)
666 vertices[i] = vertices_[i];
667
668 try
669 {
670 switch (dim)
671 {
672 case 1:
673 {
674 // formula not subject to any issues in 1d
675 if (spacedim == 1)
677 vertices, p);
678 else
679 break;
680 }
681
682 case 2:
683 {
684 const Point<dim> point =
686 p);
687
688 // formula not guaranteed to work for points outside of
689 // the cell. only take the computed point if it lies
690 // inside the reference cell
691 const double eps = 1e-15;
692 if (-eps <= point(1) && point(1) <= 1 + eps &&
693 -eps <= point(0) && point(0) <= 1 + eps)
694 {
695 return point;
696 }
697 else
698 break;
699 }
700
701 default:
702 {
703 // we should get here, based on the if-condition at the top
704 Assert(false, ExcInternalError());
705 }
706 }
707 }
708 catch (
710 {
711 // simply fall through and continue on to the standard Newton code
712 }
713 }
714 else
715 {
716 // we can't use an explicit formula,
717 }
718
719
720 // Find the initial value for the Newton iteration by a normal
721 // projection to the least square plane determined by the vertices
722 // of the cell
723 Point<dim> initial_p_unit;
724 if (this->preserves_vertex_locations())
725 {
726 initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
727 // in 1d with spacedim > 1 the affine approximation is exact
728 if (dim == 1 && polynomial_degree == 1)
729 return initial_p_unit;
730 }
731 else
732 {
733 // else, we simply use the mid point
734 for (unsigned int d = 0; d < dim; ++d)
735 initial_p_unit[d] = 0.5;
736 }
737
738 // perform the Newton iteration and return the result. note that this
739 // statement may throw an exception, which we simply pass up to the caller
740 const Point<dim> p_unit =
741 this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
742 if (p_unit[0] == std::numeric_limits<double>::infinity())
743 AssertThrow(false,
745 return p_unit;
746}
747
748
749
750template <int dim, int spacedim>
751void
754 const ArrayView<const Point<spacedim>> & real_points,
755 const ArrayView<Point<dim>> & unit_points) const
756{
757 // Go to base class functions for dim < spacedim because it is not yet
758 // implemented with optimized code.
759 if (dim < spacedim)
760 {
762 real_points,
763 unit_points);
764 return;
765 }
766
767 AssertDimension(real_points.size(), unit_points.size());
768 const std::vector<Point<spacedim>> support_points =
769 this->compute_mapping_support_points(cell);
770
771 // From the given (high-order) support points, now only pick the first
772 // 2^dim points and construct an affine approximation from those.
774 inverse_approximation(support_points, unit_cell_support_points);
775
776 const unsigned int n_points = real_points.size();
777 const unsigned int n_lanes = VectorizedArray<double>::size();
778
779 // Use the more heavy VectorizedArray code path if there is more than
780 // one point left to compute
781 for (unsigned int i = 0; i < n_points; i += n_lanes)
782 if (n_points - i > 1)
783 {
785 for (unsigned int j = 0; j < n_lanes; ++j)
786 if (i + j < n_points)
787 for (unsigned int d = 0; d < spacedim; ++d)
788 p_vec[d][j] = real_points[i + j][d];
789 else
790 for (unsigned int d = 0; d < spacedim; ++d)
791 p_vec[d][j] = real_points[i][d];
792
794 internal::MappingQImplementation::
795 do_transform_real_to_unit_cell_internal<dim, spacedim>(
796 p_vec,
797 inverse_approximation.compute(p_vec),
798 support_points,
799 polynomials_1d,
800 renumber_lexicographic_to_hierarchic);
801
802 // If the vectorized computation failed, it could be that only some of
803 // the lanes failed but others would have succeeded if we had let them
804 // compute alone without interference (like negative Jacobian
805 // determinants) from other SIMD lanes. Repeat the computation in this
806 // unlikely case with scalar arguments.
807 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
808 if (unit_point[0][j] == std::numeric_limits<double>::infinity())
809 unit_points[i + j] = internal::MappingQImplementation::
810 do_transform_real_to_unit_cell_internal<dim, spacedim>(
811 real_points[i + j],
812 inverse_approximation.compute(real_points[i + j]),
813 support_points,
814 polynomials_1d,
815 renumber_lexicographic_to_hierarchic);
816 else
817 for (unsigned int d = 0; d < dim; ++d)
818 unit_points[i + j][d] = unit_point[d][j];
819 }
820 else
821 unit_points[i] = internal::MappingQImplementation::
822 do_transform_real_to_unit_cell_internal<dim, spacedim>(
823 real_points[i],
824 inverse_approximation.compute(real_points[i]),
825 support_points,
826 polynomials_1d,
827 renumber_lexicographic_to_hierarchic);
828}
829
830
831
832template <int dim, int spacedim>
835{
836 // add flags if the respective quantities are necessary to compute
837 // what we need. note that some flags appear in both the conditions
838 // and in subsequent set operations. this leads to some circular
839 // logic. the only way to treat this is to iterate. since there are
840 // 5 if-clauses in the loop, it will take at most 5 iterations to
841 // converge. do them:
842 UpdateFlags out = in;
843 for (unsigned int i = 0; i < 5; ++i)
844 {
845 // The following is a little incorrect:
846 // If not applied on a face,
847 // update_boundary_forms does not
848 // make sense. On the other hand,
849 // it is necessary on a
850 // face. Currently,
851 // update_boundary_forms is simply
852 // ignored for the interior of a
853 // cell.
856
861
862 if (out &
867
868 // The contravariant transformation is used in the Piola
869 // transformation, which requires the determinant of the Jacobi
870 // matrix of the transformation. Because we have no way of
871 // knowing here whether the finite element wants to use the
872 // contravariant or the Piola transforms, we add the JxW values
873 // to the list of flags to be updated for each cell.
876
877 // the same is true when computing normal vectors: they require
878 // the determinant of the Jacobian
879 if (out & update_normal_vectors)
881 }
882
883 return out;
884}
885
886
887
888template <int dim, int spacedim>
889std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
891 const Quadrature<dim> &q) const
892{
893 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
894 std::make_unique<InternalData>(polynomial_degree);
895 auto &data = dynamic_cast<InternalData &>(*data_ptr);
896 data.initialize(this->requires_update_flags(update_flags), q, q.size());
897
898 return data_ptr;
899}
900
901
902
903template <int dim, int spacedim>
904std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
906 const UpdateFlags update_flags,
907 const hp::QCollection<dim - 1> &quadrature) const
908{
909 AssertDimension(quadrature.size(), 1);
910
911 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
912 std::make_unique<InternalData>(polynomial_degree);
913 auto &data = dynamic_cast<InternalData &>(*data_ptr);
914 data.initialize_face(this->requires_update_flags(update_flags),
916 ReferenceCells::get_hypercube<dim>(), quadrature[0]),
917 quadrature[0].size());
918
919 return data_ptr;
920}
921
922
923
924template <int dim, int spacedim>
925std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
927 const UpdateFlags update_flags,
928 const Quadrature<dim - 1> &quadrature) const
929{
930 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
931 std::make_unique<InternalData>(polynomial_degree);
932 auto &data = dynamic_cast<InternalData &>(*data_ptr);
933 data.initialize_face(this->requires_update_flags(update_flags),
935 ReferenceCells::get_hypercube<dim>(), quadrature),
936 quadrature.size());
937
938 return data_ptr;
939}
940
941
942
943template <int dim, int spacedim>
947 const CellSimilarity::Similarity cell_similarity,
948 const Quadrature<dim> & quadrature,
949 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
951 &output_data) const
952{
953 // ensure that the following static_cast is really correct:
954 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
956 const InternalData &data = static_cast<const InternalData &>(internal_data);
957
958 const unsigned int n_q_points = quadrature.size();
959
960 // recompute the support points of the transformation of this
961 // cell. we tried to be clever here in an earlier version of the
962 // library by checking whether the cell is the same as the one we
963 // had visited last, but it turns out to be difficult to determine
964 // that because a cell for the purposes of a mapping is
965 // characterized not just by its (triangulation, level, index)
966 // triple, but also by the locations of its vertices, the manifold
967 // object attached to the cell and all of its bounding faces/edges,
968 // etc. to reliably test that the "cell" we are on is, therefore,
969 // not easily done
970 data.mapping_support_points = this->compute_mapping_support_points(cell);
972
973 // if the order of the mapping is greater than 1, then do not reuse any cell
974 // similarity information. This is necessary because the cell similarity
975 // value is computed with just cell vertices and does not take into account
976 // cell curvature.
977 const CellSimilarity::Similarity computed_cell_similarity =
978 (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
979
980 if (dim > 1 && data.tensor_product_quadrature)
981 {
982 internal::MappingQImplementation::
983 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
984 computed_cell_similarity,
985 data,
986 output_data.quadrature_points,
987 output_data.jacobian_grads);
988 }
989 else
990 {
991 internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
993 data,
994 output_data.quadrature_points);
995
996 internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
997 computed_cell_similarity,
999 data);
1000
1002 spacedim>(
1003 computed_cell_similarity,
1005 data,
1006 output_data.jacobian_grads);
1007 }
1008
1010 dim,
1011 spacedim>(computed_cell_similarity,
1013 data,
1014 output_data.jacobian_pushed_forward_grads);
1015
1017 dim,
1018 spacedim>(computed_cell_similarity,
1020 data,
1021 output_data.jacobian_2nd_derivatives);
1022
1023 internal::MappingQImplementation::
1024 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1025 computed_cell_similarity,
1027 data,
1029
1031 dim,
1032 spacedim>(computed_cell_similarity,
1034 data,
1035 output_data.jacobian_3rd_derivatives);
1036
1037 internal::MappingQImplementation::
1038 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1039 computed_cell_similarity,
1041 data,
1043
1044 const UpdateFlags update_flags = data.update_each;
1045 const std::vector<double> &weights = quadrature.get_weights();
1046
1047 // Multiply quadrature weights by absolute value of Jacobian determinants or
1048 // the area element g=sqrt(DX^t DX) in case of codim > 0
1049
1050 if (update_flags & (update_normal_vectors | update_JxW_values))
1051 {
1052 AssertDimension(output_data.JxW_values.size(), n_q_points);
1053
1054 Assert(!(update_flags & update_normal_vectors) ||
1055 (output_data.normal_vectors.size() == n_q_points),
1056 ExcDimensionMismatch(output_data.normal_vectors.size(),
1057 n_q_points));
1058
1059
1060 if (computed_cell_similarity != CellSimilarity::translation)
1061 for (unsigned int point = 0; point < n_q_points; ++point)
1062 {
1063 if (dim == spacedim)
1064 {
1065 const double det = data.contravariant[point].determinant();
1066
1067 // check for distorted cells.
1068
1069 // TODO: this allows for anisotropies of up to 1e6 in 3D and
1070 // 1e12 in 2D. might want to find a finer
1071 // (dimension-independent) criterion
1072 Assert(det >
1073 1e-12 * Utilities::fixed_power<dim>(
1074 cell->diameter() / std::sqrt(double(dim))),
1076 cell->center(), det, point)));
1077
1078 output_data.JxW_values[point] = weights[point] * det;
1079 }
1080 // if dim==spacedim, then there is no cell normal to
1081 // compute. since this is for FEValues (and not FEFaceValues),
1082 // there are also no face normals to compute
1083 else // codim>0 case
1084 {
1085 Tensor<1, spacedim> DX_t[dim];
1086 for (unsigned int i = 0; i < spacedim; ++i)
1087 for (unsigned int j = 0; j < dim; ++j)
1088 DX_t[j][i] = data.contravariant[point][i][j];
1089
1090 Tensor<2, dim> G; // First fundamental form
1091 for (unsigned int i = 0; i < dim; ++i)
1092 for (unsigned int j = 0; j < dim; ++j)
1093 G[i][j] = DX_t[i] * DX_t[j];
1094
1095 output_data.JxW_values[point] =
1096 std::sqrt(determinant(G)) * weights[point];
1097
1098 if (computed_cell_similarity ==
1100 {
1101 // we only need to flip the normal
1102 if (update_flags & update_normal_vectors)
1103 output_data.normal_vectors[point] *= -1.;
1104 }
1105 else
1106 {
1107 if (update_flags & update_normal_vectors)
1108 {
1109 Assert(spacedim == dim + 1,
1110 ExcMessage(
1111 "There is no (unique) cell normal for " +
1113 "-dimensional cells in " +
1114 Utilities::int_to_string(spacedim) +
1115 "-dimensional space. This only works if the "
1116 "space dimension is one greater than the "
1117 "dimensionality of the mesh cells."));
1118
1119 if (dim == 1)
1120 output_data.normal_vectors[point] =
1121 cross_product_2d(-DX_t[0]);
1122 else // dim == 2
1123 output_data.normal_vectors[point] =
1124 cross_product_3d(DX_t[0], DX_t[1]);
1125
1126 output_data.normal_vectors[point] /=
1127 output_data.normal_vectors[point].norm();
1128
1129 if (cell->direction_flag() == false)
1130 output_data.normal_vectors[point] *= -1.;
1131 }
1132 }
1133 } // codim>0 case
1134 }
1135 }
1136
1137
1138
1139 // copy values from InternalData to vector given by reference
1140 if (update_flags & update_jacobians)
1141 {
1142 AssertDimension(output_data.jacobians.size(), n_q_points);
1143 if (computed_cell_similarity != CellSimilarity::translation)
1144 for (unsigned int point = 0; point < n_q_points; ++point)
1145 output_data.jacobians[point] = data.contravariant[point];
1146 }
1147
1148 // copy values from InternalData to vector given by reference
1149 if (update_flags & update_inverse_jacobians)
1150 {
1151 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1152 if (computed_cell_similarity != CellSimilarity::translation)
1153 for (unsigned int point = 0; point < n_q_points; ++point)
1154 output_data.inverse_jacobians[point] =
1155 data.covariant[point].transpose();
1156 }
1157
1158 return computed_cell_similarity;
1159}
1160
1161
1162
1163template <int dim, int spacedim>
1164void
1167 const unsigned int face_no,
1168 const hp::QCollection<dim - 1> & quadrature,
1169 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1171 &output_data) const
1172{
1173 AssertDimension(quadrature.size(), 1);
1174
1175 // ensure that the following cast is really correct:
1176 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1178 const InternalData &data = static_cast<const InternalData &>(internal_data);
1179
1180 // if necessary, recompute the support points of the transformation of this
1181 // cell (note that we need to first check the triangulation pointer, since
1182 // otherwise the second test might trigger an exception if the triangulations
1183 // are not the same)
1184 if ((data.mapping_support_points.size() == 0) ||
1185 (&cell->get_triangulation() !=
1187 (cell != data.cell_of_current_support_points))
1188 {
1189 data.mapping_support_points = this->compute_mapping_support_points(cell);
1191 }
1192
1194 *this,
1195 cell,
1196 face_no,
1199 ReferenceCells::get_hypercube<dim>(),
1200 face_no,
1201 cell->face_orientation(face_no),
1202 cell->face_flip(face_no),
1203 cell->face_rotation(face_no),
1204 quadrature[0].size()),
1205 quadrature[0],
1206 data,
1207 output_data);
1208}
1209
1210
1211
1212template <int dim, int spacedim>
1213void
1216 const unsigned int face_no,
1217 const unsigned int subface_no,
1218 const Quadrature<dim - 1> & quadrature,
1219 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1221 &output_data) const
1222{
1223 // ensure that the following cast is really correct:
1224 Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1226 const InternalData &data = static_cast<const InternalData &>(internal_data);
1227
1228 // if necessary, recompute the support points of the transformation of this
1229 // cell (note that we need to first check the triangulation pointer, since
1230 // otherwise the second test might trigger an exception if the triangulations
1231 // are not the same)
1232 if ((data.mapping_support_points.size() == 0) ||
1233 (&cell->get_triangulation() !=
1235 (cell != data.cell_of_current_support_points))
1236 {
1237 data.mapping_support_points = this->compute_mapping_support_points(cell);
1239 }
1240
1242 *this,
1243 cell,
1244 face_no,
1245 subface_no,
1247 ReferenceCells::get_hypercube<dim>(),
1248 face_no,
1249 subface_no,
1250 cell->face_orientation(face_no),
1251 cell->face_flip(face_no),
1252 cell->face_rotation(face_no),
1253 quadrature.size(),
1254 cell->subface_case(face_no)),
1255 quadrature,
1256 data,
1257 output_data);
1258}
1259
1260
1261
1262template <int dim, int spacedim>
1263void
1267 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1269 &output_data) const
1270{
1271 Assert(dim == spacedim, ExcNotImplemented());
1272
1273 // ensure that the following static_cast is really correct:
1274 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1276 const InternalData &data = static_cast<const InternalData &>(internal_data);
1277
1278 const unsigned int n_q_points = quadrature.size();
1279
1280 data.mapping_support_points = this->compute_mapping_support_points(cell);
1282
1283 internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
1285 data,
1286 output_data.quadrature_points);
1287
1288 internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
1290
1291 internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
1294 data,
1295 output_data.jacobian_grads);
1296
1298 dim,
1299 spacedim>(CellSimilarity::none,
1301 data,
1302 output_data.jacobian_pushed_forward_grads);
1303
1305 dim,
1306 spacedim>(CellSimilarity::none,
1308 data,
1309 output_data.jacobian_2nd_derivatives);
1310
1311 internal::MappingQImplementation::
1312 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
1315 data,
1317
1319 dim,
1320 spacedim>(CellSimilarity::none,
1322 data,
1323 output_data.jacobian_3rd_derivatives);
1324
1325 internal::MappingQImplementation::
1326 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
1329 data,
1331
1332 const UpdateFlags update_flags = data.update_each;
1333 const std::vector<double> &weights = quadrature.get_weights();
1334
1335 if ((update_flags & (update_normal_vectors | update_JxW_values)) != 0u)
1336 {
1337 AssertDimension(output_data.JxW_values.size(), n_q_points);
1338
1339 Assert(!(update_flags & update_normal_vectors) ||
1340 (output_data.normal_vectors.size() == n_q_points),
1341 ExcDimensionMismatch(output_data.normal_vectors.size(),
1342 n_q_points));
1343
1344
1345 for (unsigned int point = 0; point < n_q_points; ++point)
1346 {
1347 const double det = data.contravariant[point].determinant();
1348
1349 // check for distorted cells.
1350
1351 // TODO: this allows for anisotropies of up to 1e6 in 3D and
1352 // 1e12 in 2D. might want to find a finer
1353 // (dimension-independent) criterion
1354 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1355 cell->diameter() / std::sqrt(double(dim))),
1357 cell->center(), det, point)));
1358
1359 // The normals are n = J^{-T} * \hat{n} before normalizing.
1360 Tensor<1, spacedim> normal;
1361 for (unsigned int d = 0; d < spacedim; d++)
1362 normal[d] =
1363 data.covariant[point][d] * quadrature.normal_vector(point);
1364
1365 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1366
1367 if ((update_flags & update_normal_vectors) != 0u)
1368 {
1369 normal /= normal.norm();
1370 output_data.normal_vectors[point] = normal;
1371 }
1372 }
1373 }
1374
1375 // copy values from InternalData to vector given by reference
1376 if ((update_flags & update_jacobians) != 0u)
1377 {
1378 AssertDimension(output_data.jacobians.size(), n_q_points);
1379 for (unsigned int point = 0; point < n_q_points; ++point)
1380 output_data.jacobians[point] = data.contravariant[point];
1381 }
1382
1383 // copy values from InternalData to vector given by reference
1384 if ((update_flags & update_inverse_jacobians) != 0u)
1385 {
1386 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1387 for (unsigned int point = 0; point < n_q_points; ++point)
1388 output_data.inverse_jacobians[point] =
1389 data.covariant[point].transpose();
1390 }
1391}
1392
1393
1394
1395template <int dim, int spacedim>
1396void
1399 const ArrayView<const Point<dim>> & unit_points,
1400 const UpdateFlags update_flags,
1402 &output_data) const
1403{
1404 if (update_flags == update_default)
1405 return;
1406
1407 Assert(update_flags & update_inverse_jacobians ||
1408 update_flags & update_jacobians ||
1409 update_flags & update_quadrature_points,
1411
1412 output_data.initialize(unit_points.size(), update_flags);
1413 const std::vector<Point<spacedim>> support_points =
1414 this->compute_mapping_support_points(cell);
1415
1416 const unsigned int n_points = unit_points.size();
1417 const unsigned int n_lanes = VectorizedArray<double>::size();
1418
1419 // Use the more heavy VectorizedArray code path if there is more than
1420 // one point left to compute
1421 for (unsigned int i = 0; i < n_points; i += n_lanes)
1422 if (n_points - i > 1)
1423 {
1425 for (unsigned int j = 0; j < n_lanes; ++j)
1426 if (i + j < n_points)
1427 for (unsigned int d = 0; d < dim; ++d)
1428 p_vec[d][j] = unit_points[i + j][d];
1429 else
1430 for (unsigned int d = 0; d < dim; ++d)
1431 p_vec[d][j] = unit_points[i][d];
1432
1433 const auto result =
1435 polynomials_1d,
1436 support_points,
1437 p_vec,
1438 polynomial_degree == 1,
1439 renumber_lexicographic_to_hierarchic);
1440
1441 if (update_flags & update_quadrature_points)
1442 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1443 for (unsigned int d = 0; d < spacedim; ++d)
1444 output_data.quadrature_points[i + j][d] = result.first[d][j];
1445
1446 if (update_flags & update_jacobians)
1447 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1448 for (unsigned int d = 0; d < spacedim; ++d)
1449 for (unsigned int e = 0; e < dim; ++e)
1450 output_data.jacobians[i + j][d][e] = result.second[e][d][j];
1451
1452 if (update_flags & update_inverse_jacobians)
1453 {
1455 result.second);
1457 inv_jac = jac.covariant_form();
1458 for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
1459 for (unsigned int d = 0; d < dim; ++d)
1460 for (unsigned int e = 0; e < spacedim; ++e)
1461 output_data.inverse_jacobians[i + j][d][e] = inv_jac[d][e][j];
1462 }
1463 }
1464 else
1465 {
1466 const auto result =
1468 polynomials_1d,
1469 support_points,
1470 unit_points[i],
1471 polynomial_degree == 1,
1472 renumber_lexicographic_to_hierarchic);
1473
1474 if (update_flags & update_quadrature_points)
1475 output_data.quadrature_points[i] = result.first;
1476
1477 if (update_flags & update_jacobians)
1478 {
1479 DerivativeForm<1, spacedim, dim> jac = result.second;
1480 output_data.jacobians[i] = jac.transpose();
1481 }
1482
1483 if (update_flags & update_inverse_jacobians)
1484 {
1485 DerivativeForm<1, spacedim, dim> jac(result.second);
1487 for (unsigned int d = 0; d < dim; ++d)
1488 for (unsigned int e = 0; e < spacedim; ++e)
1489 output_data.inverse_jacobians[i][d][e] = inv_jac[d][e];
1490 }
1491 }
1492}
1493
1494
1495
1496template <int dim, int spacedim>
1497void
1499 const ArrayView<const Tensor<1, dim>> & input,
1500 const MappingKind mapping_kind,
1501 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1502 const ArrayView<Tensor<1, spacedim>> & output) const
1503{
1505 mapping_kind,
1506 mapping_data,
1507 output);
1508}
1509
1510
1511
1512template <int dim, int spacedim>
1513void
1515 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1516 const MappingKind mapping_kind,
1517 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1518 const ArrayView<Tensor<2, spacedim>> & output) const
1519{
1521 mapping_kind,
1522 mapping_data,
1523 output);
1524}
1525
1526
1527
1528template <int dim, int spacedim>
1529void
1531 const ArrayView<const Tensor<2, dim>> & input,
1532 const MappingKind mapping_kind,
1533 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1534 const ArrayView<Tensor<2, spacedim>> & output) const
1535{
1536 switch (mapping_kind)
1537 {
1540 mapping_kind,
1541 mapping_data,
1542 output);
1543 return;
1544
1549 mapping_kind,
1550 mapping_data,
1551 output);
1552 return;
1553 default:
1554 Assert(false, ExcNotImplemented());
1555 }
1556}
1557
1558
1559
1560template <int dim, int spacedim>
1561void
1563 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1564 const MappingKind mapping_kind,
1565 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1566 const ArrayView<Tensor<3, spacedim>> & output) const
1567{
1568 AssertDimension(input.size(), output.size());
1569 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1571 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1572
1573 switch (mapping_kind)
1574 {
1576 {
1579 "update_covariant_transformation"));
1580
1581 for (unsigned int q = 0; q < output.size(); ++q)
1582 for (unsigned int i = 0; i < spacedim; ++i)
1583 for (unsigned int j = 0; j < spacedim; ++j)
1584 {
1585 double tmp[dim];
1586 for (unsigned int K = 0; K < dim; ++K)
1587 {
1588 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1589 for (unsigned int J = 1; J < dim; ++J)
1590 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1591 }
1592 for (unsigned int k = 0; k < spacedim; ++k)
1593 {
1594 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1595 for (unsigned int K = 1; K < dim; ++K)
1596 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1597 }
1598 }
1599 return;
1600 }
1601
1602 default:
1603 Assert(false, ExcNotImplemented());
1604 }
1605}
1606
1607
1608
1609template <int dim, int spacedim>
1610void
1612 const ArrayView<const Tensor<3, dim>> & input,
1613 const MappingKind mapping_kind,
1614 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1615 const ArrayView<Tensor<3, spacedim>> & output) const
1616{
1617 switch (mapping_kind)
1618 {
1623 mapping_kind,
1624 mapping_data,
1625 output);
1626 return;
1627 default:
1628 Assert(false, ExcNotImplemented());
1629 }
1630}
1631
1632
1633
1634template <int dim, int spacedim>
1635void
1638 std::vector<Point<spacedim>> & a) const
1639{
1640 // if we only need the midpoint, then ask for it.
1641 if (this->polynomial_degree == 2)
1642 {
1643 for (unsigned int line_no = 0;
1644 line_no < GeometryInfo<dim>::lines_per_cell;
1645 ++line_no)
1646 {
1648 (dim == 1 ?
1649 static_cast<
1651 cell->line(line_no));
1652
1653 const Manifold<dim, spacedim> &manifold =
1654 ((line->manifold_id() == numbers::flat_manifold_id) &&
1655 (dim < spacedim) ?
1656 cell->get_manifold() :
1657 line->get_manifold());
1658 a.push_back(manifold.get_new_point_on_line(line));
1659 }
1660 }
1661 else
1662 // otherwise call the more complicated functions and ask for inner points
1663 // from the manifold description
1664 {
1665 std::vector<Point<spacedim>> tmp_points;
1666 for (unsigned int line_no = 0;
1667 line_no < GeometryInfo<dim>::lines_per_cell;
1668 ++line_no)
1669 {
1671 (dim == 1 ?
1672 static_cast<
1674 cell->line(line_no));
1675
1676 const Manifold<dim, spacedim> &manifold =
1677 ((line->manifold_id() == numbers::flat_manifold_id) &&
1678 (dim < spacedim) ?
1679 cell->get_manifold() :
1680 line->get_manifold());
1681
1682 const std::array<Point<spacedim>, 2> vertices{
1683 {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
1684 cell->vertex(
1686
1687 const std::size_t n_rows =
1688 support_point_weights_perimeter_to_interior[0].size(0);
1689 a.resize(a.size() + n_rows);
1690 auto a_view = make_array_view(a.end() - n_rows, a.end());
1691 manifold.get_new_points(
1692 make_array_view(vertices.begin(), vertices.end()),
1693 support_point_weights_perimeter_to_interior[0],
1694 a_view);
1695 }
1696 }
1697}
1698
1699
1700
1701template <>
1702void
1705 std::vector<Point<3>> & a) const
1706{
1707 const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
1708
1709 // used if face quad at boundary or entirely in the interior of the domain
1710 std::vector<Point<3>> tmp_points;
1711
1712 // loop over all faces and collect points on them
1713 for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
1714 {
1715 const Triangulation<3>::face_iterator face = cell->face(face_no);
1716
1717#ifdef DEBUG
1718 const bool face_orientation = cell->face_orientation(face_no),
1719 face_flip = cell->face_flip(face_no),
1720 face_rotation = cell->face_rotation(face_no);
1721 const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
1722 lines_per_face = GeometryInfo<3>::lines_per_face;
1723
1724 // some sanity checks up front
1725 for (unsigned int i = 0; i < vertices_per_face; ++i)
1726 Assert(face->vertex_index(i) ==
1727 cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
1728 face_no, i, face_orientation, face_flip, face_rotation)),
1730
1731 // indices of the lines that bound a face are given by GeometryInfo<3>::
1732 // face_to_cell_lines
1733 for (unsigned int i = 0; i < lines_per_face; ++i)
1734 Assert(face->line(i) ==
1736 face_no, i, face_orientation, face_flip, face_rotation)),
1738#endif
1739 // extract the points surrounding a quad from the points
1740 // already computed. First get the 4 vertices and then the points on
1741 // the four lines
1742 boost::container::small_vector<Point<3>, 200> tmp_points(
1744 GeometryInfo<2>::lines_per_cell * (polynomial_degree - 1));
1745 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1746 tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
1747 if (polynomial_degree > 1)
1748 for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1749 ++line)
1750 for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
1751 tmp_points[4 + line * (polynomial_degree - 1) + i] =
1753 (polynomial_degree - 1) *
1755 i];
1756
1757 const std::size_t n_rows =
1758 support_point_weights_perimeter_to_interior[1].size(0);
1759 a.resize(a.size() + n_rows);
1760 auto a_view = make_array_view(a.end() - n_rows, a.end());
1761 face->get_manifold().get_new_points(
1762 make_array_view(tmp_points.begin(), tmp_points.end()),
1763 support_point_weights_perimeter_to_interior[1],
1764 a_view);
1765 }
1766}
1767
1768
1769
1770template <>
1771void
1774 std::vector<Point<3>> & a) const
1775{
1776 std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
1777 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1778 vertices[i] = cell->vertex(i);
1779
1780 Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 1),
1782 for (unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
1783 for (unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
1784 {
1785 Point<2> point(line_support_points[q1 + 1][0],
1786 line_support_points[q2 + 1][0]);
1787 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
1789 }
1790
1791 const std::size_t n_rows = weights.size(0);
1792 a.resize(a.size() + n_rows);
1793 auto a_view = make_array_view(a.end() - n_rows, a.end());
1794 cell->get_manifold().get_new_points(
1795 make_array_view(vertices.begin(), vertices.end()), weights, a_view);
1796}
1797
1798
1799
1800template <int dim, int spacedim>
1801void
1804 std::vector<Point<spacedim>> &) const
1805{
1806 Assert(false, ExcInternalError());
1807}
1808
1809
1810
1811template <int dim, int spacedim>
1812std::vector<Point<spacedim>>
1814 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1815{
1816 // get the vertices first
1817 std::vector<Point<spacedim>> a;
1818 a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
1819 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1820 a.push_back(cell->vertex(i));
1821
1822 if (this->polynomial_degree > 1)
1823 {
1824 // check if all entities have the same manifold id which is when we can
1825 // simply ask the manifold for all points. the transfinite manifold can
1826 // do the interpolation better than this class, so if we detect that we
1827 // do not have to change anything here
1828 Assert(dim <= 3, ExcImpossibleInDim(dim));
1829 bool all_manifold_ids_are_equal = (dim == spacedim);
1830 if (all_manifold_ids_are_equal &&
1832 &cell->get_manifold()) == nullptr)
1833 {
1834 for (auto f : GeometryInfo<dim>::face_indices())
1835 if (&cell->face(f)->get_manifold() != &cell->get_manifold())
1836 all_manifold_ids_are_equal = false;
1837
1838 if (dim == 3)
1839 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1840 if (&cell->line(l)->get_manifold() != &cell->get_manifold())
1841 all_manifold_ids_are_equal = false;
1842 }
1843
1844 if (all_manifold_ids_are_equal)
1845 {
1846 const std::size_t n_rows = support_point_weights_cell.size(0);
1847 a.resize(a.size() + n_rows);
1848 auto a_view = make_array_view(a.end() - n_rows, a.end());
1850 a.end() - n_rows),
1851 support_point_weights_cell,
1852 a_view);
1853 }
1854 else
1855 switch (dim)
1856 {
1857 case 1:
1858 add_line_support_points(cell, a);
1859 break;
1860 case 2:
1861 // in 2d, add the points on the four bounding lines to the
1862 // exterior (outer) points
1863 add_line_support_points(cell, a);
1864
1865 // then get the interior support points
1866 if (dim != spacedim)
1867 add_quad_support_points(cell, a);
1868 else
1869 {
1870 const std::size_t n_rows =
1871 support_point_weights_perimeter_to_interior[1].size(0);
1872 a.resize(a.size() + n_rows);
1873 auto a_view = make_array_view(a.end() - n_rows, a.end());
1875 make_array_view(a.begin(), a.end() - n_rows),
1876 support_point_weights_perimeter_to_interior[1],
1877 a_view);
1878 }
1879 break;
1880
1881 case 3:
1882 // in 3d also add the points located on the boundary faces
1883 add_line_support_points(cell, a);
1884 add_quad_support_points(cell, a);
1885
1886 // then compute the interior points
1887 {
1888 const std::size_t n_rows =
1889 support_point_weights_perimeter_to_interior[2].size(0);
1890 a.resize(a.size() + n_rows);
1891 auto a_view = make_array_view(a.end() - n_rows, a.end());
1893 make_array_view(a.begin(), a.end() - n_rows),
1894 support_point_weights_perimeter_to_interior[2],
1895 a_view);
1896 }
1897 break;
1898
1899 default:
1900 Assert(false, ExcNotImplemented());
1901 break;
1902 }
1903 }
1904
1905 return a;
1906}
1907
1908
1909
1910template <int dim, int spacedim>
1913 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
1914{
1915 return BoundingBox<spacedim>(this->compute_mapping_support_points(cell));
1916}
1917
1918
1919
1920template <int dim, int spacedim>
1921bool
1923 const ReferenceCell &reference_cell) const
1924{
1925 Assert(dim == reference_cell.get_dimension(),
1926 ExcMessage("The dimension of your mapping (" +
1928 ") and the reference cell cell_type (" +
1929 Utilities::to_string(reference_cell.get_dimension()) +
1930 " ) do not agree."));
1931
1932 return reference_cell.is_hyper_cube();
1933}
1934
1935
1936
1937//--------------------------- Explicit instantiations -----------------------
1938#include "mapping_q.inst"
1939
1940
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:699
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
DerivativeForm< 1, spacedim, dim, Number > transpose() const
Definition: fe_dgq.h:111
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
AlignedVector< DerivativeForm< 1, dim, spacedim > > covariant
Definition: mapping_q.h:521
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
Definition: mapping_q.h:547
virtual std::size_t memory_consumption() const override
Definition: mapping_q.cc:64
std::vector< Point< spacedim > > mapping_support_points
Definition: mapping_q.h:541
AlignedVector< DerivativeForm< 1, dim, spacedim > > contravariant
Definition: mapping_q.h:530
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:215
InternalData(const unsigned int polynomial_degree)
Definition: mapping_q.cc:52
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
Definition: mapping_q.cc:85
void compute_shape_function_values(const std::vector< Point< dim > > &unit_points)
Definition: mapping_q.cc:271
const std::vector< unsigned int > renumber_lexicographic_to_hierarchic
Definition: mapping_q.h:651
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:945
const Table< 2, double > support_point_weights_cell
Definition: mapping_q.h:699
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping_q.cc:1813
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_q.cc:890
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
Definition: mapping_q.cc:1214
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition: mapping_q.cc:1912
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:926
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
Definition: mapping_q.cc:1498
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
Definition: mapping_q.cc:1165
const unsigned int polynomial_degree
Definition: mapping_q.h:627
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:752
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:905
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:488
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
Definition: mapping_q.cc:1397
const std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
Definition: mapping_q.h:685
const std::vector< Point< 1 > > line_support_points
Definition: mapping_q.h:637
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:453
const std::vector< Polynomials::Polynomial< double > > polynomials_1d
Definition: mapping_q.h:644
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_q.cc:435
const std::vector< Point< dim > > unit_cell_support_points
Definition: mapping_q.h:663
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:631
MappingQ(const unsigned int polynomial_degree)
Definition: mapping_q.cc:362
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim > > &a) const
Definition: mapping_q.cc:1636
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim > > &a) const
Definition: mapping_q.cc:1802
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) 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
Definition: mapping_q.cc:1264
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
Definition: mapping_q.cc:1922
unsigned int get_degree() const
Definition: mapping_q.cc:444
Abstract base class for mapping classes.
Definition: mapping.h:311
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
static DataSetDescriptor cell()
Definition: qprojector.h:563
bool is_tensor_product() const
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:325
const std::vector< double > & get_weights() const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2635
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660
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)
Definition: fe_values.cc:2704
std::vector< Tensor< 1, spacedim > > normal_vectors
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< Point< spacedim > > quadrature_points
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:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
Point< 3 > vertices[4]
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_3rd_derivatives
@ 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.
@ update_jacobian_2nd_derivatives
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
MappingKind
Definition: mapping.h:72
@ mapping_covariant_gradient
Definition: mapping.h:93
@ mapping_contravariant
Definition: mapping.h:87
@ mapping_contravariant_hessian
Definition: mapping.h:149
@ mapping_covariant_hessian
Definition: mapping.h:143
@ mapping_contravariant_gradient
Definition: mapping.h:99
@ mapping_piola_gradient
Definition: mapping.h:113
@ mapping_piola_hessian
Definition: mapping.h:155
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition: polynomial.cc:702
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
T fixed_power(const T t)
Definition: utilities.h:1123
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
inline ::Table< 2, double > compute_support_point_weights_cell(const unsigned int polynomial_degree)
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_3rd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 4, dim, spacedim > > &jacobian_3rd_derivatives)
void maybe_update_jacobian_pushed_forward_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< Tensor< 3, spacedim > > &jacobian_pushed_forward_grads)
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)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)
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 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, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data)
void maybe_update_jacobian_2nd_derivatives(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 3, dim, spacedim > > &jacobian_2nd_derivatives)
void maybe_update_jacobian_grads(const CellSimilarity::Similarity cell_similarity, const typename QProjector< dim >::DataSetDescriptor data_set, const typename ::MappingQ< dim, spacedim >::InternalData &data, std::vector< DerivativeForm< 2, dim, spacedim > > &jacobian_grads)
std::vector<::Table< 2, double > > compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, const unsigned int dim)
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(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:201
const types::manifold_id flat_manifold_id
Definition: types.h:269
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)