An adaptive finite element method for solving 3D electromagnetic volume integral equation with applications in microwave thermometry

An adaptive ﬁnite element method (AFEM) for the numerical solution of an electromagnetic volume integral equation (VIE) is presented. To solve the model VIE, the problem is formulated as an optimal control problem for minimization of Tikhonov’s regularization functional. A posteriori error estimates in the obtained ﬁnite element reconstruction and in the underlying Tikhonov’s functional are derived. Based on these estimates, adaptive ﬁnite element algorithms are formulated and numerically tested on the problem of microwave hyperthermia in cancer treatment. In this problem, the temperature change of a target in the computational domain results in the change of its dielectric properties. Numerical examples of monitoring this change show robust and qualitative three dimensional reconstructions of the target using the proposed adaptive algorithms. © 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The paper is devoted to the efficient and reliable solution of an electromagnetic volume integral equation (VIE) in three dimensions using an adaptive finite element method (AFEM).Problems of solution of VIE arise in many applications related to biomedical imaging and specially, to the microwave medical imaging.The specific goal of our work is to test several developed algorithms on the problem of monitoring of hyperthermia in order to reconstruct a three-dimensional target in real time.
Different mathematical analytical methods applied for investigation of three-dimensional vector electromagnetic problems related to VIE are presented in [18,19].Methods applied in computational electromagnetics are widely presented in [20][21][22], and the recently developed numerical methods for solution of inverse problems in electromagnetics are given in [23][24][25]33].
An adaptive finite element method for a Fredholm integral equation of the second kind was considered in [1].In the current work, we consider more complicated case of the VIE since it is an ill-posed problem: the solution of it is sensitive to the small perturbations in the data function and the inverse operator of VIE is not a compact operator, see [7][8][9].Thus, the problem of solution of VIE is formulated as an optimal control problem for minimization of Tikhonov's regularization functional.
In the paper, we derive a posteriori error estimates in the finite element solution of VIE and in the error of the underlying Tikhonov's functional.To do that, we specify results of [10,29] for the case of electromagnetic volume integral equation and derive new a posteriori error estimates for our specific problem.Another novelty of the current work is that we allow reconstructed function to be discontinuous what leads to new terms with jumps in a posteriori error estimates.We are not aware about some results of application of an adaptive FEM for solution of electromagnetic VIE in three dimensions with application to real-life experiments.We have developed software implemented in Matlab and C++/PETSC for solution of real-life problem of monitoring of hyperthermia which is provided at [38].Implementation of all algorithms proposed in this work is included in [38].
The main concept of an adaptive finite element method, which we use in this work for solution of VIE, is as follows: since it is not computationally efficient to use very fine mesh during the numerical solution of VIE, an adaptive finite element method allows to obtain good accuracy of the computed solution via local mesh refinements.To be able use this idea, we minimize Tikhonov functional on a sequence of locally refined meshes which are refined based on a posteriori error estimates derived in this work.Finally, different adaptive algorithms are formulated and numerically tested on the real-life problem of monitoring of hyperthermia.
Microwave hyperthermia [11] is used for treatment of cancer and is a good compliment to other types of traditional cancer therapies like radio-and chemo-therapy [12,13].The main goal of microwave hyperthermia is to increase the tumor temperature to the therapeutic levels of 40 − 44 • C while keeping healthy tissue at the normal temperature.During this procedure, it is very important to control thermal dose.Thus, there is a need for robust and real-time reconstruction methods which are capable to solve this problem.The problem of microwave hyperthermia is considered as the problem of differential microwave imaging and is formulated in many works, see [14][15][16].In the current work, we use similar experimental set-up as in [14][15][16], and consider an annular phased array system of 16 monopole antennas constructed to work at the ISM band of 915 MHz.This special frequency band is chosen because it provides a good compromise between the penetration depth of the EM waves into biological tissues and the imaging resolution.Moreover, the dielectric changes due to temperature change are significant enough at this frequency band to be detected through microwave imaging techniques.For instance, it is shown in [17] that 1% dielectric change per degree centigrade occurs between 28 • C and 53 • C and thus, it can be used in our experiment in the hyperthermic temperature range.
A linear least squares method combined with an adaptive finite element method has been presented in the recent work [5] for solution of microwave hyperthermia thermometry.However, the method considered in the present work is completely different from the AFEM method of [5]: in [5] AFEM was used for solution of image restoration problem in 2D which was formulated as solution of a Fredholm integral equation of the first kind with standard convolution kernel and measurements taken from the solution of least squares problem, see more details in Test 1 of Section 8. Thus, the adaptive method of [5] can be considered as the second method in a two-stage procedure: on the first stage, standard least squares method obtains initial reconstruction of the target, and on the second stage, AFEM improves this reconstruction.AFEM method of [5] has shown significant improvement of reconstruction obtained by the least squares method, via local adaptive mesh refinements.However, VIE was not solved in adaptive algorithm of [5] and thus, real measurements of scattered parameters were not used in it.
In the current work, we fill this gap by solving VIE which describes the process of microwave thermometry in 3D via adaptive algorithms.Our numerical examples of Section 8 show that an adaptive finite element method is a reliable and fast tool to obtain reconstructed images in 3D and thus, this method can be adopted in different microwave applications to obtain real-time feedback.
The paper is outlined as follows.In Section 2, the model VIE is derived from the vector wave equation for the electric field in non-magnetic media, and in Section 3, the mathematical model of VIE is presented.In Section 4, an ill-posed problem is formulated, and Section 5 presents the finite element method for minimization of Tikhonov's functional.Derivation of a posteriori error estimates is developed in Section 6. Section 7 describes different algorithms for solution of our ill-posed problem which are based on a posteriori error estimates of Section 6. Next, 3D numerical examples of reconstruction of the heated object during the process of hyperthermia are presented in Section 8. Finally, we end the work with concluding remarks in section 9.

Derivation of the volume integral equation
Let ⊂ R 3 be a convex bounded domain with the boundary ∂ ∈ C 2 .Following the conventional physical notations of [2], r = (x, y, z) is taken as the position vector in Cartesian coordinates.Our model PDE problem corresponding to volume integral equation is the following vector wave equation (see, for example, [2]) for the electric field in a non-magnetic medium with the Silver-Müller radiation condition at infinity: Here, ε r (r) = ε(r)/ε 0 and σ (r) are the relative dielectric permittivity and electric conductivity functions, respectively, ε 0 , μ 0 are dielectric permittivity and permeability of free space, respectively, and c = 1/ √ ε 0 μ 0 is the speed of light in free space.
Equation ( 1) can be obtained by taking the time-dependent Maxwell's equations for the electric field and then representing in these equations the time-dependent field as E(r, t) = E(r)e −iωt and the source function as J source = J (r)e −iωt , where ω is the angular frequency.We introduce the spatially distributed complex dielectric function ε (r): Then equation (1) transforms to the vector wave equation which can be rewritten as Here, ε b is the dielectric permittivity corresponding to the background medium.
To proceed further, we introduce the dyadic Green's function Ḡ(r, r ) for the problem in a homogeneous medium (see, for example, [2]): where I is an identity operator.Rewriting equation ( 5) according to the form of equation ( 6), we have We multiply (7) by Ḡ(r, r ) and ( 6) by E(r), correspondingly, to get: Then subtracting (9) from (8) and integrating, we obtain Here, (•, •) denotes the standard scalar product in space.
The above equation is reduced to the following equation We use integration by parts for the two terms on the left hand side of (11), see details in [6], and apply the Silver-Müller radiation condition at infinity, to get: Applying (12) in the left hand side of (11), we obtain Using the principle of linear superposition, see details in [2], we finally get from (11) the following equation The first term in the right hand side of the above equation corresponds to the incident electric field E inc , see [2].Hence, (14) becomes According to [4], the scattered field is defined as E sca = E(r) − E inc such that using the reciprocity Ḡ(r , r) = Ḡ(r, r ) and (15) transforms to which is the model volume integral equation.

The mathematical model
The main mathematical model considered in this work is given by the volume integral equation ( 16) for the scattered field E sca (r) in an inhomogeneous medium which is used in many microwave applications [2][3][4][5].More precisely (see details in the previous section), our mathematical model is given by the volume integral equation in the domain ⊂ R 3 : Here, is the object function, where ε (r) and ε b (r) are the complex dielectric permittivity of the medium under investigation and the dielectric permittivity of the background, respectively.
Using the Born approximation E(r) ≈ E inc (r) [2] in the right hand side of (17), we obtain the following volume integral equation According to [4], equation ( 18) transforms to the following volume integral equation for a bi-static pair of antennas (i, j): with a constant C = −ω 2 ε b /4iωμ.Using the reciprocity principle in the above equation, E inc,i is the incident electric field produced by the transmitter i and E inc, j is the incident electric field of the receiver j.In the left hand side of (18), the scattered field E sca (r) is replaced with the corresponding scattered S-parameter S scat j,i after introducing the appropriate scaling factor, i.e.C .Therefore, equation (19) uses the scattered S-parameter measurements obtained by a bi-static pair of antennas (i, j) instead of using the scattered electric field E sca (r) and is applied in many microwave applications [2][3][4][5] for image reconstruction.

Statement of an ill-posed problem
Let W 1 , W 2 , Q be three Hilbert spaces, Q ⊂ W 1 be a space which is compactly embedded in W 1 .We denote scalar products and norms in these spaces, correspondingly, as a bounded linear operator.The volume integral equation (19) can be written in the operator form (20) with an operator A : (21) and the vector of the right hand side d = S sca j,i .

The Problem (P
be unknown in .Determine O (r) for r ∈ assuming that C = −ω 2 ε b /4iωμ, E inc, j (r), E inc,i (r) and S sca j,i in (22) are known.
The problem (P) of finding a solution to the volume integral equation ( 22) is an ill-posed problem since it is a Fredholm integral equation of the first kind [9].To solve the problem (P), we introduce Tikhonov's functional where α > 0 is a regularization parameter.
Our goal is to find the function O ∈ Q which minimizes Tikhonov functional (23).We search for a stationary point of the above functional with respect to O satisfying the following equation for any b ∈ Q : where J α (O ) is the Fréchet derivative of the functional (23).
When the operator A in ( 20) is such that A : L 2 → L 2 , then the following Lemma is valid: For the proof see [9].
then the following Lemma is valid: with a convex growth factor b, i.e., |∇b| < b.The proof follows from Lemma 2 in [26].
The next Lemma is known for a bounded linear operator.We formulate this lemma for our specific case of the functional (23).Lemma 3. Let the operator A : W 1 → W 2 be a bounded linear operator which has the Fréchet derivative of the functional (23).Then the functional J α (O ) is strongly convex on the space Q and

The finite element method for minimization of Tikhonov functional
We discretize the computational domain ⊂ R 3 into elements K satisfying the minimal angle condition [27,28] by the mesh K h = {K }.Let h = h(r) be the mesh function such that where h K is the diameter of K which we define as the longest side of K .
We define the finite element space V h ⊂ V as where P 1 (K ) denotes the set of piecewise-linear functions on K .
Let us formulate the finite element method (FEM) for the functional (26).The FEM for the functional (23) with A : L 2 → L 2 can be formulated in a similar fashion with only one distinction in the formulation of the regularization term.
To find optimal solution O of the functional (26), the following problem should be solved: The finite element method for (30) reads: find The function O is approximated by where {ϕ i } N i=1 are the standard continuous piecewise linear functions and O i denote the unknown discrete function-values at the mesh point Substituting (32) into (31) with v = ϕ j , we get the discrete system of equations This system can be rewritten as which is equivalent to the following linear system of equations In system (35), matrices A, K are the finally assembled block matrices, corresponding to the first two terms in the left hand side of (34).O denotes the vector of nodal values of finite element approximation O h , and d is the finally assembled right hand side of (34).In (33), d h is the nodal interpolant of d.
From (33), we also get the discrete gradient which will be used in the iterative update of algorithms in Section 7: (∇ϕ i , ∇ϕ j )O i . (36)

A posteriori error estimates
In this section, we outline main derivations of a posteriori error estimates for the following errors: To make it clear how a posteriori error for the error in Tikhonov functional ( 23) can be obtained, we present general framework for it, see details in [29].
First we note that where is the remainder of the second order.We assume that the computed finite element approximation O h of O is close to the regularized solution O of (23).Thus, the term R(O h ) → 0 and we can neglect it in (37).
Let P h be the L 2 ( ) orthogonal projection operator.We define by f I h the standard nodal interpolant [30] of f into the space of continuous piecewise-linear functions on the mesh K h .Then by one of properties of the orthogonal projection, we have It follows from [32] that where C I = C I ( ) is a positive constant depending only on the domain .
We now use the splitting together with the Galerkin orthogonality principle Now by inserting (40) into (37), we get the following error representation: together with the interpolation property (39).
In the proofs of Theorems 1 and 2, we allow the function O to be discontinuous along the faces of elements K in K h and introduce notion of jump.We denote the jump of the function O h computed from the two neighboring elements K + and K − sharing the common side as We now make both error estimates more explicit.

A posteriori error estimate in the computed solution
In this section, we derive a posteriori error estimate for the error O − O h in the computed solution O h of the functional (23).Let the scalar product be given by (•, •) L 2 := (•, •), as well as norms •, • L 2 := •, • .Whenever the norm should be specified, we will write it henceforth.We denote by (45) Theorem 1.Let O h ∈ V h be a finite element approximations of the regularized solution O ∈ H 1 ( ) on the finite element mesh K h and O * be the exact solution of (20).Then there exists a constant D = Const.> 0 satisfying Lipschitz continuity property in a small neighborhood of the exact solution O * , such that the following a posteriori error estimate holds Here, the constant D can be chosen as D = J α (O h ) .

Proof.
Let O h be a finite element approximation of O and be the minimizer of Tikhonov functional (23).According to Lemma 3, the functional ( 23) is strongly convex on the space L 2 with the strong convexity constant γ .This implies, see where J α (O h ) , J α (O ) are the Fréchet derivatives of the functional (23).
Next, we use the splitting and substitute it in (48) to obtain Here, O I h is the standard interpolant of O .
Applying the Galerkin orthogonality principle (noting that J α (O ) = 0) in (49), we will get (51) The right hand side of (51) can be estimated using (46) as Substituting the above equation into (51), we obtain Using the interpolation property we obtain a posteriori error estimate with the interpolation constant C I : In (54), [O h ] is the jump of the function O h over the element K defined as in (44), h K is the diameter of the element K .In (54) we applied the estimate which is used often in finite element analysis, see for example, [31].
Substituting the above estimates into the right hand side of (53), we get The constant D can be estimated as follows.First we note that by the definition of the Frechét derivative of Tikhonov functional (23 where the remainder , is small and can be neglected.Then Comparing (56) with (46), we conclude that D = J α (O h ) .

A posteriori error estimates for Tikhonov functional
In Theorem 2, we derive a posteriori error estimates for the error in Tikhonov functional (23) obtained on the finite element mesh K h .Theorem 2. Suppose that there exists a minimizer O ∈ H 1 ( ) of Tikhonov functional (23) on the mesh K h .Suppose also that there exists a finite element approximation O h of O for J α (O ) on the mesh K h with the mesh function h.Then the following approximate a posteriori error estimate for the error e = || J α (O ) Proof.By the definition of the Frechét derivative of Tikhonov functional (23) with O and approximation O h , we have in the estimate (58).We again use the splitting (49) and the Galerkin orthogonality [30] J where O I h is a standard interpolant of O on the mesh K h [30].Taking norms in (60), we get where the term O − O I h can be estimated through the following interpolation estimate Substituting above estimate into (61), we get We can estimate h O H 1 ( ) in the right hand side of (62) similar to (54) to get the final estimate as follows Here, J α (O h ) is the Frechet derivative given by (64) for A : For A : L 2 → L 2 the Frechet derivative is derived in a similar fashion as in Section 5 and will be (ϕ i , ϕ j )O i . (64)

Algorithms for solution of the ill-posed problem
In this section, we formulate different algorithms for the solution of the ill-posed problem (P): conventional conjugate gradient algorithm (CG) and two different versions of adaptive finite element algorithms, AFEM2 and AFEM3.
The conjugate gradient Algorithm 1 (CG) can be used on the coarse as well as on the l times refined finite element mesh K l := K hl .Let O l := O hl is the finite element reconstruction obtained on l times refined mesh K l .The adaptive Algorithm 2 (AFEM2) is based on the solution of linear system of equation ( 35) obtained after finite element discretization of the functional (26) for A : H 1 → L 2 , or the functional (23) for A : L 2 → L 2 , to get the reconstruction O l on l times refined mesh K l .The Adaptive Algorithm 3 (AFEM3) first uses the discrete gradient g l on the l times adaptively refined mesh K l and then the conjugate gradient method (Algorithm 1) to update the reconstruction O l .

Conjugate gradient algorithm
In this algorithm, we denote the gradient for A : H 1 → L 2 for the functional (26) by where n is the number of iteration in the conjugate gradient Algorithm 1 and function O n h is the computed finite element solution of (31) at the iteration n of this algorithm.
Similarly, we denote the gradient of the functional (23) for A : Algorithm 1 Conjugate Gradient algorithm (CG).
1: Discretize the computational space domain using the finite element mesh K h .Choose the regularization parameter α. where Here, d 0 = −g 0 , and ψ is the step-sizes in the gradient update.4: For the chosen tolerance θ , stop computing O n h at the iteration M := n and obtain the function are stabilized.Otherwise set n := n + 1 and go to step 2.

Adaptive algorithms
In this section, we present two different adaptive algorithms for the solution of our ill-posed problem (P).The adaptive Algorithm 1 is based on solution of the system of linear equations (SLE) (35) while the adaptive Algorithm 2 uses computation of the discrete gradient (65) for the functional (26) when A : H 1 → L 2 , or the discrete gradient (66) for the functional (23) in the case when A : Both algorithms refine the finite element mesh locally using a posteriori error estimates provided by Theorems 1 and 2. These estimates allow us to improve the accuracy of the reconstruction O h of the problem P. First, using the estimate (47) in

Algorithm 2 Adaptive algorithm using solution of SLE (AFEM2).
1: Generate the initial FEM mesh K 0 in .Compute the sequence of reconstructions O l , l ≥ 0, on the refined meshes, in the following steps: 2: Determine the optimal regularization parameter α l .Compute the finite element solution O l on K l by solving the system of linear equations (35).
3: Use Theorems 1 and 2 for refinement: refine locally the mesh K l at all points where    Remark 1.In adaptive Algorithms 2 and 3, the number of elements in the refined mesh is controlled by the parameters β 1,2 ∈ (0, 1).The values of these parameters depend on the values of computed max |R(O l )| and max |g l | at Step 3 of both algorithms.The computational mesh K l will be refined in very narrow region if values of β 1,2 will be chosen close to 1.And opposite, if we will choose β 1,2 close to 0, then almost all elements in the finite element mesh K l will be refined.Numerical tests of Section 8 show that values β 1,2 ∈ [0.3, 0.5] are optimal one since these values provide refinement exactly at such places of K l where the inclusion is located.
1: Generate the initial FEM mesh K 0 in .Compute the sequence of reconstructions O l , l ≥ 0, on the refined meshes, in the following steps: 2: Determine the optimal regularization parameter α l .Compute the gradient g l as in (65) and then the finite element solution O l on K l as in the conjugate gradient Algorithm 1. 3: Refine locally the mesh K l as at the step 3 of Algorithm 2.  Remark 2. In some computations of Section 8 which results are reported in Tests 1, 2 and Table 2, the Balancing principle [34] was used at step 2 of adaptive Algorithms 2 and 3, for optimal computation of the regularization parameter α.
The balancing principle was implemented via fixed point algorithm of [5].Different algorithms can be chosen for optimal computation of regularization parameter α, see [9,34] and references therein for details of these algorithms.

Numerical results
In this section, we present numerical results for solution of the ill-posed problem (P) using all algorithms of Section 7 as well as the adaptive finite element algorithm (AFEM1) of [5].The goal of our simulations is to demonstrate performance

Table 1
The cooling process of a target is simulated where its temperature changes from 55 • C to 29 • C over a ten-minute window of time.As a result of this temperature change, there are changes in permittivity and conductivity of the target.Timeline (min) of all algorithms to reconstruct the function O in terms of location/shape and computational time.All simulations were performed in linux on Intel(R) Core(TM) i7-9700 CPU @ 3.00 GHz computer with the software package WavES [37].The software and perl scripts are provided at [38].
For the full-wave EM simulation, CST Microwave Studio [35] was used to simulate the scenario as shown in Fig. 2. We use the same antenna set-up for all experiments as in [5].Let us briefly present it now.All of the 16 monopole antennas, optimized to work at 915 MHz, are immersed in a coupling fluid which is a 20:7 alcohol:water mixture as suggested in [15].The choice of the matching liquid is based on a trade-off to achieve a balance between conductive loss, relative permittivity, and antenna matching.
An off-centered cylindrical inclusion with an irregular cross-section in the x-y plane is considered as the target, and its initial state at t = 0 is labeled as the baseline.The dielectric parameters of the target were then changed in a parametric sweep of CST according to Table 1 to imitate the dynamic process of the cooling over a 10-minute window of time.These changes in dielectric properties due to the temperature change are the extrapolated version of those reported in [15].Finally, the obtained S-parameter at each time step is subtracted from that of the baseline to feed the reconstruction algorithms with a differential input.Compared with 2D finite element reconstructions of [5], here we show performance of all adaptive algorithms in 3D.The dimensionless computational domain is set as and the mesh size h was h = 0.002 in x, y and z directions.The optimal value of the regularization parameter α in all methods was computed via balancing principle by the fixed point algorithm of [5] for an initial guess of the regularization parameter α 0 = 1.0.
All computations were performed on a workstation with one processor.However, since it uses C++/PETSc [36], then in principle, it can be executed on many processors as PETSc supports parallel implementations.We have estimated the computational time using the following formula for relative time T r T r = t (l + 1)n . (69) Here, t is the total computational time, n is the total number of nodes in all refined meshes, and l = 0, 1, . . . is the number of mesh refinements.Performance comparison of different reconstruction algorithms was computed and computational times (in seconds) are presented in Table 2 while Table 3 presents the relative computational time of the AFEM-based methods.The tolerance in the stopping criterion in all algorithms was θ = 5e-15.Fig. 3 also shows that an adaptive finite element method improves the reconstruction obtained by the least squares method via local adaptive mesh refinements.

Performance of different reconstruction algorithms for solution of the problem min
Fig. 4. Dynamics of reconstructed isosurfaces |O hl | on l times adaptively refined meshes at time t = 2 min using adaptive Algorithm 2 for solution of different problems: solution of problem (70) (blue isosurface) and solution of problem (71) (white isosurface).Both problems were computed with β 1,2 = 0.5 in the refinement criterion.Exact scatterer is outlined by the finite element mesh.See Table 3 for computational time and Table 4 for convergence of this test at time t = 2 min and β = 0.5.

Test 1 (CG)
The goal of computations in this test was to check performance of the conjugate gradient Algorithm 1 (CG) and balancing principle for optimal choice of the regularization parameter α in the functional (23) for different regularization terms: namely, for regularization terms α 2 O 2 2 and α 2 .We start running the conjugate gradient Algorithm 1 with the initial guess O 0 h = 0.The fixed point algorithm of [5] was used for optimal choice of the regularization parameter α which implemented balancing principle [34].
Convergence of the fixed point algorithm is presented in Fig. 1.As can be seen from the figure, for all test examples the optimal regularization parameter α was found on the interval α ∈ [0.4361, 0.4447] for the problem min and on the interval [1.3392, 1.342] for the problem min (71)  Table 2 shows performance of different reconstruction algorithms for solution of problems (70), (71) in terms of time.Using this table, we observe that all reconstruction algorithms with balancing principle usually requires more computational time than usual reconstruction algorithms with constant value of the regularization parameter α.Analyzing results of Fig. 1 and Table 2, we can conclude that the choice of the regularization parameter α = 1 is a good one for solution of both problems (70) and (71).Thus, all further tests to check performance of Algorithms 2 (AFEM2) and 3 (AFEM3) we will perform for the solution of the problem (70) or (71) with the constant value of the regularization parameter α = 1.See also Fig. 4.

Test 2 (AFEM1)
First we perform the three-dimensional AFEM computations using the adaptive algorithm of [5] which we will denote AFEM1 and briefly present below.Note that in [5] were performed only 2D AFEM computations.In AFEM1 algorithm, first the initial reconstruction O 0 has been obtained as minimization of the functional (23) with the operator A : L 2 → L 2 .Thus, the method of normal equations was used to get the solution O 0 = (K e T K e + αI) −1 K e T d.

(72)
Here, K e is the known kernel matrix resulted from discretization of the operator A(O ) in (21).Then the SVD-decomposition of K e = U V T was used in (72) to get the regularized least squares solution Next, the adaptive finite element method was used for improvement of the already obtained least squares solution O 0 thorough the minimization of the following functional where F (u) is the convolution operator similar to that of [10].Thus, the adaptive algorithm developed in [5] can be considered as an algorithm for image restoration where the actual values of the operator A(O ) were not used.We summarize it in the Algorithm 4.
1: First stage: generate the initial FEM mesh K 0 in .Obtain the initial reconstruction O 0 on K 0 via (73).
2: Second stage: compute the sequence of reconstructions u l , l ≥ 0, on the refined meshes, in the following steps: 3: Use the fixed point algorithm (balancing principle) for determining the optimal regularization parameter α.Compute the finite element solution u l on K l using the finite element method for solution of the problem (75) 4: Refine locally the mesh K l at all points where where β ∈ (0, 1) and residual R(u l ) := hu l is computed on the mesh K l .
5: Construct a new mesh K l+1 in .Interpolate O l to O l+1 and u l to u l+1 using barycentric interpolation.
6: For the chosen tolerance θ , stop computing u l at the iteration M := l and obtain the function u M if either R(u l ) L2 ( ) ≤ θ or norms u l L2 ( ) are stabilized.Otherwise go to step 3.
The three-dimensional computations were performed with β = 0.3, 0.5 in the mesh refinement criterion (76).Such values of β allow locally refine the mesh exactly where the inclusion is located, see the middle line of Fig. 5 for β = 0.5 and the top line of Fig. 7 for β = 0.3.
The optimal value of the regularization parameter α in these computations was found on the interval α ∈ [0.47, 0.69] for all refined meshes.The optimal reconstructed solution |u h | was found on a two times refined mesh for times t = 2 min and t = 4 min and on a three-times locally refined mesh for times t = 6, 8, 10 min, see Fig. 5.The top and middle lines of Fig. 5 present 3D surface solution (with and without mesh).The bottom line of Fig. 5 shows two-dimensional slices of the computed solution |u hl | on l times adaptively refined mesh K l passing through the point (36,21,7).See also Fig. 6.
Results in terms of computational time of AFEM2 for different parameters β 1,2 are presented in Table 3. Convergence tests for all times t = 2, 4, 6, 8, 10 are shown in Tables 4-8 for β 1,2 = 0.3, 0.5, and β 1,2 = 0.7 in the refinement criterion at step 3 of AFEM2.Using these tables we observe that the parameters β 1,2 control the number of the refined elements and thus, the number of the nodes in the refined mesh K hl : the computational mesh K hl is refined in very narrow region when β 1,2 = 0.7, and when we choose β 1,2 = 0.3 almost all elements in the finite element mesh K hl are refined.Also, it is required more mesh refinements for β 1,2 = 0.7 than for β 1,2 = 0.3, 0.5 in order to achieve desired tolerance in the stopping criterion.Analyzing convergence results of Tables 4-8 we can conclude that parameter β = 0.7 provides not only narrow mesh refinement which don't cover the domain where the inclusion is located but also increasing in the computational gradient g l 2 .We observe that values β 1,2 = 0.3, 0.5 are optimal one and provide decreasing of the computational gradient g l 2 and refinement of the mesh K hl exactly where the inclusion is located, see Figs. 8 and  9.
Analyzing results in Tables 4-8 and Figs. 8 and 9, we can conclude that the optimal reconstruction results in terms of time and convergence are obtained for β = 0.3 and β = 0.5.Using Figs. 8 and 9 we also observe that AFEM2 algorithm managed well reconstruct temperature changes in the object during the hyperthermia process.Results in terms of computational time of the adaptive Algorithm 3 for different parameters β 1,2 are also similar to the results obtained in adaptive Algorithm 2. Thus, both algorithms have the same performance.

Concluding remarks
An adaptive finite element method for the numerical solution of an electromagnetic volume integral equation was presented and applied to the real experimental scenario.The problem was formulated as an optimal control problem for minimization of Tikhonov's regularization functional.A posteriori error estimates in the obtained finite element reconstruction and error in Tikhonov's functional were derived.Based on these estimates, adaptive finite element algorithms were formulated and tested on 3D reconstructions of a target during the process of microwave thermometry.Our numerical results show that AFEM-based algorithms are a reliable and efficient tool for solution of microwave thermometry problem in real time.Our numerical tests also show that an adaptive finite element method improves the reconstruction obtained by the least squares method via local adaptive mesh refinements.
Future work can be related to the testing of the developed algorithms on large computational meshes using MPI in real time.The computational time for a large number of meshes processed in a single-core processor can be estimated by applying the results of Table 3.For example, using formula (69) and Table 3, we can estimate computational time of both algorithms, AFEM2 or AFEM3, on the coarse mesh consisting of 2 • 10 6 nodes without any adaptive refinement: the computational time will be computed for l = 0 as t = T r • nno ≈ Const.Test 3: Convergence of the adaptive Algorithm 3 (AFEM3) on a locally refined meshes K l := K hl for different β := β 1,2 at the time t = 4 min.Here, l is the number of the refinements of the mesh K l , nno is number of nodes in the mesh       We hope that the method presented in the paper can be adopted in different microwave applications to obtain real-time reconstructed images.
) and (42) functions O I h ∈ V h denote the interpolants of O .To derive the error O − O h in the reconstructed function O h , we will use the convexity property of Tikhonov functional (23) with the convexity parameter γ

Theorem 1 ,
we can conclude that the maximum of the computed error |O − O h | in the reconstructed function O h is located in neighborhoods of such points in the finite element mesh K h where the computed |hO h | + |[O h ]| achieves its maximal value.Next, the estimate (57) of Theorem 2 tells us that the maximum of the error in Tikhonov's functional | J α (O ) − J α (O h )| is bounded by the residual | J α (O h )| multiplied by weights |hO h + |[O h ]| and is located in the neighborhoods of mesh points where the product of the residual by weights attains its maximal value.Thus, the idea of the local finite element mesh refinement of both adaptive algorithms is that it should be refined all neighborhoods of all points in the mesh K h where the computed |hO h | + |[O h ]| achieves its maximum values or where | J α (O h )| achieves its maximal values.We define the minimizer of Tikhonov functional (23) as O , its approximated finite element solution O l h on l times adaptively refined mesh K l := K h l as O l and gradient computed on the mesh K l as g l .Let us define by O M := O M h and g M := g M h values of approximated function and computed gradient, respectively, obtained after M iterations of the conjugate gradient algorithm.

4 : 5 :
Construct a new mesh K l+1 in .Interpolate d l to d l+1 and O l to O l+1 .For the chosen tolerance θ , stop computing O l at the iteration M := l and obtain the function O M if either g l L2 ( ) ≤ θ or norms O l L2 ( ) are stabilized.Otherwise go to step 2.

4 : 5 :
Construct a new mesh K l+1 in .Interpolate d l to d l+1 and O l to O l+1 .For the chosen tolerance θ , stop computing O l at the iteration M := l and obtain the function O M if either g l L2 ( ) ≤ θ or norms O l L2 ( ) are stabilized.Otherwise go to step 2.

Fig. 1 .
Fig. 1.Test 1. Convergence of fixed point algorithm for computation of the optimal regularization parameter α using balancing principle for the functional (23): top row is for the operator A : L 2 → L 2 , bottom row is for the operator A : H 1 → L 2 .Tests are done for the conjugate gradient Algorithm 1 without adaptive mesh refinement.

Fig. 3 .
Fig. 3. Isosurfaces of reconstructions using different algorithms at t = 2 min.Table 2 presents the clarification for different colors of isosurfaces for the performed tests.(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Test 2: Reconstructions obtained by the adaptive algorithm of [5] with β = 0.5 in the refinement procedure.Top and middle lines: 3D surface solution |u hl |, and |u hl | together with the refined mesh K l , correspondingly.Bottom line: 2D slices of |u h l | passing through the point(36,21,7).Here, l denotes the number of the locally refined mesh K hl where the solution u hl was obtained.

Fig. 6 .
Fig.6.Test 2: Reconstructions obtained by the adaptive algorithm of[5] with β = 0.5 in the refinement procedure.We present 3D slices of the reconstruction |u hl | at different times.Here, l denotes the number of the locally refined mesh K l := K hl where the solution u hl was obtained.

Fig. 7 .
Fig. 7. Test 2: Reconstructions obtained by the adaptive algorithm of [5] with β = 0.3 in the refinement procedure.Top row: adaptively refined meshes for different times.Bottom line: 3D slices of the reconstruction |u hl | at different times.Here, l denotes the number of the locally refined mesh K l := K hl where the solution u hl was obtained.

Fig. 8 .
Fig. 8. Test 3: Reconstructions obtained by the adaptive Algorithm 2 with β = β 1 = β 2 = 0.3 in the refinement procedure of this algorithm.We show 3D slices of |O hl | obtained at different times on the l times locally refined mesh.Top row: prospect view, second row: xy-projection, third row: yz-projection.Fourth row: zoomed view of obtained |O hl |.Bottom row: surface view of reconstructions.Exact scatterer on all figures is outlined by the finite element mesh.

Fig. 9 .
Fig. 9. Test 3: Reconstructions obtained by the adaptive Algorithm 2 with β = β 1 = β 2 = 0.5 in the refinement procedure of this algorithm.We show 3D slices of |O hl | obtained at different times on the l times locally refined mesh.Top row: prospect view, second row: xy-projection, third row: yz-projection.Fourth row: zoomed view of obtained |O hl | superimposed with exact scatterer.Fifth row: surface view of reconstruction.Bottom row: adaptively refined meshes K h l corresponding to the obtained reconstructions |O hl |.Exact scatterer in the first -fourth rows is outlined by the finite element mesh.
• 10 −14 Start with the initial approximations O 0 h = O 0 and compute the sequences of O n h in the following steps: 2: Compute gradient g n via (65) for A : H 1 → L 2 or via (66) for A : L 2 → L 2 .3: Update the function O h := O n+1 h on K h via the conjugate gradient method

Table 3
Test 3: Performance of AFEM2 in terms of computational time in seconds and relative time for solution of the problem (71) with α = 1 and different parameters β := β 1,2 .Here, K l := K hl represents the finite element mesh refined l times where the final reconstruction O hl was obtained.See Tables 4-8 for number of nodes and convergence on the refined meshes K l .

Table 4 Test 3 :
Convergence of the adaptive Algorithm 3 (AFEM3) on a locally refined meshes K l := K hl for different β := β 1,2 at the time t = 2 min.Here, l is the number of the refinements of the mesh K l , nno is the number of nodes in the mesh K l .

Table 6
Test 3: Convergence of the adaptive Algorithm 3 (AFEM3) on a locally refined meshes K l := K hl for different β := β 1,2 at the time t = 6 min.Here, l is the number of the refinements of the mesh K l , nno is number of nodes in the mesh K l .

Table 7
Test 3: Convergence of the adaptive Algorithm 3 (AFEM3) on a locally refined meshes K l := K hl for different β := β 1,2 at the time t = 8 min.Here, l is the number of the refinements of the mesh K l , nno is number of nodes in the mesh K l .

Table 8
Test 3: Convergence of the adaptive Algorithm 3 (AFEM3) on a locally refined meshes K l := K hl for different β := β 1,2 at the time t = 10 min.Here, l is the number of the refinements of the mesh K l , nno is number of nodes in the mesh K l .