Finite element discretization

In SpinDoctor, the finite element mesh generation is performed using an external package called Tetgen in 3D [2], and Triangle in 2D [19]. Each finite element mesh consists of

• a list of nodes in three dimensions:   ;

• a list of tetrahedral elements (  indices referencing the nodes).

The list of nodes includes double nodes that are placed at the interfaces between compartments connected by permeable membranes. This allows for representing discontinuous magnetization fields and (or and ) on the same boundary . To distinguish between the different compartments, let   with    for  . The set contains the indices of the nodes representing compartment , including interface nodes. In the adjacent compartments, the corresponding interface nodes will have different indices, distinct from .

In SpinDoctor, the finite element space is the space of compartment-wise continuous piecewise linear functions on tetrahedral elements in three dimensions. This space has a set of basis functions whose number is exactly the number of finite element nodes (including double nodes), and that are defined on the entire domain :

Let the finite element nodes be denoted by . The basis function ,  , is a piece-wise linear function, non-zero on the tetrahedra of that touch the node , and zero on all other tetrahedra (including tetrahedra of other compartments different than that do touch ). At the interface between two compartments, the value of is set to be the value it has inside its own compartment, distinct from that of the adjacent compartment. On a tetrahedron of that touches , is equal to on and it is equal to on the other 3 vertices of the tetrahedron. This completely describes the piece-wise linear function. The index sets may then be defined by    , the set of indices of the finite element nodal functions whose supports lie entirely within .

Any function in the finite element space can be written as a linear combination of the above basis functions:

where   is the vector of coefficients and   is the vector of finite element nodal basis functions.

To discretize the Bloch-Torrey and Laplace operators, we construct the following finite element matrices:  , known in the FEM literature as the mass, stiffness, and flux matrices, respectively. These matrices matrices we need are defined as follows, for  :

the latter being defined as the sum of interface integrals:

We remind the reader that and are compartment indices, while and are finite element nodal indices. Note that the above sum formulation counts each interface twice, but only considers the contribution to one side of the boundary at the time. It also correctly accounts for corner nodes (if any) that belong to two different permeable boundaries at the same time.

The weak form of the Bloch-Torrey PDE also requires computing the first order moments of the product of pairs of finite element basis functions. We let these three matrices be denoted by , and :

where  . For a given gradient vector  , we define

In addition, we define the finite element relaxation matrix   given by

In SpinDoctor, these matrices are assembled from local element matrices and the assembly process is based on vectorized routines from [20], which replace expensive loops over elements by operations with 3-dimensional arrays. All local element matrices in the assembly of , , and are evaluated at the same time and stored in a full matrix of size   , where denotes the number of tetrahedral elements.

The matrices , , and are assembled as coordinate weighted mass matrices, where the three coordinate functions  , , act as nodal weights in the assembly process, given by , , ; the vectors of , , and coordinates of the finite element nodes.

With compartment-wise constant -relaxation times, is a block-diagonally scaled version of the mass matrix , where the weights are the inverses of the relaxation times.

Finite element solution to the Bloch-Torrey PDE

The solution to the Bloch-Torrey partial differential equation may be projected onto the finite element nodal basis , in which case it is given by

where    is the vector of finite element basis functions and the function     is the solution to the following system of ordinary differential equations (ODE):

Finite element solution to the HADC model

Similarly, given an isolated compartment and a scalar gradient  , the solution to the HADC model, , may be obtained by solving the equation

where   ,     is the unknown vector of coefficients of the solution  ,   is the vector of finite element basis functions, and   is given by

where   is the outwards unit normal. This three-column matrix represents the components of the boundary integral of the quantity , where the constant diffusivity and gradient sequence dependent part has been factored out. They can thus be assembled independently of the gradient sequence, and be reused when solving for multiple sequences or directions. They are assembled using the same routine as for .

The boundary integral is then given by

where   is computed using the mass matrix.

Finite element Matrix Formalism signal representation

Both of the above ODEs are of dimension , the number of finite element nodes. The finite element discretized Matrix Formalism representation uses a different approach, limiting the problem size to , the number of Laplace eigenfunctions (of which the choice is further explored in section. But first we need to solve an eigenvalue problem involving matrices of size  . The finite element discretization changes the continuous Laplace operator eigenvalue problem into the following discrete, generalized matrix eigenvalue problem: find    such that

of which we will retain the smallest eigenvalues and corresponding eigenvectors , with  . Note however that there are in total solutions to the problem. Moving back to the space of functions (the function space ), the eigenfunction associated to the eigenvalue is then

where the entries of the eigenvector are the coefficients of the eigenfunction in the finite element basis. Using matrix notation, this conversion can also be written  , where  ,  , and   . The integrals of the finite element discretized eigenfunctions are then given by    , where   . Similarly, the coefficients of the initial spin density in the finite element discretized eigenfunction basis are given by    , where    and   is such that  .

Then it is clear that the first order moments of the product of pairs of Laplace eigenfunctions can be written as:

These three matrices are computed using a total of six matrix-matrix multiplications. Similarly, the eigenfunction basis relaxation matrix \mathbf{T} may be obtained from the finite element nodal basis relaxation matrix \mathbf{R} by