Reference documentation for deal.II version 9.2.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_generic.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2020 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 
24 #include <deal.II/base/table.h>
26 
27 #include <deal.II/fe/fe_dgq.h>
28 #include <deal.II/fe/fe_tools.h>
29 #include <deal.II/fe/fe_values.h>
30 #include <deal.II/fe/mapping_q1.h>
32 
34 #include <deal.II/grid/tria.h>
36 
39 
44 
45 #include <algorithm>
46 #include <array>
47 #include <cmath>
48 #include <memory>
49 #include <numeric>
50 
51 
53 
54 
55 namespace internal
56 {
57  namespace MappingQ1
58  {
59  namespace
60  {
61  // These are left as templates on the spatial dimension (even though dim
62  // == spacedim must be true for them to make sense) because templates are
63  // expanded before the compiler eliminates code due to the 'if (dim ==
64  // spacedim)' statement (see the body of the general
65  // transform_real_to_unit_cell).
66  template <int spacedim>
67  Point<1>
68  transform_real_to_unit_cell(
70  & vertices,
71  const Point<spacedim> &p)
72  {
73  Assert(spacedim == 1, ExcInternalError());
74  return Point<1>((p[0] - vertices[0](0)) /
75  (vertices[1](0) - vertices[0](0)));
76  }
77 
78 
79 
80  template <int spacedim>
81  Point<2>
82  transform_real_to_unit_cell(
84  & vertices,
85  const Point<spacedim> &p)
86  {
87  Assert(spacedim == 2, ExcInternalError());
88 
89  // For accuracy reasons, we do all arithmetic in extended precision
90  // (long double). This has a noticeable effect on the hit rate for
91  // borderline cases and thus makes the algorithm more robust.
92  const long double x = p(0);
93  const long double y = p(1);
94 
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);
99 
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);
104 
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;
110 
111  const long double discriminant = b * b - 4 * a * c;
112  // exit if the point is not in the cell (this is the only case where the
113  // discriminant is negative)
114  AssertThrow(
115  discriminant > 0.0,
117 
118  long double eta1;
119  long double eta2;
120  const long double sqrt_discriminant = std::sqrt(discriminant);
121  // special case #1: if a is near-zero to make the discriminant exactly
122  // equal b, then use the linear formula
123  if (b != 0.0 && std::abs(b) == sqrt_discriminant)
124  {
125  eta1 = -c / b;
126  eta2 = -c / b;
127  }
128  // special case #2: a is zero for parallelograms and very small for
129  // near-parallelograms:
130  else if (std::abs(a) < 1e-8 * std::abs(b))
131  {
132  // if both a and c are very small then the root should be near
133  // zero: this first case will capture that
134  eta1 = 2 * c / (-b - sqrt_discriminant);
135  eta2 = 2 * c / (-b + sqrt_discriminant);
136  }
137  // finally, use the plain version:
138  else
139  {
140  eta1 = (-b - sqrt_discriminant) / (2 * a);
141  eta2 = (-b + sqrt_discriminant) / (2 * a);
142  }
143  // pick the one closer to the center of the cell.
144  const long double eta =
145  (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
146 
147  /*
148  * There are two ways to compute xi from eta, but either one may have a
149  * zero denominator.
150  */
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 =
155  std::max(std::max(std::abs(x0), std::abs(x1)),
156  std::max(std::abs(x2), std::abs(x3)));
157 
158  if (std::abs(xi_denominator0) > 1e-10 * max_x)
159  {
160  const double xi = (x + subexpr0) / xi_denominator0;
161  return {xi, static_cast<double>(eta)};
162  }
163  else
164  {
165  const long double max_y =
166  std::max(std::max(std::abs(y0), std::abs(y1)),
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) > 1e-10 * max_y)
172  {
173  const double xi = (subexpr1 + y) / xi_denominator1;
174  return {xi, static_cast<double>(eta)};
175  }
176  else // give up and try Newton iteration
177  {
178  AssertThrow(
179  false,
180  (typename Mapping<spacedim,
181  spacedim>::ExcTransformationFailed()));
182  }
183  }
184  // bogus return to placate compiler. It should not be possible to get
185  // here.
186  Assert(false, ExcInternalError());
187  return {std::numeric_limits<double>::quiet_NaN(),
188  std::numeric_limits<double>::quiet_NaN()};
189  }
190 
191 
192 
193  template <int spacedim>
194  Point<3>
195  transform_real_to_unit_cell(
197  & /*vertices*/,
198  const Point<spacedim> & /*p*/)
199  {
200  // It should not be possible to get here
201  Assert(false, ExcInternalError());
202  return Point<3>();
203  }
204 
205 
206 
207  template <int dim, int spacedim>
208  void
209  compute_shape_function_values_general(
210  const unsigned int n_shape_functions,
211  const std::vector<Point<dim>> &unit_points,
212  typename ::MappingQGeneric<dim, spacedim>::InternalData &data)
213  {
214  const unsigned int n_points = unit_points.size();
215 
216  // Construct the tensor product polynomials used as shape functions for
217  // the Qp mapping of cells at the boundary.
218  const TensorProductPolynomials<dim> tensor_pols(
220  data.line_support_points.get_points()));
221  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
222 
223  // then also construct the mapping from lexicographic to the Qp shape
224  // function numbering
225  const std::vector<unsigned int> renumber =
226  FETools::hierarchic_to_lexicographic_numbering<dim>(
227  data.polynomial_degree);
228 
229  std::vector<double> values;
230  std::vector<Tensor<1, dim>> grads;
231  if (data.shape_values.size() != 0)
232  {
233  Assert(data.shape_values.size() == n_shape_functions * n_points,
234  ExcInternalError());
235  values.resize(n_shape_functions);
236  }
237  if (data.shape_derivatives.size() != 0)
238  {
239  Assert(data.shape_derivatives.size() ==
240  n_shape_functions * n_points,
241  ExcInternalError());
242  grads.resize(n_shape_functions);
243  }
244 
245  std::vector<Tensor<2, dim>> grad2;
246  if (data.shape_second_derivatives.size() != 0)
247  {
248  Assert(data.shape_second_derivatives.size() ==
249  n_shape_functions * n_points,
250  ExcInternalError());
251  grad2.resize(n_shape_functions);
252  }
253 
254  std::vector<Tensor<3, dim>> grad3;
255  if (data.shape_third_derivatives.size() != 0)
256  {
257  Assert(data.shape_third_derivatives.size() ==
258  n_shape_functions * n_points,
259  ExcInternalError());
260  grad3.resize(n_shape_functions);
261  }
262 
263  std::vector<Tensor<4, dim>> grad4;
264  if (data.shape_fourth_derivatives.size() != 0)
265  {
266  Assert(data.shape_fourth_derivatives.size() ==
267  n_shape_functions * n_points,
268  ExcInternalError());
269  grad4.resize(n_shape_functions);
270  }
271 
272 
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)
278  for (unsigned int point = 0; point < n_points; ++point)
279  {
280  tensor_pols.evaluate(
281  unit_points[point], values, grads, grad2, grad3, grad4);
282 
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]];
286 
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]];
290 
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]];
294 
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]];
298 
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]];
302  }
303  }
304 
305 
306  void
307  compute_shape_function_values_hardcode(
308  const unsigned int n_shape_functions,
309  const std::vector<Point<1>> & unit_points,
311  {
312  (void)n_shape_functions;
313  const unsigned int n_points = unit_points.size();
314  for (unsigned int k = 0; k < n_points; ++k)
315  {
316  double x = unit_points[k](0);
317 
318  if (data.shape_values.size() != 0)
319  {
320  Assert(data.shape_values.size() == n_shape_functions * n_points,
321  ExcInternalError());
322  data.shape(k, 0) = 1. - x;
323  data.shape(k, 1) = x;
324  }
325  if (data.shape_derivatives.size() != 0)
326  {
327  Assert(data.shape_derivatives.size() ==
328  n_shape_functions * n_points,
329  ExcInternalError());
330  data.derivative(k, 0)[0] = -1.;
331  data.derivative(k, 1)[0] = 1.;
332  }
333  if (data.shape_second_derivatives.size() != 0)
334  {
335  Assert(data.shape_second_derivatives.size() ==
336  n_shape_functions * n_points,
337  ExcInternalError());
338  data.second_derivative(k, 0)[0][0] = 0;
339  data.second_derivative(k, 1)[0][0] = 0;
340  }
341  if (data.shape_third_derivatives.size() != 0)
342  {
343  Assert(data.shape_third_derivatives.size() ==
344  n_shape_functions * n_points,
345  ExcInternalError());
346 
348  data.third_derivative(k, 0) = zero;
349  data.third_derivative(k, 1) = zero;
350  }
351  if (data.shape_fourth_derivatives.size() != 0)
352  {
353  Assert(data.shape_fourth_derivatives.size() ==
354  n_shape_functions * n_points,
355  ExcInternalError());
356 
358  data.fourth_derivative(k, 0) = zero;
359  data.fourth_derivative(k, 1) = zero;
360  }
361  }
362  }
363 
364 
365  void
366  compute_shape_function_values_hardcode(
367  const unsigned int n_shape_functions,
368  const std::vector<Point<2>> & unit_points,
370  {
371  (void)n_shape_functions;
372  const unsigned int n_points = unit_points.size();
373  for (unsigned int k = 0; k < n_points; ++k)
374  {
375  double x = unit_points[k](0);
376  double y = unit_points[k](1);
377 
378  if (data.shape_values.size() != 0)
379  {
380  Assert(data.shape_values.size() == n_shape_functions * n_points,
381  ExcInternalError());
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;
386  }
387  if (data.shape_derivatives.size() != 0)
388  {
389  Assert(data.shape_derivatives.size() ==
390  n_shape_functions * n_points,
391  ExcInternalError());
392  data.derivative(k, 0)[0] = (y - 1.);
393  data.derivative(k, 1)[0] = (1. - y);
394  data.derivative(k, 2)[0] = -y;
395  data.derivative(k, 3)[0] = y;
396  data.derivative(k, 0)[1] = (x - 1.);
397  data.derivative(k, 1)[1] = -x;
398  data.derivative(k, 2)[1] = (1. - x);
399  data.derivative(k, 3)[1] = x;
400  }
401  if (data.shape_second_derivatives.size() != 0)
402  {
403  Assert(data.shape_second_derivatives.size() ==
404  n_shape_functions * n_points,
405  ExcInternalError());
406  data.second_derivative(k, 0)[0][0] = 0;
407  data.second_derivative(k, 1)[0][0] = 0;
408  data.second_derivative(k, 2)[0][0] = 0;
409  data.second_derivative(k, 3)[0][0] = 0;
410  data.second_derivative(k, 0)[0][1] = 1.;
411  data.second_derivative(k, 1)[0][1] = -1.;
412  data.second_derivative(k, 2)[0][1] = -1.;
413  data.second_derivative(k, 3)[0][1] = 1.;
414  data.second_derivative(k, 0)[1][0] = 1.;
415  data.second_derivative(k, 1)[1][0] = -1.;
416  data.second_derivative(k, 2)[1][0] = -1.;
417  data.second_derivative(k, 3)[1][0] = 1.;
418  data.second_derivative(k, 0)[1][1] = 0;
419  data.second_derivative(k, 1)[1][1] = 0;
420  data.second_derivative(k, 2)[1][1] = 0;
421  data.second_derivative(k, 3)[1][1] = 0;
422  }
423  if (data.shape_third_derivatives.size() != 0)
424  {
425  Assert(data.shape_third_derivatives.size() ==
426  n_shape_functions * n_points,
427  ExcInternalError());
428 
430  for (unsigned int i = 0; i < 4; ++i)
431  data.third_derivative(k, i) = zero;
432  }
433  if (data.shape_fourth_derivatives.size() != 0)
434  {
435  Assert(data.shape_fourth_derivatives.size() ==
436  n_shape_functions * n_points,
437  ExcInternalError());
439  for (unsigned int i = 0; i < 4; ++i)
440  data.fourth_derivative(k, i) = zero;
441  }
442  }
443  }
444 
445 
446 
447  void
448  compute_shape_function_values_hardcode(
449  const unsigned int n_shape_functions,
450  const std::vector<Point<3>> & unit_points,
452  {
453  (void)n_shape_functions;
454  const unsigned int n_points = unit_points.size();
455  for (unsigned int k = 0; k < n_points; ++k)
456  {
457  double x = unit_points[k](0);
458  double y = unit_points[k](1);
459  double z = unit_points[k](2);
460 
461  if (data.shape_values.size() != 0)
462  {
463  Assert(data.shape_values.size() == n_shape_functions * n_points,
464  ExcInternalError());
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;
473  }
474  if (data.shape_derivatives.size() != 0)
475  {
476  Assert(data.shape_derivatives.size() ==
477  n_shape_functions * n_points,
478  ExcInternalError());
479  data.derivative(k, 0)[0] = (y - 1.) * (1. - z);
480  data.derivative(k, 1)[0] = (1. - y) * (1. - z);
481  data.derivative(k, 2)[0] = -y * (1. - z);
482  data.derivative(k, 3)[0] = y * (1. - z);
483  data.derivative(k, 4)[0] = (y - 1.) * z;
484  data.derivative(k, 5)[0] = (1. - y) * z;
485  data.derivative(k, 6)[0] = -y * z;
486  data.derivative(k, 7)[0] = y * z;
487  data.derivative(k, 0)[1] = (x - 1.) * (1. - z);
488  data.derivative(k, 1)[1] = -x * (1. - z);
489  data.derivative(k, 2)[1] = (1. - x) * (1. - z);
490  data.derivative(k, 3)[1] = x * (1. - z);
491  data.derivative(k, 4)[1] = (x - 1.) * z;
492  data.derivative(k, 5)[1] = -x * z;
493  data.derivative(k, 6)[1] = (1. - x) * z;
494  data.derivative(k, 7)[1] = x * z;
495  data.derivative(k, 0)[2] = (x - 1) * (1. - y);
496  data.derivative(k, 1)[2] = x * (y - 1.);
497  data.derivative(k, 2)[2] = (x - 1.) * y;
498  data.derivative(k, 3)[2] = -x * y;
499  data.derivative(k, 4)[2] = (1. - x) * (1. - y);
500  data.derivative(k, 5)[2] = x * (1. - y);
501  data.derivative(k, 6)[2] = (1. - x) * y;
502  data.derivative(k, 7)[2] = x * y;
503  }
504  if (data.shape_second_derivatives.size() != 0)
505  {
506  Assert(data.shape_second_derivatives.size() ==
507  n_shape_functions * n_points,
508  ExcInternalError());
509  data.second_derivative(k, 0)[0][0] = 0;
510  data.second_derivative(k, 1)[0][0] = 0;
511  data.second_derivative(k, 2)[0][0] = 0;
512  data.second_derivative(k, 3)[0][0] = 0;
513  data.second_derivative(k, 4)[0][0] = 0;
514  data.second_derivative(k, 5)[0][0] = 0;
515  data.second_derivative(k, 6)[0][0] = 0;
516  data.second_derivative(k, 7)[0][0] = 0;
517  data.second_derivative(k, 0)[1][1] = 0;
518  data.second_derivative(k, 1)[1][1] = 0;
519  data.second_derivative(k, 2)[1][1] = 0;
520  data.second_derivative(k, 3)[1][1] = 0;
521  data.second_derivative(k, 4)[1][1] = 0;
522  data.second_derivative(k, 5)[1][1] = 0;
523  data.second_derivative(k, 6)[1][1] = 0;
524  data.second_derivative(k, 7)[1][1] = 0;
525  data.second_derivative(k, 0)[2][2] = 0;
526  data.second_derivative(k, 1)[2][2] = 0;
527  data.second_derivative(k, 2)[2][2] = 0;
528  data.second_derivative(k, 3)[2][2] = 0;
529  data.second_derivative(k, 4)[2][2] = 0;
530  data.second_derivative(k, 5)[2][2] = 0;
531  data.second_derivative(k, 6)[2][2] = 0;
532  data.second_derivative(k, 7)[2][2] = 0;
533 
534  data.second_derivative(k, 0)[0][1] = (1. - z);
535  data.second_derivative(k, 1)[0][1] = -(1. - z);
536  data.second_derivative(k, 2)[0][1] = -(1. - z);
537  data.second_derivative(k, 3)[0][1] = (1. - z);
538  data.second_derivative(k, 4)[0][1] = z;
539  data.second_derivative(k, 5)[0][1] = -z;
540  data.second_derivative(k, 6)[0][1] = -z;
541  data.second_derivative(k, 7)[0][1] = z;
542  data.second_derivative(k, 0)[1][0] = (1. - z);
543  data.second_derivative(k, 1)[1][0] = -(1. - z);
544  data.second_derivative(k, 2)[1][0] = -(1. - z);
545  data.second_derivative(k, 3)[1][0] = (1. - z);
546  data.second_derivative(k, 4)[1][0] = z;
547  data.second_derivative(k, 5)[1][0] = -z;
548  data.second_derivative(k, 6)[1][0] = -z;
549  data.second_derivative(k, 7)[1][0] = z;
550 
551  data.second_derivative(k, 0)[0][2] = (1. - y);
552  data.second_derivative(k, 1)[0][2] = -(1. - y);
553  data.second_derivative(k, 2)[0][2] = y;
554  data.second_derivative(k, 3)[0][2] = -y;
555  data.second_derivative(k, 4)[0][2] = -(1. - y);
556  data.second_derivative(k, 5)[0][2] = (1. - y);
557  data.second_derivative(k, 6)[0][2] = -y;
558  data.second_derivative(k, 7)[0][2] = y;
559  data.second_derivative(k, 0)[2][0] = (1. - y);
560  data.second_derivative(k, 1)[2][0] = -(1. - y);
561  data.second_derivative(k, 2)[2][0] = y;
562  data.second_derivative(k, 3)[2][0] = -y;
563  data.second_derivative(k, 4)[2][0] = -(1. - y);
564  data.second_derivative(k, 5)[2][0] = (1. - y);
565  data.second_derivative(k, 6)[2][0] = -y;
566  data.second_derivative(k, 7)[2][0] = y;
567 
568  data.second_derivative(k, 0)[1][2] = (1. - x);
569  data.second_derivative(k, 1)[1][2] = x;
570  data.second_derivative(k, 2)[1][2] = -(1. - x);
571  data.second_derivative(k, 3)[1][2] = -x;
572  data.second_derivative(k, 4)[1][2] = -(1. - x);
573  data.second_derivative(k, 5)[1][2] = -x;
574  data.second_derivative(k, 6)[1][2] = (1. - x);
575  data.second_derivative(k, 7)[1][2] = x;
576  data.second_derivative(k, 0)[2][1] = (1. - x);
577  data.second_derivative(k, 1)[2][1] = x;
578  data.second_derivative(k, 2)[2][1] = -(1. - x);
579  data.second_derivative(k, 3)[2][1] = -x;
580  data.second_derivative(k, 4)[2][1] = -(1. - x);
581  data.second_derivative(k, 5)[2][1] = -x;
582  data.second_derivative(k, 6)[2][1] = (1. - x);
583  data.second_derivative(k, 7)[2][1] = x;
584  }
585  if (data.shape_third_derivatives.size() != 0)
586  {
587  Assert(data.shape_third_derivatives.size() ==
588  n_shape_functions * n_points,
589  ExcInternalError());
590 
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))
595  {
596  for (unsigned int m = 0; m < 8; ++m)
597  data.third_derivative(k, m)[i][j][l] = 0;
598  }
599  else
600  {
601  data.third_derivative(k, 0)[i][j][l] = -1.;
602  data.third_derivative(k, 1)[i][j][l] = 1.;
603  data.third_derivative(k, 2)[i][j][l] = 1.;
604  data.third_derivative(k, 3)[i][j][l] = -1.;
605  data.third_derivative(k, 4)[i][j][l] = 1.;
606  data.third_derivative(k, 5)[i][j][l] = -1.;
607  data.third_derivative(k, 6)[i][j][l] = -1.;
608  data.third_derivative(k, 7)[i][j][l] = 1.;
609  }
610  }
611  if (data.shape_fourth_derivatives.size() != 0)
612  {
613  Assert(data.shape_fourth_derivatives.size() ==
614  n_shape_functions * n_points,
615  ExcInternalError());
617  for (unsigned int i = 0; i < 8; ++i)
618  data.fourth_derivative(k, i) = zero;
619  }
620  }
621  }
622  } // namespace
623  } // namespace MappingQ1
624 } // namespace internal
625 
626 
627 
628 template <int dim, int spacedim>
630  const unsigned int polynomial_degree)
631  : polynomial_degree(polynomial_degree)
632  , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
633  , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
634  , tensor_product_quadrature(false)
635 {}
636 
637 
638 
639 template <int dim, int spacedim>
640 std::size_t
642 {
643  return (
656 }
657 
658 
659 template <int dim, int spacedim>
660 void
662  const UpdateFlags update_flags,
663  const Quadrature<dim> &q,
664  const unsigned int n_original_q_points)
665 {
666  // store the flags in the internal data object so we can access them
667  // in fill_fe_*_values()
668  this->update_each = update_flags;
669 
670  const unsigned int n_q_points = q.size();
671 
672  const bool needs_higher_order_terms =
673  this->update_each &
678 
680  covariant.resize(n_original_q_points);
681 
683  contravariant.resize(n_original_q_points);
684 
686  volume_elements.resize(n_original_q_points);
687 
689 
690  // use of MatrixFree only for higher order elements and with more than one
691  // point where tensor products do not make sense
692  if (polynomial_degree < 2 || n_q_points == 1)
694 
695  if (dim > 1)
696  {
697  // find out if the one-dimensional formula is the same
698  // in all directions
700  {
701  const std::array<Quadrature<1>, dim> quad_array =
702  q.get_tensor_basis();
703  for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
704  {
705  if (quad_array[i - 1].size() != quad_array[i].size())
706  {
707  tensor_product_quadrature = false;
708  break;
709  }
710  else
711  {
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)
721  {
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)
724  {
725  tensor_product_quadrature = false;
726  break;
727  }
728  }
729  }
730  }
731 
732  if (tensor_product_quadrature)
733  {
734  // use a 1D FE_DGQ and adjust the hierarchic -> lexicographic
735  // numbering manually (building an FE_Q<dim> is relatively
736  // expensive due to constraints)
737  const FE_DGQ<1> fe(polynomial_degree);
738  shape_info.reinit(q.get_tensor_basis()[0], fe);
740  FETools::lexicographic_to_hierarchic_numbering<dim>(
742  shape_info.n_q_points = q.size();
745  }
746  }
747  }
748 
749  // Only fill the big arrays on demand in case we cannot use the tensor
750  // product quadrature code path
751  if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
752  {
753  // see if we need the (transformation) shape function values
754  // and/or gradients and resize the necessary arrays
756  shape_values.resize(n_shape_functions * n_q_points);
757 
758  if (this->update_each &
759  (update_covariant_transformation |
760  update_contravariant_transformation | update_JxW_values |
768  shape_derivatives.resize(n_shape_functions * n_q_points);
769 
770  if (this->update_each &
772  shape_second_derivatives.resize(n_shape_functions * n_q_points);
773 
776  shape_third_derivatives.resize(n_shape_functions * n_q_points);
777 
780  shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
781 
782  // now also fill the various fields with their correct values
784  }
785 }
786 
787 
788 
789 template <int dim, int spacedim>
790 void
792  const UpdateFlags update_flags,
793  const Quadrature<dim> &q,
794  const unsigned int n_original_q_points)
795 {
796  initialize(update_flags, q, n_original_q_points);
797 
798  if (dim > 1 && tensor_product_quadrature)
799  {
800  constexpr unsigned int facedim = dim - 1;
801  const FE_DGQ<1> fe(polynomial_degree);
802  shape_info.reinit(q.get_tensor_basis()[0], fe);
804  FETools::lexicographic_to_hierarchic_numbering<facedim>(
806  shape_info.n_q_points = n_original_q_points;
809  }
810 
811  if (dim > 1)
812  {
813  if (this->update_each &
816  {
817  aux.resize(dim - 1,
818  std::vector<Tensor<1, spacedim>>(n_original_q_points));
819 
820  // Compute tangentials to the unit cell.
821  for (const unsigned int i : GeometryInfo<dim>::face_indices())
822  {
823  unit_tangentials[i].resize(n_original_q_points);
824  std::fill(unit_tangentials[i].begin(),
825  unit_tangentials[i].end(),
827  if (dim > 2)
828  {
830  .resize(n_original_q_points);
831  std::fill(
833  .begin(),
835  .end(),
837  }
838  }
839  }
840  }
841 }
842 
843 
844 
845 template <>
846 void
848  const std::vector<Point<1>> &unit_points)
849 {
850  // if the polynomial degree is one, then we can simplify code a bit
851  // by using hard-coded shape functions.
852  if (polynomial_degree == 1)
853  internal::MappingQ1::compute_shape_function_values_hardcode(
854  n_shape_functions, unit_points, *this);
855  else
856  {
857  // otherwise ask an object that describes the polynomial space
858  internal::MappingQ1::compute_shape_function_values_general<1, 1>(
859  n_shape_functions, unit_points, *this);
860  }
861 }
862 
863 template <>
864 void
866  const std::vector<Point<2>> &unit_points)
867 {
868  // if the polynomial degree is one, then we can simplify code a bit
869  // by using hard-coded shape functions.
870  if (polynomial_degree == 1)
871  internal::MappingQ1::compute_shape_function_values_hardcode(
872  n_shape_functions, unit_points, *this);
873  else
874  {
875  // otherwise ask an object that describes the polynomial space
876  internal::MappingQ1::compute_shape_function_values_general<2, 2>(
877  n_shape_functions, unit_points, *this);
878  }
879 }
880 
881 template <>
882 void
884  const std::vector<Point<3>> &unit_points)
885 {
886  // if the polynomial degree is one, then we can simplify code a bit
887  // by using hard-coded shape functions.
888  if (polynomial_degree == 1)
889  internal::MappingQ1::compute_shape_function_values_hardcode(
890  n_shape_functions, unit_points, *this);
891  else
892  {
893  // otherwise ask an object that describes the polynomial space
894  internal::MappingQ1::compute_shape_function_values_general<3, 3>(
895  n_shape_functions, unit_points, *this);
896  }
897 }
898 
899 template <int dim, int spacedim>
900 void
902  const std::vector<Point<dim>> &unit_points)
903 {
904  // for non-matching combinations of dim and spacedim, just run the general
905  // case
906  internal::MappingQ1::compute_shape_function_values_general<dim, spacedim>(
907  n_shape_functions, unit_points, *this);
908 }
909 
910 
911 namespace internal
912 {
913  namespace MappingQGenericImplementation
914  {
915  namespace
916  {
925  compute_support_point_weights_on_quad(
926  const unsigned int polynomial_degree)
927  {
928  ::Table<2, double> loqvs;
929 
930  // we are asked to compute weights for interior support points, but
931  // there are no interior points if degree==1
932  if (polynomial_degree == 1)
933  return loqvs;
934 
935  const unsigned int M = polynomial_degree - 1;
936  const unsigned int n_inner_2d = M * M;
937  const unsigned int n_outer_2d = 4 + 4 * M;
938 
939  // set the weights of transfinite interpolation
940  loqvs.reinit(n_inner_2d, n_outer_2d);
941  QGaussLobatto<2> gl(polynomial_degree + 1);
942  for (unsigned int i = 0; i < M; ++i)
943  for (unsigned int j = 0; j < M; ++j)
944  {
945  const Point<2> p =
946  gl.point((i + 1) * (polynomial_degree + 1) + (j + 1));
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];
955  }
956 
957  // the sum of weights of the points at the outer rim should be one.
958  // check this
959  for (unsigned int unit_point = 0; unit_point < n_inner_2d; ++unit_point)
960  Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
961  loqvs[unit_point].end(),
962  0.) -
963  1) < 1e-13 * polynomial_degree,
964  ExcInternalError());
965 
966  return loqvs;
967  }
968 
969 
970 
978  compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
979  {
980  ::Table<2, double> lohvs;
981 
982  // we are asked to compute weights for interior support points, but
983  // there are no interior points if degree==1
984  if (polynomial_degree == 1)
985  return lohvs;
986 
987  const unsigned int M = polynomial_degree - 1;
988 
989  const unsigned int n_inner = Utilities::fixed_power<3>(M);
990  const unsigned int n_outer = 8 + 12 * M + 6 * M * M;
991 
992  // set the weights of transfinite interpolation
993  lohvs.reinit(n_inner, n_outer);
994  QGaussLobatto<3> gl(polynomial_degree + 1);
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)
998  {
999  const Point<3> p = gl.point((i + 1) * (M + 2) * (M + 2) +
1000  (j + 1) * (M + 2) + (k + 1));
1001  const unsigned int index_table = i * M * M + j * M + k;
1002 
1003  // vertices
1004  for (unsigned int v = 0; v < 8; ++v)
1005  lohvs(index_table, v) =
1007 
1008  // lines
1009  {
1010  constexpr std::array<unsigned int, 4> line_coordinates_y(
1011  {{0, 1, 4, 5}});
1012  const Point<2> py(p[0], p[2]);
1013  for (unsigned int l = 0; l < 4; ++l)
1014  lohvs(index_table, 8 + line_coordinates_y[l] * M + j) =
1016  }
1017 
1018  {
1019  constexpr std::array<unsigned int, 4> line_coordinates_x(
1020  {{2, 3, 6, 7}});
1021  const Point<2> px(p[1], p[2]);
1022  for (unsigned int l = 0; l < 4; ++l)
1023  lohvs(index_table, 8 + line_coordinates_x[l] * M + k) =
1025  }
1026 
1027  {
1028  constexpr std::array<unsigned int, 4> line_coordinates_z(
1029  {{8, 9, 10, 11}});
1030  const Point<2> pz(p[0], p[1]);
1031  for (unsigned int l = 0; l < 4; ++l)
1032  lohvs(index_table, 8 + line_coordinates_z[l] * M + i) =
1034  }
1035 
1036  // quads
1037  lohvs(index_table, 8 + 12 * M + 0 * M * M + i * M + j) =
1038  1. - p[0];
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) =
1041  1. - p[1];
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) =
1044  1. - p[2];
1045  lohvs(index_table, 8 + 12 * M + 5 * M * M + j * M + k) = p[2];
1046  }
1047 
1048  // the sum of weights of the points at the outer rim should be one.
1049  // check this
1050  for (unsigned int unit_point = 0; unit_point < n_inner; ++unit_point)
1051  Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
1052  lohvs[unit_point].end(),
1053  0.) -
1054  1) < 1e-13 * polynomial_degree,
1055  ExcInternalError());
1056 
1057  return lohvs;
1058  }
1059 
1060 
1061 
1066  std::vector<::Table<2, double>>
1067  compute_support_point_weights_perimeter_to_interior(
1068  const unsigned int polynomial_degree,
1069  const unsigned int dim)
1070  {
1071  Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
1072  std::vector<::Table<2, double>> output(dim);
1073  if (polynomial_degree <= 1)
1074  return output;
1075 
1076  // fill the 1D interior weights
1077  QGaussLobatto<1> quadrature(polynomial_degree + 1);
1078  output[0].reinit(polynomial_degree - 1,
1080  for (unsigned int q = 0; q < polynomial_degree - 1; ++q)
1081  for (const unsigned int i : GeometryInfo<1>::vertex_indices())
1082  output[0](q, i) =
1084  i);
1085 
1086  if (dim > 1)
1087  output[1] = compute_support_point_weights_on_quad(polynomial_degree);
1088 
1089  if (dim > 2)
1090  output[2] = compute_support_point_weights_on_hex(polynomial_degree);
1091 
1092  return output;
1093  }
1094 
1098  template <int dim>
1100  compute_support_point_weights_cell(const unsigned int polynomial_degree)
1101  {
1102  Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
1103  if (polynomial_degree <= 1)
1104  return ::Table<2, double>();
1105 
1106  QGaussLobatto<dim> quadrature(polynomial_degree + 1);
1107  const std::vector<unsigned int> h2l =
1108  FETools::hierarchic_to_lexicographic_numbering<dim>(
1110 
1111  ::Table<2, double> output(quadrature.size() -
1114  for (unsigned int q = 0; q < output.size(0); ++q)
1115  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1117  quadrature.point(h2l[q + GeometryInfo<dim>::vertices_per_cell]),
1118  i);
1119 
1120  return output;
1121  }
1122 
1123 
1124 
1132  template <int dim, int spacedim>
1134  compute_mapped_location_of_point(
1135  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1136  &data)
1137  {
1138  AssertDimension(data.shape_values.size(),
1139  data.mapping_support_points.size());
1140 
1141  // use now the InternalData to compute the point in real space.
1142  Point<spacedim> p_real;
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);
1145 
1146  return p_real;
1147  }
1148 
1149 
1150 
1154  template <int dim>
1155  Point<dim>
1156  do_transform_real_to_unit_cell_internal(
1157  const typename ::Triangulation<dim, dim>::cell_iterator &cell,
1158  const Point<dim> & p,
1159  const Point<dim> &initial_p_unit,
1160  typename ::MappingQGeneric<dim, dim>::InternalData &mdata)
1161  {
1162  const unsigned int spacedim = dim;
1163 
1164  const unsigned int n_shapes = mdata.shape_values.size();
1165  (void)n_shapes;
1166  Assert(n_shapes != 0, ExcInternalError());
1167  AssertDimension(mdata.shape_derivatives.size(), n_shapes);
1168 
1169  std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
1170  AssertDimension(points.size(), n_shapes);
1171 
1172 
1173  // Newton iteration to solve
1174  // f(x)=p(x)-p=0
1175  // where we are looking for 'x' and p(x) is the forward transformation
1176  // from unit to real cell. We solve this using a Newton iteration
1177  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
1178  // The start value is set to be the linear approximation to the cell
1179 
1180  // The shape values and derivatives of the mapping at this point are
1181  // previously computed.
1182 
1183  Point<dim> p_unit = initial_p_unit;
1184 
1185  mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
1186 
1187  Point<spacedim> p_real =
1188  compute_mapped_location_of_point<dim, spacedim>(mdata);
1189  Tensor<1, spacedim> f = p_real - p;
1190 
1191  // early out if we already have our point
1192  if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter())
1193  return p_unit;
1194 
1195  // we need to compare the position of the computed p(x) against the
1196  // given point 'p'. We will terminate the iteration and return 'x' if
1197  // they are less than eps apart. The question is how to choose eps --
1198  // or, put maybe more generally: in which norm we want these 'p' and
1199  // 'p(x)' to be eps apart.
1200  //
1201  // the question is difficult since we may have to deal with very
1202  // elongated cells where we may achieve 1e-12*h for the distance of
1203  // these two points in the 'long' direction, but achieving this
1204  // tolerance in the 'short' direction of the cell may not be possible
1205  //
1206  // what we do instead is then to terminate iterations if
1207  // \| p(x) - p \|_A < eps
1208  // where the A-norm is somehow induced by the transformation of the
1209  // cell. in particular, we want to measure distances relative to the
1210  // sizes of the cell in its principal directions.
1211  //
1212  // to define what exactly A should be, note that to first order we have
1213  // the following (assuming that x* is the solution of the problem, i.e.,
1214  // p(x*)=p):
1215  // p(x) - p = p(x) - p(x*)
1216  // = -grad p(x) * (x*-x) + higher order terms
1217  // This suggest to measure with a norm that corresponds to
1218  // A = {[grad p(x]^T [grad p(x)]}^{-1}
1219  // because then
1220  // \| p(x) - p \|_A \approx \| x - x* \|
1221  // Consequently, we will try to enforce that
1222  // \| p(x) - p \|_A = \| f \| <= eps
1223  //
1224  // Note that using this norm is a bit dangerous since the norm changes
1225  // in every iteration (A isn't fixed by depends on xk). However, if the
1226  // cell is not too deformed (it may be stretched, but not twisted) then
1227  // the mapping is almost linear and A is indeed constant or nearly so.
1228  const double eps = 1.e-11;
1229  const unsigned int newton_iteration_limit = 20;
1230 
1231  unsigned int newton_iteration = 0;
1232  double last_f_weighted_norm;
1233  do
1234  {
1235 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
1236  std::cout << "Newton iteration " << newton_iteration << std::endl;
1237 #endif
1238 
1239  // f'(x)
1241  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1242  {
1243  const Tensor<1, dim> & grad_transform = mdata.derivative(0, k);
1244  const Point<spacedim> &point = points[k];
1245 
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];
1249  }
1250 
1251  // Solve [f'(x)]d=f(x)
1252  AssertThrow(
1253  determinant(df) > 0,
1255  Tensor<2, spacedim> df_inverse = invert(df);
1256  const Tensor<1, spacedim> delta =
1257  df_inverse * static_cast<const Tensor<1, spacedim> &>(f);
1258 
1259 #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
1260  std::cout << " delta=" << delta << std::endl;
1261 #endif
1262 
1263  // do a line search
1264  double step_length = 1;
1265  do
1266  {
1267  // update of p_unit. The spacedim-th component of transformed
1268  // point is simply ignored in codimension one case. When this
1269  // component is not zero, then we are projecting the point to
1270  // the surface or curve identified by the cell.
1271  Point<dim> p_unit_trial = p_unit;
1272  for (unsigned int i = 0; i < dim; ++i)
1273  p_unit_trial[i] -= step_length * delta[i];
1274 
1275  // shape values and derivatives
1276  // at new p_unit point
1277  mdata.compute_shape_function_values(
1278  std::vector<Point<dim>>(1, p_unit_trial));
1279 
1280  // f(x)
1281  Point<spacedim> p_real_trial =
1282  internal::MappingQGenericImplementation::
1283  compute_mapped_location_of_point<dim, spacedim>(mdata);
1284  const Tensor<1, spacedim> f_trial = p_real_trial - p;
1285 
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
1290  << " ||f*||_A ="
1291  << (df_inverse * f_trial).norm() << std::endl;
1292 #endif
1293 
1294  // see if we are making progress with the current step length
1295  // and if not, reduce it by a factor of two and try again
1296  //
1297  // strictly speaking, we should probably use the same norm as we
1298  // use for the outer algorithm. in practice, line search is just
1299  // a crutch to find a "reasonable" step length, and so using the
1300  // l2 norm is probably just fine
1301  if (f_trial.norm() < f.norm())
1302  {
1303  p_real = p_real_trial;
1304  p_unit = p_unit_trial;
1305  f = f_trial;
1306  break;
1307  }
1308  else if (step_length > 0.05)
1309  step_length /= 2;
1310  else
1311  AssertThrow(
1312  false,
1313  (typename Mapping<dim,
1314  spacedim>::ExcTransformationFailed()));
1315  }
1316  while (true);
1317 
1318  ++newton_iteration;
1319  if (newton_iteration > newton_iteration_limit)
1320  AssertThrow(
1321  false,
1323  last_f_weighted_norm = (df_inverse * f).norm();
1324  }
1325  while (last_f_weighted_norm > eps);
1326 
1327  return p_unit;
1328  }
1329 
1330 
1331 
1335  template <int dim>
1336  Point<dim>
1337  do_transform_real_to_unit_cell_internal_codim1(
1338  const typename ::Triangulation<dim, dim + 1>::cell_iterator &cell,
1339  const Point<dim + 1> & p,
1340  const Point<dim> &initial_p_unit,
1341  typename ::MappingQGeneric<dim, dim + 1>::InternalData &mdata)
1342  {
1343  const unsigned int spacedim = dim + 1;
1344 
1345  const unsigned int n_shapes = mdata.shape_values.size();
1346  (void)n_shapes;
1347  Assert(n_shapes != 0, ExcInternalError());
1348  Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError());
1349  Assert(mdata.shape_second_derivatives.size() == n_shapes,
1350  ExcInternalError());
1351 
1352  std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
1353  Assert(points.size() == n_shapes, ExcInternalError());
1354 
1355  Point<spacedim> p_minus_F;
1356 
1357  Tensor<1, spacedim> DF[dim];
1358  Tensor<1, spacedim> D2F[dim][dim];
1359 
1360  Point<dim> p_unit = initial_p_unit;
1361  Point<dim> f;
1362  Tensor<2, dim> df;
1363 
1364  // Evaluate first and second derivatives
1365  mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
1366 
1367  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1368  {
1369  const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
1370  const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
1371  const Point<spacedim> &point_k = points[k];
1372 
1373  for (unsigned int j = 0; j < dim; ++j)
1374  {
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;
1378  }
1379  }
1380 
1381  p_minus_F = p;
1382  p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
1383 
1384 
1385  for (unsigned int j = 0; j < dim; ++j)
1386  f[j] = DF[j] * p_minus_F;
1387 
1388  for (unsigned int j = 0; j < dim; ++j)
1389  {
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;
1393  }
1394 
1395 
1396  const double eps = 1.e-12 * cell->diameter();
1397  const unsigned int loop_limit = 10;
1398 
1399  unsigned int loop = 0;
1400 
1401  while (f.norm() > eps && loop++ < loop_limit)
1402  {
1403  // Solve [df(x)]d=f(x)
1404  const Tensor<1, dim> d =
1405  invert(df) * static_cast<const Tensor<1, dim> &>(f);
1406  p_unit -= d;
1407 
1408  for (unsigned int j = 0; j < dim; ++j)
1409  {
1410  DF[j].clear();
1411  for (unsigned int l = 0; l < dim; ++l)
1412  D2F[j][l].clear();
1413  }
1414 
1415  mdata.compute_shape_function_values(
1416  std::vector<Point<dim>>(1, p_unit));
1417 
1418  for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
1419  {
1420  const Tensor<1, dim> &grad_phi_k = mdata.derivative(0, k);
1421  const Tensor<2, dim> &hessian_k = mdata.second_derivative(0, k);
1422  const Point<spacedim> &point_k = points[k];
1423 
1424  for (unsigned int j = 0; j < dim; ++j)
1425  {
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;
1429  }
1430  }
1431 
1432  // TODO: implement a line search here in much the same way as for
1433  // the corresponding function above that does so for dim==spacedim
1434  p_minus_F = p;
1435  p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
1436 
1437  for (unsigned int j = 0; j < dim; ++j)
1438  {
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;
1442  }
1443  }
1444 
1445 
1446  // Here we check that in the last execution of while the first
1447  // condition was already wrong, meaning the residual was below
1448  // eps. Only if the first condition failed, loop will have been
1449  // increased and tested, and thus have reached the limit.
1450  AssertThrow(
1451  loop < loop_limit,
1453 
1454  return p_unit;
1455  }
1456 
1462  template <int dim, int spacedim>
1463  void
1464  maybe_update_q_points_Jacobians_and_grads_tensor(
1465  const CellSimilarity::Similarity cell_similarity,
1466  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1467  & data,
1468  std::vector<Point<spacedim>> & quadrature_points,
1469  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1470  {
1471  const UpdateFlags update_flags = data.update_each;
1472 
1473  const unsigned int n_shape_values = data.n_shape_functions;
1474  const unsigned int n_q_points = data.shape_info.n_q_points;
1475  constexpr unsigned int n_lanes = VectorizedArray<double>::size();
1476  constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
1477  constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
1478 
1479  const bool evaluate_values = update_flags & update_quadrature_points;
1480  const bool evaluate_gradients =
1481  (cell_similarity != CellSimilarity::translation) &&
1482  (update_flags & update_contravariant_transformation);
1483  const bool evaluate_hessians =
1484  (cell_similarity != CellSimilarity::translation) &&
1485  (update_flags & update_jacobian_grads);
1486 
1487  Assert(!evaluate_values || n_q_points > 0, ExcInternalError());
1488  Assert(!evaluate_values || n_q_points == quadrature_points.size(),
1489  ExcDimensionMismatch(n_q_points, quadrature_points.size()));
1490  Assert(!evaluate_gradients || data.n_shape_functions > 0,
1491  ExcInternalError());
1492  Assert(!evaluate_gradients || n_q_points == data.contravariant.size(),
1493  ExcDimensionMismatch(n_q_points, data.contravariant.size()));
1494  Assert(!evaluate_hessians || n_q_points == jacobian_grads.size(),
1495  ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
1496 
1497  // shortcut in case we have an identity interpolation and only request
1498  // the quadrature points
1499  if (evaluate_values && !evaluate_gradients & !evaluate_hessians &&
1500  data.shape_info.element_type ==
1502  {
1503  for (unsigned int q = 0; q < n_q_points; ++q)
1504  quadrature_points[q] =
1505  data.mapping_support_points[data.shape_info
1506  .lexicographic_numbering[q]];
1507  return;
1508  }
1509 
1510  // prepare arrays
1511  if (evaluate_values || evaluate_gradients || evaluate_hessians)
1512  {
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));
1517 
1518  if (evaluate_hessians)
1519  data.hessians_quad.resize(n_comp * n_q_points * n_hessians);
1520 
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)
1525  {
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] =
1529  data
1530  .mapping_support_points[renumber_to_lexicographic[i]][d];
1531  }
1532 
1533  // do the actual tensorized evaluation
1534  SelectEvaluator<dim, -1, 0, n_comp, VectorizedArray<double>>::
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(),
1541  evaluate_values,
1542  evaluate_gradients,
1543  evaluate_hessians);
1544  }
1545 
1546  // do the postprocessing
1547  if (evaluate_values)
1548  {
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;
1554  ++in_comp)
1555  quadrature_points[i][out_comp * n_lanes + in_comp] =
1556  data.values_quad[out_comp * n_q_points + i][in_comp];
1557  }
1558 
1559  if (evaluate_gradients)
1560  {
1561  std::fill(data.contravariant.begin(),
1562  data.contravariant.end(),
1564  // We need to reinterpret the data after evaluate has been applied.
1565  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1566  for (unsigned int point = 0; point < n_q_points; ++point)
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;
1571  ++in_comp)
1572  {
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) *
1579  dim +
1580  j][in_comp];
1581  }
1582  }
1583  if (update_flags & update_covariant_transformation)
1584  if (cell_similarity != CellSimilarity::translation)
1585  for (unsigned int point = 0; point < n_q_points; ++point)
1586  data.covariant[point] =
1587  (data.contravariant[point]).covariant_form();
1588 
1589  if (update_flags & update_volume_elements)
1590  if (cell_similarity != CellSimilarity::translation)
1591  for (unsigned int point = 0; point < n_q_points; ++point)
1592  data.volume_elements[point] =
1593  data.contravariant[point].determinant();
1594 
1595  if (evaluate_hessians)
1596  {
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}};
1600 
1601  // We need to reinterpret the data after evaluate has been applied.
1602  for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
1603  for (unsigned int point = 0; point < n_q_points; ++point)
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;
1608  ++in_comp)
1609  {
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) *
1622  n_hessians +
1623  j][in_comp];
1624  jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1625  [new_hessian_comp_i][new_hessian_comp_j] =
1626  value;
1627  jacobian_grads[new_point][out_comp * n_lanes + in_comp]
1628  [new_hessian_comp_j][new_hessian_comp_i] =
1629  value;
1630  }
1631  }
1632  }
1633 
1634 
1641  template <int dim, int spacedim>
1642  void
1643  maybe_compute_q_points(
1644  const typename QProjector<dim>::DataSetDescriptor data_set,
1645  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1646  & data,
1647  std::vector<Point<spacedim>> &quadrature_points)
1648  {
1649  const UpdateFlags update_flags = data.update_each;
1650 
1651  if (update_flags & update_quadrature_points)
1652  for (unsigned int point = 0; point < quadrature_points.size();
1653  ++point)
1654  {
1655  const double * shape = &data.shape(point + data_set, 0);
1656  Point<spacedim> result =
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];
1661  quadrature_points[point] = result;
1662  }
1663  }
1664 
1665 
1666 
1675  template <int dim, int spacedim>
1676  void
1677  maybe_update_Jacobians(
1678  const CellSimilarity::Similarity cell_similarity,
1679  const typename ::QProjector<dim>::DataSetDescriptor data_set,
1680  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1681  &data)
1682  {
1683  const UpdateFlags update_flags = data.update_each;
1684 
1685  if (update_flags & update_contravariant_transformation)
1686  // if the current cell is just a
1687  // translation of the previous one, no
1688  // need to recompute jacobians...
1689  if (cell_similarity != CellSimilarity::translation)
1690  {
1691  const unsigned int n_q_points = data.contravariant.size();
1692 
1693  std::fill(data.contravariant.begin(),
1694  data.contravariant.end(),
1696 
1697  Assert(data.n_shape_functions > 0, ExcInternalError());
1698 
1699  const Tensor<1, spacedim> *supp_pts =
1700  data.mapping_support_points.data();
1701 
1702  for (unsigned int point = 0; point < n_q_points; ++point)
1703  {
1704  const Tensor<1, dim> *data_derv =
1705  &data.derivative(point + data_set, 0);
1706 
1707  double result[spacedim][dim];
1708 
1709  // peel away part of sum to avoid zeroing the
1710  // entries and adding for the first time
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];
1718 
1719  // write result into contravariant data. for
1720  // j=dim in the case dim<spacedim, there will
1721  // never be any nonzero data that arrives in
1722  // here, so it is ok anyway because it was
1723  // initialized to zero at the initialization
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];
1727  }
1728  }
1729 
1730  if (update_flags & update_covariant_transformation)
1731  if (cell_similarity != CellSimilarity::translation)
1732  {
1733  const unsigned int n_q_points = data.contravariant.size();
1734  for (unsigned int point = 0; point < n_q_points; ++point)
1735  {
1736  data.covariant[point] =
1737  (data.contravariant[point]).covariant_form();
1738  }
1739  }
1740 
1741  if (update_flags & update_volume_elements)
1742  if (cell_similarity != CellSimilarity::translation)
1743  {
1744  const unsigned int n_q_points = data.contravariant.size();
1745  for (unsigned int point = 0; point < n_q_points; ++point)
1746  data.volume_elements[point] =
1747  data.contravariant[point].determinant();
1748  }
1749  }
1750 
1757  template <int dim, int spacedim>
1758  void
1759  maybe_update_jacobian_grads(
1760  const CellSimilarity::Similarity cell_similarity,
1761  const typename QProjector<dim>::DataSetDescriptor data_set,
1762  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1763  & data,
1764  std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
1765  {
1766  const UpdateFlags update_flags = data.update_each;
1767  if (update_flags & update_jacobian_grads)
1768  {
1769  const unsigned int n_q_points = jacobian_grads.size();
1770 
1771  if (cell_similarity != CellSimilarity::translation)
1772  for (unsigned int point = 0; point < n_q_points; ++point)
1773  {
1774  const Tensor<2, dim> *second =
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)
1780  result[i][j][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)
1786  result[i][j][l] +=
1787  (second[k][j][l] *
1788  data.mapping_support_points[k][i]);
1789 
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];
1794  }
1795  }
1796  }
1797 
1804  template <int dim, int spacedim>
1805  void
1806  maybe_update_jacobian_pushed_forward_grads(
1807  const CellSimilarity::Similarity cell_similarity,
1808  const typename QProjector<dim>::DataSetDescriptor data_set,
1809  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1810  & data,
1811  std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
1812  {
1813  const UpdateFlags update_flags = data.update_each;
1814  if (update_flags & update_jacobian_pushed_forward_grads)
1815  {
1816  const unsigned int n_q_points =
1817  jacobian_pushed_forward_grads.size();
1818 
1819  if (cell_similarity != CellSimilarity::translation)
1820  {
1821  double tmp[spacedim][spacedim][spacedim];
1822  for (unsigned int point = 0; point < n_q_points; ++point)
1823  {
1824  const Tensor<2, dim> *second =
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)
1836  result[i][j][l] +=
1837  (second[k][j][l] *
1838  data.mapping_support_points[k][i]);
1839 
1840  // first push forward the j-components
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)
1844  {
1845  tmp[i][j][l] =
1846  result[i][0][l] * data.covariant[point][j][0];
1847  for (unsigned int jr = 1; jr < dim; ++jr)
1848  {
1849  tmp[i][j][l] += result[i][jr][l] *
1850  data.covariant[point][j][jr];
1851  }
1852  }
1853 
1854  // now, pushing forward the l-components
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)
1858  {
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)
1862  {
1863  jacobian_pushed_forward_grads[point][i][j][l] +=
1864  tmp[i][j][lr] * data.covariant[point][l][lr];
1865  }
1866  }
1867  }
1868  }
1869  }
1870  }
1871 
1878  template <int dim, int spacedim>
1879  void
1880  maybe_update_jacobian_2nd_derivatives(
1881  const CellSimilarity::Similarity cell_similarity,
1882  const typename QProjector<dim>::DataSetDescriptor data_set,
1883  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1884  & data,
1885  std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
1886  {
1887  const UpdateFlags update_flags = data.update_each;
1888  if (update_flags & update_jacobian_2nd_derivatives)
1889  {
1890  const unsigned int n_q_points = jacobian_2nd_derivatives.size();
1891 
1892  if (cell_similarity != CellSimilarity::translation)
1893  {
1894  for (unsigned int point = 0; point < n_q_points; ++point)
1895  {
1896  const Tensor<3, dim> *third =
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]);
1914 
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] =
1920  result[i][j][l][m];
1921  }
1922  }
1923  }
1924  }
1925 
1933  template <int dim, int spacedim>
1934  void
1935  maybe_update_jacobian_pushed_forward_2nd_derivatives(
1936  const CellSimilarity::Similarity cell_similarity,
1937  const typename QProjector<dim>::DataSetDescriptor data_set,
1938  const typename ::MappingQGeneric<dim, spacedim>::InternalData
1939  &data,
1940  std::vector<Tensor<4, spacedim>>
1941  &jacobian_pushed_forward_2nd_derivatives)
1942  {
1943  const UpdateFlags update_flags = data.update_each;
1945  {
1946  const unsigned int n_q_points =
1947  jacobian_pushed_forward_2nd_derivatives.size();
1948 
1949  if (cell_similarity != CellSimilarity::translation)
1950  {
1951  double tmp[spacedim][spacedim][spacedim][spacedim];
1952  for (unsigned int point = 0; point < n_q_points; ++point)
1953  {
1954  const Tensor<3, dim> *third =
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]);
1972 
1973  // push forward the j-coordinate
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)
1978  {
1979  jacobian_pushed_forward_2nd_derivatives
1980  [point][i][j][l][m] =
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]
1985  [i][j][l]
1986  [m] +=
1987  result[i][jr][l][m] *
1988  data.covariant[point][j][jr];
1989  }
1990 
1991  // push forward the l-coordinate
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)
1996  {
1997  tmp[i][j][l][m] =
1998  jacobian_pushed_forward_2nd_derivatives[point]
1999  [i][j][0]
2000  [m] *
2001  data.covariant[point][l][0];
2002  for (unsigned int lr = 1; lr < dim; ++lr)
2003  tmp[i][j][l][m] +=
2004  jacobian_pushed_forward_2nd_derivatives
2005  [point][i][j][lr][m] *
2006  data.covariant[point][l][lr];
2007  }
2008 
2009  // push forward the m-coordinate
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)
2014  {
2015  jacobian_pushed_forward_2nd_derivatives
2016  [point][i][j][l][m] =
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]
2020  [i][j][l]
2021  [m] +=
2022  tmp[i][j][l][mr] *
2023  data.covariant[point][m][mr];
2024  }
2025  }
2026  }
2027  }
2028  }
2029 
2036  template <int dim, int spacedim>
2037  void
2038  maybe_update_jacobian_3rd_derivatives(
2039  const CellSimilarity::Similarity cell_similarity,
2040  const typename QProjector<dim>::DataSetDescriptor data_set,
2041  const typename ::MappingQGeneric<dim, spacedim>::InternalData
2042  & data,
2043  std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
2044  {
2045  const UpdateFlags update_flags = data.update_each;
2046  if (update_flags & update_jacobian_3rd_derivatives)
2047  {
2048  const unsigned int n_q_points = jacobian_3rd_derivatives.size();
2049 
2050  if (cell_similarity != CellSimilarity::translation)
2051  {
2052  for (unsigned int point = 0; point < n_q_points; ++point)
2053  {
2054  const Tensor<4, dim> *fourth =
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]);
2074 
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];
2082  }
2083  }
2084  }
2085  }
2086 
2094  template <int dim, int spacedim>
2095  void
2096  maybe_update_jacobian_pushed_forward_3rd_derivatives(
2097  const CellSimilarity::Similarity cell_similarity,
2098  const typename QProjector<dim>::DataSetDescriptor data_set,
2099  const typename ::MappingQGeneric<dim, spacedim>::InternalData
2100  &data,
2101  std::vector<Tensor<5, spacedim>>
2102  &jacobian_pushed_forward_3rd_derivatives)
2103  {
2104  const UpdateFlags update_flags = data.update_each;
2106  {
2107  const unsigned int n_q_points =
2108  jacobian_pushed_forward_3rd_derivatives.size();
2109 
2110  if (cell_similarity != CellSimilarity::translation)
2111  {
2112  double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
2113  for (unsigned int point = 0; point < n_q_points; ++point)
2114  {
2115  const Tensor<4, dim> *fourth =
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]);
2135 
2136  // push-forward the j-coordinate
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)
2142  {
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];
2150  }
2151 
2152  // push-forward the l-coordinate
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)
2158  {
2159  jacobian_pushed_forward_3rd_derivatives
2160  [point][i][j][l][m][n] =
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];
2168  }
2169 
2170  // push-forward the m-coordinate
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)
2176  {
2177  tmp[i][j][l][m][n] =
2178  jacobian_pushed_forward_3rd_derivatives
2179  [point][i][j][l][0][n] *
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];
2186  }
2187 
2188  // push-forward the n-coordinate
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)
2194  {
2195  jacobian_pushed_forward_3rd_derivatives
2196  [point][i][j][l][m][n] =
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];
2204  }
2205  }
2206  }
2207  }
2208  }
2209  } // namespace
2210  } // namespace MappingQGenericImplementation
2211 } // namespace internal
2212 
2213 
2214 
2215 template <int dim, int spacedim>
2217  : polynomial_degree(p)
2220  internal::MappingQGenericImplementation::
2221  compute_support_point_weights_perimeter_to_interior(
2222  this->polynomial_degree,
2223  dim))
2225  internal::MappingQGenericImplementation::
2226  compute_support_point_weights_cell<dim>(this->polynomial_degree))
2227 {
2228  Assert(p >= 1,
2229  ExcMessage("It only makes sense to create polynomial mappings "
2230  "with a polynomial degree greater or equal to one."));
2231 }
2232 
2233 
2234 
2235 template <int dim, int spacedim>
2237  const MappingQGeneric<dim, spacedim> &mapping)
2243 {}
2244 
2245 
2246 
2247 template <int dim, int spacedim>
2248 std::unique_ptr<Mapping<dim, spacedim>>
2250 {
2251  return std_cxx14::make_unique<MappingQGeneric<dim, spacedim>>(*this);
2252 }
2253 
2254 
2255 
2256 template <int dim, int spacedim>
2257 unsigned int
2259 {
2260  return polynomial_degree;
2261 }
2262 
2263 
2264 
2265 template <int dim, int spacedim>
2268  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2269  const Point<dim> & p) const
2270 {
2271  // set up the polynomial space
2272  const TensorProductPolynomials<dim> tensor_pols(
2275  Assert(tensor_pols.n() == Utilities::fixed_power<dim>(polynomial_degree + 1),
2276  ExcInternalError());
2277 
2278  // then also construct the mapping from lexicographic to the Qp shape function
2279  // numbering
2280  const std::vector<unsigned int> renumber =
2281  FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
2282 
2283  const std::vector<Point<spacedim>> support_points =
2284  this->compute_mapping_support_points(cell);
2285 
2286  Point<spacedim> mapped_point;
2287  for (unsigned int i = 0; i < tensor_pols.n(); ++i)
2288  mapped_point +=
2289  support_points[i] * tensor_pols.compute_value(renumber[i], p);
2290 
2291  return mapped_point;
2292 }
2293 
2294 
2295 // In the code below, GCC tries to instantiate MappingQGeneric<3,4> when
2296 // seeing which of the overloaded versions of
2297 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
2298 // error messages and, generally, nothing very good. Avoid this by ensuring
2299 // that this class exists, but does not have an inner InternalData
2300 // type, thereby ruling out the codim-1 version of the function
2301 // below when doing overload resolution.
2302 template <>
2303 class MappingQGeneric<3, 4>
2304 {};
2305 
2306 
2307 
2308 // visual studio freaks out when trying to determine if
2309 // do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good
2310 // candidate. So instead of letting the compiler pick the correct overload, we
2311 // use template specialization to make sure we pick up the right function to
2312 // call:
2313 
2314 template <int dim, int spacedim>
2315 Point<dim>
2318  const Point<spacedim> &,
2319  const Point<dim> &) const
2320 {
2321  // default implementation (should never be called)
2322  Assert(false, ExcInternalError());
2323  return Point<dim>();
2324 }
2325 
2326 template <>
2327 Point<1>
2330  const Point<1> & p,
2331  const Point<1> & initial_p_unit) const
2332 {
2333  const int dim = 1;
2334  const int spacedim = 1;
2335 
2336  const Quadrature<dim> point_quadrature(initial_p_unit);
2337 
2339  if (spacedim > dim)
2340  update_flags |= update_jacobian_grads;
2342  get_data(update_flags, point_quadrature));
2343 
2345 
2346  // dispatch to the various specializations for spacedim=dim,
2347  // spacedim=dim+1, etc
2348  return internal::MappingQGenericImplementation::
2349  do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata);
2350 }
2351 
2352 template <>
2353 Point<2>
2356  const Point<2> & p,
2357  const Point<2> & initial_p_unit) const
2358 {
2359  const int dim = 2;
2360  const int spacedim = 2;
2361 
2362  const Quadrature<dim> point_quadrature(initial_p_unit);
2363 
2365  if (spacedim > dim)
2366  update_flags |= update_jacobian_grads;
2368  get_data(update_flags, point_quadrature));
2369 
2371 
2372  // dispatch to the various specializations for spacedim=dim,
2373  // spacedim=dim+1, etc
2374  return internal::MappingQGenericImplementation::
2375  do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata);
2376 }
2377 
2378 template <>
2379 Point<3>
2382  const Point<3> & p,
2383  const Point<3> & initial_p_unit) const
2384 {
2385  const int dim = 3;
2386  const int spacedim = 3;
2387 
2388  const Quadrature<dim> point_quadrature(initial_p_unit);
2389 
2391  if (spacedim > dim)
2392  update_flags |= update_jacobian_grads;
2394  get_data(update_flags, point_quadrature));
2395 
2397 
2398  // dispatch to the various specializations for spacedim=dim,
2399  // spacedim=dim+1, etc
2400  return internal::MappingQGenericImplementation::
2401  do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata);
2402 }
2403 
2404 
2405 
2406 template <>
2407 Point<1>
2410  const Point<2> & p,
2411  const Point<1> & initial_p_unit) const
2412 {
2413  const int dim = 1;
2414  const int spacedim = 2;
2415 
2416  const Quadrature<dim> point_quadrature(initial_p_unit);
2417 
2419  if (spacedim > dim)
2420  update_flags |= update_jacobian_grads;
2422  get_data(update_flags, point_quadrature));
2423 
2425 
2426  // dispatch to the various specializations for spacedim=dim,
2427  // spacedim=dim+1, etc
2428  return internal::MappingQGenericImplementation::
2429  do_transform_real_to_unit_cell_internal_codim1<1>(cell,
2430  p,
2431  initial_p_unit,
2432  *mdata);
2433 }
2434 
2435 
2436 
2437 template <>
2438 Point<2>
2441  const Point<3> & p,
2442  const Point<2> & initial_p_unit) const
2443 {
2444  const int dim = 2;
2445  const int spacedim = 3;
2446 
2447  const Quadrature<dim> point_quadrature(initial_p_unit);
2448 
2450  if (spacedim > dim)
2451  update_flags |= update_jacobian_grads;
2453  get_data(update_flags, point_quadrature));
2454 
2456 
2457  // dispatch to the various specializations for spacedim=dim,
2458  // spacedim=dim+1, etc
2459  return internal::MappingQGenericImplementation::
2460  do_transform_real_to_unit_cell_internal_codim1<2>(cell,
2461  p,
2462  initial_p_unit,
2463  *mdata);
2464 }
2465 
2466 template <>
2467 Point<1>
2470  const Point<3> &,
2471  const Point<1> &) const
2472 {
2473  Assert(false, ExcNotImplemented());
2474  return {};
2475 }
2476 
2477 
2478 
2479 template <int dim, int spacedim>
2480 Point<dim>
2482  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2483  const Point<spacedim> & p) const
2484 {
2485  // Use an exact formula if one is available. this is only the case
2486  // for Q1 mappings in 1d, and in 2d if dim==spacedim
2487  if (this->preserves_vertex_locations() && (polynomial_degree == 1) &&
2488  ((dim == 1) || ((dim == 2) && (dim == spacedim))))
2489  {
2490  // The dimension-dependent algorithms are much faster (about 25-45x in
2491  // 2D) but fail most of the time when the given point (p) is not in the
2492  // cell. The dimension-independent Newton algorithm given below is
2493  // slower, but more robust (though it still sometimes fails). Therefore
2494  // this function implements the following strategy based on the
2495  // p's dimension:
2496  //
2497  // * In 1D this mapping is linear, so the mapping is always invertible
2498  // (and the exact formula is known) as long as the cell has non-zero
2499  // length.
2500  // * In 2D the exact (quadratic) formula is called first. If either the
2501  // exact formula does not succeed (negative discriminant in the
2502  // quadratic formula) or succeeds but finds a solution outside of the
2503  // unit cell, then the Newton solver is called. The rationale for the
2504  // second choice is that the exact formula may provide two different
2505  // answers when mapping a point outside of the real cell, but the
2506  // Newton solver (if it converges) will only return one answer.
2507  // Otherwise the exact formula successfully found a point in the unit
2508  // cell and that value is returned.
2509  // * In 3D there is no (known to the authors) exact formula, so the Newton
2510  // algorithm is used.
2511  const std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
2512  vertices = this->get_vertices(cell);
2513  try
2514  {
2515  switch (dim)
2516  {
2517  case 1:
2518  {
2519  // formula not subject to any issues in 1d
2520  if (spacedim == 1)
2521  return internal::MappingQ1::transform_real_to_unit_cell(
2522  vertices, p);
2523  else
2524  break;
2525  }
2526 
2527  case 2:
2528  {
2529  const Point<dim> point =
2530  internal::MappingQ1::transform_real_to_unit_cell(vertices,
2531  p);
2532 
2533  // formula not guaranteed to work for points outside of
2534  // the cell. only take the computed point if it lies
2535  // inside the reference cell
2536  const double eps = 1e-15;
2537  if (-eps <= point(1) && point(1) <= 1 + eps &&
2538  -eps <= point(0) && point(0) <= 1 + eps)
2539  {
2540  return point;
2541  }
2542  else
2543  break;
2544  }
2545 
2546  default:
2547  {
2548  // we should get here, based on the if-condition at the top
2549  Assert(false, ExcInternalError());
2550  }
2551  }
2552  }
2553  catch (
2555  {
2556  // simply fall through and continue on to the standard Newton code
2557  }
2558  }
2559  else
2560  {
2561  // we can't use an explicit formula,
2562  }
2563 
2564 
2565  // Find the initial value for the Newton iteration by a normal
2566  // projection to the least square plane determined by the vertices
2567  // of the cell
2568  Point<dim> initial_p_unit;
2569  if (this->preserves_vertex_locations())
2570  initial_p_unit = cell->real_to_unit_cell_affine_approximation(p);
2571  else
2572  {
2573  // for the MappingQEulerian type classes, we want to still call the cell
2574  // iterator's affine approximation. do so by creating a dummy
2575  // triangulation with just the first vertices.
2576  //
2577  // we do this by first getting all support points, then
2578  // throwing away all but the vertices, and finally calling
2579  // the same function as above
2580  std::vector<Point<spacedim>> a =
2581  this->compute_mapping_support_points(cell);
2583  std::vector<CellData<dim>> cells(1);
2584  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2585  cells[0].vertices[i] = i;
2587  tria.create_triangulation(a, cells, SubCellData());
2588  initial_p_unit =
2589  tria.begin_active()->real_to_unit_cell_affine_approximation(p);
2590  }
2591  // in 1d with spacedim > 1 the affine approximation is exact
2592  if (dim == 1 && polynomial_degree == 1)
2593  {
2594  return initial_p_unit;
2595  }
2596  else
2597  {
2598  // in case the function above should have given us something back that
2599  // lies outside the unit cell, then project it back into the reference
2600  // cell in hopes that this gives a better starting point to the
2601  // following iteration
2602  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2603 
2604  // perform the Newton iteration and return the result. note that this
2605  // statement may throw an exception, which we simply pass up to the
2606  // caller
2607  return this->transform_real_to_unit_cell_internal(cell,
2608  p,
2609  initial_p_unit);
2610  }
2611 }
2612 
2613 
2614 
2615 template <int dim, int spacedim>
2618  const UpdateFlags in) const
2619 {
2620  // add flags if the respective quantities are necessary to compute
2621  // what we need. note that some flags appear in both the conditions
2622  // and in subsequent set operations. this leads to some circular
2623  // logic. the only way to treat this is to iterate. since there are
2624  // 5 if-clauses in the loop, it will take at most 5 iterations to
2625  // converge. do them:
2626  UpdateFlags out = in;
2627  for (unsigned int i = 0; i < 5; ++i)
2628  {
2629  // The following is a little incorrect:
2630  // If not applied on a face,
2631  // update_boundary_forms does not
2632  // make sense. On the other hand,
2633  // it is necessary on a
2634  // face. Currently,
2635  // update_boundary_forms is simply
2636  // ignored for the interior of a
2637  // cell.
2639  out |= update_boundary_forms;
2640 
2645 
2646  if (out &
2651 
2652  // The contravariant transformation is used in the Piola
2653  // transformation, which requires the determinant of the Jacobi
2654  // matrix of the transformation. Because we have no way of
2655  // knowing here whether the finite element wants to use the
2656  // contravariant or the Piola transforms, we add the JxW values
2657  // to the list of flags to be updated for each cell.
2659  out |= update_volume_elements;
2660 
2661  // the same is true when computing normal vectors: they require
2662  // the determinant of the Jacobian
2663  if (out & update_normal_vectors)
2664  out |= update_volume_elements;
2665  }
2666 
2667  return out;
2668 }
2669 
2670 
2671 
2672 template <int dim, int spacedim>
2673 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2675  const Quadrature<dim> &q) const
2676 {
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());
2681 
2682  return data_ptr;
2683 }
2684 
2685 
2686 
2687 template <int dim, int spacedim>
2688 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2690  const UpdateFlags update_flags,
2691  const Quadrature<dim - 1> &quadrature) const
2692 {
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),
2698  quadrature.size());
2699 
2700  return data_ptr;
2701 }
2702 
2703 
2704 
2705 template <int dim, int spacedim>
2706 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
2708  const UpdateFlags update_flags,
2709  const Quadrature<dim - 1> &quadrature) const
2710 {
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),
2716  quadrature.size());
2717 
2718  return data_ptr;
2719 }
2720 
2721 
2722 
2723 template <int dim, int spacedim>
2726  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2727  const CellSimilarity::Similarity cell_similarity,
2728  const Quadrature<dim> & quadrature,
2729  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
2731  &output_data) const
2732 {
2733  // ensure that the following static_cast is really correct:
2734  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
2735  ExcInternalError());
2736  const InternalData &data = static_cast<const InternalData &>(internal_data);
2737 
2738  const unsigned int n_q_points = quadrature.size();
2739 
2740  // recompute the support points of the transformation of this
2741  // cell. we tried to be clever here in an earlier version of the
2742  // library by checking whether the cell is the same as the one we
2743  // had visited last, but it turns out to be difficult to determine
2744  // that because a cell for the purposes of a mapping is
2745  // characterized not just by its (triangulation, level, index)
2746  // triple, but also by the locations of its vertices, the manifold
2747  // object attached to the cell and all of its bounding faces/edges,
2748  // etc. to reliably test that the "cell" we are on is, therefore,
2749  // not easily done
2751  data.cell_of_current_support_points = cell;
2752 
2753  // if the order of the mapping is greater than 1, then do not reuse any cell
2754  // similarity information. This is necessary because the cell similarity
2755  // value is computed with just cell vertices and does not take into account
2756  // cell curvature.
2757  const CellSimilarity::Similarity computed_cell_similarity =
2758  (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
2759 
2760  if (dim > 1 && data.tensor_product_quadrature)
2761  {
2762  internal::MappingQGenericImplementation::
2763  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
2764  computed_cell_similarity,
2765  data,
2766  output_data.quadrature_points,
2767  output_data.jacobian_grads);
2768  }
2769  else
2770  {
2771  internal::MappingQGenericImplementation::maybe_compute_q_points<dim,
2772  spacedim>(
2774  data,
2775  output_data.quadrature_points);
2776 
2777  internal::MappingQGenericImplementation::maybe_update_Jacobians<dim,
2778  spacedim>(
2779  computed_cell_similarity,
2781  data);
2782 
2783  internal::MappingQGenericImplementation::maybe_update_jacobian_grads<
2784  dim,
2785  spacedim>(computed_cell_similarity,
2787  data,
2788  output_data.jacobian_grads);
2789  }
2790 
2791  internal::MappingQGenericImplementation::
2792  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
2793  computed_cell_similarity,
2795  data,
2796  output_data.jacobian_pushed_forward_grads);
2797 
2798  internal::MappingQGenericImplementation::
2799  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
2800  computed_cell_similarity,
2802  data,
2803  output_data.jacobian_2nd_derivatives);
2804 
2805  internal::MappingQGenericImplementation::
2806  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
2807  computed_cell_similarity,
2809  data,
2811 
2812  internal::MappingQGenericImplementation::
2813  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
2814  computed_cell_similarity,
2816  data,
2817  output_data.jacobian_3rd_derivatives);
2818 
2819  internal::MappingQGenericImplementation::
2820  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
2821  computed_cell_similarity,
2823  data,
2825 
2826  const UpdateFlags update_flags = data.update_each;
2827  const std::vector<double> &weights = quadrature.get_weights();
2828 
2829  // Multiply quadrature weights by absolute value of Jacobian determinants or
2830  // the area element g=sqrt(DX^t DX) in case of codim > 0
2831 
2832  if (update_flags & (update_normal_vectors | update_JxW_values))
2833  {
2834  AssertDimension(output_data.JxW_values.size(), n_q_points);
2835 
2836  Assert(!(update_flags & update_normal_vectors) ||
2837  (output_data.normal_vectors.size() == n_q_points),
2838  ExcDimensionMismatch(output_data.normal_vectors.size(),
2839  n_q_points));
2840 
2841 
2842  if (computed_cell_similarity != CellSimilarity::translation)
2843  for (unsigned int point = 0; point < n_q_points; ++point)
2844  {
2845  if (dim == spacedim)
2846  {
2847  const double det = data.contravariant[point].determinant();
2848 
2849  // check for distorted cells.
2850 
2851  // TODO: this allows for anisotropies of up to 1e6 in 3D and
2852  // 1e12 in 2D. might want to find a finer
2853  // (dimension-independent) criterion
2854  Assert(det >
2855  1e-12 * Utilities::fixed_power<dim>(
2856  cell->diameter() / std::sqrt(double(dim))),
2858  cell->center(), det, point)));
2859 
2860  output_data.JxW_values[point] = weights[point] * det;
2861  }
2862  // if dim==spacedim, then there is no cell normal to
2863  // compute. since this is for FEValues (and not FEFaceValues),
2864  // there are also no face normals to compute
2865  else // codim>0 case
2866  {
2867  Tensor<1, spacedim> DX_t[dim];
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];
2871 
2872  Tensor<2, dim> G; // First fundamental form
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];
2876 
2877  output_data.JxW_values[point] =
2878  std::sqrt(determinant(G)) * weights[point];
2879 
2880  if (computed_cell_similarity ==
2882  {
2883  // we only need to flip the normal
2884  if (update_flags & update_normal_vectors)
2885  output_data.normal_vectors[point] *= -1.;
2886  }
2887  else
2888  {
2889  if (update_flags & update_normal_vectors)
2890  {
2891  Assert(spacedim == dim + 1,
2892  ExcMessage(
2893  "There is no (unique) cell normal for " +
2895  "-dimensional cells in " +
2896  Utilities::int_to_string(spacedim) +
2897  "-dimensional space. This only works if the "
2898  "space dimension is one greater than the "
2899  "dimensionality of the mesh cells."));
2900 
2901  if (dim == 1)
2902  output_data.normal_vectors[point] =
2903  cross_product_2d(-DX_t[0]);
2904  else // dim == 2
2905  output_data.normal_vectors[point] =
2906  cross_product_3d(DX_t[0], DX_t[1]);
2907 
2908  output_data.normal_vectors[point] /=
2909  output_data.normal_vectors[point].norm();
2910 
2911  if (cell->direction_flag() == false)
2912  output_data.normal_vectors[point] *= -1.;
2913  }
2914  }
2915  } // codim>0 case
2916  }
2917  }
2918 
2919 
2920 
2921  // copy values from InternalData to vector given by reference
2922  if (update_flags & update_jacobians)
2923  {
2924  AssertDimension(output_data.jacobians.size(), n_q_points);
2925  if (computed_cell_similarity != CellSimilarity::translation)
2926  for (unsigned int point = 0; point < n_q_points; ++point)
2927  output_data.jacobians[point] = data.contravariant[point];
2928  }
2929 
2930  // copy values from InternalData to vector given by reference
2931  if (update_flags & update_inverse_jacobians)
2932  {
2933  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
2934  if (computed_cell_similarity != CellSimilarity::translation)
2935  for (unsigned int point = 0; point < n_q_points; ++point)
2936  output_data.inverse_jacobians[point] =
2937  data.covariant[point].transpose();
2938  }
2939 
2940  return computed_cell_similarity;
2941 }
2942 
2943 
2944 
2945 namespace internal
2946 {
2947  namespace MappingQGenericImplementation
2948  {
2949  namespace
2950  {
2960  template <int dim, int spacedim>
2961  void
2962  maybe_compute_face_data(
2963  const ::MappingQGeneric<dim, spacedim> &mapping,
2964  const typename ::Triangulation<dim, spacedim>::cell_iterator
2965  & cell,
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
2971  &data,
2973  &output_data)
2974  {
2975  const UpdateFlags update_flags = data.update_each;
2976 
2977  if (update_flags &
2980  {
2981  if (update_flags & update_boundary_forms)
2982  AssertDimension(output_data.boundary_forms.size(), n_q_points);
2983  if (update_flags & update_normal_vectors)
2984  AssertDimension(output_data.normal_vectors.size(), n_q_points);
2985  if (update_flags & update_JxW_values)
2986  AssertDimension(output_data.JxW_values.size(), n_q_points);
2987 
2988  Assert(data.aux.size() + 1 >= dim, ExcInternalError());
2989 
2990  // first compute some common data that is used for evaluating
2991  // all of the flags below
2992 
2993  // map the unit tangentials to the real cell. checking for d!=dim-1
2994  // eliminates compiler warnings regarding unsigned int expressions <
2995  // 0.
2996  for (unsigned int d = 0; d != dim - 1; ++d)
2997  {
2999  data.unit_tangentials.size(),
3000  ExcInternalError());
3001  Assert(
3002  data.aux[d].size() <=
3003  data
3004  .unit_tangentials[face_no +
3006  .size(),
3007  ExcInternalError());
3008 
3009  mapping.transform(
3011  data
3012  .unit_tangentials[face_no +
3015  data,
3016  make_array_view(data.aux[d]));
3017  }
3018 
3019  if (update_flags & update_boundary_forms)
3020  {
3021  // if dim==spacedim, we can use the unit tangentials to compute
3022  // the boundary form by simply taking the cross product
3023  if (dim == spacedim)
3024  {
3025  for (unsigned int i = 0; i < n_q_points; ++i)
3026  switch (dim)
3027  {
3028  case 1:
3029  // in 1d, we don't have access to any of the
3030  // data.aux fields (because it has only dim-1
3031  // components), but we can still compute the
3032  // boundary form by simply looking at the number of
3033  // the face
3034  output_data.boundary_forms[i][0] =
3035  (face_no == 0 ? -1 : +1);
3036  break;
3037  case 2:
3038  output_data.boundary_forms[i] =
3039  cross_product_2d(data.aux[0][i]);
3040  break;
3041  case 3:
3042  output_data.boundary_forms[i] =
3043  cross_product_3d(data.aux[0][i], data.aux[1][i]);
3044  break;
3045  default:
3046  Assert(false, ExcNotImplemented());
3047  }
3048  }
3049  else //(dim < spacedim)
3050  {
3051  // in the codim-one case, the boundary form results from the
3052  // cross product of all the face tangential vectors and the
3053  // cell normal vector
3054  //
3055  // to compute the cell normal, use the same method used in
3056  // fill_fe_values for cells above
3057  AssertDimension(data.contravariant.size(), n_q_points);
3058 
3059  for (unsigned int point = 0; point < n_q_points; ++point)
3060  {
3061  if (dim == 1)
3062  {
3063  // J is a tangent vector
3064  output_data.boundary_forms[point] =
3065  data.contravariant[point].transpose()[0];
3066  output_data.boundary_forms[point] /=
3067  (face_no == 0 ? -1. : +1.) *
3068  output_data.boundary_forms[point].norm();
3069  }
3070 
3071  if (dim == 2)
3072  {
3074  data.contravariant[point].transpose();
3075 
3076  Tensor<1, spacedim> cell_normal =
3077  cross_product_3d(DX_t[0], DX_t[1]);
3078  cell_normal /= cell_normal.norm();
3079 
3080  // then compute the face normal from the face
3081  // tangent and the cell normal:
3082  output_data.boundary_forms[point] =
3083  cross_product_3d(data.aux[0][point], cell_normal);
3084  }
3085  }
3086  }
3087  }
3088 
3089  if (update_flags & update_JxW_values)
3090  for (unsigned int i = 0; i < output_data.boundary_forms.size();
3091  ++i)
3092  {
3093  output_data.JxW_values[i] =
3094  output_data.boundary_forms[i].norm() * weights[i];
3095 
3096  if (subface_no != numbers::invalid_unsigned_int)
3097  {
3098  const double area_ratio =
3100  cell->subface_case(face_no), subface_no);
3101  output_data.JxW_values[i] *= area_ratio;
3102  }
3103  }
3104 
3105  if (update_flags & update_normal_vectors)
3106  for (unsigned int i = 0; i < output_data.normal_vectors.size();
3107  ++i)
3108  output_data.normal_vectors[i] =
3109  Point<spacedim>(output_data.boundary_forms[i] /
3110  output_data.boundary_forms[i].norm());
3111 
3112  if (update_flags & update_jacobians)
3113  for (unsigned int point = 0; point < n_q_points; ++point)
3114  output_data.jacobians[point] = data.contravariant[point];
3115 
3116  if (update_flags & update_inverse_jacobians)
3117  for (unsigned int point = 0; point < n_q_points; ++point)
3118  output_data.inverse_jacobians[point] =
3119  data.covariant[point].transpose();
3120  }
3121  }
3122 
3123 
3130  template <int dim, int spacedim>
3131  void
3132  do_fill_fe_face_values(
3133  const ::MappingQGeneric<dim, spacedim> &mapping,
3134  const typename ::Triangulation<dim, spacedim>::cell_iterator
3135  & cell,
3136  const unsigned int face_no,
3137  const unsigned int subface_no,
3138  const typename QProjector<dim>::DataSetDescriptor data_set,
3139  const Quadrature<dim - 1> & quadrature,
3140  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3141  &data,
3143  &output_data)
3144  {
3145  if (dim > 1 && data.tensor_product_quadrature)
3146  {
3147  maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
3149  data,
3150  output_data.quadrature_points,
3151  output_data.jacobian_grads);
3152  }
3153  else
3154  {
3155  maybe_compute_q_points<dim, spacedim>(
3156  data_set, data, output_data.quadrature_points);
3157  maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
3158  data_set,
3159  data);
3160  maybe_update_jacobian_grads<dim, spacedim>(
3161  CellSimilarity::none, data_set, data, output_data.jacobian_grads);
3162  }
3163  maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
3165  data_set,
3166  data,
3167  output_data.jacobian_pushed_forward_grads);
3168  maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
3170  data_set,
3171  data,
3172  output_data.jacobian_2nd_derivatives);
3173  maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
3175  data_set,
3176  data,
3178  maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
3180  data_set,
3181  data,
3182  output_data.jacobian_3rd_derivatives);
3183  maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
3185  data_set,
3186  data,
3188 
3189  maybe_compute_face_data(mapping,
3190  cell,
3191  face_no,
3192  subface_no,
3193  quadrature.size(),
3194  quadrature.get_weights(),
3195  data,
3196  output_data);
3197  }
3198  } // namespace
3199  } // namespace MappingQGenericImplementation
3200 } // namespace internal
3201 
3202 
3203 
3204 template <int dim, int spacedim>
3205 void
3207  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
3208  const unsigned int face_no,
3209  const Quadrature<dim - 1> & quadrature,
3210  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
3212  &output_data) const
3213 {
3214  // ensure that the following cast is really correct:
3215  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
3216  ExcInternalError());
3217  const InternalData &data = static_cast<const InternalData &>(internal_data);
3218 
3219  // if necessary, recompute the support points of the transformation of this
3220  // cell (note that we need to first check the triangulation pointer, since
3221  // otherwise the second test might trigger an exception if the triangulations
3222  // are not the same)
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))
3227  {
3229  data.cell_of_current_support_points = cell;
3230  }
3231 
3232  internal::MappingQGenericImplementation::do_fill_fe_face_values(
3233  *this,
3234  cell,
3235  face_no,
3238  cell->face_orientation(face_no),
3239  cell->face_flip(face_no),
3240  cell->face_rotation(face_no),
3241  quadrature.size()),
3242  quadrature,
3243  data,
3244  output_data);
3245 }
3246 
3247 
3248 
3249 template <int dim, int spacedim>
3250 void
3252  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
3253  const unsigned int face_no,
3254  const unsigned int subface_no,
3255  const Quadrature<dim - 1> & quadrature,
3256  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
3258  &output_data) const
3259 {
3260  // ensure that the following cast is really correct:
3261  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
3262  ExcInternalError());
3263  const InternalData &data = static_cast<const InternalData &>(internal_data);
3264 
3265  // if necessary, recompute the support points of the transformation of this
3266  // cell (note that we need to first check the triangulation pointer, since
3267  // otherwise the second test might trigger an exception if the triangulations
3268  // are not the same)
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))
3273  {
3275  data.cell_of_current_support_points = cell;
3276  }
3277 
3278  internal::MappingQGenericImplementation::do_fill_fe_face_values(
3279  *this,
3280  cell,
3281  face_no,
3282  subface_no,
3284  subface_no,
3285  cell->face_orientation(face_no),
3286  cell->face_flip(face_no),
3287  cell->face_rotation(face_no),
3288  quadrature.size(),
3289  cell->subface_case(face_no)),
3290  quadrature,
3291  data,
3292  output_data);
3293 }
3294 
3295 
3296 
3297 namespace internal
3298 {
3299  namespace MappingQGenericImplementation
3300  {
3301  namespace
3302  {
3303  template <int dim, int spacedim, int rank>
3304  void
3305  transform_fields(
3306  const ArrayView<const Tensor<rank, dim>> & input,
3307  const MappingKind mapping_kind,
3308  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3309  const ArrayView<Tensor<rank, spacedim>> & output)
3310  {
3311  AssertDimension(input.size(), output.size());
3312  Assert((dynamic_cast<const typename ::
3313  MappingQGeneric<dim, spacedim>::InternalData *>(
3314  &mapping_data) != nullptr),
3315  ExcInternalError());
3316  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3317  &data =
3318  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3319  InternalData &>(mapping_data);
3320 
3321  switch (mapping_kind)
3322  {
3323  case mapping_contravariant:
3324  {
3325  Assert(
3326  data.update_each & update_contravariant_transformation,
3328  "update_contravariant_transformation"));
3329 
3330  for (unsigned int i = 0; i < output.size(); ++i)
3331  output[i] =
3332  apply_transformation(data.contravariant[i], input[i]);
3333 
3334  return;
3335  }
3336 
3337  case mapping_piola:
3338  {
3339  Assert(
3340  data.update_each & update_contravariant_transformation,
3342  "update_contravariant_transformation"));
3343  Assert(
3344  data.update_each & update_volume_elements,
3346  "update_volume_elements"));
3347  Assert(rank == 1, ExcMessage("Only for rank 1"));
3348  if (rank != 1)
3349  return;
3350 
3351  for (unsigned int i = 0; i < output.size(); ++i)
3352  {
3353  output[i] =
3354  apply_transformation(data.contravariant[i], input[i]);
3355  output[i] /= data.volume_elements[i];
3356  }
3357  return;
3358  }
3359  // We still allow this operation as in the
3360  // reference cell Derivatives are Tensor
3361  // rather than DerivativeForm
3362  case mapping_covariant:
3363  {
3364  Assert(
3365  data.update_each & update_contravariant_transformation,
3367  "update_covariant_transformation"));
3368 
3369  for (unsigned int i = 0; i < output.size(); ++i)
3370  output[i] = apply_transformation(data.covariant[i], input[i]);
3371 
3372  return;
3373  }
3374 
3375  default:
3376  Assert(false, ExcNotImplemented());
3377  }
3378  }
3379 
3380 
3381  template <int dim, int spacedim, int rank>
3382  void
3383  transform_gradients(
3384  const ArrayView<const Tensor<rank, dim>> & input,
3385  const MappingKind mapping_kind,
3386  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3387  const ArrayView<Tensor<rank, spacedim>> & output)
3388  {
3389  AssertDimension(input.size(), output.size());
3390  Assert((dynamic_cast<const typename ::
3391  MappingQGeneric<dim, spacedim>::InternalData *>(
3392  &mapping_data) != nullptr),
3393  ExcInternalError());
3394  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3395  &data =
3396  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3397  InternalData &>(mapping_data);
3398 
3399  switch (mapping_kind)
3400  {
3402  {
3403  Assert(
3404  data.update_each & update_covariant_transformation,
3406  "update_covariant_transformation"));
3407  Assert(
3408  data.update_each & update_contravariant_transformation,
3410  "update_contravariant_transformation"));
3411  Assert(rank == 2, ExcMessage("Only for rank 2"));
3412 
3413  for (unsigned int i = 0; i < output.size(); ++i)
3414  {
3416  apply_transformation(data.contravariant[i],
3417  transpose(input[i]));
3418  output[i] =
3419  apply_transformation(data.covariant[i], A.transpose());
3420  }
3421 
3422  return;
3423  }
3424 
3426  {
3427  Assert(
3428  data.update_each & update_covariant_transformation,
3430  "update_covariant_transformation"));
3431  Assert(rank == 2, ExcMessage("Only for rank 2"));
3432 
3433  for (unsigned int i = 0; i < output.size(); ++i)
3434  {
3436  apply_transformation(data.covariant[i],
3437  transpose(input[i]));
3438  output[i] =
3439  apply_transformation(data.covariant[i], A.transpose());
3440  }
3441 
3442  return;
3443  }
3444 
3446  {
3447  Assert(
3448  data.update_each & update_covariant_transformation,
3450  "update_covariant_transformation"));
3451  Assert(
3452  data.update_each & update_contravariant_transformation,
3454  "update_contravariant_transformation"));
3455  Assert(
3456  data.update_each & update_volume_elements,
3458  "update_volume_elements"));
3459  Assert(rank == 2, ExcMessage("Only for rank 2"));
3460 
3461  for (unsigned int i = 0; i < output.size(); ++i)
3462  {
3464  apply_transformation(data.covariant[i], input[i]);
3465  const Tensor<2, spacedim> T =
3466  apply_transformation(data.contravariant[i],
3467  A.transpose());
3468 
3469  output[i] = transpose(T);
3470  output[i] /= data.volume_elements[i];
3471  }
3472 
3473  return;
3474  }
3475 
3476  default:
3477  Assert(false, ExcNotImplemented());
3478  }
3479  }
3480 
3481 
3482 
3483  template <int dim, int spacedim>
3484  void
3485  transform_hessians(
3486  const ArrayView<const Tensor<3, dim>> & input,
3487  const MappingKind mapping_kind,
3488  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3489  const ArrayView<Tensor<3, spacedim>> & output)
3490  {
3491  AssertDimension(input.size(), output.size());
3492  Assert((dynamic_cast<const typename ::
3493  MappingQGeneric<dim, spacedim>::InternalData *>(
3494  &mapping_data) != nullptr),
3495  ExcInternalError());
3496  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3497  &data =
3498  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3499  InternalData &>(mapping_data);
3500 
3501  switch (mapping_kind)
3502  {
3504  {
3505  Assert(
3506  data.update_each & update_covariant_transformation,
3508  "update_covariant_transformation"));
3509  Assert(
3510  data.update_each & update_contravariant_transformation,
3512  "update_contravariant_transformation"));
3513 
3514  for (unsigned int q = 0; q < output.size(); ++q)
3515  for (unsigned int i = 0; i < spacedim; ++i)
3516  {
3517  double tmp1[dim][dim];
3518  for (unsigned int J = 0; J < dim; ++J)
3519  for (unsigned int K = 0; K < dim; ++K)
3520  {
3521  tmp1[J][K] =
3522  data.contravariant[q][i][0] * input[q][0][J][K];
3523  for (unsigned int I = 1; I < dim; ++I)
3524  tmp1[J][K] +=
3525  data.contravariant[q][i][I] * input[q][I][J][K];
3526  }
3527  for (unsigned int j = 0; j < spacedim; ++j)
3528  {
3529  double tmp2[dim];
3530  for (unsigned int K = 0; K < dim; ++K)
3531  {
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];
3535  }
3536  for (unsigned int k = 0; k < spacedim; ++k)
3537  {
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];
3543  }
3544  }
3545  }
3546  return;
3547  }
3548 
3550  {
3551  Assert(
3552  data.update_each & update_covariant_transformation,
3554  "update_covariant_transformation"));
3555 
3556  for (unsigned int q = 0; q < output.size(); ++q)
3557  for (unsigned int i = 0; i < spacedim; ++i)
3558  {
3559  double tmp1[dim][dim];
3560  for (unsigned int J = 0; J < dim; ++J)
3561  for (unsigned int K = 0; K < dim; ++K)
3562  {
3563  tmp1[J][K] =
3564  data.covariant[q][i][0] * input[q][0][J][K];
3565  for (unsigned int I = 1; I < dim; ++I)
3566  tmp1[J][K] +=
3567  data.covariant[q][i][I] * input[q][I][J][K];
3568  }
3569  for (unsigned int j = 0; j < spacedim; ++j)
3570  {
3571  double tmp2[dim];
3572  for (unsigned int K = 0; K < dim; ++K)
3573  {
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];
3577  }
3578  for (unsigned int k = 0; k < spacedim; ++k)
3579  {
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];
3585  }
3586  }
3587  }
3588 
3589  return;
3590  }
3591 
3592  case mapping_piola_hessian:
3593  {
3594  Assert(
3595  data.update_each & update_covariant_transformation,
3597  "update_covariant_transformation"));
3598  Assert(
3599  data.update_each & update_contravariant_transformation,
3601  "update_contravariant_transformation"));
3602  Assert(
3603  data.update_each & update_volume_elements,
3605  "update_volume_elements"));
3606 
3607  for (unsigned int q = 0; q < output.size(); ++q)
3608  for (unsigned int i = 0; i < spacedim; ++i)
3609  {
3610  double factor[dim];
3611  for (unsigned int I = 0; I < dim; ++I)
3612  factor[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)
3617  {
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];
3621  }
3622  for (unsigned int j = 0; j < spacedim; ++j)
3623  {
3624  double tmp2[dim];
3625  for (unsigned int K = 0; K < dim; ++K)
3626  {
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];
3630  }
3631  for (unsigned int k = 0; k < spacedim; ++k)
3632  {
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];
3638  }
3639  }
3640  }
3641 
3642  return;
3643  }
3644 
3645  default:
3646  Assert(false, ExcNotImplemented());
3647  }
3648  }
3649 
3650 
3651 
3652  template <int dim, int spacedim, int rank>
3653  void
3654  transform_differential_forms(
3655  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
3656  const MappingKind mapping_kind,
3657  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3658  const ArrayView<Tensor<rank + 1, spacedim>> & output)
3659  {
3660  AssertDimension(input.size(), output.size());
3661  Assert((dynamic_cast<const typename ::
3662  MappingQGeneric<dim, spacedim>::InternalData *>(
3663  &mapping_data) != nullptr),
3664  ExcInternalError());
3665  const typename ::MappingQGeneric<dim, spacedim>::InternalData
3666  &data =
3667  static_cast<const typename ::MappingQGeneric<dim, spacedim>::
3668  InternalData &>(mapping_data);
3669 
3670  switch (mapping_kind)
3671  {
3672  case mapping_covariant:
3673  {
3674  Assert(
3675  data.update_each & update_contravariant_transformation,
3677  "update_covariant_transformation"));
3678 
3679  for (unsigned int i = 0; i < output.size(); ++i)
3680  output[i] = apply_transformation(data.covariant[i], input[i]);
3681 
3682  return;
3683  }
3684  default:
3685  Assert(false, ExcNotImplemented());
3686  }
3687  }
3688  } // namespace
3689  } // namespace MappingQGenericImplementation
3690 } // namespace internal
3691 
3692 
3693 
3694 template <int dim, int spacedim>
3695 void
3697  const ArrayView<const Tensor<1, dim>> & input,
3698  const MappingKind mapping_kind,
3699  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3700  const ArrayView<Tensor<1, spacedim>> & output) const
3701 {
3702  internal::MappingQGenericImplementation::transform_fields(input,
3703  mapping_kind,
3704  mapping_data,
3705  output);
3706 }
3707 
3708 
3709 
3710 template <int dim, int spacedim>
3711 void
3713  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
3714  const MappingKind mapping_kind,
3715  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3716  const ArrayView<Tensor<2, spacedim>> & output) const
3717 {
3718  internal::MappingQGenericImplementation::transform_differential_forms(
3719  input, mapping_kind, mapping_data, output);
3720 }
3721 
3722 
3723 
3724 template <int dim, int spacedim>
3725 void
3727  const ArrayView<const Tensor<2, dim>> & input,
3728  const MappingKind mapping_kind,
3729  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3730  const ArrayView<Tensor<2, spacedim>> & output) const
3731 {
3732  switch (mapping_kind)
3733  {
3734  case mapping_contravariant:
3735  internal::MappingQGenericImplementation::transform_fields(input,
3736  mapping_kind,
3737  mapping_data,
3738  output);
3739  return;
3740 
3744  internal::MappingQGenericImplementation::transform_gradients(
3745  input, mapping_kind, mapping_data, output);
3746  return;
3747  default:
3748  Assert(false, ExcNotImplemented());
3749  }
3750 }
3751 
3752 
3753 
3754 template <int dim, int spacedim>
3755 void
3757  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
3758  const MappingKind mapping_kind,
3759  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3760  const ArrayView<Tensor<3, spacedim>> & output) const
3761 {
3762  AssertDimension(input.size(), output.size());
3763  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
3764  ExcInternalError());
3765  const InternalData &data = static_cast<const InternalData &>(mapping_data);
3766 
3767  switch (mapping_kind)
3768  {
3770  {
3773  "update_covariant_transformation"));
3774 
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)
3778  {
3779  double tmp[dim];
3780  for (unsigned int K = 0; K < dim; ++K)
3781  {
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];
3785  }
3786  for (unsigned int k = 0; k < spacedim; ++k)
3787  {
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];
3791  }
3792  }
3793  return;
3794  }
3795 
3796  default:
3797  Assert(false, ExcNotImplemented());
3798  }
3799 }
3800 
3801 
3802 
3803 template <int dim, int spacedim>
3804 void
3806  const ArrayView<const Tensor<3, dim>> & input,
3807  const MappingKind mapping_kind,
3808  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
3809  const ArrayView<Tensor<3, spacedim>> & output) const
3810 {
3811  switch (mapping_kind)
3812  {
3813  case mapping_piola_hessian:
3816  internal::MappingQGenericImplementation::transform_hessians(
3817  input, mapping_kind, mapping_data, output);
3818  return;
3819  default:
3820  Assert(false, ExcNotImplemented());
3821  }
3822 }
3823 
3824 
3825 
3826 template <int dim, int spacedim>
3827 void
3829  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
3830  std::vector<Point<spacedim>> & a) const
3831 {
3832  // if we only need the midpoint, then ask for it.
3833  if (this->polynomial_degree == 2)
3834  {
3835  for (unsigned int line_no = 0;
3836  line_no < GeometryInfo<dim>::lines_per_cell;
3837  ++line_no)
3838  {
3839  const typename Triangulation<dim, spacedim>::line_iterator line =
3840  (dim == 1 ?
3841  static_cast<
3843  cell->line(line_no));
3844 
3845  const Manifold<dim, spacedim> &manifold =
3846  ((line->manifold_id() == numbers::flat_manifold_id) &&
3847  (dim < spacedim) ?
3848  cell->get_manifold() :
3849  line->get_manifold());
3850  a.push_back(manifold.get_new_point_on_line(line));
3851  }
3852  }
3853  else
3854  // otherwise call the more complicated functions and ask for inner points
3855  // from the manifold description
3856  {
3857  std::vector<Point<spacedim>> tmp_points;
3858  for (unsigned int line_no = 0;
3859  line_no < GeometryInfo<dim>::lines_per_cell;
3860  ++line_no)
3861  {
3862  const typename Triangulation<dim, spacedim>::line_iterator line =
3863  (dim == 1 ?
3864  static_cast<
3866  cell->line(line_no));
3867 
3868  const Manifold<dim, spacedim> &manifold =
3869  ((line->manifold_id() == numbers::flat_manifold_id) &&
3870  (dim < spacedim) ?
3871  cell->get_manifold() :
3872  line->get_manifold());
3873 
3874  const std::array<Point<spacedim>, 2> vertices{
3875  {cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
3876  cell->vertex(
3878 
3879  const std::size_t n_rows =
3881  a.resize(a.size() + n_rows);
3882  auto a_view = make_array_view(a.end() - n_rows, a.end());
3883  manifold.get_new_points(
3884  make_array_view(vertices.begin(), vertices.end()),
3886  a_view);
3887  }
3888  }
3889 }
3890 
3891 
3892 
3893 template <>
3894 void
3897  std::vector<Point<3>> & a) const
3898 {
3899  const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
3900 
3901  // used if face quad at boundary or entirely in the interior of the domain
3902  std::vector<Point<3>> tmp_points;
3903 
3904  // loop over all faces and collect points on them
3905  for (unsigned int face_no = 0; face_no < faces_per_cell; ++face_no)
3906  {
3907  const Triangulation<3>::face_iterator face = cell->face(face_no);
3908 
3909 #ifdef DEBUG
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);
3913  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face,
3914  lines_per_face = GeometryInfo<3>::lines_per_face;
3915 
3916  // some sanity checks up front
3917  for (unsigned int i = 0; i < vertices_per_face; ++i)
3918  Assert(face->vertex_index(i) ==
3919  cell->vertex_index(GeometryInfo<3>::face_to_cell_vertices(
3920  face_no, i, face_orientation, face_flip, face_rotation)),
3921  ExcInternalError());
3922 
3923  // indices of the lines that bound a face are given by GeometryInfo<3>::
3924  // face_to_cell_lines
3925  for (unsigned int i = 0; i < lines_per_face; ++i)
3926  Assert(face->line(i) ==
3928  face_no, i, face_orientation, face_flip, face_rotation)),
3929  ExcInternalError());
3930 #endif
3931  // extract the points surrounding a quad from the points
3932  // already computed. First get the 4 vertices and then the points on
3933  // the four lines
3934  boost::container::small_vector<Point<3>, 200> tmp_points(
3937  for (const unsigned int v : GeometryInfo<2>::vertex_indices())
3938  tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no, v)];
3939  if (polynomial_degree > 1)
3940  for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
3941  ++line)
3942  for (unsigned int i = 0; i < polynomial_degree - 1; ++i)
3943  tmp_points[4 + line * (polynomial_degree - 1) + i] =
3945  (polynomial_degree - 1) *
3946  GeometryInfo<3>::face_to_cell_lines(face_no, line) +
3947  i];
3948 
3949  const std::size_t n_rows =
3951  a.resize(a.size() + n_rows);
3952  auto a_view = make_array_view(a.end() - n_rows, a.end());
3953  face->get_manifold().get_new_points(
3954  make_array_view(tmp_points.begin(), tmp_points.end()),
3956  a_view);
3957  }
3958 }
3959 
3960 
3961 
3962 template <>
3963 void
3966  std::vector<Point<3>> & a) const
3967 {
3968  std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
3969  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
3970  vertices[i] = cell->vertex(i);
3971 
3972  Table<2, double> weights(Utilities::fixed_power<2>(polynomial_degree - 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)
3976  {
3978  line_support_points.point(q2 + 1)[0]);
3979  for (const unsigned int i : GeometryInfo<2>::vertex_indices())
3980  weights(q, i) = GeometryInfo<2>::d_linear_shape_function(point, i);
3981  }
3982 
3983  const std::size_t n_rows = weights.size(0);
3984  a.resize(a.size() + n_rows);
3985  auto a_view = make_array_view(a.end() - n_rows, a.end());
3986  cell->get_manifold().get_new_points(
3987  make_array_view(vertices.begin(), vertices.end()), weights, a_view);
3988 }
3989 
3990 
3991 
3992 template <int dim, int spacedim>
3993 void
3996  std::vector<Point<spacedim>> &) const
3997 {
3998  Assert(false, ExcInternalError());
3999 }
4000 
4001 
4002 
4003 template <int dim, int spacedim>
4004 std::vector<Point<spacedim>>
4006  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
4007 {
4008  // get the vertices first
4009  std::vector<Point<spacedim>> a;
4010  a.reserve(Utilities::fixed_power<dim>(polynomial_degree + 1));
4011  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4012  a.push_back(cell->vertex(i));
4013 
4014  if (this->polynomial_degree > 1)
4015  {
4016  // check if all entities have the same manifold id which is when we can
4017  // simply ask the manifold for all points. the transfinite manifold can
4018  // do the interpolation better than this class, so if we detect that we
4019  // do not have to change anything here
4020  Assert(dim <= 3, ExcImpossibleInDim(dim));
4021  bool all_manifold_ids_are_equal = (dim == spacedim);
4022  if (all_manifold_ids_are_equal &&
4024  &cell->get_manifold()) == nullptr)
4025  {
4026  for (auto f : GeometryInfo<dim>::face_indices())
4027  if (&cell->face(f)->get_manifold() != &cell->get_manifold())
4028  all_manifold_ids_are_equal = false;
4029 
4030  if (dim == 3)
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;
4034  }
4035 
4036  if (all_manifold_ids_are_equal)
4037  {
4038  const std::size_t n_rows = support_point_weights_cell.size(0);
4039  a.resize(a.size() + n_rows);
4040  auto a_view = make_array_view(a.end() - n_rows, a.end());
4041  cell->get_manifold().get_new_points(make_array_view(a.begin(),
4042  a.end() - n_rows),
4044  a_view);
4045  }
4046  else
4047  switch (dim)
4048  {
4049  case 1:
4050  add_line_support_points(cell, a);
4051  break;
4052  case 2:
4053  // in 2d, add the points on the four bounding lines to the
4054  // exterior (outer) points
4055  add_line_support_points(cell, a);
4056 
4057  // then get the interior support points
4058  if (dim != spacedim)
4059  add_quad_support_points(cell, a);
4060  else
4061  {
4062  const std::size_t n_rows =
4064  a.resize(a.size() + n_rows);
4065  auto a_view = make_array_view(a.end() - n_rows, a.end());
4066  cell->get_manifold().get_new_points(
4067  make_array_view(a.begin(), a.end() - n_rows),
4069  a_view);
4070  }
4071  break;
4072 
4073  case 3:
4074  // in 3d also add the points located on the boundary faces
4075  add_line_support_points(cell, a);
4076  add_quad_support_points(cell, a);
4077 
4078  // then compute the interior points
4079  {
4080  const std::size_t n_rows =
4082  a.resize(a.size() + n_rows);
4083  auto a_view = make_array_view(a.end() - n_rows, a.end());
4084  cell->get_manifold().get_new_points(
4085  make_array_view(a.begin(), a.end() - n_rows),
4087  a_view);
4088  }
4089  break;
4090 
4091  default:
4092  Assert(false, ExcNotImplemented());
4093  break;
4094  }
4095  }
4096 
4097  return a;
4098 }
4099 
4100 
4101 
4102 //--------------------------- Explicit instantiations -----------------------
4103 #include "mapping_q_generic.inst"
4104 
4105 
Transformed quadrature weights.
std::vector< Tensor< 2, dim > > shape_second_derivatives
static ::ExceptionBase & ExcTransformationFailed()
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
const types::manifold_id flat_manifold_id
Definition: types.h:273
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Tensor< 1, spacedim, Number > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number > &grad_F, const Tensor< 1, dim, Number > &d_x)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
const unsigned int polynomial_degree
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1424
Contravariant transformation.
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)
Table< 2, double > support_point_weights_cell
const std::vector< Point< dim > > & get_points() const
void reinit(const Quadrature< 1 > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
const std::vector< double > & get_weights() const
virtual void add_quad_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual void add_line_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, std::vector< Point< spacedim >> &a) const
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
Definition: fe_dgq.h:112
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping.cc:28
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:442
constexpr void clear()
const Point< dim > & point(const unsigned int i) const
QGaussLobatto< 1 > line_support_points
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Determinant of the Jacobian.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12013
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:318
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
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
MappingQGeneric(const unsigned int polynomial_degree)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
Point< 2 > second
Definition: grid_out.cc:4353
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:346
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Definition: point.h:111
InternalData(const unsigned int polynomial_degree)
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:2407
const unsigned int polynomial_degree
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:779
QGaussLobatto< 1 > line_support_points
T fixed_power(const T t)
Definition: utilities.h:1072
std::vector< Table< 2, double > > support_point_weights_perimeter_to_interior
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:316
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static const char T
void compute_shape_function_values(const std::vector< Point< dim >> &unit_points)
std::vector< std::vector< Tensor< 1, spacedim > > > aux
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:10521
#define Assert(cond, exc)
Definition: exceptions.h:1419
UpdateFlags
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim >> &grads, std::vector< Tensor< 2, dim >> &grad_grads, std::vector< Tensor< 3, dim >> &third_derivatives, std::vector< Tensor< 4, dim >> &fourth_derivatives) const override
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
Abstract base class for mapping classes.
Definition: mapping.h:302
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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:607
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::vector< Point< spacedim > > mapping_support_points
std::vector< Tensor< 3, dim > > shape_third_derivatives
VectorType::value_type * end(VectorType &V)
Point< 3 > vertices[4]
DerivativeForm< 1, spacedim, dim, Number > transpose() const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Expression fabs(const Expression &x)
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< double > volume_elements
Gradient of volume element.
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int n_shape_functions
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
std::vector< Tensor< 1, dim > > shape_derivatives
unsigned int size() const
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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
std::vector< Point< spacedim > > quadrature_points
static const char A
static Point< dim > project_to_unit_cell(const Point< dim > &p)
unsigned int get_degree() const
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const
Definition: manifold.cc:123
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
Definition: cuda.h:31
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
size_type size(const unsigned int i) const
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
Definition: tensor.h:450
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
virtual bool preserves_vertex_locations() const override
VectorType::value_type * begin(VectorType &V)
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Normal vectors.
Triangulation< dim, spacedim >::cell_iterator cell_of_current_support_points
virtual std::size_t memory_consumption() const override
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
static ::ExceptionBase & ExcNotImplemented()
static const types::blas_int zero
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
numbers::NumberTraits< Number >::real_type norm() const
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:349
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< Tensor< 4, dim > > shape_fourth_derivatives
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
double compute_value(const unsigned int i, const Point< dim > &p) const
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2381
std::vector< double > shape_values
static const bool value
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
bool is_tensor_product() const
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:10361
T max(const T &t, const MPI_Comm &mpi_communicator)
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
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:834
Covariant transformation.
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
std::vector< Tensor< 1, spacedim > > normal_vectors
internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< double > > shape_info
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
UpdateFlags update_each
Definition: mapping.h:630
static ::ExceptionBase & ExcInternalError()