SUGISAKAMASANORI

Journal of Advances in Artificial Life Robotics
Vol. 1(1); June (2020), pp. 27–34
ON LINE ISSN 2435-8061; ISSN-L 2435-8061
https://alife-robotics.org/jallr.html

Research Article

Parallel Full-Wave Electromagnetic
Anatomical Human Body Models

Field

AnalysisUsing

Amane Takei, Makoto Sakamoto

Department of Electrical and System Engineering, University of Miyazaki, 1-1 Gakuen Kibanadai-Nishi, Miyazaki, 889-2192, Japan

ARTICLE INFO
Article History
Received 24 November 2019
Accepted 10 July 2020
Keywords
Iterative domain decomposition method
Parallel finite element method
Berenger’s PML
Large-scale analysis
Microwave analysis

ABSTRACTThe authors have investigated a high-accuracy analysis for the electromagnetic field of anatomical
human body models using the parallel finite element method. In the full-wave electromagnetic
field analysis by the parallel finite element analysis with anatomical human models, efficientconversion from voxels to finite element meshes and speed-up of an analysis become important
topics. In this paper, we propose parallel mesh generation method with mesh smoothing, and show
to be able to generate large-scale mesh very effective. Then, we evaluated the performance of two
kinds of iterative method by numerical experiments using large-scale human body model. As theresults, the COCR method is very useful as the interface solver in the DDM for high-frequency
electromagnetic field analysis using anatomical human model.

© 2022 The Author. Published by Sugisaka Masanori at ALife Robotics Corporation Ltd.
This is an open access article distributed under the CC BY-NC 4.0 license(http://creativecommons.org/licenses/by-nc/4.0/).

complex symmetric, and we have employed the
conjugate orthogonal conjugate gradient (COCG)
method for the interface problem. However, because of
the ill-conditioned coefficient matrix of the large-scale
analysis, it suffers from low convergence rate. To solve
this issue, we propose a domain decomposition method(DDM) algorithm based on the conjugate orthogonal
conjugate residual (COCR) method [4]. As shown in thispaper, in order to verify the accuracy of the developed
analytical solver, we calculate TEAM Workshop #29,
which is a benchmark problem, and evaluate the accuracyof the analytical solver. Moreover, an improvement of theperformance by using COCR method with huge-scale
anatomical human model that has 1,700 million meshes
is revolutionary in huge-scale analyses.

1. Introduction
Recently, a numerical electromagnetic analysis using
anatomical human model from CT/MRI is required
strongly[1],[2]. Important topics in numericalelectromagnetic field analysis using numerical humanbody model are efficient conversion from voxels to finite
element meshes and speed-up of analyses.
In our research group, we are investigating anddeveloping a full-wave electromagnetic analysis method
for large-scale analysis that has several hundred million
meshes based on the parallel finite element method [3].
To apply the parallel finite element analysis code in thecomputation using the anatomical human model, aparallel mesh generation algorithm and a meshsmoothing function are verified in the viewpoint of
accuracy and performances using a super computer
system.
On the other hand, the high-frequency
electromagnetic problem with our formulation is

2. Formulations

Corresponding author’s E-mail: takei@cc.miyazaki-u.ac.jp, URL: www.miyazaki-u.ac.jp/english/27

2.1. Governing equations

where K denotes the coefficient matrix, u the unknown
vector, and f the known right-hand side vector. As shownin the following equation, the domain Ω is decomposed
into N pieces so that there is no overlap in the boundary
between subdomains, namely

Let Ω be a domain with a boundary 𝜕𝜕Ω. The vector
wave equations [5] that describe an electromagnetic fieldwith a single angular frequency ω (rad/s) are derived
from Maxwell’s equations containing the displacement
current. The vector wave equations describing an electricfield E (V/m) are given by Eqs. (1a)-(1c) for a current
density J (A/m2) with j as an imaginary unit:
1curl � curl 𝐄𝐄� − 𝜔𝜔2 𝜀𝜀𝐄𝐄 − 𝑗𝑗𝑗𝑗𝑗𝑗𝐄𝐄 = 𝑗𝑗𝑗𝑗 𝐉𝐉 in Ω ,𝜇𝜇

𝐄𝐄 × 𝐧𝐧 = 𝟎𝟎 on ∂Ω ,

Ω = � Ω(𝑖𝑖) .

After domain decomposition, (4) is rewritten as follows:

(1b)

(1)

⎡ 𝐾𝐾𝐼𝐼𝐼𝐼
00
⎢
⎢
⋮
⎢
0
⎢
⎢ (1)𝑇𝑇 (1)𝑇𝑇⎢𝑅𝑅𝐵𝐵 𝐾𝐾𝐼𝐼𝐼𝐼⎣

Permittivity and permeability are given by ε (F/m) andμ (H/m), respectively. By solving Eq. (1a) under the
essential boundary condition in (1b), the electric field E
is calculated.
Magnetic fields are calculated by a post-process using
Faraday's law of induction, which is expressed by Eq. (2).
(2)

Here, we describe the finite element discretization.
The electric field E is approximated with Nédélec
elements (edge elements). The finite element
approximation is performed as follows.
Find Eh such thatΩ= 𝑗𝑗𝑗𝑗 �Ω(3)(𝑖𝑖) (𝑖𝑖)Here, Jh is an electric current density approximated by
the conventional piecewise linear tetrahedral elements.

⋱
0

⋯

⋯

(𝑁𝑁)𝑇𝑇

⋯

(1)

(𝑖𝑖)

𝐾𝐾𝐼𝐼𝐼𝐼 𝑢𝑢𝐼𝐼 = 𝑓𝑓𝐼𝐼𝑅𝑅𝐵𝐵

𝐾𝐾𝐼𝐼𝐼𝐼 𝑅𝑅𝐵𝐵⋮⋮
(𝑁𝑁) (𝑁𝑁)
𝐾𝐾𝐼𝐼𝐼𝐼 𝑅𝑅𝐵𝐵

⎤
⎥
⎥
⎥
⎥
𝑁𝑁
⎥
(𝑖𝑖)𝑇𝑇 (𝑖𝑖) (𝑖𝑖)� 𝑅𝑅𝐵𝐵 𝐾𝐾𝐵𝐵𝐵𝐵 𝑅𝑅𝐵𝐵 ⎥⎦𝑖𝑖=1(𝑁𝑁)𝑇𝑇

𝐾𝐾𝐼𝐼𝐼𝐼(6)(𝑖𝑖)

3.1. Domain decomposition method𝑆𝑆

The domain decomposition method (DDM) is
introduced to high-frequency problems using the E
method. Let us put the finite element equations of (3) inmatrix form, as follows:

(𝑖𝑖)

(𝑖𝑖)

− 𝐾𝐾𝐼𝐼𝐼𝐼 𝑅𝑅𝐵𝐵 𝑢𝑢𝐵𝐵 , 𝑖𝑖 = 1, … , 𝑁𝑁 ,𝑁𝑁

3. Parallel algorithms for solving large linear
system𝐾𝐾𝐾𝐾 = 𝑓𝑓,⋱(1) (1)

0⋮
0(𝑁𝑁)
𝐾𝐾𝐼𝐼𝐼𝐼

where the subscripts I and B correspond to unknowns in
the interior of subdomains and on the interface boundary,(𝑖𝑖)
respectively. 𝑅𝑅𝐵𝐵 maps the global degrees of freedom(DOF) of the interface onto the local DOF of the
subdomain interface. Equation (6) leads to linear systems
as follows:� �1�μ�rot𝐄𝐄𝐡𝐡 ∙ rot𝐄𝐄𝐡𝐡∗ 𝑑𝑑𝑑𝑑 − 𝜔𝜔2 � 𝜀𝜀𝐄𝐄𝐡𝐡 ∙ 𝐄𝐄𝐡𝐡∗ 𝑑𝑑𝑑𝑑
𝐉𝐉 𝐡𝐡 ∙ 𝐄𝐄𝐡𝐡∗ 𝑑𝑑𝑑𝑑.

⋯(1)⎡ 𝑓𝑓𝐼𝐼 ⎤⎡ 𝑢𝑢𝐼𝐼 ⎤
⎢ ⋮ ⎥
⎢ ⋮ ⎥
⎢ ⋮ ⎥ =⎢ ⋮ ⎥⎢ (𝑁𝑁) ⎥⎢ (𝑁𝑁) ⎥⎢𝑢𝑢𝐼𝐼 ⎥
⎢𝑓𝑓𝐼𝐼 ⎥
⎣ 𝑢𝑢𝐵𝐵 ⎦⎣ 𝑓𝑓𝐵𝐵 ⎦2.2. Finite element formulation

Ω

(5)

𝑖𝑖=1

(1a)

rot𝐄𝐄 − 𝑗𝑗𝑗𝑗𝜇𝜇𝐇𝐇 = 𝟎𝟎.

𝑁𝑁

𝑆𝑆 𝑢𝑢𝐵𝐵 = 𝑔𝑔 ,(𝑖𝑖)𝑇𝑇

𝑖𝑖=1(𝑖𝑖)(𝑖𝑖)𝑇𝑇(𝑖𝑖) −1 (𝑖𝑖)
= 𝐾𝐾𝐵𝐵𝐵𝐵 − 𝐾𝐾𝐼𝐼𝐼𝐼 �𝐾𝐾𝐼𝐼𝐼𝐼 � 𝐾𝐾𝐼𝐼𝐼𝐼 ,𝑖𝑖=1(𝑖𝑖)𝑇𝑇(𝑖𝑖) −1 (𝑖𝑖)𝑔𝑔(𝑖𝑖) = 𝑓𝑓𝐵𝐵 − 𝐾𝐾𝐼𝐼𝐼𝐼 �𝐾𝐾𝐼𝐼𝐼𝐼 � 𝑓𝑓𝐼𝐼

(4)

(9)

(𝑖𝑖)𝑇𝑇𝑔𝑔 = � 𝑅𝑅𝐵𝐵 𝑔𝑔(𝑖𝑖) ,
(𝑖𝑖)

(8)

(𝑖𝑖)

𝑆𝑆 = � 𝑅𝑅𝐵𝐵 𝑆𝑆 (𝑖𝑖) 𝑅𝑅𝐵𝐵 ,
𝑁𝑁

(7)

.

(10)

We call (7) the subdomain problems, (8) the interface
problem, 𝑆𝑆 the Schur complement matrix. In this paper,
28

𝑢𝑢𝐵𝐵0 ∶ an initial guess
𝑟𝑟 0 = 𝑔𝑔 − 𝑆𝑆𝑆𝑆𝐵𝐵0
𝑡𝑡 0 = 𝑀𝑀−1 𝑟𝑟 0
𝑝𝑝−1 = 0, 𝑞𝑞−1 = 0, 𝛽𝛽 −1 = 0
𝑠𝑠 0 = 𝑆𝑆𝑡𝑡 0
𝑓𝑓𝑓𝑓𝑓𝑓 𝑛𝑛 = 0,1, …
𝑝𝑝𝑛𝑛 = 𝑡𝑡 𝑛𝑛 + 𝛽𝛽 𝑛𝑛−1 𝑝𝑝𝑛𝑛−1
𝑞𝑞𝑛𝑛 = 𝑆𝑆𝑝𝑝𝑛𝑛 = 𝑠𝑠 𝑛𝑛 + 𝛽𝛽 𝑛𝑛−1 𝑞𝑞𝑛𝑛−1
𝑤𝑤 𝑛𝑛 = 𝑀𝑀−1 𝑞𝑞𝑛𝑛
𝛼𝛼 𝑛𝑛 = (𝑡𝑡̅𝑛𝑛 , 𝑠𝑠 𝑛𝑛 )⁄(𝑞𝑞�𝑛𝑛 , 𝑤𝑤 𝑛𝑛 )𝑢𝑢𝐵𝐵𝑛𝑛+1 = 𝑢𝑢𝐵𝐵𝑛𝑛 + 𝛼𝛼 𝑛𝑛 𝑝𝑝𝑛𝑛
𝑡𝑡 𝑛𝑛+1 = 𝑡𝑡 𝑛𝑛 − 𝛼𝛼 𝑛𝑛 𝑤𝑤 𝑛𝑛𝑖𝑖𝑖𝑖 (‖𝑟𝑟 𝑛𝑛+1 ‖⁄‖𝑟𝑟 0 ‖ < 𝛿𝛿 ) 𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏;
𝑠𝑠 𝑛𝑛+1 = 𝑆𝑆𝑆𝑆 𝑛𝑛+1
𝛽𝛽 𝑛𝑛 = (𝑡𝑡̅𝑛𝑛+1 , 𝑠𝑠 𝑛𝑛+1 )⁄(𝑡𝑡̅𝑛𝑛 , 𝑠𝑠 𝑛𝑛 )(7) are solved by the direct method based on Gaussianelimination, on the other hand, (8) is solved by theiterative methods. Here, the Krylov subspace methods
require the coefficient matrix-vector multiplication ineach iterative procedure, however, such calculation for(8) can be replaced by likely solving subdomain
problems. Therefore, in our DDM algorithm, subdomain
problems are also solved at each iteration for solving theinterface problem.
3.2. COCG method

Fig. 2. Algorithm of preconditioned COCR method.

In this paper, we introduce two kinds of iterative
methods for the subdomain problems and the interface
problem. Firstly, Fig. 1 shows the preconditioned COCG
method. In Fig.1, 𝑀𝑀 is the preconditioning matrix,(𝑥𝑥, 𝑦𝑦) = ∑mi=1 𝑥𝑥̅𝑖𝑖 𝑦𝑦𝑖𝑖 an inner product of complex vectors 𝑥𝑥
and 𝑦𝑦, ‖∙‖ the Euclidean norm, and 𝛿𝛿 a positive constantas the convergence criterion.

4. Accuracy verification
4.1. TEAM29
We analyze TEAM Workshop #29[8] in order toverify the accuracy of the developed analytical solver.2
TEAM 29 is a benchmark problem and involves aresonator model.

𝑢𝑢𝐵𝐵0 ∶ an initial guess
𝑟𝑟 0 = 𝑔𝑔 − 𝑆𝑆𝑆𝑆𝐵𝐵00−1 0𝑧𝑧 = 𝑀𝑀 𝑟𝑟
𝑝𝑝0 = 𝑧𝑧 0
𝑓𝑓𝑓𝑓𝑓𝑓 𝑛𝑛 = 0,1, …
𝑞𝑞𝑛𝑛 = 𝑆𝑆𝑆𝑆𝑛𝑛
𝛼𝛼 𝑛𝑛 = (𝑟𝑟̅ 𝑛𝑛 , 𝑧𝑧 𝑛𝑛 )⁄(𝑝𝑝̅ 𝑛𝑛 , 𝑞𝑞𝑛𝑛 )𝑢𝑢𝐵𝐵𝑛𝑛+1 = 𝑢𝑢𝐵𝐵𝑛𝑛 + 𝛼𝛼 𝑛𝑛 𝑝𝑝𝑛𝑛
𝑟𝑟 𝑛𝑛+1 = 𝑟𝑟 𝑛𝑛 − 𝛼𝛼 𝑛𝑛 𝑞𝑞𝑛𝑛
𝑖𝑖𝑖𝑖 (‖𝑟𝑟 𝑛𝑛+1 ‖⁄‖𝑟𝑟 0 ‖ < 𝛿𝛿 ) 𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏;
𝑧𝑧 𝑛𝑛+1 = 𝑀𝑀−1 𝑟𝑟 𝑛𝑛+1
𝛽𝛽 𝑛𝑛 = (𝑟𝑟̅ 𝑛𝑛+1 , 𝑧𝑧 𝑛𝑛+1 )⁄(𝑟𝑟̅ 𝑛𝑛 , 𝑧𝑧 𝑛𝑛 )𝑝𝑝𝑛𝑛+1 = 𝑧𝑧 𝑛𝑛+1 + 𝛽𝛽 𝑛𝑛 𝑝𝑝𝑛𝑛

Fig. 1. Algorithm of preconditioned COCG method.
3.3. COCR method
The preconditioned COCR method [4] is shown in Fig.
2. [6] The COCR method requires more or less the sameamount of computational cost and working memory periteration as the COCG method, and is expected to getsmooth convergence behavior compared with the COCG
method. In particular, the COCR method for (8) has more
effect than the COCG method in a large-scale analysis.

Fig. 3. TEAM 29 cavity resonator model.
Table 1. TEAM 29 model data.
No. of Elements

DOFs

No. of Subdomains

121,277

149,668

100×7

The resonator is cylindrical and has a diameter of 1.9[m] and a height of 1.45 [m]. In the analysis, a dielectric
phantom having a relative permittivity εr of 80 and anelectric conductivity σ of 0.52 [S/m] is positioned, andthe resonance state is investigated. The analysis domain

29

Table 3. Resonance frequencies. (Units: [MHz]. The error rate[%] between measured data and the numerical solution is shown inparentheses.)boundary is a perfect conductor. The analysis model isshown in Fig. 3. Table 1 shows the specifications of the
TEAM 29 model. The highest calculation efficiency isachieved when the number of elements contained in one
subdomain is approximately 170, and the number ofpartial domains is determined such that the number of
elements contained in one subdomain is equal to 170.
4.2. Frequency response analysis
The results of the analysis of TEAM 29 are subjected
to frequency response analysis in order to confirm theaccuracy of the developed analytical solver. In order todetect the resonant frequency and compare solutions withactual measurements, the frequencies of some range are
analyzed. The frequency band of 60 [MHz] to 140 [MHz]is calculated for 2-MHz steps, and the response for every
frequency step is investigated. In addition, calculations
near the resonance frequency are performed in 0.4-MHz
intervals. The computing environment used in the present
study is a 25-PC cluster equipped with Intel Core i72600K multi-core CPUs (total: 100 cores) and 32 GB
memory is used. The compiler used is the gcc. In addition,
Message Passing Interface (MPI) is used for theparallelization library. The average calculation time perfrequency step and the averagely used memory are shownin Table 2. Fig. 4 shows the frequency response of themagnetic field. The measured and calculated values areshown in Table 3.

DOFs

No. of Subdomains

121,277

149,668

100×7

Measureddata

FDTD 25mm mesh

1st

68.6

67 (2.33)

2nd

110

110

3rd

134

134

Result65.6(4.37)
109.0(0.91)
134.4(0.30)As shown in Fig. 4, a solution resonating around theresonance frequency of the actually measured value isobtained. In the comparison of the measured andcalculated values, the error rate is 4.37 [%] in the 1stmode, 0.91 [%] in the 2nd mode, and 0.30 [%] in the 3rdmode. As the mode increases, the error rate decreases.
However, it is the same tendency as the analysis result
obtained by the FDTD method. Moreover, the error rateis less than 5 [%], and a solution with high accuracy is
obtained.
Therefore, the solution obtained by the developed
solver is proven to have a sufficiently high accuracy.
Moreover, in the analysis of the dipole antenna applying
the PML described in the following sections, the error
tolerance index is defined as 5 [%] in order to evaluate
the accuracy.
5. Pre-processing

Table 2. Specifications of the TEAM 29 model.
No. of Elements

Resonancemode

5.1. Mesh generation and smoothing
In many numerical human body models employ abinary data format wherein types of organs (including airarea) are encoded using voxels. In NICT model1, the sizeof the adult male model is 320 voxels wide, 160 voxels
deep, and 866 voxels high. Represented as the “char”
type, it occupies 44,339,200 bytes of memory. In thisresearch, each voxel is divided into five tetrahedra.
Fig. 5 shows a parallel procedure for generating thedomain decomposed mesh. Firstly, the voxel mesh of thenumerical human body model is read. Secondly, using
ParMETIS [7], the input voxel data is divided into a
number of parts. After part decomposition, all processes
can perform remaining procedures independently
without communication.

Fig. 4. Numerical and measured frequency response
analysis results.

30

frequency step is investigated. Computations are
performed using 1,152 cores (72 nodes) of the FX10
supercomputer. A comparison between measured
resonant frequencies [8] is shown in Table 4. Resonances
do not occur in the frequency band of 60–70 MHz using
the original voxel mesh and the one-level smoothing
mesh. On the other hand, the obtained solution is in very
good agreement with that obtained using the two-levelsmoothing. The maximum error rate between theobtained solution and the measurement is 6.10% in themode.

Fig. 5. Procedure for generation of
the domain decomposed mesh.

Table 4. Resonant frequencies in MHz.( ( ): error rate
between Measured Data (68.6 MHz) and Numerical Solutions (%))
Smoothing type
ResultNone (Original voxel model)
N/A
One-level smoothing
N/A
Two-level smoothing
64.4 (6.10)

A mesh smoothing technique based on the marching
cube-like algorithm is employed in part of thetransformation from voxels to tetrahedras of our meshgeneration procedure. Fig. 6 shows the schematic of themesh smoothing. The dotted line (A) in Fig. 6 shows theoriginal boundary between different materials. Firstly,
the boundary is smoothed by triangular prisms, as shownby the dashed line (B). This is referred to as one-levelsmoothing. Next, by searching for angles of 45 degrees,
the points on the angles are moved so that these angles
become 0 degrees, as shown by (C). The solid line (D)shows the smoothed boundary obtained by these twolevel smoothing operations.

Fig. 7 shows conditions of the analysis model. The
dipole antenna is set above the breast. The distance
between the breast and antenna H is 0.06 m. The antenna
length L is 0.5 m, which is same as half wavelength at
300 MHz. The lengths of the model’s sides are x = 0.64m, y = 0.32 m, and z = 1.732 m. [6]Fig. 7. Conditions of the analysis model.
Table 5 show the generation conditions and elapsed
times of the HDD mesh using 12-nodes (96-cores) on the
FX10. In this paper, full body models of tetrahedral
elements are generated from 4mm voxel models. Models“1, 2, 3, 4 and 5” shown in Table 5 with 27.71 million
tetrahedral elements are obtained from the 4mm voxel
model.
The mesh decomposition tool used in this research isdesigned to add the voxel partition function to
ADVENTURE Metis [9]. As shown in Table 5, thelongest calculation time is spent on the part andsubdomain decompositions. It is difficult to reduce thecalculation cost of the part and subdomain
decompositions with largest calculation cost in theFig. 6. Schematic of the mesh smoothing.
5.2. Performance evaluation
First, accuracy verification for mesh smoothing is
performed using the voxel mesh model of TEAM 29. An
accuracy verification is performed using this smoothedmodel. The processing time necessary to smoothing isless than 1 second. To detect the resonant frequency andcompare solutions with actual measurements, theresonance state is investigated. The frequency band of
60–70 MHz is calculated for 0.4 MHz steps aroundresonant frequencies, and the response for every

31

processes x 8 threads) of the Oakleaf-FX10. The iteration
count and calculation time for each frequency comparing
COCG method and COCR method are shown in Table 6.
The number of iterations required by the DDM islargest for the lowest frequency (1 MHz) in each iterative
method. Generally, the matrix equation becomes semidefinite in the full-wave electromagnetic analysis basedon the vector wave equation. However, when thefrequency becomes lower, the matrix equation becomes
non-definite. Both results indicate this thing. Fig. 8
shows convergence histories in case of 1 (MHz) using
2mm model. As a result, the iterative calculation is
improved drastically by using COCR method as interface
solver in DDM. In the previous research by Sogabe andChou, COCR method solves the matrix which size is2,534 x 2,534[4]. In this result, the number of iteration is
reduced from the result by using COCG method. On theother hand, more conspicuous improvement is seen byusing COCR method in the interface problem with theanatomical human model.
Then, to confirm performances of the huge-scale
analysis, the analyses in the smallest frequency are
performed using 1mm model that is refined Model1 twiceas the most difficult case of the computation. The modelhas 1,760 million tetrahedral elements. Computations are
performed using 600-nodes (1,200 processes x 8 threads)of the Oakleaf-FX10. The iteration count and calculation
time in 1 MHz comparing COCG method and COCR
method are shown in Table 7. Fig. 9 shows convergencehistories. We succeed calculation using huge-scale modelthat has 1,760 million tetrahedral elements. The iterative
calculation is also improved by using COCR method inthis case.

original ADVENTURE Metis. However, sincecommunication among nodes in the subdomain
decomposition is not required, the calculation time canbe reduced by increasing the number of nodes. On theother hand, the calculation times for the transformationof the voxel data to a tetrahedral mesh and thecommunication table generation represent only a smallpart of the overall calculation time and the large-scaleinitial mesh generation is not required in advance. Thesefigures of merit demonstrate that the developed function
is very effective.
Table 5. Generation conditions and elapsed time for the
HDD mesh of 4mm.
Mode
l

Num.
of
parts

1
2
3
4
5

300
600
1200
2400
4800

Num.
of
doms
/part
1000
500
250
125
62

Total
(sec)
50.12
40.64
40.21
41.35
44.21

Part
decomp-osition
(sec)
20.02
21.77
23.78
26.20
30.23

Tetra &
Com.
table
(sec)
0.09
0.08
0.08
0.09
0.11

Subdomaindecomposition (sec)
15.25
8.85
8.37
5.85
5.83

6. Numerical experiments
6.1. Performance evaluation
The electromagnetic field in the human body isanalyzed via refined meshes of Model 1 in Tables 3. All
computations are performed using Oakleaf-FX10
supercomputer at the University of Tokyo. We analyzed
frequencies of 1 (MHz), 8 (MHz), 70 (MHz) and 300(MHz) using the full-wave electromagnetic field analysis
based on the DDM. The electromagnetic field source is adipole antenna assumed to have a current source of 0.8(A). In the analyses using NICT numeric human bodymodels, it is necessary for the user to provide the electric
conductivity and the permittivity of the material ID foreach internal organ. Electric conductivities and relative
permittivities given to internal organs. These physicalvalues are based on the human organs property databasefor computational biomechanics open to the public by theInstitute of Physical and Chemical Research (RIKEN) in
Japan [10].

Table 6. Numerical result using 2mm Model.

Frequency
(MHz)
1
8
70
300

6.2. Performance evaluation
Full-wave electromagnetic analyses [11] areperformed using the 2mm model that is refined Model1
once. The model has 220 million tetrahedral elements.
Computations are performed using 75-nodes (15032

Iter. method
COCG
COCR
COCG
COCR
COCG
COCR
COCG
COCR

Num.
of Iter.
13,798
477
1,383
221
793
128
747
86

Time
(sec)
10,330
465
1,107
278
675
197
628
160

Peak
(%)
4.97
3.86
4.65
3.02
4.38
2.50
4.43
2.09

Acknowledgment
This work was supported by JSPS KAKENHI Grant
Number JP17H03256, JP17H02829, JP18KK0278,
JP19H01098
References
1. NICT EMC group home page:
http://emc.nict.go.jp/bio/index.html
2. E. Conil, A. Hadjem, F. Lacroux, M. F. Wong, and J.
Wiart,“Variability analysis of SAR from 20 MHz to 2.4
GHz for different adult and child models using finitedifference time-domain,”Phys. Med. Biol., vol.53, no.6,
pp.1511–1525 (2008)
3. A. Takei, K. Murotani, S. Sugimoto, M. Ogino and H.
Kawai, “High-accuracy electromagnetic field simulation
using numerical human body models”, IEEE Trans. Magn.,
vol. 52, no. 3 (2016)
4. T. Sogabe and S.L. Zhang, “A COCR method for solving
complex symmetric linear systems”s, J. Comput. Appl.
Math., vol.199, no.2, pp.297-303 (2007)
5. A. Takei, K. Murotani, S. Sugimoto, M. Ogino, T. Yamada,
S. Yoshimura “Performance evaluation of parallel finite
element electromagnetic field analysis using numericalhuman models,” JASSE, vol.1, no.1, pp.127-140 (2014)
6. A. Takei, M. Ogino, S. Sugimoto,” High-frequency
electromagnetic field analysis by COCR method using
anatomical human body models,” IEEE Transactions on
Magnetics, vol. 54, no. 3 (2018)
7. G. Karypis, V. Kumar, ”A fast and highly quality
multilevel scheme for partitioning irregular graphs”,
SIAM Journal on Scientific Computing, vol.20, no.1, pp.
359-392 (1998)
8. Y. Kanai, “Description of TEAM workshop problem 29:
while body cavity resonator”, Tech. Rep. TEAM
Workshop in Tucson (1998)
9. K. Murotani, S. Sugimoto, H. Kawai and S. Yoshimura,“Hierarchical domain decomposer with parallel meshrefinement function for billions of DOF scales finite
element analysis analyses,” International Journal ofComputational Methods, vol. 11, no.4 (2014)
10. The human organs property database for computational
biomechanics in Institute of Physical and Chemical
Research (RIKEN) Home page：
http://cfd-duo.riken.go.jp/cbms-mp/index.htm
11. M. Ogino, A. Takei, N. Notsu, S. Sugimoto and S.
Yoshimura, “Finite element analysis of high frequency
electromagnetic fields using a domain decomposition
method based on the COCR method”, Theor. Appl. Mech.,
vol.61, pp.173-181 (2013)

Fig. 8. Convergence histories in case of 1 MHz using
2mm model from once refined MODEL1.
Table 7. Numerical result using 1mm model.
Frequency
(MHz)
1

Iter.
method
COCG
COCR

Num.
of Iter.
7,348
1,189

Time
(sec)
7,071
1,319

Peak
(%)
4.37
3.82

Fig. 9. Convergence histories in case of 1 MHz using
1mm model from twice refined MODEL1.
7. Conclusion
In this paper, we propose parallel mesh generation
method with mesh smoothing, and show to be able togenerate large-scale mesh very effective. Then, we
evaluated the performance of two kinds of iterative
method by numerical experiments using large-scale
human body model. As the results, selecting iterative
method is very important for high-frequency
electromagnetic field analysis using anatomical humanmodel. The COCR method is expected to reduce both
iteration counts and computation time compared to the
COCG method in the analysis using huge-scale humanbody model.

33

Authors Introduction
Prof. Amane Takei

Amane Takei is working as
Associate
ProfessorforDepartment of Electrical andsystems Engineering, Universityof Miyazaki, Japan. His researchinterest
includeshigh
performance computing forcomputational
electromagnetism,
iterative
methods for the solution ofsparse linear systems, domain
decomposition methods for large-scale problems. Prof.
Takei is a member of IEEE, an expert advisor of The
Institute of Electronics, Information and Communication
Engineers (IEICE), a delegate of the Kyushu branch of
Institute of Electrical Engineers of Japan (IEEJ), a directorof Japan Society for Simulation Technology (JSST).
Prof. Makoto Sakamoto

Makoto Sakamoto received the
Ph.D. degree in computerscience and systems engineeringfrom Yamaguchi University. He
is presently an associate
professor in the Faculty ofEngineering, University ofMiyazaki. He is a theoreticalcomputer scientist, and his
current main research interests
are automata theory, languages
and computation. He is also interested in digital geometry,
digital image processing, computer vision, computer
graphics, virtual reality, augmented reality, entertainment
computing, complex systems and so on.

34