Reference documentation for deal.II version Git 10f3285 2019-01-18 22:14:04 -0700
Classes | Enumerations
Matrix-free infrastructure


struct  internal::MatrixFreeFunctions::DoFInfo
class  FEEvaluationBase< dim, n_components_, Number, is_face >
class  FEEvaluationAccess< dim, n_components_, Number, is_face >
class  FEEvaluationAccess< dim, 1, Number, is_face >
class  FEEvaluationAccess< dim, dim, Number, is_face >
class  FEEvaluationAccess< 1, 1, Number, is_face >
class  FEEvaluation< dim, fe_degree, n_q_points_1d, n_components_, Number >
struct  internal::MatrixFreeFunctions::MappingInfoStorage< structdim, spacedim, Number >
struct  internal::MatrixFreeFunctions::MappingInfo< dim, Number >
class  MatrixFree< dim, Number >
struct  MatrixFree< dim, Number >::DoFHandlers
struct  MatrixFree< dim, Number >::AdditionalData
struct  internal::MatrixFreeFunctions::ShapeInfo< Number >


enum  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants : unsigned char {
  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::full, internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved, internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous, internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved_contiguous,
  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved_contiguous_strided, internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved_contiguous_mixed_strides
enum  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex : unsigned char { internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior = 0, internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior = 1, internal::MatrixFreeFunctions::DoFInfo::dof_access_cell = 2 }
enum  internal::MatrixFreeFunctions::GeometryType : unsigned char { internal::MatrixFreeFunctions::cartesian = 0, internal::MatrixFreeFunctions::affine = 1, internal::MatrixFreeFunctions::flat_faces = 2, internal::MatrixFreeFunctions::general = 3 }
enum  MatrixFree< dim, Number >::AdditionalData::TasksParallelScheme { MatrixFree< dim, Number >::AdditionalData::none = internal::MatrixFreeFunctions::TaskInfo::none, MatrixFree< dim, Number >::AdditionalData::partition_partition, MatrixFree< dim, Number >::AdditionalData::partition_color, MatrixFree< dim, Number >::AdditionalData::color = internal::MatrixFreeFunctions::TaskInfo::color }
enum  MatrixFree< dim, Number >::DataAccessOnFaces { MatrixFree< dim, Number >::DataAccessOnFaces::none, MatrixFree< dim, Number >::DataAccessOnFaces::values, MatrixFree< dim, Number >::DataAccessOnFaces::gradients, MatrixFree< dim, Number >::DataAccessOnFaces::unspecified }
enum  MatrixFree< dim, Number >::DoFHandlers::ActiveDoFHandler { MatrixFree< dim, Number >::DoFHandlers::usual, MatrixFree< dim, Number >::DoFHandlers::hp }
enum  internal::MatrixFreeFunctions::ElementType {
  internal::MatrixFreeFunctions::tensor_symmetric_collocation = 0, internal::MatrixFreeFunctions::tensor_symmetric_hermite = 1, internal::MatrixFreeFunctions::tensor_symmetric = 2, internal::MatrixFreeFunctions::tensor_general = 3,
  internal::MatrixFreeFunctions::truncated_tensor = 4, internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0 = 5

Detailed Description

This module describes the matrix-free infrastructure in deal.II. In essence, the framework provided by the FEEvaluation class on top of the data storage in MatrixFree is a specialized operator evaluation framework. It is currently only compatible with a subset of the elements provided by the library which have a special structure, namely those where the basis can be described as a tensor product of one-dimensional polynomials. This opens for efficient transformation between vector entries and values or gradients in quadrature points with a technique that is called sum factorization. This technique has its origin in the spectral element community, started by the work of Orszag in 1980. While this technique is initially nothing else than a particular technique for assembling vectors (or matrices) that is faster than the general-purpose vehicle FEValues, its efficiency makes it possible to use these integration facilities to directly evaluate the matrix-vector products in iterative solvers, rather than first assembling a matrix and then using that matrix for doing matrix-vector products. This step is initially non-intuitive and goes against what many people were taught in their mathematics and computer science education, including most of the deal.II developers, because it appears to be wasteful to re-compute integrals over and over again, instead of using precomputed data. However, as the tutorial programs step-37, step-48, and step-59 show, these concepts usually outperform traditional algorithms on modern computer architectures.

The two main reasons that favor matrix-free computations are the following:

  1. Matrix-free methods skip the storage of big global sparse matrices and compute the underlying weak forms on the fly. Since the memory transfer, i.e., the speed at which the data can be read from RAM memory, is the bottleneck for matrix-based computations rather than the actual arithmetic done using this data, a matrix-free evaluation that reads less data can be advantageous even if it does more computations. This concept is building upon a trend in computer architecture which is best described by the term memory wall, saying that compute performance has increased more rapidly than the memory performance. Thus, a certain degree of arithmetic operations is essentially for free, and this share has become larger during the last twenty years. It has enabled this radical algorithm switch going from a matrix-based to a matrix-free implementation of matrix-vector products for iterative solvers, besides their classical use in explicit time integration. Of course, the implementation must be efficient and there cannot be an excess in computations to make it a win in total. The deal.II library uses SIMD vectorization and highly optimized kernels based on templates of the polynomial degree to achieve this goal. To give a perspective, a sparse matrix-vector product for quadratic elements FE_Q used to be equally fast as the matrix-free implementation on processors designed around 2005-2007 (e.g. Pentium 4 or AMD Opteron Barcelona with 2-4 cores per chip). By 2018, the matrix-free evaluation is around eight times as fast (measured on Intel Skylake Server, 14 cores).
  2. Matrix-free methods have a better complexity per degree of freedom as the degree is increased, due to sum factorization. The work per degree of freedom increases as \(\mathcal O(k)\) in the degree \(k\) for matrix-free schemes, whereas it increases as \(\mathcal O(k^d)\) for matrix-based methods. This gives higher order schemes an edge. A particularly nice feature in matrix-free evaluation is that the \(\mathcal O(1)\) terms often dominate, so it appears that higher order methods are as fast in terms of evaluation time as low order ones, when they have the same number of degrees of freedom. For the implementation in deal.II, best throughput is typically achieved for polynomial degrees between three and six.

To summarize, matrix-free computations are the way to go for higher order elements (where higher order means everything except linear shape functions) and use in explicit time stepping (step-48) or iterative solvers where also preconditioning can be done in a matrix-free way, as demonstrated in the step-37 and step-59 tutorial programs.

The matrix-free evaluation infrastructure

The top level interface is provided by the FEEvaluation class, which also contains an extensive description of different use cases.

The FEEvaluation class hierarchy

The class FEEvaluation is derived from the class FEEvaluationAccess, which in turn inherits from FEEvaluationBase. The FEEvaluation class itself is templated not only on the dimension, the number of components, and the number type (e.g. double or float), but also on the polynomial degree and on the number of quadrature points per spatial direction. This information is used to pass the loop lengths in sum factorization to the respective kernels (see tensor_product_kernels.h and evaluation_kernels.h) and ensure optimal efficiency. All methods that access the vectors or provide access into the data fields on an individual quadrature point are inherited from FEEvaluationAccess.

The motivation for the FEEvaluationAccess classes is to allow for specializations of the value and gradient access of interpolated solution fields depending on the number of components. Whereas the base class FEEvaluationBase returns the gradient as a Tensor<1,n_components,Tensor<1,dim,VectorizedArray<Number>>>, with the outer tensor going over the components and the inner tensor going through the dim components of the gradient. For a scalar field, i.e., n_components=1, we can skip the outer tensor and simply use Tensor<1,dim,VectorizedArray<Number>> as the gradient type. Likewise, for a system with n_components=dim, the appropriate format for the gradient is Tensor<2,dim,VectorizedArray<Number>>.

The FEFaceEvaluation class

Face integrals, like for inhomogeneous Neumann conditions in continuous FEM or for the large class of discontinuous Galerkin schemes, require the evaluation of quantities on the quadrature point of a face, besides the cell integrals. The facilities for face evaluation are mostly shared with FEEvaluation, in the sense that FEFaceEvaluation also inherits from FEEvaluationAccess. All data fields regarding the degrees of freedom and shape functions can be reused, the latter because all information consists of 1D shape data anyway. With respect to the mapping data, however, a specialization is used because the data is of structdim=dim-1. As a consequence, the FEEvaluationAccess and FEEvaluationBase are given a template argument is_face to hold pointers to the cell and face mapping information, respectively. Besides access to the function values with FEEvaluationAccess::get_value() or gradients with FEEvaluationAccess::get_gradient(), the face evaluator also enables the access to the normal vector by FEEvaluationAccess::get_normal_vector() and a specialized field FEEvaluationAccess::get_normal_derivative(), which returns the derivative of the solution field normal to the face. This quantity is computed as the gradient (in real space) multiplied by the normal vector. The combination of the gradient and normal vector is typical of many (simple) second-order elliptic equations, such as the discretization of the Laplacian with the interior penalty method. If the gradient alone is not needed, the combined operation significantly reduces the data access, because only dim data entries for normal * Jacobian per quadrature point are necessary, as opposed to dim^2 fields for the Jacobian and dim fields for the normal when accessing them individually.

An important optimization for the computation of face integrals is to think about the amount of vector data that must be accessed to evaluate the integrals on a face. Think for example of the case of FE_DGQ, i.e., Lagrange polynomials that have some of their nodes on the element boundary. For evaluation of the function values, only \((k+1)^{d-1}\) degrees of freedom contribute via a non-zero basis function, whereas the rest of the \((k+1)^d\) basis functions evaluate to zero on that boundary. Since vector access is one of the bottlenecks in matrix-free computations, the access to the vector should be restricted to the interesting entries. To enable this setup, the method FEFaceEvaluation::gather_evaluate() (and FEFaceEvaluation::integrate_scatter() for the integration equivalent) combines the vector access with the interpolation to the quadrature points. There exist two specializations, including the aforementioned "non-zero" value case, which is stored as the field internal::MatrixFreeFunctions::ShapeInfo::nodal_at_cell_boundaries. A similar property is also possible for the case where only the value and the first derivative of a selected number of basis functions evaluate to nonzero on a face. The associated element type is FE_DGQHermite and the decision is stored on the property internal::MatrixFreeFunctions::tensor_symmetric_hermite. The decision on whether such an optimized kernel can be used is made automatically inside FEFaceEvaluation::gather_evaluate() and FEFaceEvaluation::integrate_scatter(). It might seem inefficient to do this decision for every integration task, but in the end this is a single if statement (conditional jump) that is easily predicable for a modern CPU as the decision is always the same inside an integration loop. (One only pays by somewhat increased compile times because the compiler needs to generate code for all paths, though).

The data storage through the MatrixFree class

The tasks performed by FEEvaluation and FEFaceEvaluation can be split into the three categories: index access into vectors, evaluation and integration on the unit cell, and operation on quadrature points including the geometry evaluation. This split is reflected by the major data fields contained by MatrixFree, using internal::MatrixFreeFunctions::DoFInfo, internal::MatrixFreeFunctions::ShapeInfo, and internal::MatrixFreeFunctions::MappingInfo for each these three categories, respectively. Their design principles and internal layout is described in the following subsections.

The main interface all these data structure adhere to is that integration tasks are broken down into a range of cells or faces that one can index into by a single integer index. The information about an integer range for the cell integrals, inner face integrals, and boundary integrals is provided by the class internal::MatrixFreeFunctions::TaskInfo, using the data fields cell_partition_data, face_partition_data, and boundary_partition_data. This class also contains information about subranges of indices for scheduling tasks in parallel using threads, and a grouping of the index range within {cell,face,boundary}_partition_data for interleaving cell and face integrals such that the access to vector entries for cell and face integrals re-uses data already in caches.

Index storage: the internal::MatrixFreeFunctions::DoFInfo struct

The main purpose of the DoFInfo class is to provide the indices consumed by the vector access functions FEEvaluationBase::read_dof_values() and FEEvaluationBase::distribute_local_to_global(). The indices are laid out as follows:

  1. Indices are stored in MPI-local index space to enable direct array access, rather than translating a global index into a local one. The latter would be absolutely detrimental to performance.
  2. The indices are stored in a field called internal::MatrixFreeFunctions::DoFInfo::dof_indices, which is a long index array. The access granularity in terms of a cell index is controlled by the auxiliary field internal::MatrixFreeFunctions::DoFInfo::row_starts that is similar to the rowstart index in a compressed matrix storage. The scheme supports variable lengths because we support hp adaptivity and index indirections due to constraints that are contained in the main index array. Due to vectorization over cells, the access granularity would initially be in terms of cell batches. However, we must be able to access also an individual cell, for example for face integrals with the batches of faces that are in general different from the cell batches and access is thus not linear. Furthermore, the support for multi-component systems becomes transparent if we provide a start index to every single component separately. Thus, the row_starts field is of length n_cell_batches()*VectorizedArray<Number>::n_array_elements*n_components.
  3. The translation between components within a system of multiple base elements is controlled by the four variables

    1. std::vector<unsigned int> n_components (components per base element),
    2. std::vector<unsigned int> start_components (translation from the base element to the unique component number),
    3. std::vector<unsigned int> component_to_base_index (translation from unique component number to base index), and
    4. std::vector<std::vector<unsigned int>> component_dof_indices_offset (offset of the particular component's range of degrees of freedom within the full list of degrees of freedom on a cell).

  4. Information to extract the FE index in hp adaptive computations.
  5. Information about the 'first access' into a particular vector entry that is used to zero out the entries in a destination vectors within the MatrixFree::loop shortly before accessing them the first time. This is used to avoid writing zeros to the whole vector which destroys data locality.

The setup of the data structures in internal::MatrixFreeFunctions::DoFInfo is done in internal::MatrixFreeFunctions::DoFInfo::read_dof_indices, where we first assume a very general finite element layout, be it continuous or discontinuous elements, and where we resolve the constraints due to hanging nodes. This initial step is done in the original ordering of cells. In a later stage, these cells will in general be rearranged to reflect the order by which we go through the cells in the final loop, and we also look for patterns in the DoF indices that can be utilized, such as contiguous index ranges within a cell. This reordering is done to enable overlap of communication and computation with MPI (if enabled) and to form better group of batches with vectorization over cells. The data storage of indices is linear in this final order, and arranged in internal::MatrixFreeFunctions::DoFInfo::reorder_cells.

Since the amount of data to store indices is not negligible, it is worthwhile to reduce the amount of data for special configurations that carry more structure. One example is the case of FE_DGQ where a single index per cell is enough to describe all its degrees of freedom, with the others coming in consecutive order. The class internal::MatrixFreeFunctions::DoFInfo contains a special array of vectors internal::MatrixFreeFunctions::DoFInfo::dof_indices_contiguous that contains a single number per cell. Since both cell and face integrals use different access patterns and the data in this special case is small, we are better off storing 3 such vectors, one for the faces decorated as interior (index 0), one for the faces decorated as exterior (index 1), and one for the cells (index 2), rather than using the indirection through internal::MatrixFreeFunctions::FaceInfo. There is a series of additional special storage formats available in DoFInfo. We refer to the documentation of the struct internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants for the options implemented in deal.II and their motivation.

Finally, the DoFInfo class also holds a shared pointer describing the parallel partitioning of the vectors. Due to the restriction of Utilities::MPI::Partitioner, the indices within an individual DoFHandler object passed to the MatrixFree::reinit() function must be contiguous within each MPI process, i.e., the local range must consist of at most one chunk. Besides the basic partitioner, the class also provides a set of tighter index sets involving only a subset of all ghost indices that are added to the vectors' ghost range. These exchange patterns are designed to be combined with the reduced index access via the internal::MatrixFreeFunctions::ShapeInfo::nodal_at_cell_boundaries for example.

The MatrixFree class supports multiple DoFHandler objects to be passed to the MatrixFree::reinit() function. For each of these DoFHandler objects, a separate internal::MatrixFreeFunctions::DoFInfo object is created. In MatrixFree, we store a std::vector of internal::MatrixFreeFunctions::DoFInfo objects to account for this fact.

The internal::MatrixFreeFunctions::ShapeInfo structure

The evaluation of one-dimensional shape functions on one-dimensional quadrature points is stored in the class internal::MatrixFreeFunctions::ShapeInfo. More precisely, we hold all function values, gradients, and hessians. Furthermore, the values and derivatives of shape functions on the faces, i.e., the points 0 and 1 of the unit interval, are also stored. For face integrals on hanging nodes, the coarser of the two adjacent cells must interpolate the values not to the full quadrature but to a subface only (evaluation points either scaled to [0, 1/2] or [1/2, 1]). This case is handled by the data fields values_within_subface, gradients_within_subface, and hessians_within_subface. This data structure also checks for symmetry in the shape functions with respect to the center of the reference cell (in which case the so-called even-odd transformation is applied, further reducing computations).

The internal::MatrixFreeFunctions::MappingInfo structure

The evaluated geometry information is stored in the class internal::MatrixFreeFunctions::MappingInfo. Similarly to the internal::MatrixFreeFunctions::DoFInfo class, multiple variants are possible within a single MatrixFree instance, in this case based on multiple quadrature formulas. Furthermore, separate data for both cells and faces is stored. Since there is more logic involved and there are synergies between the fields, the std::vector of fields is kept within internal::MatrixFreeFunctions::MappingInfo. The individual field is of type internal::MatrixFreeFunctions::MappingInfoStorage and holds arrays with the inverse Jacobians, the JxW values, normal vectors, normal vectors times inverse Jacobians (for FEEvaluationAccess::get_normal_derivative()), quadrature points in real space, and quadrature points on the reference element. We use an auxiliary index array that points to the start of the data for each cell, namely the data_index_offsets field for the Jacobians, JxW values, and normal vectors, and quadrature_point_offsets for the quadrature points. This offset enables hp adaptivity with variable lengths of fields similar to what is done for DoFInfo, but it also enables something we call geometry compression. In order to reduce the data access, we detect simple geometries of cells where Jacobians are constant within a cell or also across cells, using internal::MatrixFreeFunctions::MappingInfo::cell_type:

  1. Cartesian cells are cells where the Jacobian is diagonal and the same on every quadrature point of the cell. Only a single field needs to be stored per cell. Due to the similarity within the cell, we also check for other cell batches with the same Jacobian for all cells on the current processor. This can further reduce the memory access. Since the JxW values in the general case store the Jacobian times the quadrature weight, but we only want to keep a single field for a Cartesian cell, we misuse the name JxW in the Cartesian case and only store the determinant of the Jacobian, without the quadrature weight. As a consequence, we need to be careful in FEEvaluationBase::submit_value() and similar for this case as we must still multiply by the weight.
  2. Affine cells have constant Jacobian within the whole cell, so only a single field needs to be stored per cell. Due to the similarity within the cell, we also check for other cell batches with the same Jacobian for all cells on the current processor. Since the JxW values in the general case store the Jacobian times the quadrature weight, but we only want to keep a single field for an affine cell, we misuse the name JxW in the affine case, just as in the Cartesian case, and only store the determinant of the Jacobian, without the quadrature weight. As a consequence, we need to be careful in FEEvaluationBase::submit_value() and similar for this case as we must still multiply by the weight.
  3. On faces, we can have the special case that the normal vector is the same in all quadrature points also when the JxW values are different. This is the case for faces which are flat. To reduce the data access, we keep this as a third option of compressed indices in internal::MatrixFreeFunctions::GeometryType. As opposed to the Cartesian and affine case where only a single field is reserved in the arrays, flat faces keep a separate entry for all quadrature points (to keep a single index field data_index_offsets), but only access the first one.
  4. The general type indices a cell or face where no compression was found. In this case, we also do not look for opportunities to find the same pattern on more than one cell, even though such cases might exist such as for extruded meshes. This search operation, which is based on inserting data into a std::map using a custom floating point comparator FPArrayComparator, is efficient enough when a single data field per cell is used. However, it would be pretty expensive if done for all quadrature points of all cells (with many different cases).

The implementation of internal::MatrixFreeFunctions::MappingInfo is split into cell and face parts, so the two components can be easily held apart. What makes the code a bit awkward to read is the fact that we need to batch several objects together from the original scalar evaluation done in an FEValues object, that we need to identify data fields that are repetitive, and that we need to define the compression over several cells with a std::map for the Cartesian and affine cases.

The data computation part of internal::MatrixFreeFunctions::MappingInfo is parallelized by tasks besides the obvious MPI parallelization. Each processor computes the information on a subrange, before the data is eventually copied into a single combined data field.

Identification and parallelization of face integrals

The current scheme for face integrals in MatrixFree builds an independent list of tasks for all of the faces, rather than going through the 2*dim faces of a cell explicitly. This has the advantage that all information on a face is processed only once. Typical DG methods compute numerical fluxes that are conservative, i.e., that look the same from both sides of the face and whatever information leaves one cell must exactly enter the neighbor again. With this scheme, they must only be computed once. Also, this ensures that the geometry information must only be loaded once, too. (A possible disadvantage is that a face-based approach with independent numbering makes thread-based parallelism much more complicated than a cell-based approach where only the information of the current cell is written into and neighbors are only read.)

Since faces are independent of cells, they get their own layout of vectorization. It is the nature of faces that whatever is a contiguous batch of cells gets intertwined when seen from a batch of faces (where we only keep faces together that have the same face index within a cell and so on). The setup of the face loop, which is done in the file face_setup_internal.h, tries to provide face batches that at least partly resemble the cell patches, to increase the data locality. Along these lines, the face work is also interleaved with cell work in the typical MatrixFree::loop context, i.e., the cell_range and face_range arguments returned to the function calls are usually pretty short.

Since all integrals from both sides are performed at once, the question arises which one of the two processors at subdomain boundaries is assigned a face. The authors of this module have performed extensive experiments and found out that the scheme that is applied for the degree of freedom storage, namely to assign all items with possible overlap to a single processor, is pretty imbalanced with up to 20% difference in the number of faces. For better performance, a balanced scheme is implemented in face_setup_internal.h that splits all interfaces between each pair of processors into two chunks, one being done by one processor and one by the other. Even though this increases the number of messages to be sent over MPI, this is worth it because the load gets more balanced. Also, messages are rather big at around 5-50kB when the local problem size is 100,000 DoFs in 3D. At this message size, the latency is typically less than the throughput anyway.

Face data is not initialized by default, but must be triggered by the face update flags in MatrixFree::AdditionalData, namely mapping_update_flags_inner_faces or mapping_update_flags_boundary_faces set to a value different from update_default.

Invoking MatrixFree::loop

The MatrixFree class supports two types of loops over the entities. The first one, which has been available on the deal.II master branch since 2012, is to only perform cell integrals, using one of the three cell_loop functions that takes a function pointer to the cell operation. The second setup, introduced in 2018, is a loop where also face and/or boundary integrals can be performed, called simply loop. This takes three function pointers, addressing the cell work, inner face work, and boundary face work, respectively.

Besides scheduling the work in an appropriate way, the loop performs two more tasks:

  1. Data exchange on the src and dst vector, calling update_ghost_values() and compress(VectorOperation::add), respectively. The exchange can be done in an asynchronous fashion overlapping the communication with work on cells that do not need data from remote processors, if the respective flag MatrixFree::AdditionalData::overlap_communication_computation is set to true (the default).
  2. Zero the dst vector using the respective flag. The advantage of doing this inside the loop is that the loop knows which entries in the vectors are (first) touched by some of the subranges in the cell and face loops. Thus, it can zero the vector piece by piece to ensure that we do not need to access the vector entries twice (once for zeroing, once for adding contributions). This might seem like a tiny optimization, but indeed the operator evaluation can be so quick that simply zeroing a vector can take around 20% of the operator evaluation time, so it is really worth the effort! Since there is some experimentation to this parameter, the DoFInfo class keeps a static variable internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector where this can be adjusted (if someone thinks that something else would be better, for example because future computers look different than they did in 2018 when this was introduced).

Finally, the MatrixFree::loop functions also take an argument to pass the type of data access on face integrals, described by the struct MatrixFree::DataAccessOnFaces, to reduce the amount of data that needs to be exchanged between processors. Unfortunately, there is currently no way of communicating this information, that gets available inside FEFaceEvaluation by the combination of the type of evaluation (values and/or gradients) and the underlying shape functions, to the MatrixFree::loop for avoiding to manually set this kind of information at a second spot.

Enumeration Type Documentation

Enum for various storage variants of the indices. This storage format is used to implement more efficient indexing schemes in case the underlying data structures allow for them, and to inform the access functions in FEEvaluationBase::read_write_operation() on which array to get the data from. One example of more efficient storage is the enum value IndexStorageVariants::contiguous, which means that one can get the indices to all degrees of freedom of a cell by reading only the first index for each cell, whereas all subsequent indices are merely an offset from the first index.


This value indicates that no index compression was found and the only valid storage is to access all indices present on the cell, possibly including constraints. For a cell/face of this index type, the data access in FEEvaluationBase is directed to the array dof_indices with the index row_starts[cell_index*n_vectorization*n_components].first.


This value indicates that the indices are interleaved for access with vectorized gather and scatter operation. This storage variant is possible in case there are no constraints on the cell and the indices in the batch of cells are not pointing to the same global index in different slots of a vectorized array (in order to support scatter operations). For a cell/face of this index type, the data access in FEEvaluationBase is directed to the array dof_indices_interleaved with the index row_starts[cell_index*n_vectorization*n_components].first.


This value indicates that the indices within a cell are all contiguous, and one can get the index to the cell by reading that single value for each of the cells in the cell batch. For a cell/face of this index type, the data access in FEEvaluationBase is directed to the array dof_indices_contiguous with the index cell_index*n_vectorization*n_components.


This value indicates that the indices with a cell are contiguous and interleaved for vectorization, i.e., the first DoF index on a cell to the four or eight cells in the vectorization batch come first, than the second DoF index, and so on. Furthermore, the interleaving between cells implies that only the batches for vectorization can be accessed efficiently, whereas there is a strided access for getting only some of the entries.

The two additional categories interleaved_contiguous_strided and interleaved_contiguous_mixed_strides are a consequence of this storage type. The former is for faces where at least one of the two adjacent sides will break with the interleaved storage. We then have to make a strided access as described in the next category. The last category interleaved_contiguous_mixed_strides appears in the ghost layer, see the more detailed description of that category below. Again, this is something that cannot be avoided in general once we interleave the indices between cells.

For a cell/face of this index type, the data access in FEEvaluationBase is directed to the array dof_indices_contiguous with the index cell_index*n_vectorization*n_components.


Similar to interleaved_contiguous storage, but for the case when the interleaved indices are only contiguous within the degrees of freedom, but not also over the components of a vectorized array. This happens typically on faces with DG where the cells have interleaved_contiguous storage but the faces' numbering is not the same as the cell's numbering. For a cell/face of this index type, the data access in FEEvaluationBase is directed to the array dof_indices_contiguous with the index cell_index*n_vectorization*n_components.


Similar to interleaved_contiguous_separate storage, but for the case when the interleaved indices are not n_vectorization apart. This happens typically within the ghost layer of DG where the remote owner has applied an interleaved storage and the current processor only sees some of the cells. For a cell/face of this index type, the data access in FEEvaluationBase is directed to the array dof_indices_contiguous with the index cell_index*n_vectorization*n_components, including the array dof_indices_interleave_strides for the information about the actual stride.

Definition at line 237 of file dof_info.h.

Enum used to distinguish the data arrays for the vectorization type in cells and faces.


The data index for the faces designated as interior


The data index for the faces designated as exterior


The data index for the cells

Definition at line 324 of file dof_info.h.

An enum to identify various types of cells and faces. The most general type is what we typically compute in the FEValues context but for many geometries we can save significant storage.


The cell or face is Cartesian.


The cell or face can be described with an affine mapping.


The face is flat, i.e., the normal factor on a face is the same on all quadrature points. This type is not assigned for cells.


There is no special information available for compressing the representation of the object under consideration.

Definition at line 50 of file mapping_info.h.

template<int dim, typename Number = double>
enum MatrixFree::AdditionalData::TasksParallelScheme

Collects options for task parallelism. See the documentation of the member variable MatrixFree::AdditionalData::tasks_parallel_scheme for a thorough description.


Perform application in serial.


Partition the cells into two levels and afterwards form chunks.


Partition on the global level and color cells within the partitions.


Use the traditional coloring algorithm: this is like TasksParallelScheme::partition_color, but only uses one partition.

Definition at line 179 of file matrix_free.h.

template<int dim, typename Number = double>
enum MatrixFree::DataAccessOnFaces

This class defines the type of data access for face integrals in loop () that is passed on to the update_ghost_values and compress functions of the parallel vectors, with the purpose of being able to reduce the amount of data that must be exchanged. The data exchange is a real bottleneck in particular for high-degree DG methods, therefore a more restrictive way of exchange is clearly beneficial. Note that this selection applies to FEFaceEvaluation objects assigned to the exterior side of cells accessing FaceToCellTopology::exterior_cells only; all interior objects are available in any case.


The loop does not involve any FEFaceEvaluation access into neighbors, as is the case with only boundary integrals (but no interior face integrals) or when doing mass matrices in a MatrixFree::cell_loop() like setup.


The loop does only involve FEFaceEvaluation access into neighbors by function values, such as FEFaceEvaluation::gather_evaluate(src, true, false), but no access to shape function derivatives (which typically need to access more data). For FiniteElement types where only some of the shape functions have support on a face, such as an FE_DGQ element with Lagrange polynomials with nodes on the element surface, the data exchange is reduced from (k+1)^dim to (k+1)^(dim-1).


The loop does involve FEFaceEvaluation access into neighbors by function values and gradients, but no second derivatives, such as FEFaceEvaluation::gather_evaluate(src, true, true). For FiniteElement types where only some of the shape functions have non-zero value and first derivative on a face, such as an FE_DGQHermite element, the data exchange is reduced, e.g. from (k+1)^dim to 2(k+1)^(dim-1). Note that for bases that do not have this special property, the full neighboring data is sent anyway.


General setup where the user does not want to make a restriction. This is typically more expensive than the other options, but also the most conservative one because the full data of elements behind the faces to be computed locally will be exchanged.

Definition at line 635 of file matrix_free.h.

template<int dim, typename Number = double>
enum MatrixFree::DoFHandlers::ActiveDoFHandler

Use DoFHandler.


Use hp::DoFHandler.

Definition at line 1657 of file matrix_free.h.

An enum that encodes the type of element detected during initialization. FEEvaluation will select the most efficient algorithm based on the given element type.


Tensor product shape function where the shape value evaluation in the quadrature point corresponds to the identity operation and no interpolation needs to be performed (collocation approach, also called spectral evaluation). This is for example the case for an element with nodes in the Gauss-Lobatto support points and integration in the Gauss-Lobatto quadrature points of the same order.


Symmetric tensor product shape functions fulfilling a Hermite identity with values and first derivatives zero at the element end points in 1D.


Usual tensor product shape functions whose shape values and quadrature points are symmetric about the midpoint of the unit interval 0.5


Tensor product shape functions without further particular properties


Polynomials of complete degree rather than tensor degree which can be described by a truncated tensor product


Tensor product shape functions that are symmetric about the midpoint of the unit interval 0.5 that additionally add a constant shape function according to FE_Q_DG0.

Definition at line 41 of file shape_info.h.