FDTD Voxels-in-Cell Method With Debye Media

The Voxels-in-Cell (VIC) method was recently introduced to reduce the computational cost of the finite-difference time-domain (FDTD) method with objects composed of dielectric voxels. It relies on using an FDTD cell larger than the voxels, with eight or more voxels in each VIC cell. With the objective of using it in bio-electromagnetics applications, this article extends the VIC method to voxels filled with Debye media. Besides the theory and the algorithm of the extended VIC method, several numerical experiments are reported with a canonical object and with human body phantoms composed of voxels. The experiments show that the accuracy of the method is preserved while large reductions of the computational requirements can be achieved, especially the computational time can be reduced by about one order of magnitude.



Abstract-The Voxels-in-Cell (VIC) method was recently introduced for reducing the computational cost of the finitedifference time-domain (FDTD) method with objects composed with dielectric voxels.It relies on using a FDTD cell larger than the voxels, with eight or more voxels in each VIC cell.With the objective of using it in bio-electromagnetics applications, this paper extends the VIC method to voxels filled with Debye media.Beside the theory and the algorithm of the extended VIC method, several numerical experiments are reported with a canonical object and with human body phantoms composed with voxels.The experiments show that the accuracy of the method is preserved while large reductions of the computational requirements can be achieved, especially the computational time can be reduced by about one order of magnitude.

I. INTRODUCTION
he Voxels-in-cell (VIC) method was introduced in [1] with the objective of reducing the computational cost of FDTD simulations, when the object of interest is given as a set of voxels.The idea behind the method consists in using a FDTD cell larger than the voxels, with for instance eight voxels in each cell.The framework of the VIC method [1] is summarized in Fig. 1.As compared to using cells equal in size to the voxels size, the overall number of cells is widely reduced while the FDTD time step can be enlarged, resulting in a large reduction of the computational burden.
Initially presented with pure dielectric objects [1], the VIC method is now extended to dispersive Debye media.This allows the method to be used in bio-electromagnetics, where the human phantoms are composed with voxels assumed as Debye media [2].Such an achievement was not possible with methods [3,4,5,6], since they are limited to a single interface splitting the FDTD cell into two parts, one filled with vacuum or dielectric, the other with a dispersive medium.Inversely, Fig. 1.The summary of the VIC method, as with dielectric media in [1].With Debye media, in addition to E field computed in the slices and H field in the cells, the electric flux density D is computed in each voxel.(a) The Yee cell as usually defined, with E nodes on its edges, is filled with voxels.(b) The Yee cell shifted with x/2 and y/2 with respect to (a), with Ez node at its center, denoted as Ez(ig, jg, kg).(c) The Ez field in each slice, denoted as Ez(kv), is uniform, but discontinuous at the interfaces between slices, thus it is updated in each slice using the Maxwell-Ampere law that yields (12) in Debye media.From which Dz in each voxel of the slices is updated using equation (10).(d) Ez(ig, jg, kg) at the center of the cell in (b) or (c), to be used later for updating the neighbouring H components, is updated using equation (13) and Ez(kv) in the slices.Note that there is no need of computing Dz at the Ez (ig, jg, kg) node.The other electric field components, Ex and Ey, are treated using the same logic.
the VIC method can handle the several non-parallel interfaces present in the FDTD cells when they are filled with voxels.
The modifications of the theory for dealing with Debye media are described and experiments are reported with a canonical object and with the head of a human phantom, in both cases with VIC cells two times larger than the voxels (eight voxels in each VIC cell).The experiments show that the accuracy of the method is roughly similar to that with dielectric media [1].A critical issue is addressed, namely the reduction of the computational requirements by using the VIC method in place of a fine grid with cells size equal to the voxels size.It is shown that reductions of the computational time of one order of magnitude can be achieved.Finally, simulations with two whole human phantoms illustrate what can be done with the VIC method using a personal computer.

B.1. A homogeneous FDTD cell:
Let us assume that the electric flux density Dz is centred in the cell of sizes x, y, z.(same arrangement as in Fig. 1 (b), with Dz in place of Ez).The magnetic field nodes Hy-and Hy+ are situated x/2 from Dz, whereas the nodes Hx-and Hx+ are y/2 away from Dz.With the FDTD method, these components are used to update Dz by discretizing the integral form of the Maxwell-Ampere equation where t is the time, x y is the surface perpendicular to Dz and Ch is the integral where L is the contour path that surrounds surface x y and holds the four H nodes. Discretization of ( ) at time n+1/2 reads Assuming that the flux Dz is homogeneous in the cell, discretizing the derivative on time in (1) yields where ∆ / is the difference between Dz values at times n and n+1 The electric flux Dz and field Ez share the same FDTD nodes and the relationship that connects them for the one-pole Debye model reads where is the optical permittivity is the static permittivity, is the vacuum permittivity, is the static conductivity, is the relaxation time and is the angular frequency.Component Ez at time n+1 from (6) can be obtained in function of Ez and ∆Dz at previous time steps by using the discretized auxiliary equation similar to [[7], eq.( 16)].Using (5), flux Dz can be eliminated from the auxiliary equation that becomes We note that we have replaced the update of flux Dz with the update of its increment Dz over the time step, inversely to what is done in [7].This permits the needed storage to be reduced, only Dz at one time step instead of Dz at two time steps.In addition, the algorithm is slightly more rapid.In human body applications, in general flux Dz is not needed, in case it is needed at some nodes, it can be obtained from Dz.The same substitution of Dz with Dz is used in the following with non-homogeneous cells which also reduces the needed memory.

B.2. A non-homogeneous FDTD cell:
Let us now address the case in Fig. 1 where the cell is non-homogeneous and split into Mi Mj Mk voxels, each with its own medium and indexes (iv, jv, kv), where iv ϵ (1, Mi), jv ϵ (1, Mj), kv ϵ (1, Mk).In a given slice, the component Ez is tangential to the interfaces between the voxels of the slice, and thus it is continuous through those interfaces.Since the cell size in FDTD is small compared to the wavelength, we can assume that Ez is uniform in the slice, i.e. it is the same in all the voxels of the slice.Inversely, Dz is not continuous at the interfaces, so that each voxel of the slice has its own Dz.The Maxwell-Ampere (1) law holds everywhere in the physical medium, including the slices composed of voxels.This allows us to derive the FDTD update of Ez and Dz as with a homogeneous cell, with only slightly more complicated mathematics, but without any further assumption on the Physics.For slice kv, (1) yields the discrete integral at time t that can be discretized as where ∆ / ( , , ) is the increment of Dz(iv, jv, kv) in each voxel, similar to (5).In each voxel the same auxiliary equation (7) as in a homogeneous cell holds, which provides us with Mi Mj equations that can be rewritten as where coefficients , , , and , , depend on the physical parameters in voxel ( , , ).Equations ( 9) and (10) form a set of MiMj +1 equations in slice kv for the unknown ( ) and the MiMj unknowns ∆ / ( , , ).The set can be easily solved as the set of three equations in [ [7], eq.( 17)] by using the MiMj equations (10), so as to replace the Mi Mj unknowns ∆ / ( , , ) in (9).This yields The sole unknown in (11), ( ), can then be obtained as with / from (3) and Once ( ) has been computed, the MiMj unknowns ∆ / ( , , ) can be obtained by using the auxiliary equation (10) that holds in each voxel.
After updating Ez(kv) and its x and y counterparts in all the slices of all the FDTD cells, the H components should be updated to complete the FDTD iteration.This can be done as with dielectric voxels in [1].The standard FDTD algorithm can be left unchanged, we just need to place the average value < > over z at each original Ez , , node of the FDTD cells, that is Finally, with cells filled with voxels, the FDTD advance of Ez and ∆Dz components can be summarized as follows 1. Compute ( ) in each slice using (12).2. Compute ∆ ( , , ) in each voxel using (10).3. Compute the average of ( ) in the cell (13), to be used later for the advance of H components with the standard FDTD equation.
It can be shown that the above VIC method reduces to the previously presented method [8] in the special case of the transverse-magnetic (TM) polarization in two dimensional problems.

III. THE METHODS WITH WHICH VIC IS COMPARED
In the next sections, the VIC calculations are compared with similar calculations as in [1]: 1/ a reference calculation using a FDTD fine grid, where the cell size is identical to the voxels size x, y, z, with each cell filled with one voxel.The E nodes are located on the edges of the fine cells, that are also the edges of the voxels, and the updates of the E components comply with the Maxwell-Ampere law.With dielectric media [1], this just consists in averaging the permittivity's in the four voxels surrounding the considered E node.With Debye media, the E field remains the same in the volume of the cell where the E node is centered, because it is continuous at the interfaces, but the D fields are different in the four media of that volume composed with four quarters of voxels.The problem is the same as in a slice of the VIC cells.This means that we can get one equation like (9) and four equations like (10).From which E can be updated with an equation like (12) and then Dz with four equations like (10).This solution is rigorous.A simpler implementation of Debye voxels has been sometime used with FDTD, it consists in shifting by half a cell the voxels in two directions, so as the E field to be updated is at the center of one voxel.Then, ( 4) and ( 7) can be applied.However, the shifts depend on the component, due to the FDTD staggered grid, which means that the medium is not rigorously the same for the three E components.This method is only an approximation.We performed some comparisons of the two methods.A significant difference was observed at some nodes, in the proximity of high contrast at interfaces or corners.In this paper, all the reported fine grid results were computed with the rigorous method.
2/ a calculation denoted as AVG, using a coarse cell that has the same sizes and location as the VIC cell.The medium parameters used to update E equal the average values in the volume with dimensions of x, y, z and whose center is located at the E node.This reads, for component Ez(ig, jg, kg) in Fig. 1 ꭙ , , = where ꭙ is either , , , or .The medium is thus assumed as homogeneous in the AVG cell, so that both E and D can be updated using ( 4) and (7).

IV. EXPERIMENTS WITH A COMPLEX CANONICAL OBJECT
In this section, we report experiments with the 60x60x60-VIC-cell canonical object mimicking the topology of the human body.Here, the three dielectric media in [1] are replaced with three human body media, skin, fat and bone.The parameters of the media are shown in Fig. 2. The computational settings are the same as in [1], i.e. voxel size 1 mm 3 , VIC cell size (2 mm) 3 , VIC time step 3.66 ps, the object surrounded with a vacuum and a PML, and the incident plane wave propagating in y direction and polarized in z direction, with waveform where = 134.5 ps.The situations of the corners of the object are defined as in [1] where they are represented and numbered in the 2D case in Fig. 6 and numbered in the 3D case in Fig. 16.In 3D the corner represented in Fig. 16, in grey, is the situation 1, while for instance situation 3 in 3D corresponds to a corner situated at (i0-1/2, j0+1/2, k0), i.e. shifted to the right with y/2 with respect with the corner in situation 1. Figs 3 reports the Ey and Ez fields at the inner corner denoted as Co-3 (Fig. 2) which is a point the three media have in common, while Fig. 4 shows E fields at the inner edge Ed-2 between skin and fat and Ed-3 between fat and bone.The exact locations A, C, R, S, R', S', T of the plotted E are the same as in ([1], Figs.16 and 22).The same situations of the object are plotted, so that we can compare Figs. 3 and 4  Results in Fig. 3 are roughly similar to the ones in ( [1], Fig. 27).However, the magnitude of the fields, especially the peak values, are lower with the Debye media than with the pure dielectric media.And the oscillations of the fields are significantly reduced, which is due to the presence of loss terms with the Debye media, while the pure dielectric media are lossless.Similar reductions of the magnitudes are observed in the vicinity of the edges (Fig. 4), but with little reduction of the oscillations of the fields.not singular (node T in Fig. 4), is superimposed to the fine grid component, as with pure dielectric media.
The observations on Figs. 3 and 4 are general.At any outputted node the agreement of VIC with fine grid calculation is similar to that observed with pure dielectric media, in both cases there is a significant difference only for the components that are singular at corners or edges.Concerning the magnitude of the field within the object, it depends on the situation of the observation.From the outer surface of the object to its interior, the field decreases faster than in the pure dielectric case, because of the propagation in the lossy Debye media.Fig. 5 reports the average of the VIC and AVG errors on the 60x60x60-cell object surrounded with a layer of vacuum 3 cells in thickness, computed as in [1].Comparing Fig. 5 with ([1], Fig. 30), the errors are quite similar at early time, but at late time the errors in Fig. 5 are significantly smaller because of the losses in the human body media.

V. SIMULATIONS WITH A HUMAN PHANTOM HEAD
We performed simulations with the head of the Voxel-Phantom provided by IT'IS Foundation (Zurich, Switzerland) [9].The voxel size is 1 mm 3 .The simulation scenario is depicted in Fig. 6.In both the fine grid and the VIC grid, we imported a part of the Phantom composed with 192x252x224 voxels, that can be viewed as composed with 224 slices perpendicular to direction z and one voxel in thickness (Fig. 6).Each slice has 192x252 voxels in (x, y) directions.Note that the voxels of each slice are filled with either head media or vacuum.The fine grid was a FDTD domain of 240x300x272 cells, surrounded with a 12-cell PML.The head (Fig. 6) was imported in the centre of the domain, which means that there were a 24-cell layer of vacuum between the slices in Fig. 6 and the PML.In the VIC grid, the cell was (2 mm) 3 in size.The imported head (Fig. 6) was composed with VIC cells, each one with eight voxels inside, which means 96 VIC cells in x direction, 126 in y direction, and 112 in z direction.The VIC counterpart of the 192x252x224 voxels or fine grid cells was thus a VIC grid of size 96x126x112.It was surrounded with a layer of vacuum 12 cells in thickness, having the same physical thickness as the 24-cell layer of vacuum in the fine grid.For the convenience of encoding, the 12-cell layer of vacuum was the addition of a 2-cell VIC layer and a 10-cell layer of usual FDTD cells, in the manner of ([1], Fig. 4).A 12-cell PML was placed outside the vacuum.In total, the overall FDTD domains were 264x324x296 for the fine grid calculation and 144x174x160 for the VIC and AVG calculations.The Debye parameters for all human tissues in the phantom are presented in [10].The incident plane wave and the FDTD steps were the same as in Sec.IV.Inversely to the previous experiments, we did not vary the situation of the phantom in the VIC grid, we considered only the situation resulting from the import of the 192x252x224-voxel part of the phantom that just fits a 96x126x112-VIC grid.

Fig. 6:
The 3D human head model was imported from the human phantom [9] in the form of 224 slices of 192x252 voxels.Each tissue was associated with its own one-pole Debye parameters [10].The observation points P1, P2, and P3, were located behind the left eyeball, in the cranial portion of skull, and in the white matter, respectively.Fig. 7 reports the three field components at their nodes closest to points P1, P2, and P3 in Fig. 6.We can see that the VIC results agree very well with the reference fine grid solutions, at P1, P2, and P3.Inversely, the AVG results depart from the reference solutions at P1 and P2, especially at P1 where Ex and Ez are strongly erroneous.The interpretation of the results is complex in such a strongly non-homogeneous structure, where errors propagating from neighbouring nodes may be added to the local error.However, beside the excellent agreement of VIC results with the references, the discordance of AVG results in Fig. 7 seems in accordance with the situations of locations P1, P2, and P3.At P1, the eight voxels that surround the plotted Ez VIC node (see the situations of the components w.r. the voxels in [1], Fig. 15) are filled with different media, some with fat, some with muscle, whereas the voxels that surround Ex and Ey are all filled with fat.From this, Ez is discontinuous in the volume composed with the eight voxels [1], which results in the large error of AVG in Fig. 7a.At P2, this is Ey that is discontinuous, the surrounding voxels are filled either with skull or cerebrospinal fluid, whereas Ex and Ez VIC nodes are in homogeneous skull medium.The lower AVG errors in Fig. 7b are, at least in part, due to the lower contrast of the media at P2.At P1 the permittivity's of fat and muscle are about 5.5 and 57, respectively, while at P2 the permittivity's of skull and cerebrospinal Fluid are about 14 and 70.At P3, VIC, AVG, and fine grid results are almost superimposed.This is because P3 is in homogeneous region (white matter) from which no local error is produced by AVG which is equivalent to VIC.Fig. 8 reports the average of the VIC and AVG errors on the phantom head surrounded with 2 VIC cells of vacuum, computed as in [1], except that the nodes where |E| is lower than 1 V/m are excluded.The errors are smaller than with the canonical object in Fig. 5, because 1/the object is larger so that the part of the object where the incident field is small is larger, and 2/the propagation of the field in the lossy Debye media is longer.The ratio of AVG error to VIC error is slightly smaller than in Fig. 5 which may be, at least in part, due to the lower contrast between media in the head which may reduce the AVG error.As in the dielectric case [1], the VIC method is stable, we have not found any instability with more than 5 000 VIC time steps in all the simulations.

VI. THE OPTIMIZED ALGORITHMS FOR REDUCTION OF THE CPU TIME
With dielectric media and VIC cells two times larger than the voxels, the CPU time is in theory reduced with a factor of 16, as compared to a simulation using a cell equal in size to the voxels [1].With Debye media, things are different, because of the more complex relationship between E and D vectors, which results in more complex updates.More specifically, E must be updated in each slice with (12), while D must be updated in each voxel with (10), which means 10 update equations per cell.From this, there is no such simple theoretical reduction of CPU time as with dielectric media.Estimates of the ratio of CPU times by counting the numbers of operations suggest that the reduction of CPU times is lower than with dielectric media.This has been confirmed by experiments.In the above as well as in [1], it is implicitly assumed that all the VIC cells surrounding components Ex, Ey, Ez, are nonhomogeneous, i.e. as if all the eight voxels were different.Actually, this is not true in most applications, especially in bioelectromagnetics, where the phantoms are composed with homogenous organs, separated with interfaces.So that parts of the VIC cells are homogeneous, with eight identical voxels, and the remaining VIC cells are non-homogeneous, usually with less than eight different voxels.As an example, in the phantom head in Fig. 6, approximately 60% of VIC cells in the head are homogeneous as indicated in Table I, where, for instance, VIC-Ex denotes the VIC cell having Ex at its centre.This suggests that the CPU time could be reduced by using simplified updates in some VIC cells instead of the general  12) and (10), especially in the homogeneous ones where the update reduces to (4) and (7), which means 2 update equations instead of 10.We refer to the general case were all the eight voxels are different as the "brute force" update.The reduction of CPU time can be achieved gradually by initially focusing on the update process for the simplest VIC cells.To begin with, assuming the brute force update is already implemented, the first step is to incorporate the treatment of homogeneous VIC cells.Assuming that each VIC cell, e.g.Ez VIC cell, has been provided with a flag before loop on time which equals 1 (iflag=1) if the VIC cell is homogeneous, then the block diagram of the update of the VIC cells would be as in Fig. 9.With just this simple addition of treatment of homogeneous cells, the CPU time for updating the Debye object will be widely reduced, since about 60% of the cells in the head will have iflag=1 and thus will be updated with ( 4) and (7).Concerning the VIC cells outside the Debye object and filled with vacuum, if any, as the ones in the slices of the imported head (Fig. 6), they can be updated as usual FDTD cells in vacuum, so that one can consider that a branch with for instance iflag=0 could be added to the diagram in Fig. 9 to update these cells.The next step in the optimization process involves addressing VIC cells that have at least one homogeneous slice.This comprises three cases, 1/the two slices are homogeneous but with different media, 2/slice 1 is homogeneous while slice 2 is non-homogeneous, and 3/the inverse case.Table II provides with the number of non-homogeneous cells in the three cases, for the same head phantom as in Table I.We can introduce additional flags, before the time loop, corresponding to the three cases.Specifically, we assign "11" if both slices are homogeneous but with different media, "12" if the first slice is homogeneous and the second slice is nonhomogeneous, "21" if the first slice is non-homogeneous and the second slice is homogeneous.This optimization is visualized using a block diagram as Fig. 10, similar to Fig. 9. From Table II, about 60 % of the non-homogeneous cells in Table I belong to one of the three cases where at least one slice is homogeneous, which means about 25 % of all the VIC cells where equation ( 12) and 4 equations ( 10) can be replaced with equations ( 4) and ( 7) at least in one slice.The remaining cells to be updated with brute force only represent about 17% of the total number of VIC cells (about 100 000 out of 580 000).Table III reports the CPU times for 100 VIC iterations (11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz), with the head phantom and brute force, algorithm 1, and algorithm 2 calculations, along with the CPU times for 200 iterations of the two fine grid calculations (section III).All the CPU times are the ones devoted to the Debye cells, excluding the surrounding cells of vacuum in Fig. 6.With the rigorous fine grid method, equation (12) and four equations (10) are used to update E and the four D's in the non-homogeneous cells, while (4) and (7) are used in homogeneous cells.This is the same strategy as with VIC i.e. an algorithm like Algorithm 1 of VIC can be applied.With the approximate fine grid method, all the cells are homogeneous, and updates ( 4) and ( 7) are used everywhere.III (ratios of VIC times over fine grid times).With respect with the rigorous fine grid method, most of the reduction of CPU times is achieved with the brute force updates, further reductions with algorithms 1 and 2 are relatively small.This is because we applied to the fine grid the equivalent of algorithm 1 which reduces the CPU time (81.3 s v.s.102.6).Inversely, with the approximate fine grid method, no reduction of CPU time is possible, so that when using algorithms 1 or 2 with VIC, the difference with the brute force VIC is larger.
Finally, with the head phantom, the reduction of CPU time is smaller than with dielectric media [1], but it remains large, close to an order of magnitude, especially with the refinements of algorithms 1 or 2. This was also observed with the canonical Debye object in section IV with which the reduction factors were higher than with the head phantom because the proportion of homogeneous VIC cells was larger.
Concerning the memory requirements, the VIC method needs storage of the 3 H components, the 3 E components in the two slices and at two time steps (12), and the 3 D components in each voxel.For a number of VIC cells Nvic, this results in 39 Nvic real numbers.Similarly, for a number of fine grid cells Nfine = 8 Nvic, the rigorous fine grid method needs storage of 3 H, 6 E, 12 D, in total 21 Nfine = 168 Nvic.For the approximate fine grid method, 12 Nfine = 96 Nvic.From this, the VIC method memory needs are reduced with factors about 4.3 and 2.4 with respect with the fine grid methods.In actual computations, additional arrays may be used to store coefficients of the update equations before the loop on time.Especially c1, c2, c3 in (12) for the 3 E components, in 2 slices, which needs 18 Nvic memories, to be added to the 39 Nvic for storing the field components.Even with in addition storage of some other data, as the flags used with algorithms 1 and 2, the memory needs of the VIC method are always lower than those of a fine grid calculation, at least a factor of two with the rigorous one.As an example, with the head phantom, whose reduction of CPU times are given in Table IV, the memory needs were 1.1 GB with VIC, 3.1 GB with the rigorous fine grid, and 1.8 GB with the approximate fine grid.

VII. VIC-SIMULATIONS WITH WHOLE HUMAN BODY PHANTOMS
This section illustrates what can be done using the VIC method on a personal computer (PC).We present the simulation results with two full-body phantoms.Firstly, with the full body of the head phantom that we used in section V, referred to as the "male phantom".Secondly, with the "female phantom" in [9].The male phantom consists of 610x310x1840 voxels, while the female phantom is composed of 530x300x1650 voxels, in both cases with a voxel size of 1 mm 3 .The incident wave is the same as in Sec.V. Fig. 11 presents the simulation results at two locations, within the skin of the outer abdominal region, and the heart muscle.Note that the plotted E fields significantly depend on the considered phantom, an influential parameter is probably the size.The sizes of the FDTD domains of the VIC computations are presented in Table V.They include a 12-cell layer of vacuum surrounding the phantoms and a 12-cell PML.As can be seen, both the memory needs and the CPU times are compatible with use of a PC, even without parallel computing.Using the same PC, the corresponding fine grid calculation would need a memory larger than 64 GB, which is usually not possible (most PCs are limited to 64 GB).

VIII. CONCLUSION
The Voxels-in-cell (VIC) method presented in [1] in the case of pure dielectric objects, has been extended to dispersive Debye media, in view of applications to bioelectromagnetics.Several experiments have shown that the accuracy and the limitations of the method are roughly similar to the ones observed with dielectric objects.The VIC method is by far more accurate than the averaging method, a significant error occurs only in the close vicinity of the singularities of the E field, which is inherent to the larger size of the VIC cell.The limitations are obviously the same as in [1], mainly the upper frequency cut-off which is lowered with the ratio of the VIC cell to the voxels size.
The reduction of the computational requirements is lesser than with dielectric objects, due to the complex constitutive relation of the Debye media, but it remains high, especially for the CPU time which is reduced by about one order of magnitude with VIC cells two times the voxels size.The VIC method is stable with Debye media, its stability condition remains the CFL condition of the FDTD method, as with dielectric media in [1].Finally, the VIC method could also be applied to such other dispersive media as the Drude or Lorentz media using their own auxiliary equations in place of the Debye media one.

Fig. 3 :
Fig. 3: Results at the corner Co-3.Comparison of fine grid E with VIC E (a) at node A and (b) at node C. The errors at both cases are so small in all the situations that any situation can be considered as almost exact.

Fig. 4 :
Fig. 4: Results with situation 4, (a) at edge Ed-2, where node T results are multiplied with 0.25 and (b) at edge Ed-3, where node T results are multiplied with 0.5.Concerning the difference between the VIC results and the fine grid results, i.e. the VIC error, in Figs.3 and 4it is quite small as in the pure dielectric cases ([1], Figs.27 and 28).Especially, the VIC component tangential to the edge, that is

Fig. 5 :
Fig. 5: The overall average error of the electric field Ez in situations 5 and 8 of the object (the largest errors of VIC).It includes all nodes in the object and 3 VIC cells of its surrounding vacuum.

Fig. 7 :
Fig. 7: Results with Fine Grid, VIC and AVG methods (a) at P1, (b) at P2, and (c) at P3, within the head of the phantom.

Fig. 8 .
Fig. 8.The overall average error of the electric field components.It includes all the nodes within the head of the phantom and 2 VIC cells of its surrounding vacuum, but the nodes where |E| < 1.0 V/m were excluded.
This article has been accepted for publication in IEEE Transactions on Antennas and Propagation.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TAP.2024.3378847This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/updates (

Fig. 9 .
Fig. 9. Algorithm 1: simplified updates in the homogeneous VIC cells and brute force update in the inhomogeneous ones, using just one "if" tests.

Fig. 10 :
Fig. 10: Algorithm 2: with treatment of VIC cells with at least one homogeneous slice.

Fig. 11 :
Fig. 11: The observation of Ez within skin and heart muscle of human phantoms.
This article has been accepted for publication in IEEE Transactions on Antennas and Propagation.This is the author's version which has not been fully edited and This article has been accepted for publication in IEEE Transactions on Antennas and Propagation.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TAP.2024.3378847 This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE II NUMBER
OF VIC-E CELLS HAVING AT LEAST ONE HOMOGENEOUS SLICE IN THE PHANTOM HEAD

TABLE III CPU
TIME (s) WITH VIC AND FINE GRID FOR THE PHANTOM HEAD

TABLE IV REDUCTION
FACTOR OF CPU TIME ACHIEVED WITH VIC FOR THE This article has been accepted for publication in IEEE Transactions on Antennas and Propagation.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TAP.2024.3378847Thiswork is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/approximatefine grid (62.5 s).Algorithm 1 further reduces the CPU time significantly (8.6 s), but Algorithm 2 achieves little additional reduction (8.3 s).Table IV presents the reduction factors of CPU times deduced from Table One can see in TableIIIthat even with the brute force VIC updates, the CPU time of VIC is widely smaller than with the rigorous fine grid updates (11.9 s instead of 102.6) or with the

TABLE V MEMORY
AND CPU TIME MEASUREMENTS FOR THE HUMAN BODIES