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 #ifndef dealii_mapping_h
17 #define dealii_mapping_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/array_view.h>
23 #include <deal.II/base/derivative_form.h>
24 
25 #include <deal.II/fe/fe_update_flags.h>
26 
27 #include <deal.II/grid/tria.h>
28 
29 #include <array>
30 #include <cmath>
31 #include <memory>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 template <typename ElementType, typename MemorySpaceType>
36 class ArrayView;
37 template <int dim>
38 class Quadrature;
39 template <int dim, int spacedim>
40 class FEValues;
41 template <int dim, int spacedim>
42 class FEValuesBase;
43 template <int dim, int spacedim>
44 class FEValues;
45 template <int dim, int spacedim>
46 class FEFaceValues;
47 template <int dim, int spacedim>
48 class FESubfaceValues;
49 
50 
51 /**
52  * The transformation kind used for the Mapping::transform() functions.
53  *
54  * Special finite elements may need special Mapping from the reference cell to
55  * the actual mesh cell. In order to be most flexible, this enum provides an
56  * extensible interface for arbitrary transformations. Nevertheless, these
57  * must be implemented in the transform() functions of inheriting classes in
58  * order to work.
59  *
60  * @ingroup mapping
61  */
62 enum MappingKind
63 {
64   /**
65    * No mapping, i.e., shape functions are not mapped from a reference cell
66    * but instead are defined right on the real-space cell.
67    */
68   mapping_none = 0x0000,
69 
70   /**
71    * Covariant mapping (see Mapping::transform() for details).
72    */
73   mapping_covariant = 0x0001,
74 
75   /**
76    * Contravariant mapping (see Mapping::transform() for details).
77    */
78   mapping_contravariant = 0x0002,
79 
80   /**
81    * Mapping of the gradient of a covariant vector field (see
82    * Mapping::transform() for details).
83    */
84   mapping_covariant_gradient = 0x0003,
85 
86   /**
87    * Mapping of the gradient of a contravariant vector field (see
88    * Mapping::transform() for details).
89    */
90   mapping_contravariant_gradient = 0x0004,
91 
92   /**
93    * The Piola transform usually used for Hdiv elements. Piola transform is
94    * the standard transformation of vector valued elements in H<sup>div</sup>.
95    * It amounts to a contravariant transformation scaled by the inverse of the
96    * volume element.
97    */
98   mapping_piola = 0x0100,
99 
100   /**
101    * Transformation for the gradient of a vector field corresponding to a
102    * mapping_piola transformation (see Mapping::transform() for details).
103    */
104   mapping_piola_gradient = 0x0101,
105 
106   /**
107    * The mapping used for Nedelec elements.
108    *
109    * Curl-conforming elements are mapped as covariant vectors. Nevertheless,
110    * we introduce a separate mapping kind, such that we can use the same flag
111    * for the vector and its gradient (see Mapping::transform() for details).
112    */
113   mapping_nedelec = 0x0200,
114 
115   /**
116    * The mapping used for Raviart-Thomas elements.
117    */
118   mapping_raviart_thomas = 0x0300,
119 
120   /**
121    * The mapping used for BDM elements.
122    */
123   mapping_bdm = mapping_raviart_thomas,
124 
125   /**
126    * The mappings for 2-forms and third order tensors.
127    *
128    * These are mappings typpically applied to hessians transformed to the
129    * reference cell.
130    *
131    * Mapping of the hessian of a covariant vector field (see
132    * Mapping::transform() for details).
133    */
134   mapping_covariant_hessian,
135 
136   /**
137    * Mapping of the hessian of a contravariant vector field (see
138    * Mapping::transform() for details).
139    */
140   mapping_contravariant_hessian,
141 
142   /**
143    * Mapping of the hessian of a piola vector field (see Mapping::transform()
144    * for details).
145    */
146   mapping_piola_hessian
147 };
148 
149 
150 /**
151  * @short Abstract base class for mapping classes.
152  *
153  * This class declares the interface for the functionality to describe
154  * mappings from the reference (unit) cell to a cell in real space, as well as
155  * for filling the information necessary to use the FEValues, FEFaceValues,
156  * and FESubfaceValues classes. Concrete implementations of these interfaces
157  * are provided in derived classes.
158  *
159  * <h3>Mathematics of the mapping</h3>
160  *
161  * The mapping is a transformation $\mathbf x = \mathbf F_K(\hat{\mathbf  x})$
162  * which maps points $\hat{\mathbf x}$ in the reference cell
163  * $[0,1]^\text{dim}$ to points $\mathbf x$ in the actual grid cell
164  * $K\subset{\mathbb R}^\text{spacedim}$. Many of the applications of such
165  * mappings require the Jacobian of this mapping, $J(\hat{\mathbf x}) =
166  * \hat\nabla {\mathbf F}_K(\hat{\mathbf  x})$. For instance, if
167  * dim=spacedim=2, we have
168  * @f[
169  * J(\hat{\mathbf  x}) = \left(\begin{matrix}
170  * \frac{\partial x}{\partial \hat x} & \frac{\partial x}{\partial \hat y}
171  * \\
172  * \frac{\partial y}{\partial \hat x} & \frac{\partial y}{\partial \hat y}
173  * \end{matrix}\right)
174  * @f]
175  *
176  * <h4>%Mapping of scalar functions</h4>
177  *
178  * The shape functions of scalar finite elements are typically defined on a
179  * reference cell and are then simply mapped according to the rule
180  * @f[
181  * \varphi(\mathbf x) = \varphi\bigl(\mathbf F_K(\hat{\mathbf  x})\bigr)
182  * = \hat \varphi(\hat{\mathbf  x}).
183  * @f]
184  *
185  *
186  * <h4>%Mapping of integrals</h4>
187  *
188  * Using simply a change of variables, integrals of scalar functions over a
189  * cell $K$ can be expressed as an integral over the reference cell $\hat K$.
190  * Specifically, The volume form $d\hat x$ is transformed so that
191  * @f[
192  *  \int_K u(\mathbf x)\,dx = \int_{\hat K} \hat
193  * u(\hat{\mathbf  x}) \left|\text{det}J(\hat{\mathbf  x})\right|
194  * \,d\hat x.
195  * @f]
196  *
197  * In expressions where such integrals are approximated by quadrature, this
198  * then leads to terms of the form
199  * @f[
200  *  \int_K u(\mathbf x)\,dx
201  *  \approx
202  *  \sum_{q}
203  *  \hat u(\hat{\mathbf  x}_q)
204  *  \underbrace{\left|\text{det}J(\hat{\mathbf  x}_q)\right| w_q}_{=:
205  * \text{JxW}_q}.
206  * @f]
207  * Here, the weights $\text{JxW}_q$ of each quadrature point (where <i>JxW</i>
208  * mnemonically stands for <i>Jacobian times Quadrature Weights</i>) take the
209  * role of the $dx$ in the original integral. Consequently, they appear in all
210  * code that computes integrals approximated by quadrature, and are accessed
211  * by FEValues::JxW().
212  *
213  * @todo Document what happens in the codimension-1 case.
214  *
215  *
216  * <h4>%Mapping of vector fields, differential forms and gradients of vector
217  * fields</h4>
218  *
219  * The transformation of vector fields or differential forms (gradients of
220  * scalar functions) $\mathbf v$, and gradients of vector fields $\mathbf T$
221  * follows the general form
222  *
223  * @f[
224  * \mathbf v(\mathbf x) = \mathbf A(\hat{\mathbf  x})
225  * \hat{\mathbf  v}(\hat{\mathbf  x}),
226  * \qquad
227  * \mathbf T(\mathbf x) = \mathbf A(\hat{\mathbf  x})
228  * \hat{\mathbf  T}(\hat{\mathbf  x}) \mathbf B(\hat{\mathbf  x}).
229  * @f]
230  * The differential forms <b>A</b> and <b>B</b> are determined by the kind of
231  * object being transformed. These transformations are performed through the
232  * transform() functions, and the type of object being transformed is
233  * specified by their MappingKind argument. See the documentation there for
234  * possible choices.
235  *
236  * <h4>Derivatives of the mapping</h4>
237  *
238  * Some applications require the derivatives of the mapping, of which the
239  * first order derivative is the mapping Jacobian, $J_{iJ}(\hat{\mathbf
240  * x})=\frac{\partial x_i}{\partial \hat x_J}$, described above. Higher order
241  * derivatives of the mapping are similarly defined, for example the Jacobian
242  * derivative, $\hat H_{iJK}(\hat{\mathbf  x}) = \frac{\partial^2
243  * x_i}{\partial \hat x_J \partial \hat x_K}$, and the Jacobian second
244  * derivative, $\hat K_{iJKL}(\hat{\mathbf  x}) = \frac{\partial^3
245  * x_i}{\partial \hat x_J \partial \hat x_K \partial \hat x_L}$. It is also
246  * useful to define the "pushed-forward" versions of the higher order
247  * derivatives: the Jacobian pushed-forward derivative, $H_{ijk}(\hat{\mathbf
248  * x}) = \frac{\partial^2 x_i}{\partial \hat x_J \partial \hat
249  * x_K}(J_{jJ})^{-1}(J_{kK})^{-1}$, and the Jacobian pushed-forward second
250  * derivative, $K_{ijkl}(\hat{\mathbf  x}) = \frac{\partial^3 x_i}{\partial
251  * \hat x_J \partial \hat x_K \partial \hat
252  * x_L}(J_{jJ})^{-1}(J_{kK})^{-1}(J_{lL})^{-1}$. These pushed-forward versions
253  * can be used to compute the higher order derivatives of functions defined on
254  * the reference cell with respect to the real cell coordinates. For instance,
255  * the Jacobian derivative with respect to the real cell coordinates is given
256  * by:
257  *
258  * @f[
259  * \frac{\partial}{\partial x_j}\left[J_{iJ}(\hat{\mathbf  x})\right] =
260  * H_{ikn}(\hat{\mathbf  x})J_{nJ}(\hat{\mathbf  x}),
261  * @f]
262  * and the derivative of the Jacobian inverse with respect to the real cell
263  * coordinates is similarly given by:
264  * @f[
265  * \frac{\partial}{\partial x_j}\left[\left(J_{iJ}(\hat{\mathbf
266  * x})\right)^{-1}\right] = -H_{nik}(\hat{\mathbf  x})\left(J_{nJ}(\hat{\mathbf
267  * x})\right)^{-1}.
268  * @f]
269  *
270  * In a similar fashion, higher order derivatives, with respect to the real
271  * cell coordinates, of functions defined on the reference cell can be defined
272  * using the Jacobian pushed-forward higher-order derivatives. For example,
273  * the derivative, with respect to the real cell coordinates, of the Jacobian
274  * pushed-forward derivative is given by:
275  *
276  * @f[
277  * \frac{\partial}{\partial x_l}\left[H_{ijk}(\hat{\mathbf  x})\right] =
278  * K_{ijkl}(\hat{\mathbf  x}) -H_{mjl}(\hat{\mathbf  x})H_{imk}(\hat{\mathbf
279  * x})-H_{mkl}(\hat{\mathbf  x})H_{imj}(\hat{\mathbf  x}).
280  * @f]
281  *
282  * <h3>References</h3>
283  *
284  * A general publication on differential geometry and finite elements is the
285  * survey
286  * <ul>
287  * <li>Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. <i>Finite
288  * element exterior calculus: from Hodge theory to numerical stability.</i>
289  * Bull. Amer. Math. Soc. (N.S.), 47:281-354, 2010. <a
290  * href="http://dx.doi.org/10.1090/S0273-0979-10-01278-4">DOI:
291  * 10.1090/S0273-0979-10-01278-4</a>.
292  * </ul>
293  *
294  * The description of the Piola transform has been taken from the <a
295  * href="http://www.math.uh.edu/~rohop/spring_05/downloads/">lecture notes</a>
296  * by Ronald H. W. Hoppe, University of Houston, Chapter 7.
297  *
298  * @ingroup mapping
299  */
300 template <int dim, int spacedim = dim>
301 class Mapping : public Subscriptor
302 {
303 public:
304   /**
305    * Virtual destructor.
306    */
307   virtual ~Mapping() override = default;
308 
309   /**
310    * Return a pointer to a copy of the present object. The caller of this copy
311    * then assumes ownership of it.
312    *
313    * The function is declared abstract virtual in this base class, and derived
314    * classes will have to implement it.
315    *
316    * This function is mainly used by the hp::MappingCollection class.
317    */
318   virtual std::unique_ptr<Mapping<dim, spacedim>>
319   clone() const = 0;
320 
321   /**
322    * Return the mapped vertices of a cell.
323    *
324    * Most of the time, these values will simply be the coordinates of the
325    * vertices of a cell as returned by <code>cell-@>vertex(v)</code> for
326    * vertex <code>v</code>, i.e., information stored by the triangulation.
327    * However, there are also mappings that add displacements or choose
328    * completely different locations, e.g., MappingQEulerian,
329    * MappingQ1Eulerian, or MappingFEField.
330    *
331    * The default implementation of this function simply returns the
332    * information stored by the triangulation, i.e.,
333    * <code>cell-@>vertex(v)</code>.
334    */
335   virtual std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
336   get_vertices(
337     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
338 
339   /**
340    * Return the mapped center of a cell.
341    *
342    * If you are using a (bi-,tri-)linear mapping that preserves vertex
343    * locations, this function simply returns the value also produced by
344    * `cell->center()`. However, there are also mappings that add displacements
345    * or choose completely different locations, e.g., MappingQEulerian,
346    * MappingQ1Eulerian, or MappingFEField, and mappings based on high order
347    * polynomials, for which the center may not coincide with the average of
348    * the vertex locations.
349    *
350    * By default, this function returns the push forward of the center of the
351    * reference cell. If the parameter
352    * @p map_center_of_reference_cell is set to false, than the return value
353    * will be the average of the vertex locations, as returned by the
354    * get_vertices() method.
355    *
356    * @param[in] cell The cell for which you want to compute the center
357    * @param[in] map_center_of_reference_cell A flag that switches the algorithm
358    * for the computation of the cell center from
359    * transform_unit_to_real_cell() applied to the center of the reference cell
360    * to computing the vertex averages.
361    */
362   virtual Point<spacedim>
363   get_center(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
364              const bool map_center_of_reference_cell = true) const;
365 
366   /**
367    * Return the bounding box of a mapped cell.
368    *
369    * If you are using a (bi-,tri-)linear mapping that preserves vertex
370    * locations, this function simply returns the value also produced by
371    * `cell->bounding_box()`. However, there are also mappings that add
372    * displacements or choose completely different locations, e.g.,
373    * MappingQEulerian, MappingQ1Eulerian, or MappingFEField.
374    *
375    * This function returns the bounding box containing all the vertices of the
376    * cell, as returned by the get_vertices() method. Beware of the fact that
377    * for higher order mappings this bounding box is only an approximation of the
378    * true bounding box, since it does not take into account curved faces, and it
379    * may be smaller than the true bounding box.
380    *
381    * @param[in] cell The cell for which you want to compute the bounding box
382    */
383   virtual BoundingBox<spacedim>
384   get_bounding_box(
385     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
386 
387   /**
388    * Return whether the mapping preserves vertex locations. In other words,
389    * this function returns whether the mapped location of the reference cell
390    * vertices (given by GeometryInfo::unit_cell_vertex()) equals the result of
391    * <code>cell-@>vertex()</code> (i.e., information stored by the
392    * triangulation).
393    *
394    * For example, implementations in derived classes return @p true for
395    * MappingQ, MappingQGeneric, MappingCartesian, but @p false for
396    * MappingQEulerian, MappingQ1Eulerian, and MappingFEField.
397    */
398   virtual bool
399   preserves_vertex_locations() const = 0;
400 
401   /**
402    * @name Mapping points between reference and real cells
403    * @{
404    */
405 
406   /**
407    * Map the point @p p on the unit cell to the corresponding point on the
408    * real cell @p cell.
409    *
410    * @param cell Iterator to the cell that will be used to define the mapping.
411    * @param p Location of a point on the reference cell.
412    * @return The location of the reference point mapped to real space using
413    * the mapping defined by the class derived from the current one that
414    * implements the mapping, and the coordinates of the cell identified by the
415    * first argument.
416    */
417   virtual Point<spacedim>
418   transform_unit_to_real_cell(
419     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
420     const Point<dim> &                                          p) const = 0;
421 
422   /**
423    * Map the point @p p on the real @p cell to the corresponding point on the
424    * unit cell, and return its coordinates. This function provides the inverse
425    * of the mapping provided by transform_unit_to_real_cell().
426    *
427    * In the codimension one case, this function returns the normal projection
428    * of the real point @p p on the curve or surface identified by the @p cell.
429    *
430    * @note Polynomial mappings from the reference (unit) cell coordinates to
431    * the coordinate system of a real cell are not always invertible if the
432    * point for which the inverse mapping is to be computed lies outside the
433    * cell's boundaries. In such cases, the current function may fail to
434    * compute a point on the reference cell whose image under the mapping
435    * equals the given point @p p.  If this is the case then this function
436    * throws an exception of type Mapping::ExcTransformationFailed . Whether
437    * the given point @p p lies outside the cell can therefore be determined by
438    * checking whether the returned reference coordinates lie inside or outside
439    * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell()) or
440    * whether the exception mentioned above has been thrown.
441    *
442    * @param cell Iterator to the cell that will be used to define the mapping.
443    * @param p Location of a point on the given cell.
444    * @return The reference cell location of the point that when mapped to real
445    * space equals the coordinates given by the second argument. This mapping
446    * uses the mapping defined by the class derived from the current one that
447    * implements the mapping, and the coordinates of the cell identified by the
448    * first argument.
449    */
450   virtual Point<dim>
451   transform_real_to_unit_cell(
452     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
453     const Point<spacedim> &                                     p) const = 0;
454 
455   /**
456    * Transform the point @p p on the real @p cell to the corresponding point
457    * on the reference cell, and then project this point to a (dim-1)-dimensional
458    * point in the coordinate system of the face with
459    * the given face number @p face_no. Ideally the point @p p is near the face
460    * @p face_no, but any point in the cell can technically be projected.
461    *
462    * This function does not make physical sense when dim=1, so it throws an
463    * exception in this case.
464    */
465   Point<dim - 1>
466   project_real_point_to_unit_point_on_face(
467     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
468     const unsigned int                                          face_no,
469     const Point<spacedim> &                                     p) const;
470 
471   /**
472    * @}
473    */
474 
475 
476   /**
477    * @name Exceptions
478    * @{
479    */
480 
481   /**
482    * Exception
483    */
484   DeclException0(ExcInvalidData);
485 
486 
487   /**
488    * Computing the mapping between a real space point and a point in reference
489    * space failed, typically because the given point lies outside the cell
490    * where the inverse mapping is not unique.
491    *
492    * @ingroup Exceptions
493    */
494   DeclExceptionMsg(
495     ExcTransformationFailed,
496     "Computing the mapping between a real space point and a point in reference "
497     "space failed, typically because the given point lies outside the cell "
498     "where the inverse mapping is not unique.");
499 
500   /**
501    * deal.II assumes the Jacobian determinant to be positive. When the cell
502    * geometry is distorted under the image of the mapping, the mapping becomes
503    * invalid and this exception is thrown.
504    *
505    * @ingroup Exceptions
506    */
507   DeclException3(ExcDistortedMappedCell,
508                  Point<spacedim>,
509                  double,
510                  int,
511                  << "The image of the mapping applied to cell with center ["
512                  << arg1 << "] is distorted. The cell geometry or the "
513                  << "mapping are invalid, giving a non-positive volume "
514                  << "fraction of " << arg2 << " in quadrature point " << arg3
515                  << ".");
516 
517   /**
518    * @}
519    */
520 
521   /**
522    * @name Interface with FEValues
523    * @{
524    */
525 
526 public:
527   /**
528    * Base class for internal data of mapping objects. The internal mechanism
529    * is that upon construction of a FEValues object, it asks the mapping and
530    * finite element classes that are to be used to allocate memory for their
531    * own purpose in which they may store data that only needs to be computed
532    * once. For example, most finite elements will store the values of the
533    * shape functions at the quadrature points in this object, since they do
534    * not change from cell to cell and only need to be computed once. The same
535    * may be true for Mapping classes that want to only evaluate the shape
536    * functions used for mapping once at the quadrature points.
537    *
538    * Since different FEValues objects using different quadrature rules might
539    * access the same mapping object at the same time, it is necessary to
540    * create one such object per FEValues object. FEValues does this by calling
541    * Mapping::get_data(), or in reality the implementation of the
542    * corresponding function in derived classes. Ownership of the object
543    * created by Mapping::get_data() is then transferred to the FEValues
544    * object, but a reference to this object is passed to the mapping object
545    * every time it is asked to compute information on a concrete cell. This
546    * happens when FEValues::reinit() (or the corresponding classes in
547    * FEFaceValues and FESubfaceValues) call Mapping::fill_fe_values() (and
548    * similarly via Mapping::fill_fe_face_values() and
549    * Mapping::fill_fe_subface_values()).
550    *
551    * The purpose of this class is for mapping objects to store information
552    * that can be computed once at the beginning, on the reference cell, and to
553    * access it later when computing information on a concrete cell. As such,
554    * the object handed to Mapping::fill_fe_values() is marked as
555    * <code>const</code>, because the assumption is that at the time this
556    * information is used, it will not need to modified again. However, classes
557    * derived from Mapping can also use such objects for two other purposes:
558    *
559    * - To provide scratch space for computations that are done in
560    * Mapping::fill_fe_values() and similar functions. Some of the derived
561    * classes would like to use scratch arrays and it would be a waste of time
562    * to allocate these arrays every time this function is called, just to de-
563    * allocate it again at the end of the function. Rather, one could allocate
564    * this memory once as a member variable of the current class, and simply
565    * use it in Mapping::fill_fe_values().
566    * - After calling Mapping::fill_fe_values(), FEValues::reinit()
567    * calls FiniteElement::fill_fe_values() where the finite element computes
568    * values, gradients, etc of the shape functions using both information
569    * computed once at the beginning using a mechanism similar to the one
570    * described here (see FiniteElement::InternalDataBase) as well as the data
571    * already computed by Mapping::fill_fe_values(). As part of its work, some
572    * implementations of FiniteElement::fill_fe_values() need to transform
573    * shape function data, and they do so by calling Mapping::transform(). The
574    * call to the latter function also receives a reference to the
575    * Mapping::InternalDataBase object. Since Mapping::transform() may be
576    * called many times on each cell, it is sometimes worth for derived classes
577    * to compute some information only once in Mapping::fill_fe_values() and
578    * reuse it in Mapping::transform(). This information can also be stored in
579    * the classes that derived mapping classes derive from InternalDataBase.
580    *
581    * In both of these cases, the InternalDataBase object being passed around
582    * is "morally const", i.e., no external observer can tell whether a scratch
583    * array or some intermediate data for Mapping::transform() is being
584    * modified by Mapping::fill_fe_values() or not. Consequently, the
585    * InternalDataBase objects are always passed around as <code>const</code>
586    * objects. Derived classes that would like to make use of the two
587    * additional uses outlined above therefore need to mark the member
588    * variables they want to use for these purposes as <code>mutable</code> to
589    * allow for their modification despite the fact that the surrounding object
590    * is marked as <code>const</code>.
591    */
592   class InternalDataBase
593   {
594   public:
595     /**
596      * Constructor. Sets update_flags to @p update_default and @p first_cell
597      * to @p true.
598      */
599     InternalDataBase();
600 
601     /**
602      * Copy construction is forbidden.
603      */
604     InternalDataBase(const InternalDataBase &) = delete;
605 
606     /**
607      * Virtual destructor for derived classes
608      */
609     virtual ~InternalDataBase() = default;
610 
611     /**
612      * A set of update flags specifying the kind of information that an
613      * implementation of the Mapping interface needs to compute on each cell
614      * or face, i.e., in Mapping::fill_fe_values() and friends.
615      *
616      * This set of flags is stored here by implementations of
617      * Mapping::get_data(), Mapping::get_face_data(), or
618      * Mapping::get_subface_data(), and is that subset of the update flags
619      * passed to those functions that require re-computation on every cell.
620      * (The subset of the flags corresponding to information that can be
621      * computed once and for all already at the time of the call to
622      * Mapping::get_data() -- or an implementation of that interface -- need
623      * not be stored here because it has already been taken care of.)
624      */
625     UpdateFlags update_each;
626 
627     /**
628      * Return an estimate (in bytes) for the memory consumption of this object.
629      */
630     virtual std::size_t
631     memory_consumption() const;
632   };
633 
634 
635 protected:
636   /**
637    * Given a set of update flags, compute which other quantities <i>also</i>
638    * need to be computed in order to satisfy the request by the given flags.
639    * Then return the combination of the original set of flags and those just
640    * computed.
641    *
642    * As an example, if @p update_flags contains update_JxW_values (i.e., the
643    * product of the determinant of the Jacobian and the weights provided by
644    * the quadrature formula), a mapping may require the computation of the
645    * full Jacobian matrix in order to compute its determinant. They would then
646    * return not just update_JxW_values, but also update_jacobians. (This is
647    * not how it is actually done internally in the derived classes that
648    * compute the JxW values -- they set update_contravariant_transformation
649    * instead, from which the determinant can also be computed -- but this does
650    * not take away from the instructiveness of the example.)
651    *
652    * An extensive discussion of the interaction between this function and
653    * FEValues can be found in the
654    * @ref FE_vs_Mapping_vs_FEValues
655    * documentation module.
656    *
657    * @see UpdateFlags
658    */
659   virtual UpdateFlags
660   requires_update_flags(const UpdateFlags update_flags) const = 0;
661 
662   /**
663    * Create and return a pointer to an object into which mappings can store
664    * data that only needs to be computed once but that can then be used
665    * whenever the mapping is applied to a concrete cell (e.g., in the various
666    * transform() functions, as well as in the fill_fe_values(),
667    * fill_fe_face_values() and fill_fe_subface_values() that form the
668    * interface of mappings with the FEValues class).
669    *
670    * Derived classes will return pointers to objects of a type derived from
671    * Mapping::InternalDataBase (see there for more information) and may pre-
672    * compute some information already (in accordance with what will be asked
673    * of the mapping in the future, as specified by the update flags) and for
674    * the given quadrature object. Subsequent calls to transform() or
675    * fill_fe_values() and friends will then receive back the object created
676    * here (with the same set of update flags and for the same quadrature
677    * object). Derived classes can therefore pre-compute some information in
678    * their get_data() function and store it in the internal data object.
679    *
680    * The mapping classes do not keep track of the objects created by this
681    * function. Ownership will therefore rest with the caller.
682    *
683    * An extensive discussion of the interaction between this function and
684    * FEValues can be found in the
685    * @ref FE_vs_Mapping_vs_FEValues
686    * documentation module.
687    *
688    * @param update_flags A set of flags that define what is expected of the
689    * mapping class in future calls to transform() or the fill_fe_values()
690    * group of functions. This set of flags may contain flags that mappings do
691    * not know how to deal with (e.g., for information that is in fact computed
692    * by the finite element classes, such as UpdateFlags::update_values).
693    * Derived classes will need to store these flags, or at least that subset
694    * of flags that will require the mapping to perform any actions in
695    * fill_fe_values(), in InternalDataBase::update_each.
696    * @param quadrature The quadrature object for which mapping information
697    * will have to be computed. This includes the locations and weights of
698    * quadrature points.
699    * @return A pointer to a newly created object of type InternalDataBase (or
700    * a derived class). Ownership of this object passes to the calling
701    * function.
702    *
703    * @note C++ allows that virtual functions in derived classes may return
704    * pointers to objects not of type InternalDataBase but in fact pointers to
705    * objects of classes <i>derived</i> from InternalDataBase. (This feature is
706    * called "covariant return types".) This is useful in some contexts where
707    * the calling is within the derived class and will immediately make use of
708    * the returned object, knowing its real (derived) type.
709    */
710   virtual std::unique_ptr<InternalDataBase>
711   get_data(const UpdateFlags      update_flags,
712            const Quadrature<dim> &quadrature) const = 0;
713 
714   /**
715    * Like get_data(), but in preparation for later calls to transform() or
716    * fill_fe_face_values() that will need information about mappings from the
717    * reference face to a face of a concrete cell.
718    *
719    * @param update_flags A set of flags that define what is expected of the
720    * mapping class in future calls to transform() or the fill_fe_values()
721    * group of functions. This set of flags may contain flags that mappings do
722    * not know how to deal with (e.g., for information that is in fact computed
723    * by the finite element classes, such as UpdateFlags::update_values).
724    * Derived classes will need to store these flags, or at least that subset
725    * of flags that will require the mapping to perform any actions in
726    * fill_fe_values(), in InternalDataBase::update_each.
727    * @param quadrature The quadrature object for which mapping information
728    * will have to be computed. This includes the locations and weights of
729    * quadrature points.
730    * @return A pointer to a newly created object of type InternalDataBase (or
731    * a derived class). Ownership of this object passes to the calling
732    * function.
733    *
734    * @note C++ allows that virtual functions in derived classes may return
735    * pointers to objects not of type InternalDataBase but in fact pointers to
736    * objects of classes <i>derived</i> from InternalDataBase. (This feature is
737    * called "covariant return types".) This is useful in some contexts where
738    * the calling is within the derived class and will immediately make use of
739    * the returned object, knowing its real (derived) type.
740    */
741   virtual std::unique_ptr<InternalDataBase>
742   get_face_data(const UpdateFlags          update_flags,
743                 const Quadrature<dim - 1> &quadrature) const = 0;
744 
745   /**
746    * Like get_data() and get_face_data(), but in preparation for later calls
747    * to transform() or fill_fe_subface_values() that will need information
748    * about mappings from the reference face to a child of a face (i.e.,
749    * subface) of a concrete cell.
750    *
751    * @param update_flags A set of flags that define what is expected of the
752    * mapping class in future calls to transform() or the fill_fe_values()
753    * group of functions. This set of flags may contain flags that mappings do
754    * not know how to deal with (e.g., for information that is in fact computed
755    * by the finite element classes, such as UpdateFlags::update_values).
756    * Derived classes will need to store these flags, or at least that subset
757    * of flags that will require the mapping to perform any actions in
758    * fill_fe_values(), in InternalDataBase::update_each.
759    * @param quadrature The quadrature object for which mapping information
760    * will have to be computed. This includes the locations and weights of
761    * quadrature points.
762    * @return A pointer to a newly created object of type InternalDataBase (or
763    * a derived class). Ownership of this object passes to the calling
764    * function.
765    *
766    * @note C++ allows that virtual functions in derived classes may return
767    * pointers to objects not of type InternalDataBase but in fact pointers to
768    * objects of classes <i>derived</i> from InternalDataBase. (This feature is
769    * called "covariant return types".) This is useful in some contexts where
770    * the calling is within the derived class and will immediately make use of
771    * the returned object, knowing its real (derived) type.
772    */
773   virtual std::unique_ptr<InternalDataBase>
774   get_subface_data(const UpdateFlags          update_flags,
775                    const Quadrature<dim - 1> &quadrature) const = 0;
776 
777   /**
778    * Compute information about the mapping from the reference cell to the real
779    * cell indicated by the first argument to this function. Derived classes
780    * will have to implement this function based on the kind of mapping they
781    * represent. It is called by FEValues::reinit().
782    *
783    * Conceptually, this function's represents the application of the mapping
784    * $\mathbf x=\mathbf F_K(\hat {\mathbf x})$ from reference coordinates
785    * $\mathbf\in [0,1]^d$ to real space coordinates $\mathbf x$ for a given
786    * cell $K$. Its purpose is to compute the following kinds of data:
787    *
788    * - Data that results from the application of the mapping itself, e.g.,
789    * computing the location $\mathbf x_q = \mathbf F_K(\hat{\mathbf x}_q)$ of
790    * quadrature points on the real cell, and that is directly useful to users
791    * of FEValues, for example during assembly.
792    * - Data that is necessary for finite element implementations to compute
793    * their shape functions on the real cell. To this end, the
794    * FEValues::reinit() function calls FiniteElement::fill_fe_values() after
795    * the current function, and the output of this function serves as input to
796    * FiniteElement::fill_fe_values(). Examples of information that needs to be
797    * computed here for use by the finite element classes is the Jacobian of
798    * the mapping, $\hat\nabla \mathbf F_K(\hat{\mathbf x})$ or its inverse,
799    * for example to transform the gradients of shape functions on the
800    * reference cell to the gradients of shape functions on the real cell.
801    *
802    * The information computed by this function is used to fill the various
803    * member variables of the output argument of this function. Which of the
804    * member variables of that structure should be filled is determined by the
805    * update flags stored in the Mapping::InternalDataBase object passed to
806    * this function.
807    *
808    * An extensive discussion of the interaction between this function and
809    * FEValues can be found in the
810    * @ref FE_vs_Mapping_vs_FEValues
811    * documentation module.
812    *
813    * @param[in] cell The cell of the triangulation for which this function is
814    * to compute a mapping from the reference cell to.
815    * @param[in] cell_similarity Whether or not the cell given as first
816    * argument is simply a translation, rotation, etc of the cell for which
817    * this function was called the most recent time. This information is
818    * computed simply by matching the vertices (as stored by the Triangulation)
819    * between the previous and the current cell. The value passed here may be
820    * modified by implementations of this function and should then be returned
821    * (see the discussion of the return value of this function).
822    * @param[in] quadrature A reference to the quadrature formula in use for
823    * the current evaluation. This quadrature object is the same as the one
824    * used when creating the @p internal_data object. The object is used both
825    * to map the location of quadrature points, as well as to compute the JxW
826    * values for each quadrature point (which involves the quadrature weights).
827    * @param[in] internal_data A reference to an object previously created by
828    * get_data() and that may be used to store information the mapping can
829    * compute once on the reference cell. See the documentation of the
830    * Mapping::InternalDataBase class for an extensive description of the
831    * purpose of these objects.
832    * @param[out] output_data A reference to an object whose member variables
833    * should be computed. Not all of the members of this argument need to be
834    * filled; which ones need to be filled is determined by the update flags
835    * stored inside the @p internal_data object.
836    * @return An updated value of the @p cell_similarity argument to this
837    * function. The returned value will be used for the corresponding argument
838    * when FEValues::reinit() calls FiniteElement::fill_fe_values(). In most
839    * cases, derived classes will simply want to return the value passed for @p
840    * cell_similarity. However, implementations of this function may downgrade
841    * the level of cell similarity. This is, for example, the case for classes
842    * that take not only into account the locations of the vertices of a cell
843    * (as reported by the Triangulation), but also other information specific
844    * to the mapping. The purpose is that FEValues::reinit() can compute
845    * whether a cell is similar to the previous one only based on the cell's
846    * vertices, whereas the mapping may also consider displacement fields
847    * (e.g., in the MappingQ1Eulerian and MappingFEField classes). In such
848    * cases, the mapping may conclude that the previously computed cell
849    * similarity is too optimistic, and invalidate it for subsequent use in
850    * FiniteElement::fill_fe_values() by returning a less optimistic cell
851    * similarity value.
852    *
853    * @note FEValues ensures that this function is always called with the same
854    * pair of @p internal_data and @p output_data objects. In other words, if
855    * an implementation of this function knows that it has written a piece of
856    * data into the output argument in a previous call, then there is no need
857    * to copy it there again in a later call if the implementation knows that
858    * this is the same value.
859    */
860   virtual CellSimilarity::Similarity
861   fill_fe_values(
862     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
863     const CellSimilarity::Similarity                            cell_similarity,
864     const Quadrature<dim> &                                     quadrature,
865     const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
866     dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
867       &output_data) const = 0;
868 
869   /**
870    * This function is the equivalent to Mapping::fill_fe_values(), but for
871    * faces of cells. See there for an extensive discussion of its purpose. It
872    * is called by FEFaceValues::reinit().
873    *
874    * @param[in] cell The cell of the triangulation for which this function is
875    * to compute a mapping from the reference cell to.
876    * @param[in] face_no The number of the face of the given cell for which
877    * information is requested.
878    * @param[in] quadrature A reference to the quadrature formula in use for
879    * the current evaluation. This quadrature object is the same as the one
880    * used when creating the @p internal_data object. The object is used both
881    * to map the location of quadrature points, as well as to compute the JxW
882    * values for each quadrature point (which involves the quadrature weights).
883    * @param[in] internal_data A reference to an object previously created by
884    * get_data() and that may be used to store information the mapping can
885    * compute once on the reference cell. See the documentation of the
886    * Mapping::InternalDataBase class for an extensive description of the
887    * purpose of these objects.
888    * @param[out] output_data A reference to an object whose member variables
889    * should be computed. Not all of the members of this argument need to be
890    * filled; which ones need to be filled is determined by the update flags
891    * stored inside the @p internal_data object.
892    */
893   virtual void
894   fill_fe_face_values(
895     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
896     const unsigned int                                          face_no,
897     const Quadrature<dim - 1> &                                 quadrature,
898     const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
899     dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
900       &output_data) const = 0;
901 
902   /**
903    * This function is the equivalent to Mapping::fill_fe_values(), but for
904    * subfaces (i.e., children of faces) of cells. See there for an extensive
905    * discussion of its purpose. It is called by FESubfaceValues::reinit().
906    *
907    * @param[in] cell The cell of the triangulation for which this function is
908    * to compute a mapping from the reference cell to.
909    * @param[in] face_no The number of the face of the given cell for which
910    * information is requested.
911    * @param[in] subface_no The number of the child of a face of the given cell
912    * for which information is requested.
913    * @param[in] quadrature A reference to the quadrature formula in use for
914    * the current evaluation. This quadrature object is the same as the one
915    * used when creating the @p internal_data object. The object is used both
916    * to map the location of quadrature points, as well as to compute the JxW
917    * values for each quadrature point (which involves the quadrature weights).
918    * @param[in] internal_data A reference to an object previously created by
919    * get_data() and that may be used to store information the mapping can
920    * compute once on the reference cell. See the documentation of the
921    * Mapping::InternalDataBase class for an extensive description of the
922    * purpose of these objects.
923    * @param[out] output_data A reference to an object whose member variables
924    * should be computed. Not all of the members of this argument need to be
925    * filled; which ones need to be filled is determined by the update flags
926    * stored inside the @p internal_data object.
927    */
928   virtual void
929   fill_fe_subface_values(
930     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
931     const unsigned int                                          face_no,
932     const unsigned int                                          subface_no,
933     const Quadrature<dim - 1> &                                 quadrature,
934     const typename Mapping<dim, spacedim>::InternalDataBase &   internal_data,
935     dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
936       &output_data) const = 0;
937 
938   /**
939    * @}
940    */
941 
942 public:
943   /**
944    * @name Functions to transform tensors from reference to real coordinates
945    * @{
946    */
947 
948   /**
949    * Transform a field of vectors or 1-differential forms according to the
950    * selected MappingKind.
951    *
952    * @note Normally, this function is called by a finite element, filling
953    * FEValues objects. For this finite element, there should be an alias
954    * MappingKind like @p mapping_bdm, @p mapping_nedelec, etc. This alias
955    * should be preferred to using the kinds below.
956    *
957    * The mapping kinds currently implemented by derived classes are:
958    * <ul>
959    * <li> @p mapping_contravariant: maps a vector field on the reference cell
960    * to the physical cell through the Jacobian:
961    * @f[
962    * \mathbf u(\mathbf x) = J(\hat{\mathbf  x})\hat{\mathbf  u}(\hat{\mathbf
963    * x}).
964    * @f]
965    * In physics, this is usually referred to as the contravariant
966    * transformation. Mathematically, it is the push forward of a vector field.
967    *
968    * <li> @p mapping_covariant: maps a field of one-forms on the reference
969    * cell to a field of one-forms on the physical cell. (Theoretically this
970    * would refer to a DerivativeForm<1,dim,1> but we canonically identify this
971    * type with a Tensor<1,dim>). Mathematically, it is the pull back of the
972    * differential form
973    * @f[
974    * \mathbf u(\mathbf x) = J(\hat{\mathbf  x})(J(\hat{\mathbf  x})^{T}
975    * J(\hat{\mathbf  x}))^{-1}\hat{\mathbf u}(\hat{\mathbf  x}).
976    * @f]
977    * Gradients of scalar differentiable functions are transformed this way.
978    *
979    * In the case when dim=spacedim the previous formula reduces to
980    * @f[
981    * \mathbf u(\mathbf x) = J(\hat{\mathbf  x})^{-T}\hat{\mathbf
982    * u}(\hat{\mathbf  x})
983    * @f]
984    * because we assume that the mapping $\mathbf F_K$ is always invertible,
985    * and consequently its Jacobian $J$ is an invertible matrix.
986    *
987    * <li> @p mapping_piola: A field of <i>dim-1</i>-forms on the reference
988    * cell is also represented by a vector field, but again transforms
989    * differently, namely by the Piola transform
990    * @f[
991    *  \mathbf u(\mathbf x) = \frac{1}{\text{det}\;J(\hat{\mathbf x})}
992    * J(\hat{\mathbf x}) \hat{\mathbf  u}(\hat{\mathbf x}).
993    * @f]
994    * </ul>
995    *
996    * @param[in] input An array (or part of an array) of input objects that
997    * should be mapped.
998    * @param[in] kind The kind of mapping to be applied.
999    * @param[in] internal A pointer to an object of type
1000    * Mapping::InternalDataBase that contains information previously stored by
1001    * the mapping. The object pointed to was created by the get_data(),
1002    * get_face_data(), or get_subface_data() function, and will have been
1003    * updated as part of a call to fill_fe_values(), fill_fe_face_values(), or
1004    * fill_fe_subface_values() for the current cell, before calling the current
1005    * function. In other words, this object also represents with respect to
1006    * which cell the transformation should be applied to.
1007    * @param[out] output An array (or part of an array) into which the
1008    * transformed objects should be placed. (Note that the array view is @p
1009    * const, but the tensors it points to are not.)
1010    */
1011   virtual void
1012   transform(const ArrayView<const Tensor<1, dim>> &                  input,
1013             const MappingKind                                        kind,
1014             const typename Mapping<dim, spacedim>::InternalDataBase &internal,
1015             const ArrayView<Tensor<1, spacedim>> &output) const = 0;
1016 
1017   /**
1018    * Transform a field of differential forms from the reference cell to the
1019    * physical cell.  It is useful to think of $\mathbf{T} = \nabla \mathbf u$
1020    * and $\hat{\mathbf  T} = \hat \nabla \hat{\mathbf  u}$, with $\mathbf u$ a
1021    * vector field.  The mapping kinds currently implemented by derived classes
1022    * are:
1023    * <ul>
1024    * <li> @p mapping_covariant: maps a field of forms on the reference cell to
1025    * a field of forms on the physical cell. Mathematically, it is the pull
1026    * back of the differential form
1027    * @f[
1028    * \mathbf T(\mathbf x) = \hat{\mathbf  T}(\hat{\mathbf  x})
1029    *                        J(\hat{\mathbf  x})(J(\hat{\mathbf  x})^{T}
1030    * J(\hat{\mathbf  x}))^{-1}.
1031    * @f]
1032    * Jacobians of spacedim-vector valued differentiable functions are
1033    * transformed this way.
1034    *
1035    * In the case when dim=spacedim the previous formula reduces to
1036    * @f[
1037    * \mathbf T(\mathbf x) = \hat{\mathbf  u}(\hat{\mathbf  x})
1038    *                        J(\hat{\mathbf  x})^{-1}.
1039    * @f]
1040    * </ul>
1041    *
1042    * @note It would have been more reasonable to make this transform a
1043    * template function with the rank in <code>DerivativeForm@<1, dim,
1044    * rank@></code>. Unfortunately C++ does not allow templatized virtual
1045    * functions. This is why we identify <code>DerivativeForm@<1, dim,
1046    * 1@></code> with a <code>Tensor@<1,dim@></code> when using
1047    * mapping_covariant() in the function transform() above this one.
1048    *
1049    * @param[in] input An array (or part of an array) of input objects that
1050    * should be mapped.
1051    * @param[in] kind The kind of mapping to be applied.
1052    * @param[in] internal A pointer to an object of type
1053    * Mapping::InternalDataBase that contains information previously stored by
1054    * the mapping. The object pointed to was created by the get_data(),
1055    * get_face_data(), or get_subface_data() function, and will have been
1056    * updated as part of a call to fill_fe_values(), fill_fe_face_values(), or
1057    * fill_fe_subface_values() for the current cell, before calling the current
1058    * function. In other words, this object also represents with respect to
1059    * which cell the transformation should be applied to.
1060    * @param[out] output An array (or part of an array) into which the
1061    * transformed objects should be placed. (Note that the array view is @p
1062    * const, but the tensors it points to are not.)
1063    */
1064   virtual void
1065   transform(const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1066             const MappingKind                                        kind,
1067             const typename Mapping<dim, spacedim>::InternalDataBase &internal,
1068             const ArrayView<Tensor<2, spacedim>> &output) const = 0;
1069 
1070   /**
1071    * Transform a tensor field from the reference cell to the physical cell.
1072    * These tensors are usually the Jacobians in the reference cell of vector
1073    * fields that have been pulled back from the physical cell.  The mapping
1074    * kinds currently implemented by derived classes are:
1075    * <ul>
1076    * <li> @p mapping_contravariant_gradient: it assumes $\mathbf u(\mathbf x)
1077    * = J \hat{\mathbf  u}$ so that
1078    * @f[
1079    * \mathbf T(\mathbf x) =
1080    * J(\hat{\mathbf  x}) \hat{\mathbf  T}(\hat{\mathbf  x})
1081    * J(\hat{\mathbf  x})^{-1}.
1082    * @f]
1083    * <li> @p mapping_covariant_gradient: it assumes $\mathbf u(\mathbf x) =
1084    * J^{-T} \hat{\mathbf  u}$ so that
1085    * @f[
1086    * \mathbf T(\mathbf x) =
1087    * J(\hat{\mathbf  x})^{-T} \hat{\mathbf  T}(\hat{\mathbf  x})
1088    * J(\hat{\mathbf  x})^{-1}.
1089    * @f]
1090    * <li> @p mapping_piola_gradient: it assumes $\mathbf u(\mathbf x) =
1091    * \frac{1}{\text{det}\;J(\hat{\mathbf x})} J(\hat{\mathbf x}) \hat{\mathbf
1092    * u}(\hat{\mathbf x})$ so that
1093    * @f[
1094    * \mathbf T(\mathbf x) =
1095    * \frac{1}{\text{det}\;J(\hat{\mathbf x})}
1096    * J(\hat{\mathbf  x}) \hat{\mathbf  T}(\hat{\mathbf  x})
1097    * J(\hat{\mathbf  x})^{-1}.
1098    * @f]
1099    * </ul>
1100    *
1101    * @todo The formulas for mapping_covariant_gradient,
1102    * mapping_contravariant_gradient and mapping_piola_gradient are only true
1103    * as stated for linear mappings. If, for example, the mapping is bilinear
1104    * (or has a higher order polynomial degree) then there is a missing term
1105    * associated with the derivative of $J$.
1106    *
1107    * @param[in] input An array (or part of an array) of input objects that
1108    * should be mapped.
1109    * @param[in] kind The kind of mapping to be applied.
1110    * @param[in] internal A pointer to an object of type
1111    * Mapping::InternalDataBase that contains information previously stored by
1112    * the mapping. The object pointed to was created by the get_data(),
1113    * get_face_data(), or get_subface_data() function, and will have been
1114    * updated as part of a call to fill_fe_values(), fill_fe_face_values(), or
1115    * fill_fe_subface_values() for the current cell, before calling the current
1116    * function. In other words, this object also represents with respect to
1117    * which cell the transformation should be applied to.
1118    * @param[out] output An array (or part of an array) into which the
1119    * transformed objects should be placed. (Note that the array view is @p
1120    * const, but the tensors it points to are not.)
1121    */
1122   virtual void
1123   transform(const ArrayView<const Tensor<2, dim>> &                  input,
1124             const MappingKind                                        kind,
1125             const typename Mapping<dim, spacedim>::InternalDataBase &internal,
1126             const ArrayView<Tensor<2, spacedim>> &output) const = 0;
1127 
1128   /**
1129    * Transform a tensor field from the reference cell to the physical cell.
1130    * This tensors are most of times the hessians in the reference cell of
1131    * vector fields that have been pulled back from the physical cell.
1132    *
1133    * The mapping kinds currently implemented by derived classes are:
1134    * <ul>
1135    * <li> @p mapping_covariant_gradient: maps a field of forms on the
1136    * reference cell to a field of forms on the physical cell. Mathematically,
1137    * it is the pull back of the differential form
1138    * @f[
1139    * \mathbf T_{ijk}(\mathbf x) = \hat{\mathbf  T}_{iJK}(\hat{\mathbf  x})
1140    * J_{jJ}^{\dagger} J_{kK}^{\dagger}@f],
1141    *
1142    * where @f[ J^{\dagger} = J(\hat{\mathbf  x})(J(\hat{\mathbf  x})^{T}
1143    * J(\hat{\mathbf  x}))^{-1}.
1144    * @f]
1145    * </ul>
1146    *
1147    * Hessians of spacedim-vector valued differentiable functions are
1148    * transformed this way (After subtraction of the product of the derivative
1149    * with the Jacobian gradient).
1150    *
1151    * In the case when dim=spacedim the previous formula reduces to
1152    * @f[J^{\dagger} = J^{-1}@f]
1153    *
1154    * @param[in] input An array (or part of an array) of input objects that
1155    * should be mapped.
1156    * @param[in] kind The kind of mapping to be applied.
1157    * @param[in] internal A pointer to an object of type
1158    * Mapping::InternalDataBase that contains information previously stored by
1159    * the mapping. The object pointed to was created by the get_data(),
1160    * get_face_data(), or get_subface_data() function, and will have been
1161    * updated as part of a call to fill_fe_values(), fill_fe_face_values(), or
1162    * fill_fe_subface_values() for the current cell, before calling the current
1163    * function. In other words, this object also represents with respect to
1164    * which cell the transformation should be applied to.
1165    * @param[out] output An array (or part of an array) into which the
1166    * transformed objects should be placed. (Note that the array view is @p
1167    * const, but the tensors it points to are not.)
1168    */
1169   virtual void
1170   transform(const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1171             const MappingKind                                        kind,
1172             const typename Mapping<dim, spacedim>::InternalDataBase &internal,
1173             const ArrayView<Tensor<3, spacedim>> &output) const = 0;
1174 
1175   /**
1176    * Transform a field of 3-differential forms from the reference cell to the
1177    * physical cell.  It is useful to think of $\mathbf{T}_{ijk} = D^2_{jk}
1178    * \mathbf u_i$ and $\mathbf{\hat T}_{IJK} = \hat D^2_{JK} \mathbf{\hat
1179    * u}_I$, with $\mathbf u_i$ a vector field.
1180    *
1181    * The mapping kinds currently implemented by derived classes are:
1182    * <ul>
1183    * <li> @p mapping_contravariant_hessian: it assumes $\mathbf u_i(\mathbf x)
1184    * = J_{iI} \hat{\mathbf  u}_I$ so that
1185    * @f[
1186    * \mathbf T_{ijk}(\mathbf x) =
1187    * J_{iI}(\hat{\mathbf  x}) \hat{\mathbf  T}_{IJK}(\hat{\mathbf  x})
1188    * J_{jJ}(\hat{\mathbf  x})^{-1} J_{kK}(\hat{\mathbf  x})^{-1}.
1189    * @f]
1190    * <li> @p mapping_covariant_hessian: it assumes $\mathbf u_i(\mathbf x) =
1191    * J_{iI}^{-T} \hat{\mathbf  u}_I$ so that
1192    * @f[
1193    * \mathbf T_{ijk}(\mathbf x) =
1194    * J_iI(\hat{\mathbf  x})^{-1} \hat{\mathbf  T}_{IJK}(\hat{\mathbf  x})
1195    * J_{jJ}(\hat{\mathbf  x})^{-1} J_{kK}(\hat{\mathbf  x})^{-1}.
1196    * @f]
1197    * <li> @p mapping_piola_hessian: it assumes $\mathbf u_i(\mathbf x) =
1198    * \frac{1}{\text{det}\;J(\hat{\mathbf x})} J_{iI}(\hat{\mathbf x})
1199    * \hat{\mathbf u}(\hat{\mathbf x})$ so that
1200    * @f[
1201    * \mathbf T_{ijk}(\mathbf x) =
1202    * \frac{1}{\text{det}\;J(\hat{\mathbf x})}
1203    * J_{iI}(\hat{\mathbf  x}) \hat{\mathbf  T}_{IJK}(\hat{\mathbf  x})
1204    * J_{jJ}(\hat{\mathbf  x})^{-1} J_{kK}(\hat{\mathbf  x})^{-1}.
1205    * @f]
1206    * </ul>
1207    *
1208    * @param[in] input An array (or part of an array) of input objects that
1209    * should be mapped.
1210    * @param[in] kind The kind of mapping to be applied.
1211    * @param[in] internal A pointer to an object of type
1212    * Mapping::InternalDataBase that contains information previously stored by
1213    * the mapping. The object pointed to was created by the get_data(),
1214    * get_face_data(), or get_subface_data() function, and will have been
1215    * updated as part of a call to fill_fe_values(), fill_fe_face_values(), or
1216    * fill_fe_subface_values() for the current cell, before calling the current
1217    * function. In other words, this object also represents with respect to
1218    * which cell the transformation should be applied to.
1219    * @param[out] output An array (or part of an array) into which the
1220    * transformed objects should be placed.
1221    */
1222   virtual void
1223   transform(const ArrayView<const Tensor<3, dim>> &                  input,
1224             const MappingKind                                        kind,
1225             const typename Mapping<dim, spacedim>::InternalDataBase &internal,
1226             const ArrayView<Tensor<3, spacedim>> &output) const = 0;
1227 
1228   /**
1229    * @}
1230    */
1231 
1232 
1233   // Give class @p FEValues access to the private <tt>get_...data</tt> and
1234   // <tt>fill_fe_...values</tt> functions.
1235   friend class FEValuesBase<dim, spacedim>;
1236   friend class FEValues<dim, spacedim>;
1237   friend class FEFaceValues<dim, spacedim>;
1238   friend class FESubfaceValues<dim, spacedim>;
1239 };
1240 
1241 
1242 DEAL_II_NAMESPACE_CLOSE
1243 
1244 #endif
1245