66 template <
int spacedim>
68 transform_real_to_unit_cell(
80 template <
int spacedim>
82 transform_real_to_unit_cell(
92 const long double x = p(0);
93 const long double y = p(1);
95 const long double x0 =
vertices[0](0);
96 const long double x1 =
vertices[1](0);
97 const long double x2 =
vertices[2](0);
98 const long double x3 =
vertices[3](0);
100 const long double y0 =
vertices[0](1);
101 const long double y1 =
vertices[1](1);
102 const long double y2 =
vertices[2](1);
103 const long double y3 =
vertices[3](1);
105 const long double a = (x1 - x3) * (y0 - y2) - (x0 - x2) * (y1 - y3);
106 const long double b = -(x0 - x1 - x2 + x3) * y +
107 (x - 2 * x1 + x3) * y0 - (x - 2 * x0 + x2) * y1 -
108 (x - x1) * y2 + (x - x0) * y3;
109 const long double c = (x0 - x1) * y - (x - x1) * y0 + (x - x0) * y1;
111 const long double discriminant =
b *
b - 4 * a * c;
120 const long double sqrt_discriminant = std::sqrt(discriminant);
123 if (
b != 0.0 && std::abs(
b) == sqrt_discriminant)
130 else if (std::abs(a) < 1
e-8 * std::abs(
b))
134 eta1 = 2 * c / (-
b - sqrt_discriminant);
135 eta2 = 2 * c / (-
b + sqrt_discriminant);
140 eta1 = (-
b - sqrt_discriminant) / (2 * a);
141 eta2 = (-
b + sqrt_discriminant) / (2 * a);
144 const long double eta =
145 (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
151 const long double subexpr0 = -eta * x2 + x0 * (eta - 1);
152 const long double xi_denominator0 =
153 eta * x3 - x1 * (eta - 1) + subexpr0;
154 const long double max_x =
156 std::max(std::abs(x2), std::abs(x3)));
158 if (std::abs(xi_denominator0) > 1
e-10 * max_x)
160 const double xi = (x + subexpr0) / xi_denominator0;
161 return {xi,
static_cast<double>(eta)};
165 const long double max_y =
167 std::max(std::abs(y2), std::abs(y3)));
168 const long double subexpr1 = -eta * y2 + y0 * (eta - 1);
169 const long double xi_denominator1 =
170 eta * y3 - y1 * (eta - 1) + subexpr1;
171 if (std::abs(xi_denominator1) > 1
e-10 * max_y)
173 const double xi = (subexpr1 + y) / xi_denominator1;
174 return {xi,
static_cast<double>(eta)};
181 spacedim>::ExcTransformationFailed()));
187 return {std::numeric_limits<double>::quiet_NaN(),
188 std::numeric_limits<double>::quiet_NaN()};
193 template <
int spacedim>
195 transform_real_to_unit_cell(
207 template <
int dim,
int spacedim>
209 compute_shape_function_values_general(
210 const unsigned int n_shape_functions,
212 typename ::MappingQGeneric<dim, spacedim>::InternalData &data)
214 const unsigned int n_points = unit_points.size();
220 data.line_support_points.get_points()));
225 const std::vector<unsigned int> renumber =
226 FETools::hierarchic_to_lexicographic_numbering<dim>(
227 data.polynomial_degree);
229 std::vector<double> values;
230 std::vector<Tensor<1, dim>> grads;
231 if (data.shape_values.size() != 0)
233 Assert(data.shape_values.size() == n_shape_functions * n_points,
235 values.resize(n_shape_functions);
237 if (data.shape_derivatives.size() != 0)
239 Assert(data.shape_derivatives.size() ==
240 n_shape_functions * n_points,
242 grads.resize(n_shape_functions);
245 std::vector<Tensor<2, dim>> grad2;
246 if (data.shape_second_derivatives.size() != 0)
248 Assert(data.shape_second_derivatives.size() ==
249 n_shape_functions * n_points,
251 grad2.resize(n_shape_functions);
254 std::vector<Tensor<3, dim>> grad3;
255 if (data.shape_third_derivatives.size() != 0)
257 Assert(data.shape_third_derivatives.size() ==
258 n_shape_functions * n_points,
260 grad3.resize(n_shape_functions);
263 std::vector<Tensor<4, dim>> grad4;
264 if (data.shape_fourth_derivatives.size() != 0)
266 Assert(data.shape_fourth_derivatives.size() ==
267 n_shape_functions * n_points,
269 grad4.resize(n_shape_functions);
273 if (data.shape_values.size() != 0 ||
274 data.shape_derivatives.size() != 0 ||
275 data.shape_second_derivatives.size() != 0 ||
276 data.shape_third_derivatives.size() != 0 ||
277 data.shape_fourth_derivatives.size() != 0)
281 unit_points[
point], values, grads, grad2, grad3, grad4);
283 if (data.shape_values.size() != 0)
284 for (
unsigned int i = 0; i < n_shape_functions; ++i)
285 data.shape(
point, i) = values[renumber[i]];
287 if (data.shape_derivatives.size() != 0)
288 for (
unsigned int i = 0; i < n_shape_functions; ++i)
289 data.derivative(
point, i) = grads[renumber[i]];
291 if (data.shape_second_derivatives.size() != 0)
292 for (
unsigned int i = 0; i < n_shape_functions; ++i)
293 data.second_derivative(
point, i) = grad2[renumber[i]];
295 if (data.shape_third_derivatives.size() != 0)
296 for (
unsigned int i = 0; i < n_shape_functions; ++i)
297 data.third_derivative(
point, i) = grad3[renumber[i]];
299 if (data.shape_fourth_derivatives.size() != 0)
300 for (
unsigned int i = 0; i < n_shape_functions; ++i)
301 data.fourth_derivative(
point, i) = grad4[renumber[i]];
307 compute_shape_function_values_hardcode(
308 const unsigned int n_shape_functions,
309 const std::vector<
Point<1>> & unit_points,
312 (void)n_shape_functions;
313 const unsigned int n_points = unit_points.size();
314 for (
unsigned int k = 0; k < n_points; ++k)
316 double x = unit_points[k](0);
322 data.
shape(k, 0) = 1. - x;
323 data.
shape(k, 1) = x;
328 n_shape_functions * n_points,
336 n_shape_functions * n_points,
344 n_shape_functions * n_points,
354 n_shape_functions * n_points,
366 compute_shape_function_values_hardcode(
367 const unsigned int n_shape_functions,
368 const std::vector<
Point<2>> & unit_points,
371 (void)n_shape_functions;
372 const unsigned int n_points = unit_points.size();
373 for (
unsigned int k = 0; k < n_points; ++k)
375 double x = unit_points[k](0);
376 double y = unit_points[k](1);
382 data.
shape(k, 0) = (1. - x) * (1. - y);
383 data.
shape(k, 1) = x * (1. - y);
384 data.
shape(k, 2) = (1. - x) * y;
385 data.
shape(k, 3) = x * y;
390 n_shape_functions * n_points,
404 n_shape_functions * n_points,
426 n_shape_functions * n_points,
430 for (
unsigned int i = 0; i < 4; ++i)
436 n_shape_functions * n_points,
439 for (
unsigned int i = 0; i < 4; ++i)
448 compute_shape_function_values_hardcode(
449 const unsigned int n_shape_functions,
450 const std::vector<
Point<3>> & unit_points,
453 (void)n_shape_functions;
454 const unsigned int n_points = unit_points.size();
455 for (
unsigned int k = 0; k < n_points; ++k)
457 double x = unit_points[k](0);
458 double y = unit_points[k](1);
459 double z = unit_points[k](2);
465 data.
shape(k, 0) = (1. - x) * (1. - y) * (1. - z);
466 data.
shape(k, 1) = x * (1. - y) * (1. - z);
467 data.
shape(k, 2) = (1. - x) * y * (1. - z);
468 data.
shape(k, 3) = x * y * (1. - z);
469 data.
shape(k, 4) = (1. - x) * (1. - y) * z;
470 data.
shape(k, 5) = x * (1. - y) * z;
471 data.
shape(k, 6) = (1. - x) * y * z;
472 data.
shape(k, 7) = x * y * z;
477 n_shape_functions * n_points,
479 data.
derivative(k, 0)[0] = (y - 1.) * (1. - z);
480 data.
derivative(k, 1)[0] = (1. - y) * (1. - z);
487 data.
derivative(k, 0)[1] = (x - 1.) * (1. - z);
489 data.
derivative(k, 2)[1] = (1. - x) * (1. - z);
495 data.
derivative(k, 0)[2] = (x - 1) * (1. - y);
499 data.
derivative(k, 4)[2] = (1. - x) * (1. - y);
507 n_shape_functions * n_points,
588 n_shape_functions * n_points,
591 for (
unsigned int i = 0; i < 3; ++i)
592 for (
unsigned int j = 0; j < 3; ++j)
593 for (
unsigned int l = 0;
l < 3; ++
l)
594 if ((i == j) || (j ==
l) || (
l == i))
596 for (
unsigned int m = 0; m < 8; ++m)
614 n_shape_functions * n_points,
617 for (
unsigned int i = 0; i < 8; ++i)
628 template <
int dim,
int spacedim>
630 const unsigned int polynomial_degree)
631 : polynomial_degree(polynomial_degree)
633 , line_support_points(
QGaussLobatto<1>(polynomial_degree + 1))
634 , tensor_product_quadrature(false)
639 template <
int dim,
int spacedim>
659 template <
int dim,
int spacedim>
664 const unsigned int n_original_q_points)
668 this->update_each = update_flags;
670 const unsigned int n_q_points = q.
size();
672 const bool needs_higher_order_terms =
680 covariant.resize(n_original_q_points);
683 contravariant.resize(n_original_q_points);
686 volume_elements.resize(n_original_q_points);
693 tensor_product_quadrature =
false;
699 if (tensor_product_quadrature)
701 const std::array<Quadrature<1>, dim> quad_array =
703 for (
unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
705 if (quad_array[i - 1].size() != quad_array[i].size())
707 tensor_product_quadrature =
false;
712 const std::vector<Point<1>> &points_1 =
713 quad_array[i - 1].get_points();
714 const std::vector<Point<1>> &points_2 =
715 quad_array[i].get_points();
716 const std::vector<double> &weights_1 =
717 quad_array[i - 1].get_weights();
718 const std::vector<double> &weights_2 =
719 quad_array[i].get_weights();
720 for (
unsigned int j = 0; j < quad_array[i].size(); ++j)
722 if (std::abs(points_1[j][0] - points_2[j][0]) > 1.e-10 ||
723 std::abs(weights_1[j] - weights_2[j]) > 1.e-10)
725 tensor_product_quadrature =
false;
732 if (tensor_product_quadrature)
739 shape_info.lexicographic_numbering =
740 FETools::lexicographic_to_hierarchic_numbering<dim>(
742 shape_info.n_q_points = q.
size();
743 shape_info.dofs_per_component_on_cell =
751 if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
756 shape_values.resize(n_shape_functions * n_q_points);
758 if (this->update_each &
768 shape_derivatives.resize(n_shape_functions * n_q_points);
770 if (this->update_each &
772 shape_second_derivatives.resize(n_shape_functions * n_q_points);
776 shape_third_derivatives.resize(n_shape_functions * n_q_points);
780 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
783 compute_shape_function_values(q.
get_points());
789 template <
int dim,
int spacedim>
794 const unsigned int n_original_q_points)
796 initialize(update_flags, q, n_original_q_points);
798 if (dim > 1 && tensor_product_quadrature)
800 constexpr
unsigned int facedim = dim - 1;
803 shape_info.lexicographic_numbering =
804 FETools::lexicographic_to_hierarchic_numbering<facedim>(
806 shape_info.n_q_points = n_original_q_points;
807 shape_info.dofs_per_component_on_cell =
813 if (this->update_each &
823 unit_tangentials[i].resize(n_original_q_points);
824 std::fill(unit_tangentials[i].
begin(),
825 unit_tangentials[i].
end(),
830 .resize(n_original_q_points);
848 const std::vector<
Point<1>> &unit_points)
853 internal::MappingQ1::compute_shape_function_values_hardcode(
854 n_shape_functions, unit_points, *
this);
858 internal::MappingQ1::compute_shape_function_values_general<1, 1>(
859 n_shape_functions, unit_points, *
this);
866 const std::vector<
Point<2>> &unit_points)
871 internal::MappingQ1::compute_shape_function_values_hardcode(
872 n_shape_functions, unit_points, *
this);
876 internal::MappingQ1::compute_shape_function_values_general<2, 2>(
877 n_shape_functions, unit_points, *
this);
884 const std::vector<
Point<3>> &unit_points)
889 internal::MappingQ1::compute_shape_function_values_hardcode(
890 n_shape_functions, unit_points, *
this);
894 internal::MappingQ1::compute_shape_function_values_general<3, 3>(
895 n_shape_functions, unit_points, *
this);
899 template <
int dim,
int spacedim>
906 internal::MappingQ1::compute_shape_function_values_general<dim, spacedim>(
907 n_shape_functions, unit_points, *
this);
913 namespace MappingQGenericImplementation
925 compute_support_point_weights_on_quad(
936 const unsigned int n_inner_2d = M * M;
937 const unsigned int n_outer_2d = 4 + 4 * M;
940 loqvs.
reinit(n_inner_2d, n_outer_2d);
942 for (
unsigned int i = 0; i < M; ++i)
943 for (
unsigned int j = 0; j < M; ++j)
947 const unsigned int index_table = i * M + j;
948 for (
unsigned int v = 0; v < 4; ++v)
949 loqvs(index_table, v) =
951 loqvs(index_table, 4 + i) = 1. - p[0];
952 loqvs(index_table, 4 + i + M) = p[0];
953 loqvs(index_table, 4 + j + 2 * M) = 1. - p[1];
954 loqvs(index_table, 4 + j + 3 * M) = p[1];
959 for (
unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point)
961 loqvs[unit_point].
end(),
989 const unsigned int n_inner = Utilities::fixed_power<3>(M);
990 const unsigned int n_outer = 8 + 12 * M + 6 * M * M;
993 lohvs.
reinit(n_inner, n_outer);
995 for (
unsigned int i = 0; i < M; ++i)
996 for (
unsigned int j = 0; j < M; ++j)
997 for (
unsigned int k = 0; k < M; ++k)
1000 (j + 1) * (M + 2) + (k + 1));
1001 const unsigned int index_table = i * M * M + j * M + k;
1004 for (
unsigned int v = 0; v < 8; ++v)
1005 lohvs(index_table, v) =
1010 constexpr std::array<unsigned int, 4> line_coordinates_y(
1013 for (
unsigned int l = 0;
l < 4; ++
l)
1014 lohvs(index_table, 8 + line_coordinates_y[
l] * M + j) =
1019 constexpr std::array<unsigned int, 4> line_coordinates_x(
1022 for (
unsigned int l = 0;
l < 4; ++
l)
1023 lohvs(index_table, 8 + line_coordinates_x[
l] * M + k) =
1028 constexpr std::array<unsigned int, 4> line_coordinates_z(
1031 for (
unsigned int l = 0;
l < 4; ++
l)
1032 lohvs(index_table, 8 + line_coordinates_z[
l] * M + i) =
1037 lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) =
1039 lohvs(index_table, 8 + 12 * M + 1 * M * M + i * M + j) = p[0];
1040 lohvs(index_table, 8 + 12 * M + 2 * M * M + k * M + i) =
1042 lohvs(index_table, 8 + 12 * M + 3 * M * M + k * M + i) = p[1];
1043 lohvs(index_table, 8 + 12 * M + 4 * M * M + j * M + k) =
1045 lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2];
1050 for (
unsigned int unit_point = 0; unit_point < n_inner; ++unit_point)
1052 lohvs[unit_point].
end(),
1066 std::vector<::Table<2, double>>
1067 compute_support_point_weights_perimeter_to_interior(
1069 const unsigned int dim)
1072 std::vector<::Table<2, double>> output(dim);
1104 return ::Table<2, double>();
1107 const std::vector<unsigned int> h2l =
1108 FETools::hierarchic_to_lexicographic_numbering<dim>(
1114 for (
unsigned int q = 0; q < output.
size(0); ++q)
1132 template <
int dim,
int spacedim>
1134 compute_mapped_location_of_point(
1135 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1139 data.mapping_support_points.size());
1143 for (
unsigned int i = 0; i < data.mapping_support_points.size(); ++i)
1144 p_real += data.mapping_support_points[i] * data.shape(0, i);
1156 do_transform_real_to_unit_cell_internal(
1157 const typename ::Triangulation<dim, dim>::cell_iterator &cell,
1160 typename ::MappingQGeneric<dim, dim>::InternalData &mdata)
1162 const unsigned int spacedim = dim;
1164 const unsigned int n_shapes = mdata.shape_values.size();
1169 std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
1185 mdata.compute_shape_function_values(std::vector<
Point<dim>>(1, p_unit));
1188 compute_mapped_location_of_point<dim, spacedim>(mdata);
1192 if (f.
norm_square() < 1
e-24 * cell->diameter() * cell->diameter())
1228 const double eps = 1.e-11;
1229 const unsigned int newton_iteration_limit = 20;
1231 unsigned int newton_iteration = 0;
1232 double last_f_weighted_norm;
1235 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
1236 std::cout <<
"Newton iteration " << newton_iteration << std::endl;
1241 for (
unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1246 for (
unsigned int i = 0; i < spacedim; ++i)
1247 for (
unsigned int j = 0; j < dim; ++j)
1248 df[i][j] +=
point[i] * grad_transform[j];
1259 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
1260 std::cout <<
" delta=" << delta << std::endl;
1264 double step_length = 1;
1272 for (
unsigned int i = 0; i < dim; ++i)
1273 p_unit_trial[i] -= step_length * delta[i];
1277 mdata.compute_shape_function_values(
1282 internal::MappingQGenericImplementation::
1283 compute_mapped_location_of_point<dim, spacedim>(mdata);
1286 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
1287 std::cout <<
" step_length=" << step_length << std::endl
1288 <<
" ||f || =" << f.
norm() << std::endl
1289 <<
" ||f*|| =" << f_trial.
norm() << std::endl
1291 << (df_inverse * f_trial).
norm() << std::endl;
1303 p_real = p_real_trial;
1304 p_unit = p_unit_trial;
1308 else if (step_length > 0.05)
1319 if (newton_iteration > newton_iteration_limit)
1323 last_f_weighted_norm = (df_inverse * f).
norm();
1325 while (last_f_weighted_norm >
eps);
1337 do_transform_real_to_unit_cell_internal_codim1(
1338 const typename ::Triangulation<dim, dim + 1>::cell_iterator &cell,
1341 typename ::MappingQGeneric<dim, dim + 1>::InternalData &mdata)
1343 const unsigned int spacedim = dim + 1;
1345 const unsigned int n_shapes = mdata.shape_values.size();
1349 Assert(mdata.shape_second_derivatives.size() == n_shapes,
1352 std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
1365 mdata.compute_shape_function_values(std::vector<
Point<dim>>(1, p_unit));
1367 for (
unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1370 const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
1373 for (
unsigned int j = 0; j < dim; ++j)
1375 DF[j] += grad_phi_k[j] * point_k;
1376 for (
unsigned int l = 0;
l < dim; ++
l)
1377 D2F[j][
l] += hessian_k[j][
l] * point_k;
1382 p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
1385 for (
unsigned int j = 0; j < dim; ++j)
1386 f[j] = DF[j] * p_minus_F;
1388 for (
unsigned int j = 0; j < dim; ++j)
1390 f[j] = DF[j] * p_minus_F;
1391 for (
unsigned int l = 0;
l < dim; ++
l)
1392 df[j][
l] = -DF[j] * DF[
l] + D2F[j][
l] * p_minus_F;
1396 const double eps = 1.e-12 * cell->diameter();
1397 const unsigned int loop_limit = 10;
1399 unsigned int loop = 0;
1408 for (
unsigned int j = 0; j < dim; ++j)
1411 for (
unsigned int l = 0;
l < dim; ++
l)
1415 mdata.compute_shape_function_values(
1418 for (
unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1424 for (
unsigned int j = 0; j < dim; ++j)
1426 DF[j] += grad_phi_k[j] * point_k;
1427 for (
unsigned int l = 0;
l < dim; ++
l)
1428 D2F[j][
l] += hessian_k[j][
l] * point_k;
1435 p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
1437 for (
unsigned int j = 0; j < dim; ++j)
1439 f[j] = DF[j] * p_minus_F;
1440 for (
unsigned int l = 0;
l < dim; ++
l)
1441 df[j][
l] = -DF[j] * DF[
l] + D2F[j][
l] * p_minus_F;
1462 template <
int dim,
int spacedim>
1464 maybe_update_q_points_Jacobians_and_grads_tensor(
1466 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1471 const UpdateFlags update_flags = data.update_each;
1473 const unsigned int n_shape_values = data.n_shape_functions;
1474 const unsigned int n_q_points = data.shape_info.n_q_points;
1476 constexpr
unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1477 constexpr
unsigned int n_hessians = (dim * (dim + 1)) / 2;
1480 const bool evaluate_gradients =
1483 const bool evaluate_hessians =
1490 Assert(!evaluate_gradients || data.n_shape_functions > 0,
1492 Assert(!evaluate_gradients || n_q_points == data.contravariant.size(),
1494 Assert(!evaluate_hessians || n_q_points == jacobian_grads.size(),
1499 if (evaluate_values && !evaluate_gradients & !evaluate_hessians &&
1500 data.shape_info.element_type ==
1503 for (
unsigned int q = 0; q < n_q_points; ++q)
1505 data.mapping_support_points[data.shape_info
1506 .lexicographic_numbering[q]];
1511 if (evaluate_values || evaluate_gradients || evaluate_hessians)
1513 data.values_dofs.resize(n_comp * n_shape_values);
1514 data.values_quad.resize(n_comp * n_q_points);
1515 data.gradients_quad.resize(n_comp * n_q_points * dim);
1516 data.scratch.resize(2 *
std::max(n_q_points, n_shape_values));
1518 if (evaluate_hessians)
1519 data.hessians_quad.resize(n_comp * n_q_points * n_hessians);
1521 const std::vector<unsigned int> &renumber_to_lexicographic =
1522 data.shape_info.lexicographic_numbering;
1523 for (
unsigned int i = 0; i < n_shape_values; ++i)
1524 for (
unsigned int d = 0;
d < spacedim; ++
d)
1526 const unsigned int in_comp =
d % n_lanes;
1527 const unsigned int out_comp =
d / n_lanes;
1528 data.values_dofs[out_comp * n_shape_values + i][in_comp] =
1530 .mapping_support_points[renumber_to_lexicographic[i]][
d];
1535 evaluate(data.shape_info,
1536 data.values_dofs.begin(),
1537 data.values_quad.begin(),
1538 data.gradients_quad.begin(),
1539 data.hessians_quad.begin(),
1540 data.scratch.begin(),
1547 if (evaluate_values)
1549 for (
unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1550 for (
unsigned int i = 0; i < n_q_points; ++i)
1551 for (
unsigned int in_comp = 0;
1552 in_comp < n_lanes &&
1553 in_comp < spacedim - out_comp * n_lanes;
1556 data.values_quad[out_comp * n_q_points + i][in_comp];
1559 if (evaluate_gradients)
1561 std::fill(data.contravariant.begin(),
1562 data.contravariant.end(),
1565 for (
unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1567 for (
unsigned int j = 0; j < dim; ++j)
1568 for (
unsigned int in_comp = 0;
1569 in_comp < n_lanes &&
1570 in_comp < spacedim - out_comp * n_lanes;
1573 const unsigned int total_number =
point * dim + j;
1574 const unsigned int new_comp = total_number / n_q_points;
1575 const unsigned int new_point = total_number % n_q_points;
1576 data.contravariant[new_point][out_comp * n_lanes +
1577 in_comp][new_comp] =
1578 data.gradients_quad[(out_comp * n_q_points +
point) *
1586 data.covariant[
point] =
1587 (data.contravariant[
point]).covariant_form();
1592 data.volume_elements[
point] =
1593 data.contravariant[
point].determinant();
1595 if (evaluate_hessians)
1597 constexpr
int desymmetrize_3d[6][2] = {
1598 {0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
1599 constexpr
int desymmetrize_2d[3][2] = {{0, 0}, {1, 1}, {0, 1}};
1602 for (
unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1604 for (
unsigned int j = 0; j < n_hessians; ++j)
1605 for (
unsigned int in_comp = 0;
1606 in_comp < n_lanes &&
1607 in_comp < spacedim - out_comp * n_lanes;
1610 const unsigned int total_number =
point * n_hessians + j;
1611 const unsigned int new_point = total_number % n_q_points;
1612 const unsigned int new_hessian_comp =
1613 total_number / n_q_points;
1614 const unsigned int new_hessian_comp_i =
1615 dim == 2 ? desymmetrize_2d[new_hessian_comp][0] :
1616 desymmetrize_3d[new_hessian_comp][0];
1617 const unsigned int new_hessian_comp_j =
1618 dim == 2 ? desymmetrize_2d[new_hessian_comp][1] :
1619 desymmetrize_3d[new_hessian_comp][1];
1620 const double value =
1621 data.hessians_quad[(out_comp * n_q_points +
point) *
1624 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1625 [new_hessian_comp_i][new_hessian_comp_j] =
1627 jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1628 [new_hessian_comp_j][new_hessian_comp_i] =
1641 template <
int dim,
int spacedim>
1643 maybe_compute_q_points(
1645 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1649 const UpdateFlags update_flags = data.update_each;
1655 const double * shape = &data.shape(
point + data_set, 0);
1657 (shape[0] * data.mapping_support_points[0]);
1658 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
1659 for (
unsigned int i = 0; i < spacedim; ++i)
1660 result[i] += shape[k] * data.mapping_support_points[k][i];
1675 template <
int dim,
int spacedim>
1677 maybe_update_Jacobians(
1679 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1680 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1683 const UpdateFlags update_flags = data.update_each;
1691 const unsigned int n_q_points = data.contravariant.size();
1693 std::fill(data.contravariant.begin(),
1694 data.contravariant.end(),
1700 data.mapping_support_points.data();
1705 &data.derivative(
point + data_set, 0);
1707 double result[spacedim][dim];
1711 for (
unsigned int i = 0; i < spacedim; ++i)
1712 for (
unsigned int j = 0; j < dim; ++j)
1713 result[i][j] = data_derv[0][j] * supp_pts[0][i];
1714 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
1715 for (
unsigned int i = 0; i < spacedim; ++i)
1716 for (
unsigned int j = 0; j < dim; ++j)
1717 result[i][j] += data_derv[k][j] * supp_pts[k][i];
1724 for (
unsigned int i = 0; i < spacedim; ++i)
1725 for (
unsigned int j = 0; j < dim; ++j)
1726 data.contravariant[
point][i][j] = result[i][j];
1733 const unsigned int n_q_points = data.contravariant.size();
1736 data.covariant[
point] =
1737 (data.contravariant[
point]).covariant_form();
1744 const unsigned int n_q_points = data.contravariant.size();
1746 data.volume_elements[
point] =
1747 data.contravariant[
point].determinant();
1757 template <
int dim,
int spacedim>
1759 maybe_update_jacobian_grads(
1762 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1766 const UpdateFlags update_flags = data.update_each;
1769 const unsigned int n_q_points = jacobian_grads.size();
1775 &data.second_derivative(
point + data_set, 0);
1776 double result[spacedim][dim][dim];
1777 for (
unsigned int i = 0; i < spacedim; ++i)
1778 for (
unsigned int j = 0; j < dim; ++j)
1779 for (
unsigned int l = 0;
l < dim; ++
l)
1781 (
second[0][j][
l] * data.mapping_support_points[0][i]);
1782 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
1783 for (
unsigned int i = 0; i < spacedim; ++i)
1784 for (
unsigned int j = 0; j < dim; ++j)
1785 for (
unsigned int l = 0;
l < dim; ++
l)
1788 data.mapping_support_points[k][i]);
1790 for (
unsigned int i = 0; i < spacedim; ++i)
1791 for (
unsigned int j = 0; j < dim; ++j)
1792 for (
unsigned int l = 0;
l < dim; ++
l)
1793 jacobian_grads[
point][i][j][
l] = result[i][j][
l];
1804 template <
int dim,
int spacedim>
1806 maybe_update_jacobian_pushed_forward_grads(
1809 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1813 const UpdateFlags update_flags = data.update_each;
1816 const unsigned int n_q_points =
1817 jacobian_pushed_forward_grads.size();
1821 double tmp[spacedim][spacedim][spacedim];
1825 &data.second_derivative(
point + data_set, 0);
1826 double result[spacedim][dim][dim];
1827 for (
unsigned int i = 0; i < spacedim; ++i)
1828 for (
unsigned int j = 0; j < dim; ++j)
1829 for (
unsigned int l = 0;
l < dim; ++
l)
1830 result[i][j][
l] = (
second[0][j][
l] *
1831 data.mapping_support_points[0][i]);
1832 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
1833 for (
unsigned int i = 0; i < spacedim; ++i)
1834 for (
unsigned int j = 0; j < dim; ++j)
1835 for (
unsigned int l = 0;
l < dim; ++
l)
1838 data.mapping_support_points[k][i]);
1841 for (
unsigned int i = 0; i < spacedim; ++i)
1842 for (
unsigned int j = 0; j < spacedim; ++j)
1843 for (
unsigned int l = 0;
l < dim; ++
l)
1846 result[i][0][
l] * data.covariant[
point][j][0];
1847 for (
unsigned int jr = 1; jr < dim; ++jr)
1849 tmp[i][j][
l] += result[i][jr][
l] *
1850 data.covariant[
point][j][jr];
1855 for (
unsigned int i = 0; i < spacedim; ++i)
1856 for (
unsigned int j = 0; j < spacedim; ++j)
1857 for (
unsigned int l = 0;
l < spacedim; ++
l)
1859 jacobian_pushed_forward_grads[
point][i][j][
l] =
1860 tmp[i][j][0] * data.covariant[
point][
l][0];
1861 for (
unsigned int lr = 1; lr < dim; ++lr)
1863 jacobian_pushed_forward_grads[
point][i][j][
l] +=
1864 tmp[i][j][lr] * data.covariant[
point][
l][lr];
1878 template <
int dim,
int spacedim>
1880 maybe_update_jacobian_2nd_derivatives(
1883 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1887 const UpdateFlags update_flags = data.update_each;
1890 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1897 &data.third_derivative(
point + data_set, 0);
1898 double result[spacedim][dim][dim][dim];
1899 for (
unsigned int i = 0; i < spacedim; ++i)
1900 for (
unsigned int j = 0; j < dim; ++j)
1901 for (
unsigned int l = 0;
l < dim; ++
l)
1902 for (
unsigned int m = 0; m < dim; ++m)
1903 result[i][j][
l][m] =
1904 (third[0][j][
l][m] *
1905 data.mapping_support_points[0][i]);
1906 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
1907 for (
unsigned int i = 0; i < spacedim; ++i)
1908 for (
unsigned int j = 0; j < dim; ++j)
1909 for (
unsigned int l = 0;
l < dim; ++
l)
1910 for (
unsigned int m = 0; m < dim; ++m)
1911 result[i][j][
l][m] +=
1912 (third[k][j][
l][m] *
1913 data.mapping_support_points[k][i]);
1915 for (
unsigned int i = 0; i < spacedim; ++i)
1916 for (
unsigned int j = 0; j < dim; ++j)
1917 for (
unsigned int l = 0;
l < dim; ++
l)
1918 for (
unsigned int m = 0; m < dim; ++m)
1919 jacobian_2nd_derivatives[
point][i][j][
l][m] =
1933 template <
int dim,
int spacedim>
1935 maybe_update_jacobian_pushed_forward_2nd_derivatives(
1938 const typename ::MappingQGeneric<dim, spacedim>::InternalData
1941 &jacobian_pushed_forward_2nd_derivatives)
1943 const UpdateFlags update_flags = data.update_each;
1946 const unsigned int n_q_points =
1947 jacobian_pushed_forward_2nd_derivatives.size();
1951 double tmp[spacedim][spacedim][spacedim][spacedim];
1955 &data.third_derivative(
point + data_set, 0);
1956 double result[spacedim][dim][dim][dim];
1957 for (
unsigned int i = 0; i < spacedim; ++i)
1958 for (
unsigned int j = 0; j < dim; ++j)
1959 for (
unsigned int l = 0;
l < dim; ++
l)
1960 for (
unsigned int m = 0; m < dim; ++m)
1961 result[i][j][
l][m] =
1962 (third[0][j][
l][m] *
1963 data.mapping_support_points[0][i]);
1964 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
1965 for (
unsigned int i = 0; i < spacedim; ++i)
1966 for (
unsigned int j = 0; j < dim; ++j)
1967 for (
unsigned int l = 0;
l < dim; ++
l)
1968 for (
unsigned int m = 0; m < dim; ++m)
1969 result[i][j][
l][m] +=
1970 (third[k][j][
l][m] *
1971 data.mapping_support_points[k][i]);
1974 for (
unsigned int i = 0; i < spacedim; ++i)
1975 for (
unsigned int j = 0; j < spacedim; ++j)
1976 for (
unsigned int l = 0;
l < dim; ++
l)
1977 for (
unsigned int m = 0; m < dim; ++m)
1979 jacobian_pushed_forward_2nd_derivatives
1981 result[i][0][
l][m] *
1982 data.covariant[
point][j][0];
1983 for (
unsigned int jr = 1; jr < dim; ++jr)
1984 jacobian_pushed_forward_2nd_derivatives[
point]
1987 result[i][jr][
l][m] *
1988 data.covariant[
point][j][jr];
1992 for (
unsigned int i = 0; i < spacedim; ++i)
1993 for (
unsigned int j = 0; j < spacedim; ++j)
1994 for (
unsigned int l = 0;
l < spacedim; ++
l)
1995 for (
unsigned int m = 0; m < dim; ++m)
1998 jacobian_pushed_forward_2nd_derivatives[
point]
2001 data.covariant[
point][
l][0];
2002 for (
unsigned int lr = 1; lr < dim; ++lr)
2004 jacobian_pushed_forward_2nd_derivatives
2005 [
point][i][j][lr][m] *
2006 data.covariant[
point][
l][lr];
2010 for (
unsigned int i = 0; i < spacedim; ++i)
2011 for (
unsigned int j = 0; j < spacedim; ++j)
2012 for (
unsigned int l = 0;
l < spacedim; ++
l)
2013 for (
unsigned int m = 0; m < spacedim; ++m)
2015 jacobian_pushed_forward_2nd_derivatives
2017 tmp[i][j][
l][0] * data.covariant[
point][m][0];
2018 for (
unsigned int mr = 1; mr < dim; ++mr)
2019 jacobian_pushed_forward_2nd_derivatives[
point]
2023 data.covariant[
point][m][mr];
2036 template <
int dim,
int spacedim>
2038 maybe_update_jacobian_3rd_derivatives(
2041 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2045 const UpdateFlags update_flags = data.update_each;
2048 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
2055 &data.fourth_derivative(
point + data_set, 0);
2056 double result[spacedim][dim][dim][dim][dim];
2057 for (
unsigned int i = 0; i < spacedim; ++i)
2058 for (
unsigned int j = 0; j < dim; ++j)
2059 for (
unsigned int l = 0;
l < dim; ++
l)
2060 for (
unsigned int m = 0; m < dim; ++m)
2061 for (
unsigned int n = 0; n < dim; ++n)
2062 result[i][j][
l][m][n] =
2063 (fourth[0][j][
l][m][n] *
2064 data.mapping_support_points[0][i]);
2065 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
2066 for (
unsigned int i = 0; i < spacedim; ++i)
2067 for (
unsigned int j = 0; j < dim; ++j)
2068 for (
unsigned int l = 0;
l < dim; ++
l)
2069 for (
unsigned int m = 0; m < dim; ++m)
2070 for (
unsigned int n = 0; n < dim; ++n)
2071 result[i][j][
l][m][n] +=
2072 (fourth[k][j][
l][m][n] *
2073 data.mapping_support_points[k][i]);
2075 for (
unsigned int i = 0; i < spacedim; ++i)
2076 for (
unsigned int j = 0; j < dim; ++j)
2077 for (
unsigned int l = 0;
l < dim; ++
l)
2078 for (
unsigned int m = 0; m < dim; ++m)
2079 for (
unsigned int n = 0; n < dim; ++n)
2080 jacobian_3rd_derivatives[
point][i][j][
l][m][n] =
2081 result[i][j][
l][m][n];
2094 template <
int dim,
int spacedim>
2096 maybe_update_jacobian_pushed_forward_3rd_derivatives(
2099 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2102 &jacobian_pushed_forward_3rd_derivatives)
2104 const UpdateFlags update_flags = data.update_each;
2107 const unsigned int n_q_points =
2108 jacobian_pushed_forward_3rd_derivatives.size();
2112 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
2116 &data.fourth_derivative(
point + data_set, 0);
2117 double result[spacedim][dim][dim][dim][dim];
2118 for (
unsigned int i = 0; i < spacedim; ++i)
2119 for (
unsigned int j = 0; j < dim; ++j)
2120 for (
unsigned int l = 0;
l < dim; ++
l)
2121 for (
unsigned int m = 0; m < dim; ++m)
2122 for (
unsigned int n = 0; n < dim; ++n)
2123 result[i][j][
l][m][n] =
2124 (fourth[0][j][
l][m][n] *
2125 data.mapping_support_points[0][i]);
2126 for (
unsigned int k = 1; k < data.n_shape_functions; ++k)
2127 for (
unsigned int i = 0; i < spacedim; ++i)
2128 for (
unsigned int j = 0; j < dim; ++j)
2129 for (
unsigned int l = 0;
l < dim; ++
l)
2130 for (
unsigned int m = 0; m < dim; ++m)
2131 for (
unsigned int n = 0; n < dim; ++n)
2132 result[i][j][
l][m][n] +=
2133 (fourth[k][j][
l][m][n] *
2134 data.mapping_support_points[k][i]);
2137 for (
unsigned int i = 0; i < spacedim; ++i)
2138 for (
unsigned int j = 0; j < spacedim; ++j)
2139 for (
unsigned int l = 0;
l < dim; ++
l)
2140 for (
unsigned int m = 0; m < dim; ++m)
2141 for (
unsigned int n = 0; n < dim; ++n)
2143 tmp[i][j][
l][m][n] =
2144 result[i][0][
l][m][n] *
2145 data.covariant[
point][j][0];
2146 for (
unsigned int jr = 1; jr < dim; ++jr)
2147 tmp[i][j][
l][m][n] +=
2148 result[i][jr][
l][m][n] *
2149 data.covariant[
point][j][jr];
2153 for (
unsigned int i = 0; i < spacedim; ++i)
2154 for (
unsigned int j = 0; j < spacedim; ++j)
2155 for (
unsigned int l = 0;
l < spacedim; ++
l)
2156 for (
unsigned int m = 0; m < dim; ++m)
2157 for (
unsigned int n = 0; n < dim; ++n)
2159 jacobian_pushed_forward_3rd_derivatives
2161 tmp[i][j][0][m][n] *
2162 data.covariant[
point][
l][0];
2163 for (
unsigned int lr = 1; lr < dim; ++lr)
2164 jacobian_pushed_forward_3rd_derivatives
2165 [
point][i][j][
l][m][n] +=
2166 tmp[i][j][lr][m][n] *
2167 data.covariant[
point][
l][lr];
2171 for (
unsigned int i = 0; i < spacedim; ++i)
2172 for (
unsigned int j = 0; j < spacedim; ++j)
2173 for (
unsigned int l = 0;
l < spacedim; ++
l)
2174 for (
unsigned int m = 0; m < spacedim; ++m)
2175 for (
unsigned int n = 0; n < dim; ++n)
2177 tmp[i][j][
l][m][n] =
2178 jacobian_pushed_forward_3rd_derivatives
2180 data.covariant[
point][m][0];
2181 for (
unsigned int mr = 1; mr < dim; ++mr)
2182 tmp[i][j][
l][m][n] +=
2183 jacobian_pushed_forward_3rd_derivatives
2184 [
point][i][j][
l][mr][n] *
2185 data.covariant[
point][m][mr];
2189 for (
unsigned int i = 0; i < spacedim; ++i)
2190 for (
unsigned int j = 0; j < spacedim; ++j)
2191 for (
unsigned int l = 0;
l < spacedim; ++
l)
2192 for (
unsigned int m = 0; m < spacedim; ++m)
2193 for (
unsigned int n = 0; n < spacedim; ++n)
2195 jacobian_pushed_forward_3rd_derivatives
2197 tmp[i][j][
l][m][0] *
2198 data.covariant[
point][n][0];
2199 for (
unsigned int nr = 1; nr < dim; ++nr)
2200 jacobian_pushed_forward_3rd_derivatives
2201 [
point][i][j][
l][m][n] +=
2202 tmp[i][j][
l][m][nr] *
2203 data.covariant[
point][n][nr];
2215 template <
int dim,
int spacedim>
2220 internal::MappingQGenericImplementation::
2221 compute_support_point_weights_perimeter_to_interior(
2225 internal::MappingQGenericImplementation::
2229 ExcMessage(
"It only makes sense to create polynomial mappings "
2230 "with a polynomial degree greater or equal to one."));
2235 template <
int dim,
int spacedim>
2238 : polynomial_degree(mapping.polynomial_degree)
2239 , line_support_points(mapping.line_support_points)
2240 , support_point_weights_perimeter_to_interior(
2241 mapping.support_point_weights_perimeter_to_interior)
2242 , support_point_weights_cell(mapping.support_point_weights_cell)
2247 template <
int dim,
int spacedim>
2248 std::unique_ptr<Mapping<dim, spacedim>>
2251 return std_cxx14::make_unique<MappingQGeneric<dim, spacedim>>(*this);
2256 template <
int dim,
int spacedim>
2260 return polynomial_degree;
2265 template <
int dim,
int spacedim>
2274 line_support_points.get_points()));
2275 Assert(tensor_pols.
n() == Utilities::fixed_power<dim>(polynomial_degree + 1),
2280 const std::vector<unsigned int> renumber =
2281 FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
2283 const std::vector<Point<spacedim>> support_points =
2284 this->compute_mapping_support_points(cell);
2287 for (
unsigned int i = 0; i < tensor_pols.
n(); ++i)
2289 support_points[i] * tensor_pols.
compute_value(renumber[i], p);
2291 return mapped_point;
2314 template <
int dim,
int spacedim>
2331 const Point<1> & initial_p_unit)
const
2334 const int spacedim = 1;
2341 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
2342 get_data(update_flags, point_quadrature));
2344 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
2348 return internal::MappingQGenericImplementation::
2349 do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
2357 const Point<2> & initial_p_unit)
const
2360 const int spacedim = 2;
2367 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
2368 get_data(update_flags, point_quadrature));
2370 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
2374 return internal::MappingQGenericImplementation::
2375 do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
2383 const Point<3> & initial_p_unit)
const
2386 const int spacedim = 3;
2393 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
2394 get_data(update_flags, point_quadrature));
2396 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
2400 return internal::MappingQGenericImplementation::
2401 do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
2411 const Point<1> & initial_p_unit)
const
2414 const int spacedim = 2;
2421 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
2422 get_data(update_flags, point_quadrature));
2424 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
2428 return internal::MappingQGenericImplementation::
2429 do_transform_real_to_unit_cell_internal_codim1<1>(cell,
2442 const Point<2> & initial_p_unit)
const
2445 const int spacedim = 3;
2452 auto mdata = Utilities::dynamic_unique_cast<InternalData>(
2453 get_data(update_flags, point_quadrature));
2455 mdata->mapping_support_points = this->compute_mapping_support_points(cell);
2459 return internal::MappingQGenericImplementation::
2460 do_transform_real_to_unit_cell_internal_codim1<2>(cell,
2479 template <
int dim,
int spacedim>
2487 if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
2488 ((dim == 1) || ((dim == 2) && (dim == spacedim))))
2512 vertices = this->get_vertices(cell);
2521 return internal::MappingQ1::transform_real_to_unit_cell(
2530 internal::MappingQ1::transform_real_to_unit_cell(
vertices,
2536 const double eps = 1
e-15;
2569 if (this->preserves_vertex_locations())
2570 initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
2580 std::vector<Point<spacedim>> a =
2581 this->compute_mapping_support_points(cell);
2583 std::vector<CellData<dim>> cells(1);
2585 cells[0].vertices[i] = i;
2589 tria.
begin_active()->real_to_unit_cell_affine_approximation(p);
2592 if (dim == 1 && polynomial_degree == 1)
2594 return initial_p_unit;
2607 return this->transform_real_to_unit_cell_internal(cell,
2615 template <
int dim,
int spacedim>
2627 for (
unsigned int i = 0; i < 5; ++i)
2672 template <
int dim,
int spacedim>
2673 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2677 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
2678 std_cxx14::make_unique<InternalData>(polynomial_degree);
2679 auto &data =
dynamic_cast<InternalData &
>(*data_ptr);
2680 data.initialize(this->requires_update_flags(update_flags), q, q.
size());
2687 template <
int dim,
int spacedim>
2688 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2693 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
2694 std_cxx14::make_unique<InternalData>(polynomial_degree);
2695 auto &data =
dynamic_cast<InternalData &
>(*data_ptr);
2696 data.initialize_face(this->requires_update_flags(update_flags),
2705 template <
int dim,
int spacedim>
2706 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2711 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
2712 std_cxx14::make_unique<InternalData>(polynomial_degree);
2713 auto &data =
dynamic_cast<InternalData &
>(*data_ptr);
2714 data.initialize_face(this->requires_update_flags(update_flags),
2723 template <
int dim,
int spacedim>
2734 Assert(
dynamic_cast<const InternalData *
>(&internal_data) !=
nullptr,
2736 const InternalData &data =
static_cast<const InternalData &
>(internal_data);
2738 const unsigned int n_q_points = quadrature.
size();
2750 data.mapping_support_points = this->compute_mapping_support_points(cell);
2751 data.cell_of_current_support_points = cell;
2760 if (dim > 1 && data.tensor_product_quadrature)
2762 internal::MappingQGenericImplementation::
2763 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2764 computed_cell_similarity,
2771 internal::MappingQGenericImplementation::maybe_compute_q_points<dim,
2777 internal::MappingQGenericImplementation::maybe_update_Jacobians<dim,
2779 computed_cell_similarity,
2783 internal::MappingQGenericImplementation::maybe_update_jacobian_grads<
2785 spacedim>(computed_cell_similarity,
2791 internal::MappingQGenericImplementation::
2792 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2793 computed_cell_similarity,
2798 internal::MappingQGenericImplementation::
2799 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2800 computed_cell_similarity,
2805 internal::MappingQGenericImplementation::
2806 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2807 computed_cell_similarity,
2812 internal::MappingQGenericImplementation::
2813 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2814 computed_cell_similarity,
2819 internal::MappingQGenericImplementation::
2820 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2821 computed_cell_similarity,
2826 const UpdateFlags update_flags = data.update_each;
2827 const std::vector<double> &weights = quadrature.
get_weights();
2845 if (dim == spacedim)
2847 const double det = data.contravariant[
point].determinant();
2855 1
e-12 * Utilities::fixed_power<dim>(
2856 cell->diameter() / std::sqrt(
double(dim))),
2858 cell->center(), det,
point)));
2868 for (
unsigned int i = 0; i < spacedim; ++i)
2869 for (
unsigned int j = 0; j < dim; ++j)
2870 DX_t[j][i] = data.contravariant[
point][i][j];
2873 for (
unsigned int i = 0; i < dim; ++i)
2874 for (
unsigned int j = 0; j < dim; ++j)
2875 G[i][j] = DX_t[i] * DX_t[j];
2880 if (computed_cell_similarity ==
2891 Assert(spacedim == dim + 1,
2893 "There is no (unique) cell normal for " +
2895 "-dimensional cells in " +
2897 "-dimensional space. This only works if the "
2898 "space dimension is one greater than the "
2899 "dimensionality of the mesh cells."));
2911 if (cell->direction_flag() ==
false)
2937 data.covariant[
point].transpose();
2940 return computed_cell_similarity;
2947 namespace MappingQGenericImplementation
2960 template <
int dim,
int spacedim>
2962 maybe_compute_face_data(
2963 const ::MappingQGeneric<dim, spacedim> &mapping,
2964 const typename ::Triangulation<dim, spacedim>::cell_iterator
2966 const unsigned int face_no,
2967 const unsigned int subface_no,
2968 const unsigned int n_q_points,
2969 const std::vector<double> &weights,
2970 const typename ::MappingQGeneric<dim, spacedim>::InternalData
2975 const UpdateFlags update_flags = data.update_each;
2996 for (
unsigned int d = 0;
d != dim - 1; ++
d)
2999 data.unit_tangentials.size(),
3002 data.aux[
d].size() <=
3004 .unit_tangentials[face_no +
3012 .unit_tangentials[face_no +
3023 if (dim == spacedim)
3025 for (
unsigned int i = 0; i < n_q_points; ++i)
3035 (face_no == 0 ? -1 : +1);
3065 data.contravariant[
point].transpose()[0];
3067 (face_no == 0 ? -1. : +1.) *
3078 cell_normal /= cell_normal.
norm();
3098 const double area_ratio =
3100 cell->subface_case(face_no), subface_no);
3119 data.covariant[
point].transpose();
3130 template <
int dim,
int spacedim>
3132 do_fill_fe_face_values(
3133 const ::MappingQGeneric<dim, spacedim> &mapping,
3134 const typename ::Triangulation<dim, spacedim>::cell_iterator
3136 const unsigned int face_no,
3137 const unsigned int subface_no,
3140 const typename ::MappingQGeneric<dim, spacedim>::InternalData
3145 if (dim > 1 && data.tensor_product_quadrature)
3147 maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
3155 maybe_compute_q_points<dim, spacedim>(
3160 maybe_update_jacobian_grads<dim, spacedim>(
3163 maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
3168 maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
3173 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
3178 maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
3183 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
3189 maybe_compute_face_data(mapping,
3204 template <
int dim,
int spacedim>
3208 const unsigned int face_no,
3215 Assert((
dynamic_cast<const InternalData *
>(&internal_data) !=
nullptr),
3217 const InternalData &data =
static_cast<const InternalData &
>(internal_data);
3223 if ((data.mapping_support_points.size() == 0) ||
3224 (&cell->get_triangulation() !=
3225 &data.cell_of_current_support_points->get_triangulation()) ||
3226 (cell != data.cell_of_current_support_points))
3228 data.mapping_support_points = this->compute_mapping_support_points(cell);
3229 data.cell_of_current_support_points = cell;
3232 internal::MappingQGenericImplementation::do_fill_fe_face_values(
3238 cell->face_orientation(face_no),
3239 cell->face_flip(face_no),
3240 cell->face_rotation(face_no),
3249 template <
int dim,
int spacedim>
3253 const unsigned int face_no,
3254 const unsigned int subface_no,
3261 Assert((
dynamic_cast<const InternalData *
>(&internal_data) !=
nullptr),
3263 const InternalData &data =
static_cast<const InternalData &
>(internal_data);
3269 if ((data.mapping_support_points.size() == 0) ||
3270 (&cell->get_triangulation() !=
3271 &data.cell_of_current_support_points->get_triangulation()) ||
3272 (cell != data.cell_of_current_support_points))
3274 data.mapping_support_points = this->compute_mapping_support_points(cell);
3275 data.cell_of_current_support_points = cell;
3278 internal::MappingQGenericImplementation::do_fill_fe_face_values(
3285 cell->face_orientation(face_no),
3286 cell->face_flip(face_no),
3287 cell->face_rotation(face_no),
3289 cell->subface_case(face_no)),
3299 namespace MappingQGenericImplementation
3303 template <
int dim,
int spacedim,
int rank>
3312 Assert((
dynamic_cast<const typename ::
3313 MappingQGeneric<dim, spacedim>::InternalData *
>(
3314 &mapping_data) !=
nullptr),
3316 const typename ::MappingQGeneric<dim, spacedim>::InternalData
3318 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3319 InternalData &
>(mapping_data);
3321 switch (mapping_kind)
3328 "update_contravariant_transformation"));
3330 for (
unsigned int i = 0; i < output.size(); ++i)
3342 "update_contravariant_transformation"));
3346 "update_volume_elements"));
3351 for (
unsigned int i = 0; i < output.size(); ++i)
3355 output[i] /= data.volume_elements[i];
3367 "update_covariant_transformation"));
3369 for (
unsigned int i = 0; i < output.size(); ++i)
3381 template <
int dim,
int spacedim,
int rank>
3383 transform_gradients(
3390 Assert((
dynamic_cast<const typename ::
3391 MappingQGeneric<dim, spacedim>::InternalData *
>(
3392 &mapping_data) !=
nullptr),
3394 const typename ::MappingQGeneric<dim, spacedim>::InternalData
3396 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3397 InternalData &
>(mapping_data);
3399 switch (mapping_kind)
3406 "update_covariant_transformation"));
3410 "update_contravariant_transformation"));
3413 for (
unsigned int i = 0; i < output.size(); ++i)
3430 "update_covariant_transformation"));
3433 for (
unsigned int i = 0; i < output.size(); ++i)
3450 "update_covariant_transformation"));
3454 "update_contravariant_transformation"));
3458 "update_volume_elements"));
3461 for (
unsigned int i = 0; i < output.size(); ++i)
3470 output[i] /= data.volume_elements[i];
3483 template <
int dim,
int spacedim>
3492 Assert((
dynamic_cast<const typename ::
3493 MappingQGeneric<dim, spacedim>::InternalData *
>(
3494 &mapping_data) !=
nullptr),
3496 const typename ::MappingQGeneric<dim, spacedim>::InternalData
3498 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3499 InternalData &
>(mapping_data);
3501 switch (mapping_kind)
3508 "update_covariant_transformation"));
3512 "update_contravariant_transformation"));
3514 for (
unsigned int q = 0; q < output.size(); ++q)
3515 for (
unsigned int i = 0; i < spacedim; ++i)
3517 double tmp1[dim][dim];
3518 for (
unsigned int J = 0; J < dim; ++J)
3519 for (
unsigned int K = 0; K < dim; ++K)
3522 data.contravariant[q][i][0] * input[q][0][J][K];
3523 for (
unsigned int I = 1; I < dim; ++I)
3525 data.contravariant[q][i][I] * input[q][I][J][K];
3527 for (
unsigned int j = 0; j < spacedim; ++j)
3530 for (
unsigned int K = 0; K < dim; ++K)
3532 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
3533 for (
unsigned int J = 1; J < dim; ++J)
3534 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
3536 for (
unsigned int k = 0; k < spacedim; ++k)
3538 output[q][i][j][k] =
3539 data.covariant[q][k][0] * tmp2[0];
3540 for (
unsigned int K = 1; K < dim; ++K)
3541 output[q][i][j][k] +=
3542 data.covariant[q][k][K] * tmp2[K];
3554 "update_covariant_transformation"));
3556 for (
unsigned int q = 0; q < output.size(); ++q)
3557 for (
unsigned int i = 0; i < spacedim; ++i)
3559 double tmp1[dim][dim];
3560 for (
unsigned int J = 0; J < dim; ++J)
3561 for (
unsigned int K = 0; K < dim; ++K)
3564 data.covariant[q][i][0] * input[q][0][J][K];
3565 for (
unsigned int I = 1; I < dim; ++I)
3567 data.covariant[q][i][I] * input[q][I][J][K];
3569 for (
unsigned int j = 0; j < spacedim; ++j)
3572 for (
unsigned int K = 0; K < dim; ++K)
3574 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
3575 for (
unsigned int J = 1; J < dim; ++J)
3576 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
3578 for (
unsigned int k = 0; k < spacedim; ++k)
3580 output[q][i][j][k] =
3581 data.covariant[q][k][0] * tmp2[0];
3582 for (
unsigned int K = 1; K < dim; ++K)
3583 output[q][i][j][k] +=
3584 data.covariant[q][k][K] * tmp2[K];
3597 "update_covariant_transformation"));
3601 "update_contravariant_transformation"));
3605 "update_volume_elements"));
3607 for (
unsigned int q = 0; q < output.size(); ++q)
3608 for (
unsigned int i = 0; i < spacedim; ++i)
3611 for (
unsigned int I = 0; I < dim; ++I)
3613 data.contravariant[q][i][I] / data.volume_elements[q];
3614 double tmp1[dim][dim];
3615 for (
unsigned int J = 0; J < dim; ++J)
3616 for (
unsigned int K = 0; K < dim; ++K)
3618 tmp1[J][K] = factor[0] * input[q][0][J][K];
3619 for (
unsigned int I = 1; I < dim; ++I)
3620 tmp1[J][K] += factor[I] * input[q][I][J][K];
3622 for (
unsigned int j = 0; j < spacedim; ++j)
3625 for (
unsigned int K = 0; K < dim; ++K)
3627 tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
3628 for (
unsigned int J = 1; J < dim; ++J)
3629 tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
3631 for (
unsigned int k = 0; k < spacedim; ++k)
3633 output[q][i][j][k] =
3634 data.covariant[q][k][0] * tmp2[0];
3635 for (
unsigned int K = 1; K < dim; ++K)
3636 output[q][i][j][k] +=
3637 data.covariant[q][k][K] * tmp2[K];
3652 template <
int dim,
int spacedim,
int rank>
3654 transform_differential_forms(
3661 Assert((
dynamic_cast<const typename ::
3662 MappingQGeneric<dim, spacedim>::InternalData *
>(
3663 &mapping_data) !=
nullptr),
3665 const typename ::MappingQGeneric<dim, spacedim>::InternalData
3667 static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3668 InternalData &
>(mapping_data);
3670 switch (mapping_kind)
3677 "update_covariant_transformation"));
3679 for (
unsigned int i = 0; i < output.size(); ++i)
3694 template <
int dim,
int spacedim>
3702 internal::MappingQGenericImplementation::transform_fields(input,
3710 template <
int dim,
int spacedim>
3718 internal::MappingQGenericImplementation::transform_differential_forms(
3719 input, mapping_kind, mapping_data, output);
3724 template <
int dim,
int spacedim>
3732 switch (mapping_kind)
3735 internal::MappingQGenericImplementation::transform_fields(input,
3744 internal::MappingQGenericImplementation::transform_gradients(
3745 input, mapping_kind, mapping_data, output);
3754 template <
int dim,
int spacedim>
3763 Assert(
dynamic_cast<const InternalData *
>(&mapping_data) !=
nullptr,
3765 const InternalData &data =
static_cast<const InternalData &
>(mapping_data);
3767 switch (mapping_kind)
3773 "update_covariant_transformation"));
3775 for (
unsigned int q = 0; q < output.size(); ++q)
3776 for (
unsigned int i = 0; i < spacedim; ++i)
3777 for (
unsigned int j = 0; j < spacedim; ++j)
3780 for (
unsigned int K = 0; K < dim; ++K)
3782 tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
3783 for (
unsigned int J = 1; J < dim; ++J)
3784 tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
3786 for (
unsigned int k = 0; k < spacedim; ++k)
3788 output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
3789 for (
unsigned int K = 1; K < dim; ++K)
3790 output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
3803 template <
int dim,
int spacedim>
3811 switch (mapping_kind)
3816 internal::MappingQGenericImplementation::transform_hessians(
3817 input, mapping_kind, mapping_data, output);
3826 template <
int dim,
int spacedim>
3833 if (this->polynomial_degree == 2)
3835 for (
unsigned int line_no = 0;
3836 line_no < GeometryInfo<dim>::lines_per_cell;
3843 cell->line(line_no));
3848 cell->get_manifold() :
3857 std::vector<Point<spacedim>> tmp_points;
3858 for (
unsigned int line_no = 0;
3859 line_no < GeometryInfo<dim>::lines_per_cell;
3866 cell->line(line_no));
3871 cell->get_manifold() :
3874 const std::array<Point<spacedim>, 2>
vertices{
3879 const std::size_t n_rows =
3880 support_point_weights_perimeter_to_interior[0].size(0);
3881 a.resize(a.size() + n_rows);
3885 support_point_weights_perimeter_to_interior[0],
3902 std::vector<Point<3>> tmp_points;
3905 for (
unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
3910 const bool face_orientation = cell->face_orientation(face_no),
3911 face_flip = cell->face_flip(face_no),
3912 face_rotation = cell->face_rotation(face_no);
3917 for (
unsigned int i = 0; i < vertices_per_face; ++i)
3918 Assert(face->vertex_index(i) ==
3920 face_no, i, face_orientation, face_flip, face_rotation)),
3925 for (
unsigned int i = 0; i < lines_per_face; ++i)
3928 face_no, i, face_orientation, face_flip, face_rotation)),
3934 boost::container::small_vector<Point<3>, 200> tmp_points(
3939 if (polynomial_degree > 1)
3940 for (
unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
3942 for (
unsigned int i = 0; i < polynomial_degree - 1; ++i)
3943 tmp_points[4 + line * (polynomial_degree - 1) + i] =
3945 (polynomial_degree - 1) *
3949 const std::size_t n_rows =
3950 support_point_weights_perimeter_to_interior[1].size(0);
3951 a.resize(a.size() + n_rows);
3953 face->get_manifold().get_new_points(
3955 support_point_weights_perimeter_to_interior[1],
3974 for (
unsigned int q = 0, q2 = 0; q2 < polynomial_degree - 1; ++q2)
3975 for (
unsigned int q1 = 0; q1 < polynomial_degree - 1; ++q1, ++q)
3978 line_support_points.point(q2 + 1)[0]);
3983 const std::size_t n_rows = weights.
size(0);
3984 a.resize(a.size() + n_rows);
3986 cell->get_manifold().get_new_points(
3992 template <
int dim,
int spacedim>
4003 template <
int dim,
int spacedim>
4004 std::vector<Point<spacedim>>
4009 std::vector<Point<spacedim>> a;
4010 a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
4012 a.push_back(cell->vertex(i));
4014 if (this->polynomial_degree > 1)
4021 bool all_manifold_ids_are_equal = (dim == spacedim);
4022 if (all_manifold_ids_are_equal &&
4024 &cell->get_manifold()) ==
nullptr)
4027 if (&cell->face(f)->get_manifold() != &cell->get_manifold())
4028 all_manifold_ids_are_equal =
false;
4031 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
4032 if (&cell->line(
l)->get_manifold() != &cell->get_manifold())
4033 all_manifold_ids_are_equal =
false;
4036 if (all_manifold_ids_are_equal)
4038 const std::size_t n_rows = support_point_weights_cell.size(0);
4039 a.resize(a.size() + n_rows);
4043 support_point_weights_cell,
4050 add_line_support_points(cell, a);
4055 add_line_support_points(cell, a);
4058 if (dim != spacedim)
4059 add_quad_support_points(cell, a);
4062 const std::size_t n_rows =
4063 support_point_weights_perimeter_to_interior[1].size(0);
4064 a.resize(a.size() + n_rows);
4066 cell->get_manifold().get_new_points(
4068 support_point_weights_perimeter_to_interior[1],
4075 add_line_support_points(cell, a);
4076 add_quad_support_points(cell, a);
4080 const std::size_t n_rows =
4081 support_point_weights_perimeter_to_interior[2].size(0);
4082 a.resize(a.size() + n_rows);
4084 cell->get_manifold().get_new_points(
4086 support_point_weights_perimeter_to_interior[2],
4103 #include "mapping_q_generic.inst"