Collaboration diagram for The interplay of UpdateFlags, Mapping, and FiniteElement in FEValues:

Introduction

In order to compute local contributions of an individual cell to the global matrix and right hand side, we usually employ two techniques:

First, the integral is transformed from the actual cell \(K\) to the unit/reference cell \(\hat K\). For example, for the Laplace equation, we transform
\[ A^K_{ij} = \int_K \nabla \varphi_i(\bf x) \cdot \nabla \varphi_j(\bf x) \; dx \]
into
\[ A^K_{ij} = \int_{\hat K} \left[ J^{-1}(\hat{\bf x}) \hat \nabla \varphi_i(\hat{\bf x}) \right] \cdot \left[ J^{-1}(\hat{\bf x}) \hat \nabla \varphi_j(\hat{\bf x}) \right] \; |\textrm{det}\; J(\hat{\bf x})| \;\; d\hat x, \]
where a hat indicates reference coordinates, and \(J(\hat{\bf x}_q)\) is the Jacobian \(\frac{\partial \bf F_K(\hat{\bf x})}{\partial\bf \hat x}\) of the mapping \(\bf x = \bf F_K(\hat{\bf x})\).
Second, this integral is then approximated through quadrature. This yields the formula
\[ A^K_{ij} = \sum_{q}J^{-1}\left[(\hat{\bf x}_q) \hat \nabla \varphi_i(\hat{\bf x}_q)\right] \cdot \left[J^{-1}(\hat{\bf x}_q) \hat \nabla \varphi_j(\hat{\bf x}_q)\right]\ |\textrm{det}\ J(\hat{\bf x}_q)| w_q, \]
where \(q\) indicates the index of the quadrature point, \(\hat{\bf x}_q\) its location on the reference cell, and \(w_q\) its weight.

In order to evaluate such an expression in an application code, we have to access three different kinds of objects: a quadrature object that describes locations \(\hat{\bf x}_q\) and weights \(w_q\) of quadrature points on the reference cell; a finite element object that describes the gradients \(\hat\nabla \varphi_i(\hat{\bf x}_q)\) of shape functions on the unit cell; and a mapping object that provides the Jacobian as well as its determinant. Dealing with all these objects would be cumbersome and error prone.

On the other hand, these three kinds of objects almost always appear together, and it is in fact very rare for deal.II application codes to do anything with quadrature, finite element, or mapping objects besides using them together. For this reason, deal.II uses the FEValues abstraction combining information on the shape functions, the geometry of the actual mesh cell and a quadrature rule on a reference cell. Upon construction it takes one object of each of the three mentioned categories. Later, it can be "re-initialized" for a concrete grid cell and then provides mapped quadrature points and weights, mapped shape function values and derivatives as well as some properties of the transformation from the reference cell to the actual mesh cell.

Since computation of any of these values is potentially expensive (for example when using high order mappings with high order elements), the FEValues class only computes what it is explicitly asked for. To this end, it takes a list of flags of type UpdateFlags at construction time specifying which quantities should be updated each time a cell is visited. In the case above, you want the gradients of the shape functions on the real cell, which is encoded by the flag update_gradients, as well as the product of the determinant of the Jacobian times the quadrature weight, which is mnemonically encoded using the term JxW and encoded in the flag update_JxW_values. Because these flags are represented by single bits in integer numbers, producing a set of flags amounts to setting multiple bits in an integer, which is facilitated using the operation update_gradients | update_JxW_values (in other words, and maybe slightly confusingly so, the operation "this operation and that operation" is represented by the expression "single-bit-in-an-integer-for-this-operation binary-or single-bit-in-an-integer-for-that-operation"). To make operations cheaper, FEValues and the mapping and finite element objects it depends on really only compute those pieces of information that you have specified in the update flags (plus some information necessary to compute what has been specified, see below), and not everything that could possibly be computed on a cell. This optimization makes it much cheaper to iterate over cells for assembly, but it also means that one should take care to provide the minimal set of flags possible.

In addition, once you pass a set of flags for what you want, the functions filling the data fields of FEValues are able to distinguish between values that have to be recomputed on each cell (for example mapped gradients) and quantities that do not change from cell to cell (for example the values of shape functions of the usual \(Q_p\) finite elements at the same quadrature points on different cells; this property does not hold for the shape functions of Raviart-Thomas elements, however, which must be rotated with the local cell). This allows further optimization of the computations underlying assembly.

Tracking dependencies

Let's say you want to compute the Laplace matrix as shown above. In that case, you need to specify the update_gradients flag (for \(\nabla\varphi_i(\bf x_q)\)) and the update_JxW_values flag (for computing \(|\textrm{det}\; J(\bf x_q)|w_q\)). Internally, however, the finite element requires the computation of the inverse of the full Jacobian matrix, \(J^{-1}(\bf x_q)\) (and not just the determinant of the matrix), and to compute the inverse of the Jacobian, it is also necessary to compute the Jacobian matrix first.

Since these are requirements that are not important to the user, it is not necessary to specify this in user code. Rather, given a set of update flags, the FEValues object first asks the finite element object what information it needs to compute in order to satisfy the user's request provided in the update flags. The finite element object may therefore add other flags to the update flags (e.g., in the example above, an FE_Q object would add update_covariant_transformation to the list, since that is the necessary transformation from \(\hat\nabla\hat\varphi_i(\hat{\bf x}_q)\) to \(\nabla\varphi_i(\bf x_q)\)). With these updated flags, FEValues then asks the mapping whether it also wants to add more flags to the list to satisfy the needs of both the user and the finite element object, by calling Mapping::requires_update_flags(). (This procedure of first asking the finite element and then the mapping does not have to be iterated because mappings never require information computed by the finite element classes, while finite element classes typically need information computed by mappings.) Using this final list, the FEValues object then asks both the finite element object and mapping object to create temporary structures into which they can store some temporary information that can be computed once and for all, and these flags will be used when re-computing data on each cell we will visit later on.

Update once or each

As outlined above, we have now determined the final set of things that are necessary to satisfy a user's desired pieces of information as conveyed by the update flags she provided. This information will then typically be queried on every cell the user code visits in a subsequent integration loop.

Given that many of the things mappings or finite element classes compute are potentially expensive, FEValues employs a system whereby mappings and finite element objects are encouraged to pre-compute information that can be computed once without reference to a concrete cell, and make use of this when asked to visit a particular cell of the mesh. An example is that the values of the shape functions of the common FE_Q element are defined on the reference cell, and the values on the actual cell are simply exactly the values on the reference cell – there is consequently no need to evaluate shape functions on every cell, but it is sufficient to do this once at the beginning, store the values somewhere, and when visiting a concrete cell simply copying these values from their temporary location to the output structure. (Note, however, that this is specific to the FE_Q element: this is not the case if we used a FE_RaviartThomas element, since there, computing the values of the shape functions on a cell involves knowing the Jacobian of the mapping which depends on the geometry of the cell we visit; thus, for this element, simply copying pre-computed information is not sufficient to evaluate the values of shape functions on a particular cell.)

To accommodate this structure, both mappings and finite element classes may internally split the update flags into two sets commonly referenced as update_once and update_each (though these names do not appear in any public interfaces). The former contains all those pieces of information that can be pre-computed once at the time the FEValues object starts to interact with a mapping or finite element, whereas the latter contains those flags corresponding to things that need to be computed on every cell. For example, if update_flags=update_values, then the FE_Q class will set update_once=update_values and update_each=0, whereas the Raviart-Thomas element will do it the other way around.

These sets of flags are intended to be mutually exclusive. There is, on the other hand, nothing that ever provides this decomposition to anything outside the mapping or finite element classes – it is a purely internal decomposition.

Generation of the actual data

As outlined above, data is computed at two different times: once at the beginning on the reference cell, and once whenever we move to an actual cell. The functions involved in each of these steps are discussed next:

Initialization

Computing data on the reference cell before we even visit the first real cell is a two-step process. First, the constructor of FEValues, FEFaceValues and FESubfaceValues, respectively, need to allow the Mapping and FiniteElement objects to set up internal data structures. These structures are internal in the following sense: the FEValues object asks the finite element and mapping objects to create an object of type FiniteElement::InternalDataBase and Mapping::InternalDataBase each; the actual finite element and mapping class may in fact create objects of a derived type if they wish to store some data beyond what these base classes already provide. The functions involved in this are

The FEValues object then takes over ownership of these objects and will destroy them at the end of the FEValues object's lifetime. After this, the FEValues object asks the FiniteElement and Mapping objects to fill these InternalDataBase objects with the data that pertains to what can and needs to be computed on the reference cell. This is done in these functions:

Reinitialization for a mesh cell

Once initialization is over and we call FEValues::reinit, FEFaceValues::reinit or FESubfaceValues::reinit to move to a concrete cell or face, we need to calculate the "update_each" kinds of data. This is done in the following functions:

FEValues::reinit() calls Mapping::fill_fe_values(), then FiniteElement::fill_fe_values()
FEFaceValues::reinit() calls Mapping::fill_fe_face_values(), then FiniteElement::fill_fe_face_values()
FESubfaceValues::reinit() calls Mapping::fill_fe_subface_values(), then FiniteElement::fill_fe_subface_values()

This is where the actual data fields for FEValues, stored in internal::FEValues::MappingRelatedData and internal::FEValues::FiniteElementRelatedData objects, are computed. These functions call the function in Mapping first, such that all the mapping data required by the finite element is available. Then, the FiniteElement function is called.