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