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