start white paper on more efficient quadrature
[dealii.wiki.git] / Tensor-product-structures-for-polynomials-and-quadrature.md
1 # Rationale
2
3 Most of the local objects in dealii have a tensor product
4 structure. Quadrature formulas are tensor products of one-dimensional
5 formulas. Shape functions are tensor products of polynomials in
6 different coordinate directions. Even spaces like Nédéléc and
7 Raviart-Thomas have such structures in each component. Katharina
8 Kormann and Martin Kronbichler have shown that exploiting these
9 structures can boost performance of local integrals
10 considerably. Therefore, adapting finite element and quadrature
11 classes to their concepts is a valuable task.
12
13 An important ingredient required for the optimization of local
14 integrals is that the length of all loops is known at compile
15 time. Therefore, the number of quadrature points and the number of
16 shape functions must become template arguments.
17
18 The class `FEValues` was originally introduced to allow for
19 optimizations based on the internal structure of finite element shape
20 function spaces, but this was never exploited. Additionally, `FEValues`
21 stores a lot of data to avoid recomputation. Data, which is highly
22 redundant. At a time where computers had a single working thread and
23 computations were expensive, this was in part justified. Now, that
24 most algorithms become memory-bound, we must consider retiring
25 `FEValues`.
26
27 ## Discussion
28
29 - Implementation of *hp* methods
30 - Efficient implementation of methods where we do not have tensor
31   product structure
32   - XFEM
33   - Overlapping nonmatching grids
34
35 # Tensor product quadrature
36
37 It is suggested to introduce a new class template inheriting from
38 `Quadrature`, say `TensorProductQuadrature`. In a compatibility mode,
39 this class fills the fields in `Quadrature` (deferred execution
40 preferred). Quadrature points are stored in a `TensorProductPoints`,
41 where the access operator `[]` returns a `Point<dim>` created on the
42 fly. Additionally, it allows direct access to the one-dimensional
43 point sets. Similarly, we can consider a class `TensorProductWeights`.
44
45 Note that while we are currently not making use of this (?),
46 `Quadrature` is not restricted to tensor products. Therefore,
47 `TensorProductQuadrature` is not going to be an equivalent
48 replacement.
49
50 ## Discussion
51
52 * We will need anisotropic tensor products. Do we want an isotropic
53   class as well?
54 * If anisotropic, are variadic templates the only option?
55 * The order of the quadrature point (do we make any assumptions yet?)
56   - Currently: lexicographic, low to high abscissa
57   - Low to high weight for numerical stability?
58   - End points of intervals first to conform to finite elements?
59 * Not inheriting from `Quadrature` would on one hand mean that we have
60   two independent quadrature sets, on the other hand we would not need
61   to link to the library, which is important for off-loading.
62
63 ## Draft of the classes
64
65 The convention for the indices shall be that `d=1,...,dim`, `di`
66 enumerates the coordinates in direction `d`, and `i` enumerates all
67 quadrature points.
68
69 ```
70 template <typename number, int dim, int...>
71 class TensorProductPoints
72 {
73   public:
74     template <typename number2>
75     void set_coordinates (const unsigned int d, const std::vector<number2>&);
76     
77     number operator() (const unsigned int d, const unsigned int di) const;
78     Point<dim> point(const unsigned int i) const;
79 };
80 ```
81
82 The class `TensorProductWeights` is similar and
83 `TensorProductQuadrature` combines the two.
84
85 # Tensor product structure of polynomials
86
87 As Katharina and Martin point out, combining tensor product quadrature
88 with tensor product polynomials can generate very efficient code. This
89 has been known to the spectral element community for a while, but
90 widely ignored by us. It is implemented in the matrix free framework,
91 but only for certain elements and with limited functionality. One goal
92 of these new structures is making efficient matrix free computations
93 available throughout the library.
94
95 ### Example: Lagrange polynomials
96
97 Given polynomials of degree `p-1` and `q` quadrature points in each
98 direction, the current `FEValues` in 3D computes and stores `p^3 q^3`
99 shape function values, three times as many derivatives and nine times
100 as many second derivatives. In (isotropic) tensor product
101 representation, we need `p q` values, derivatives, and second
102 derivatives. This means, that now even higher derivatives and 4D become
103 feasible. As for computation of these values, some polynomial spaces
104 allow for a recursive computation of the values of higher order
105 members. Thus, here is another point where we can reduce the
106 computational complexity by a factor up to `p`.
107
108 ### Example: vector valued elements
109
110 `FE_Nedelec` as well as `FE_RaviartThomas` can be implemented such
111 that each component is an anisotropic tensor product.
112
113 ### Example: `FE_DGP`
114
115 This example differs from the previous ones in that all multivariate
116 polynomials are still products of 1D ones, but that not all possible
117 combinations are admitted. This is true for serendipity elements as
118 well. For these we need truncated tensor products.
119
120 ## Discussion
121
122 - The matrix free code currently does a renumbering to tensor product
123   numbering of the shape functions, which is the inversion of a
124   similar process when we generate the functions. Should we reconsider
125   the ordering in `DoFHandler`?

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.