Material Coarsening Strategy for Structured Meshless Multigrid Method for Dosimetry in Anisotropic Human Body Models

—In this article, a novel multigrid-based method is developed for modeling of electric ﬁelds and currents in electrically anisotropic and heterogeneous media. The method is useful for numerical assessment of human exposure to low-frequency electromagnetic ﬁelds, which has been typically performed using isotropic models of the human body. The method is matrix-free and based on structured voxelized grids. Numerical tests indicate that the method has all typical beneﬁts of multigrid methods, i.e., convergence in a constant number of iterations and linear complexity in terms of the number of degrees of freedom. The developed method is applied to the dosimetry of the induced electric ﬁeld in ten anatomically realistic models of the human head exposed to uniform 50-Hz magnetic ﬁelds. The models are constructed from magnetic resonance images and account for the electrical anisotropy of the brain tissues via diffusion weighted images. The dosimetrymodelingshowsthatthetissueanisotropyleadstoasmallbutsigniﬁcantincrease(14%onaverage)ofthe99thpercentile inducedelectricﬁeldstrengthinthewhitematterofthebrain.However,isotropicmodelscanstillprovidesufﬁcientaccuracy fordosimetryifthisslightunderestimationofthe99thpercentileelectricﬁeldstrengthisaccountedfor.

on Electromagnetic Safety (ICES) [3] have established safety guidelines and standards to limit human exposure to EMF [4].In the low-frequency range (< 10 MHz), the guidelines/standards set restrictions on the induced electric field within the central and peripheral nervous systems (CNS and PNS).Assessment of the induced electric field requires detailed dosimetric evaluation using either numerical and/or measurement methods [4].To simplify the exposure assessment, the ICNIRP and IEEE-ICES have developed exposure reference levels, which are defined in terms of the external EMF strength.The ICNIRP reference levels have been derived based on numerical dosimetry using anatomically realistic models of the human body [1].However, the reference levels have a large conservative reduction factor (3) due to dosimetric uncertainty.
One of the most important sources of uncertainty is the electrical characterization of body tissues, which has been mentioned in both the IEEE-ICES research agenda [5] and the ICNIRP data gaps document [6].Until now, numerical dosimetry of human low-frequency EMF exposure has been performed in anatomical models where the conductivity of tissues has been assumed electrically isotropic.However, substantial volumes of electrically anisotropic tissues occupy the human body.In an average human being, the most prevalent tissue in the body is the skeletal muscle [7], which is anisotropic.In the CNS, which is relevant for human exposure to utility-frequency magnetic fields (50/60 Hz), the largest induced electric fields are found in the strongly anisotropic white matter of the brain [8] in two-thirds of the population [9].
Thus, the tissue anisotropy may substantially affect the induced electric fields.Characterization of the uncertainty caused by anisotropy could therefore improve the accuracy of dosimetry modeling and, if the effect is found to be small, reduce the size of dosimetric uncertainty factors in the international human exposure guidelines.
Previous studies addressing the effect of tissue anisotropy have been limited to diagonal-matrix approximations of the anisotropy tensor.Early findings indicated that the maxima of the induced current density and/or electric field may differ by more than 10% compared to isotropic models [10], [11].A more recent study found up to 40% error due to isotropic conductivity model for local magnetic-field exposure of the legs, also using diagonal approximation of the anisotropy tensor [12].
Existing methods that are applicable to the problem either require the assembly of the sparse system matrix [13], [14] This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/This article has been accepted for inclusion in a future issue of this journal.Content is final as presented, with the exception of pagination.
allowing the use of fully anisotropic material parameters and highly optimized sparse linear solvers but suffer from the high cost of actual matrix assembly, or they support only diagonal anisotropy but do not require the assembly of the system matrix [10], [15], [16].The drawback of assembling the sparse matrix is that its computational complexity is O(N ), where N is the number of degrees of freedom.This may add significantly to the total simulation time and memory consumption if the linear system is resolved also in O(N ) using, e.g., multigrid methods [17], [18].
Therefore, there is a great demand for a computational dosimetry method that is matrix-free, converges fast, and can handle anisotropic material parameters.
In this work, we extend an efficient multigrid-based method [15] aimed to solve current-quasistatic problems to handle full anisotropies.Key features of the method are that it is based on the finite-element method (FEM) and that it utilizes a regular grid of hexahedral elements.The latter feature makes the method directly applicable to voxel-based anatomically realistic models that are commonly used for dosimetry of human EMF exposure [4].In addition, in such models, anisotropic conductivity data can be readily generated using diffusion weighted magnetic resonance imaging.The rest of this article is organized as follows.Section II describes the method.We validate the method with a manufactured solution and compare its performance against commercially available finite element software (COMSOL Multiphysics 6.0) in Section III.We then apply the method for the dosimetry of magnetic field exposure in ten anatomically realistic head models and compare the induced electric fields in the CNS calculated using isotropic and anisotropic conductivities.Finally, Section IV concludes this article.

II. METHOD
The proposed method solves for the electric scalar potential ϕ in the current-quasistatic problem [19] ⎧ where J is the current density, E is the electric field, and ˙ A is the time-derivative of a given vector potential, which is assumed to separate into product of spatial and temporal parts.If the temporal part is of the form e i2πf t , where f is the frequency, we replace ˙ A with i2πf A and factor out the temporal part so that A depends only on spatial coordinates.The anisotropic electric conductivity σ is described by a symmetric positive definite matrix-valued function ( In this work, we restrict to real valued conductivity data and so we omit the imaginary unit everywhere for clarity.Equation ( 1) is discretized using FEM.We compute the elements of the associated stiffness matrix on the fly using precomputed analytical expressions.They are linear maps of local end for conductivity values because the mesh is a voxelization of the computational domain.The finite element system is solved using a geometric multigrid method with successive overrelaxation (SOR) as pre and postsmoothers.In the multigrid method, the computational mesh is split into nested meshes of various levels of geometric coarseness.In this case, each coarser level element is split to exactly 2 × 2 × 2 elements to make up the finer level elements.Next, we shall detail on the finite-element discretization in anisotropic media and describe a novel anisotropic material coarsening approach for the multigrid method.

A. Solution Using the FEM and Successive Over-Relaxation
The finite element approximation φ h ∈ S h of ϕ in (1) satisfies where D ⊂ R 3 is the voxelized computational domain and S h is the finite element space of trilinear first-order interpolating nodal basis functions.In the following, we assume that D is a cube for brevity, but that can be relaxed so that it is a rectangular cuboid whose side lengths are integer multiples of its voxel side lengths.Enumerating the basis functions ψ klm by their grid node indexes klm and pqr, 1 ≤ k, l, m, p, q, r ≤ N + 1, where N is the number of elements in one coordinate direction, we denote the stiffness matrix by A h , the right-hand side (RHS) vector by y and the solution coefficient vector by In this notation, the SOR iteration with overrelaxation parameter (ω) used as the pre and postsmoothers in the multigrid method is listed in Algorithm 1.
We observed during numerical testing that this algorithm can be parallelized over the index k with static scheduling multithreading as long as the number of threads stay moderate with respect to N .In the SOR algorithm, we are faced with computing the 27 row values A h klm,pqr given the elementwise conductivity data σ.It amounts to calculating linear combinations the six unique conductivity tensor entries of each eight voxels that intersect at the node klm making the operation that takes conductivity data to row values a 27 × 48 matrix.These row values are listed for completeness in Appendix A, ( 12), ( 14), (15), and (17).Such a finite element-based stencil has not been previously documented in literature to our knowledge.Furthermore, the RHS values y klm have a simple expression when calculated with a one point quadrature and they are also listed in Appendix A, (16).
Repeated application of this algorithm eventually converges to an approximate solution of (3).However, the convergence is impractically slow, which is why such iterators are typically used in combination with other methods, such as the multigrid method.

B. Anisotropic Material Coarsening for the Multigrid Method
In classical geometric multigrid, the coarse grid stiffness matrix A 2h is obtained from the fine-grid stiffness matrix A h by multiplication with the prolongation matrix and its transpose [18], but here we approximate A 2h using a material coarsening strategy.For isotropic conductivity, this amounts to averaging the conductivity over the fine grid voxels [15], but the equivalent of that for anisotropic conductivity is more involved and has not yet been discussed in literature.
In the following, we shall outline the multigrid method in order to justify the material coarsening strategy.To that end, let us denote the bilinear form appearing on the left-hand side of (3) by a(v, u; σ) = D ∇v • σ∇udx and the linear form on the RHS by f (v; σ) = − D ∇v • σ Adx.Furthermore, let us denote the space of fine grid finite element functions with S h and that of the coarse grid functions with S 2h .
Let now φ 0 h be initial approximation of φ h in S h and let φ 1 h = S n ω φ 0 h be the result of n-times application of the smoothing operator (Algorithm 1).The coarse grid correction to φ 1 h is denoted by φ 2h ∈ S 2h , and it is the solution to The coarse grid corrected solution denoted by with m iterations of the smoother.Typically, the coarse grid correction is obtained recursively with another application of the procedure outlined above until the coarsest grid is reached.For example, if there are levels of grid resolutions, the coarsest grid is S 2 −1 h .In multigrid literature, such an iteration is called a v-cycle.Similarly to [15], we obtain the coarsest grid correction with multiple applications of the smoother (n ≈ 100) and the solution to (3) is sought using multiple v-cycles by setting φ 0 h ← φ 3 h until a desired tolerance is achieved.In terms of coefficient vectors, the coarse grid correction φ 2h is p T x 2h , where x 2h solves and p is the matrix representation of the embedding S 2h → S h .The action p(y − A h x 1 h ) is the straightforward tabulation of this matrix representation and for p T x 2h we chose to use integral averages, as discussed in [15] due to the code being readily available.
The operation of coarse-grid matrix A 2h = pA h p T is more involved because the stiffness matrices are not constructed at any point.Also the calculation of the coarse grid matrix entries on the fly as double matrix products would not be feasible because that would void the benefit of solving smaller and smaller systems at every level in multigrid method altogether.
Instead, we construct the coarse grid matrix on the fly using coarsened material data σ 2h , which is anterpolated from the fine grid material σ h in the following manner.Let Q 2h be a coarse grid cube and L h (Q 2h ) be the set of nodal interpolating basis functions of S h restricted to Q 2h .Then the coarse grid conductivity associated to Q 2h is the minimizer (7) Because the minimization problem is a quadratic one, given 8 × 6 values of dense grid conductivities σ h Q 2h , we find the six independent coarse grid anisotropy matrix entries as a linear combination of the dense grid conductivities.The corresponding matrix is reported in Appendix B, (18).This material coarsening strategy has not been documented in literature to our knowledge.

III. NUMERICAL TESTS
We first validate the method and study its convergence properties with a manufactured solution, where the RHS of the equation is set so that the solution results in a known function.We also compare the performance of the method to COMSOL Multiphysics weak form partial differential equation (PDE) interface which assembles the sparse system matrix.Finally, we showcase the utility of the method in dosimetric computation with ten human head models exposed to uniform magnetic fields.In the following, we use the fractional anisotropy (FA) as the measure of the amount of anisotropy.The FA for tensor σ is defined by where λ i , 1 ≤ i ≤ 3, are the eigenvalues of σ and λ is their arithmetic mean.FA of 0 means that all eigenvalues are equal, i.e., the material is isotropic, and a value of 1 that only one of the eigenvalues is nonzero.
All numerical tests were performed on Intel Core i7-6850 K CPU at 3.60 GHz using six threads via OpenMP compiler pragmas.The performance critical codes were compiled with GCC version 11.3.0.

A. Validation Using a Manufactured Solution
The manufactured solution is defined in an origin-centered cube having side length of L = 200 mm by thus we set A = ∇ϕ manu in (3) and solve for φ h .We set the conductivity tensor to be where u is a unit vector angling xand z-axes at π/5 and π/3 radians, respectively.In all cases we choose σ 2 = 1 and σ 1 so that the FA is within the range of values appearing in the white matter of the human brain (0.03-0.94, see Section III-B).
This article has been accepted for inclusion in a future issue of this journal.Content is final as presented, with the exception of pagination.

TABLE I CONVERGENCE TO MANUFACTURED SOLUTION AT DIFFERENT HOMOGENEOUS FRACTIONAL ANISOTROPIES
Fig. 1.Convergence of manufactured solution in terms of V-cycle iterations at different FAs and mesh sizes (h).
The relaxation parameter ω for the SOR iterator was set to 1.2 and it was parallelized with six threads over index k.Generally 0.8 ≤ ω ≤ 1.4 appeared to work fine.The multigrid algorithm was applied over log 2 N − 1 levels, where N is the number of elements in a single coordinate direction.The coarse grid solution was obtained with 100 SOR iterations.The number of pre and postsmoother iterations were 3.
Table I and Fig. 1 report the results of the numerical experiments.The results verify the validity of the proposed method, and we can also conclude that the method converges in O(1) Fig. 2. Convergence time to tolerance (10 −9 ) of manufactured solution with respect to mesh size (h).Comparison to COMSOL Multiphysics 6.0 is made with isotropic homogeneous material and multigrid mesh coarsening strategy using SOR as pre and postsmoothers.Most of the time in COMSOL is spent on assembling the system matrix and setting up the prolongations.COMSOL converges in 8 V-cycles.
iterations in terms of degrees of freedom, and that anisotropy has no significant effect on the performance of the method when the FA values are in the biological range.Fig. 2 shows how long it takes for the method to reach the tolerance in wall-clock time including a comparison to similar problem calculated with COMSOL.The COMSOL model problem is isotropic making it less involved to solve iteratively but still the upfront cost of traversing the mesh and assembling the sparse linear system is significant.

B. Dosimetry in Anatomical Head Models Exposed to a Magnetic Field
The method was applied to calculate the electric field induced in the human head due to exposure to 50-Hz homogeneous 1-mT magnetic field directed in three orthogonal directions: anterior-posterior (AP), lateral-medial (LAT), and superior-inferior (TOP) [20].It is worth noting that the magnetic flux density was equal to the reference level of ICNIRP for occupational exposure [1].
Ten high-resolution (1 mm) anisotropic models of the head were constructed from magnetic resonance images obtained from healthy volunteers (ages: 25-38 years, 6 male, 4 female).The volunteers gave an informed consent prior to participating in the experiment.The study protocol was approved by the Aalto University Research Ethics committee (decision D/574/03.04/2022).All methods were carried out in accordance with approved institutional guidelines and regulations.
T1, T2, and diffusion weighted (DW) images of the volunteers were acquired using a 3-T magnetic resonance imaging (MRI) scanner (Magnetom Skyra; Siemens, Ltd., Erlangen, Germany).T1 and T2 weighted images were processed using a previously described procedure [21] to generate head models that consisted of multiple tissues and bodily fluids.Isotropic conductivity values of tissues and bodily fluids were obtained from [22] and were identical to those listed in our previous study [9].For anisotropic tissues, the isotropic conductivity values were converted into anisotropy tensors, as described in the following.
The DW-MRI data was acquired with 100 gradient orientations using 900, 1600, and 2500 s/mm 2 diffusion-weightings at 2 × 2 × 2 mm voxel size.Seventeen additional images without diffusion-weightings were acquired out of which five were in reverse-phase encoding direction.The diffusion MR images This article has been accepted for inclusion in a future issue of this journal.Content is final as presented, with the exception of pagination.were corrected for subject motion [23], eddy current induced distortions [24], and echo-planar imaging distortions [25] using the FIMRIB software Library [26].Anisotropic first-order diffusion tensor estimate was obtained from the DW-MRI data with two iterations of ordinary least squares method [27] implemented in the MRtrix3 software framework [28].
Fig. 3 illustrates the FA and the direction of anisotropy calculated from the diffusion tensor estimate D in a representative head model.
The volume normalized anisotropic conductivity σ was obtained from tissue specific isotropic conductivity σ iso and D using a volume normalization technique [29] as follows: The anisotropy was only applied to the parts of the anatomical model that are known to be significantly anisotropic, i.e., white matter of the cerebrum and cerebellum, and brainstem.The conductivity tensor elsewhere in the head was assumed to be isotropic.
For the calculations using the multigrid method, the finest level grid resolution was 1 mm resulting in approximately 10 7 elements, including air.Excluding the nodes in air, the number of degrees of freedom was on average 4.8 • 10 6 (standard deviation: 4 • 10 5 ).The calculations were done using six threads and they converged to the tolerance of 10 −6 on average in 18.4 (standard deviation: 1.7) v-cycles for isotropic models and 18.7 (standard deviation: 1.7) v-cycles for anisotropic models.
The resulting electric fields were spatially averaged over 2 × 2 × 2 mm cubes and the 99th percentile values of the magnitudes were obtained according to the ICNIRP guidelines [1].
The 99th percentile electric field values in various brain tissues are displayed in Fig. 4. Anisotropy increased the 99th  percentile electric field strength in anisotropic tissues (white matter and brainstem), but its impact on the electric field in other parts of the brain was negligible.Overall, the largest electric fields were induced in white matter, and also the enhancement of the field strength was the largest in white matter (sample mean and standard deviation: 14 ± 5.7%, range: 0.1-20.9%).Despite the increase compared to the case of isotropic conductivity, the highest induced electric field strengths remained well below the ICNIRP basic restriction for occupational exposure, which is 100 mV/m [1].
Fig. 5 illustrates the difference between the induced electric field calculated using isotropic and anisotropic conductivity values on axial slices of the head models.It can be observed that the largest relative differences are found in white matter, as expected.It is worth noting that the largest local differences can be much larger than the differences in the 99th percentile values.While the local values are not relevant for dosimetry of human exposure, they can be more important, e.g., for modeling of noninvasive brain stimulation of specific brain regions [14].
The convergence curves showing the number of v-cycle iterations and wall-clock time spent to obtain to reach the tolerance are shown in Fig. 6.

IV. CONCLUSION
In this work, we developed a matrix-free multigrid solver based on SOR smoother suitable for LF dosimetry and verified it with synthetic tests and realistic dosimetric calculations.The numerical tests show that the method exhibits correct approximation properties and that it reaches the tolerance in O(1) iterations, which translates to O(N dof ) complexity in terms of degrees of freedom.
The implementation of the method is straightforward as there is no need to form the stiffness matrix nor the prolongation operators rendering its setup time negligible.The anisotropic conductivity data consumes six times the memory compared to isotropic solver, but that cannot be avoided because the conductivity data must be anyhow stored.However, the memory requirements are 6/27 compared to solvers that explicitly construct the stiffness matrix (27 nonzero elements on each row).The reduced memory usage can be important for dosimetry in high-resolution (< 1 mm) whole-body human models.
The coarsening matrix in (18) can also be used to anterpolate isotropic conductivities down to a coarser grid still preserving some of the isotropic conductivity features at material boundaries.Anterpolated anisotropic conductivity data generated using such an approach might have favorable properties in terms of the staircasing effect that occurs when curved boundaries are approximated with rectilinear elements [30], [31].In fact, a recent study has suggested an approach based on effective direction-dependent conductivity values at tissue boundaries to reduce the staircase approximation error [32].Dosimetry modeling using realistic head models suggested that tissue anisotropy led to a systematic enhancement of the 99th percentile electric field strength especially in the white matter of the brain.The effect was small, on average 14%, compared to the dosimetric uncertainty factor of three assumed in the ICNIRP guidelines [1].However, it was still relatively large compared to the normal interindividual variability in the 99th percentile induced electric field that has been estimated using isotropic models [9].The isotropic modeling study found that the 99th percentile values in the brain varied remarkably little between individuals, the relative standard deviation being approximately 10% of the mean [9].Another important observation was that the enhancement effect due to anisotropy was systematic, meaning that isotropic models can still provide sufficient accuracy for dosimetry, if the slight underestimation of the 99th percentile electric field strength is accounted for.
The dosimetric analysis reported in this article focused on the induced fields in the central nervous system.Our recent modeling work using realistic models of human legs [12] suggested that the difference between isotropic and anisotropic models can be much larger (up to 40% on average) in peripheral nervous system than the differences found herein for white matter.However, diagonal anisotropy was used for the previous investigation [12].Therefore, to fully quantify the effect of anisotropy in dosimetry, further research is needed using whole-body or partial-body models with anisotropic skeletal muscles.
In conclusion, we developed and validated an effective approach for modeling induced electric fields in voxel-based anisotropic human body models.The dosimetric data produced using the method provide useful information on the effects of anisotropy on the induced electric fields within the human central nervous system and can aid the development of human EMF exposure guidelines and standards.

A. Stiffness Matrix Row Values and RHS Entries
The stiffness matrix row values are best calculated with a computer algebra system, but for completeness we shall report the terms in the chosen local enumeration of nodes and elements.The local matrix indexes are related to global indexes by Fig. 7(b).A single row in the stiffness matrix requires information of eight conductivity tensors of the voxels intersecting the corresponding node.The local to global mapping of these voxels is shown in Fig. 7(a), where the element number is associated to subscript in σ j,k i in ( 12)- (15).The diagonal entry is given by the expression This article has been accepted for inclusion in a future issue of this journal.Content is final as presented, with the exception of pagination.where the sign coefficients s j,k i are given by vectors All the subsequent terms also have the common factor h/18, which shall be omitted in the formulas for brevity.
The face terms (a 2 . . .a 7 ) are given by (w/o the h/18 term, s The edge and the corner terms (a 8 . . .a 19 , a 19 . . .a 27 ) listed in (17), shown at the bottom of this page and ( 15) are shorter but appear to lack clear simplified expressions in terms of chosen enumeration (again the term h/18 is omitted).
The RHS value obtained with one point quadrature rule is given by ) This article has been accepted for inclusion in a future issue of this journal.Content is final as presented, with the exception of pagination.
where r i is the centroid of the subcube i, A i is the value of the vector potential at it, and sign is the sign function applied componentwise.

Fig. 3 .
Fig. 3. Fractional anisotropy (and corresponding uniaxial anisotropy ratio) and the eigenvector corresponding to the largest eigenvalue of the diffusion tensor estimate at each location in a representative head model.The eigenvector is orthographically projected to the axial plane.

Fig. 4 .
Fig. 4. Sample means (± sample standard deviation) of the 99th percentile electric field strengths (mV/m per 1 mT) in various brain tissues (N = 10).The three magnetic field directions are illustrated on the bottom left.Nuclei tissue includes deep grey matter regions of the cerebrum.

Fig.
Fig. Induced electric field magnitudes in isotropic and anisotropic head models due to 50-Hz TOP-oriented 1-mT magnetic fields and the relative difference between anisotropic and isotropic models.

Fig. 6 .
Fig. 6.Residual with respect to V-cycles and wall-clock time in the dosimetry simulations.Median residual in shown in saturated color on left.

Fig. 7 .
Fig. 7. Enumeration of local to global interactions and local subvoxels.Local subvoxel numbering scheme associated to a node (a).Local stiffness matrix row numbering scheme (b).