When the main magnetic field in a MRI system, $$$B_0$$$, increases, Larmor frequency increases and the corresponding wavelength decreases. At small wavelengths, standing wave patterns are formed in human body under scanning. There are peaks and nulls in the energy/field distribution inside the body, which causes safety issues and deteriorates imaging accuracy, respectively. Thus, a quick and accurate electromagnetic simulation of the human body is crucial for predicting the temperature and specific absorption rate (SAR) distribution before a scan. Intensive attention has been paid recently on developing calculation tools to fulfill such requirement [1-6]. Although the finite difference time domain (FDTD) method [6] is popular in the MRI society, it is more promising to use volume integral equation (VIE) method rather than partial differential equation method for such purpose because no artificial boundary condition is introduced and many mature IE based fast techniques are available. Several fast solvers based on the method of moments (MoM) have been developed for MRI applications [2-5]. However, these works are based on the traditional strong form VIE, which usually require massive computational time and memory. We develop herein a more efficient solver which is based on the weak-form VIE [7] and accelerated by fast Fourier transform (FFT) method.When waves impinge an object, the total electric field $$$\vec{E}^{tot}\left(\vec{r}\right)$$$ in the domain of interest is governed by VIE and can be obtained by the superposition of the impinging waves $$$\vec{E}^{inc}\left(\vec{r}\right)$$$ and the scattered field. This relation is expressed in (1), where $$$\chi\left(\vec{r}\right)=1-\epsilon_0/\epsilon\left(\vec{r}\right)$$$ is the normalized contrast function and $$$\epsilon\left(\vec{r}\right)$$$ is the complex permittivity of the dielectric object. The object in our study is the highly inhomogeneous human body. It is discretized into small voxels, each of which has constant permittivity and conductivity. By introducing the electric flux density $$$\vec{D}$$$ and the electric contrast vector potential $$$\vec{A}$$$ expressed in (2) , the VIE formulation can be further expressed by (3). VIE is traditionally discretized in terms of $$$\vec{D}$$$, which leads to a hyper-singularity. Therefore, in conventional fast solvers, the interactions are categorized into near and far interaction groups based on the distance, and the near interactions are conducted by rigorous MoM. By doing this, the calculation efficiency is considerably sacrificed. To achieve an optimal efficiency, the weak form method originally proposed in [7] is applied. We expand both $$$\vec{D}$$$ and $$$\vec{A}$$$ as in equation (4) and (5) by using three dimensional roof-top basis functions. Testing (3) by using Galerkin method results in (6), where $$$M_1$$$ and $$$M_2$$$ are sparse matrices. The convolution integration (2) that determines the relationship between unknown coefficient vector $$$\left[a\right]$$$ and $$$\left[d\right]$$$ can be speeded up by the 3D-FFT method. The Green's function in the convolution kernel can be weakened by taking its spherical mean value. In this way, no singularity appears in the numerical implementation of the interactions, which greatly reduces both the implementation complexity and the time and memory required. The stable and robust iterative solver GMRES (30) is applied to solve the numerical matrix system. Due to the fact that VIE is a second kind Fredholm IE, diagonal block preconditioner is adopted here to improve the convergence. $$\vec{E}^{tot}\left(\vec{r}\right)=\vec{E}^{inc}\left(\vec{r}\right)-k_0^2\int_Vg\left(\vec{r},\vec{r}'\right)\cdot\epsilon\left(\vec{r}'\right)\chi\left(\vec{r}'\right)\vec{E}\left(\vec{r}'\right)d\vec{r}'\quad(1)$$$$\vec{A}(\vec{r})=\frac{1}{\epsilon_0}\int_Vg(\vec{r},\vec{r}')\chi(\vec{r}')\vec{D}(\vec{r}')d\vec{r}'\quad(2)$$$$\vec{E}^{inc}\left(\vec{r}\right)=\frac{\vec{D}\left(\vec{r}\right)}{\epsilon\left(\vec{r}\right)}-\left(k_0^2+\nabla\nabla\cdot\right)\vec{A}\left(\vec{r}\right)\quad(3)$$$$\vec{D}\left(\vec{r}\right)=\epsilon_0\sum{d_j\vec{f}_j\left(\vec{r}\right)}\quad(4)$$$$\vec{A}\left(\vec{r}\right)=\sum{a_j\vec{f}_j\left(\vec{r}\right)}\quad(5)$$$$\left[e\right]=M_1\cdot\left[a\right]+M_2\cdot\left[d\right]\quad(6)$$To verify the proposed application of the weak-form VIE method in MRI applications, we calculated the near scattered electric field (E-field) of a dielectric sphere with radius r = 84 mm, relative permittivity ε_r = 49 and conductivity σ = 0.6 S/m. Results are compared with that of analytic Mie series method. Good agreement can be seen from Fig. 1. We further calculate the E-field of a highly inhomogeneous human head of the Duke body model from the Virtual Family [8]. The head is embedded in a domain of 252×288×294 mm³ which is uniformly discretized into 126×144×147 voxels (resulting in 8,059,338 degrees of freedom for VIE methods). The model is excited by a plane wave that is linearly polarized in the x-direction, traveling in the negative z-direction and operating at 298.2 MHz (7 T). By the proposed method, the whole system is solved using 7.15 GB memory within merely 7.11 minutes. While, with the strong-form VIE method, 56.7 GB memory and more than 6 hours are required. About 5.5 hours are needed for the FDTD-based commercial software SEMCAD to achieve a stable result. As shown in Fig. 2, good agreement can also be observed in the comparison between the proposed method and SEMCAD. We can conclude that the FFT-accelerated weak-form VIE solver can provide a powerful and promising full-wave electromagnetic analysis for extremly inhomogeneous and high contrast problems for MRI applications.