Static regime imaging of locations of certain 3D
electromagnetic imperfections from a boundary
perturbation formula

Séraphin Mefire

► To cite this version:

Séraphin Mefire. Static regime imaging of locations of certain 3D electromagnetic imperfections from a boundary perturbation formula. Journal of Computational Mathematics -International Edition-, 2014, 32 (4), pp.412-441. 10.4208/jcm.1401-m4214 . hal-01286516

HAL Id: hal-01286516

https://hal.science/hal-01286516

Submitted on 10 Mar 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Static regime imaging of locations of certain 3D electromagnetic imperfections from a boundary perturbation formula

S.M. Mefire

IECL, CNRS UMR 7502, Université de Lorraine,
Faculté des Sciences et Technologies de Nancy,
54506 Vandœuvre-lès-Nancy Cedex, France.
Email: seraphin.mefire@univ-lorraine.fr

**Abstract.** We are concerned, in a static regime, with an imaging approach of the locations in a three-dimensional bounded domain of certain electromagnetic imperfections. This approach is related to Electrical Impedance Tomography and makes use of a new perturbation formula in the electric fields. We present two localization procedures, from a Current Projection method that deals with the single imperfection context and an inverse Fourier process that is devoted to multiple imperfections configurations. These procedures extend those that were described in our previous work, since operating for a broader class of settings. Namely, the localization is additionally performed for certain purely electric imperfections, as established from numerical simulations.

**Keywords.** inverse problems, Maxwell equations, electric fields, inhomogeneities, Electrical Impedance Tomography, Current Projection method, FFT, numerical boundary measurements, random noise, edge elements, least square systems, Incomplete Modified Gram-Schmidt preconditioning.

**MS Codes.** 65N21, 65N30, 78A25

## 1 Introduction

In the context of certain applications, the imaging process mainly consists of determining from boundary measurements the locations of a finite number of defects contained in a bounded domain, and a full reconstruction of each of these defects is not necessary. Diverse reasons can motivate this determination of locations only: a priori knowledge of shapes of the defects, effective reconstruction time cost, smallness of the defects, ... Several recent works have been devoted to such a determination when the common order of magnitude of the diameters of these defects is small. The approach in these works (see e.g. [1, 2, 3, 4, 5, 12, 18]), related to Electrical Impedance Tomography, consists of localizing the defects from a particular combination of an asymptotic formula for perturbations in the potential or field, in the presence of defects, with an inversion process. Such a formula makes use of polarization tensors [3], associated with defects, and the consideration of direct inversion processes is well suited. This approach has been developed in various situations: small conductivity inhomogeneities, small elastic inhomogeneities, ...
In the presence of three-dimensional settings of small electromagnetic inhomogeneities and in the frequency regime, H. Ammari, M. Vogelius and D. Volkov [6] have proposed the asymptotic formula for perturbations in the electromagnetic fields. This tool has been thus considered by M. Asch and S.M. Mefire [8, 9] for performing the numerical localization in a three-dimensional bounded domain of these inhomogeneities, in diverse experimental configurations. Of course, in the static regime, and in contrast to the frequency regime, the localization is no longer subject to a suitable choice of frequency or of a range of frequencies. However, the inhomogeneities that are characterized by a complex electric permittivity are not considered in the localization from the static regime. In [14, 15], we deal with the static regime in several experimental configurations. The two localization procedures introduced in these works are based on a limit model in electric field and a limit perturbation model, in combination with the Current Projection method (for the single inhomogeneity context) or with an inverse Fourier method (for the multiple inhomogeneities context). These limit models are respectively obtained by letting the frequency vanish in the time-harmonic Maxwell equations and in the asymptotic formula for perturbations in the electromagnetic fields proposed in [6]. It results from [14, 15] that these procedures achieve uniquely the numerical localization of inhomogeneities involving a magnetic contrast; purely electric inhomogeneities, namely inhomogeneities of real electric permittivity and without magnetic characteristics, cannot thus be localized from these procedures. This observation leads us to an essential question, that of knowing whether, from this limit model in electric field, both purely magnetic and purely electric inhomogeneities can be localized? The numerical treatment of this question constitutes here the main objective.

Unlike in [14, 15], we will deal here with an asymptotic formula for perturbations that derives directly from the limit model in electric field.

This work is subdivided into six sections. In Section 2, after recalling the limit model in electric field following [14], we introduce an asymptotic formula that allows to study perturbations in the electric fields due to the presence of small electromagnetic inhomogeneities in a three-dimensional bounded domain $\Omega$. This asymptotic formula contains information on the electric contrast as well as on the magnetic contrast relatively to $\Omega$, and constitutes a basis for some approximate inversion techniques. Since boundary measurements, initiated by electric currents applied on the boundary of $\Omega$, shall be generated in numerical simulations and from this asymptotic formula, we briefly recall in Section 3, following [14], the (direct) computation of the electric field then required in the evaluation of each measurement. This computation is achieved by a least squares approach, with the help of Nédélec's edge elements and nodal finite elements. In Section 4, we present two localization procedures resulting from the combination of the asymptotic formula with each of the following inversion processes: the Current Projection method and an inverse Fourier method. Each procedure uses boundary measurements (in a finite number) as data and is employed in the single inhomogeneity context. Only the procedure considering an inverse Fourier method is required for multiple inhomogeneities configurations. As a result, these localization procedures extend those that were introduced in [14], since dealing with the case of purely magnetic inhomogeneities in an identical way and furthermore numerically suitable for
the context of purely electric inhomogeneities. We describe in Section 5 localization results obtained from various settings and compare some results in the single inhomogeneity context. Finally, some conclusions are reported in Section 6.

# 2 Perturbation Formula in the Electric Fields

We start by considering some notation, before introducing the perturbation formula in the electric fields.

## 2.1 Preliminaries

We are here concerned with the model also presented in [14]. Let us consider a bounded open subset $\Omega$ of $\mathbb{R}^3$, that is convex. For simplicity we assume $\partial\Omega$, the boundary of $\Omega$, connected, to be $C^\infty$, but this regularity condition could be considerably weakened. Here, $\Omega$ contains a finite number $m$ of inhomogeneities, each one of the form $z_j + \alpha B_j$, where $B_j \subset \mathbb{R}^3$ is a bounded, smooth ($C^\infty$) domain containing the origin. The points $z_j \in \Omega$, $1 \le j \le m$, that determine the locations of the inhomogeneities are assumed to satisfy:

$$ \left\{ \begin{array}{ll} 0 < d_0 \le |z_j - z_k| & \forall j \ne k, \\ d_0 \le \operatorname{dist}(z_j, \partial\Omega) & \forall j. \end{array} \right. \qquad (1) $$

The parameter $\alpha > 0$, the common order of magnitude of the diameters of the inhomogeneities, is sufficiently small in such a way that these inhomogeneities are disjoint and their distance to $\mathbb{R}^3 \setminus \Omega$ is larger than $d_0/2$. As in [8], we call hereafter, an *imperfection*, each one of these small inhomogeneities. The total collection of imperfections takes the form $I_\alpha = \bigcup_{j=1}^m (z_j + \alpha B_j)$.

Let us denote by $\mu_0 > 0$ and $\varepsilon_0$ the magnetic permeability and the electric permittivity of the background medium. Let also $\mu_j > 0$ and $\varepsilon_j$ denote the permeability and the permittivity of the $j$-th imperfection $z_j + \alpha B_j$. We assume here that the conductivity vanishes everywhere in $\Omega$. The permittivities $\varepsilon_0$ and $\varepsilon_j$ are therefore real-valued and considered such that: $\varepsilon_0 > 0$, $\varepsilon_j > 0$. By assuming that all these parameters are constant, we represent as below the piecewise constant magnetic permeability and the piecewise constant electric permittivity: $\forall x \in \Omega$,

$$ \mu_{\alpha}(x) = \begin{cases} \mu_0 & \text{if } x \in \Omega \setminus \overline{I_{\alpha}}, \\ \mu_j & \text{if } x \in z_j + \alpha B_j, \end{cases} \qquad \varepsilon_{\alpha}(x) = \begin{cases} \varepsilon_0 & \text{if } x \in \Omega \setminus \overline{I_{\alpha}}, \\ \varepsilon_j & \text{if } x \in z_j + \alpha B_j, \end{cases} $$

with $1 \le j \le m$. If we allow the degenerate case $\alpha = 0$, then the function $\mu_\alpha$ equals the constant $\mu_0$ and the function $\varepsilon_\alpha$ equals the constant $\varepsilon_0$.

Let $\{\gamma_n\}_{0 \le n \le m}$, with $\gamma_n > 0$, for $0 \le n \le m$, be related to either the set $\{\mu_n\}_{0 \le n \le m}$ or the set $\{\varepsilon_n\}_{0 \le n \le m}$. For any fixed $1 \le j_0 \le m$, let $\gamma$ denote the function defined as: $\forall x \in \mathbb{R}^3$,

$$ \gamma(x) = \begin{cases} \gamma_0 & \text{if } x \in \mathbb{R}^3 \setminus \overline{B_{j_0}}, \\ \gamma_{j_0} & \text{if } x \in B_{j_0}. \end{cases} $$
Let us represent by $x_l$, $1 \le l \le 3$, the coordinates of $x \in \mathbb{R}^3$. For any fixed $1 \le l \le 3$, denote by $\phi_l$ the solution to:

$$ \left\{ \begin{array}{ll} \operatorname{div}(\gamma(x) \operatorname{grad} \phi_l(x)) = 0 & \text{for } x \in \mathbb{R}^3, \\ \phi_l(x) - x_l \to 0 & \text{as } |x| \to \infty. \end{array} \right. $$

The scalar potential $\phi_l$ depends in fact only on $\gamma_0$ and $\gamma_{j_0}$ through the ratio $c = \frac{\gamma_0}{\gamma_{j_0}}$. Here, the essential assumption is that the constant $c$ cannot be zero or a negative real number. With this aspect ratio, we define (as in [6]) the polarization tensor, $M^{j_0}(c)$, of the inhomogeneity $B_{j_0}$ as follows: $\forall 1 \le i, l \le 3$,

$$ M_{il}^{j_0}(c) = c^{-1} \int_{B_{j_0}} \frac{\partial \phi_l}{\partial x_i} dx. \qquad (2) $$

Following [6] for instance, the tensor $M^{j_0}(c)$ is symmetric, and is furthermore positive definite if $c \in \mathbb{R}_+^*$.

As in [14], we are here concerned with the following limit model. In the presence of imperfections, the electric field denoted $E_\alpha$ satisfies:

$$ \left\{ \begin{array}{ll} \operatorname{curl}(\frac{1}{\mu_\alpha} \operatorname{curl} E_\alpha) = 0 & \text{in } \Omega, \\ \operatorname{div}(\varepsilon_\alpha E_\alpha) = 0 & \text{in } \Omega, \\ E_\alpha \times \nu = g & \text{on } \partial\Omega, \end{array} \right. \qquad (3) $$

where $g$ is a prescribed datum on $\partial\Omega$, and $\nu$ represents the outward unit normal to $\Omega$, defined on $\partial\Omega$. The electric field denoted $E_0$, in the absence of all the imperfections, is such that:

$$ \left\{ \begin{array}{ll} \operatorname{curl}(\frac{1}{\mu_0} \operatorname{curl} E_0) = 0 & \text{in } \Omega, \\ \operatorname{div}(\varepsilon_0 E_0) = 0 & \text{in } \Omega, \\ E_0 \times \nu = g & \text{on } \partial\Omega. \end{array} \right. \qquad (4) $$

Let

$$ H(\operatorname{curl}; \Omega) = \{u \in (L^2(\Omega))^3; \operatorname{curl} u \in (L^2(\Omega))^3\} $$

be endowed with its usual Hermitian product denoted here by $(\cdot, \cdot)_{H(\operatorname{curl}; \Omega)}$; the corresponding norm is denoted by $\| \cdot \|_{H(\operatorname{curl}; \Omega)}$. The vector fields $E_\alpha$ and $E_0$ shall be determined in $H(\operatorname{curl}; \Omega)$. By representing the surface divergence by $\operatorname{div}_{\partial\Omega}$, let us consider the space

$$ TH^{-\frac{1}{2}}(\operatorname{div}; \partial\Omega) = \{q \in (H^{-\frac{1}{2}}(\partial\Omega))^3; \operatorname{div}_{\partial\Omega} q \in H^{-\frac{1}{2}}(\partial\Omega), q \cdot \nu = 0 \text{ on } \partial\Omega\}, $$

with its usual norm denoted here by $\| \cdot \|_{TH^{-\frac{1}{2}}(\operatorname{div}; \partial\Omega)}$.

Of course, the datum $g$ is taken in $TH^{-\frac{1}{2}}(\operatorname{div}; \partial\Omega)$, and we consider $u_g \in H(\operatorname{curl}; \Omega)$ such that (see e.g. [7]):

$$ \left\{ \begin{array}{ll} u_g \times \nu = g & \text{on } \partial\Omega, \\ \|u_g\|_{H(\operatorname{curl}; \Omega)} \le C_{\Omega} \|g\|_{TH^{-\frac{1}{2}}(\operatorname{div}; \partial\Omega)}, & \text{on } \partial\Omega. \end{array} \right. \qquad (5) $$
where $C_\Omega > 0$ is a constant depending only on $\Omega$. With this extension field, the determination of $E_\alpha$ satisfying (3) is reduced to the problem which consists of finding $\mathcal{E}_\alpha$ such that:

$$ egin{cases} \mathrm{curl}(\frac{1}{\mu_\alpha} \mathrm{curl} \mathcal{E}_\alpha) = - \mathrm{curl}(\frac{1}{\mu_\alpha} \mathrm{curl} u_g) & \text{in } \Omega, \\ \mathrm{div}(\varepsilon_\alpha \mathcal{E}_\alpha + \varepsilon_\alpha u_g) = 0 & \text{in } \Omega, \\ \mathcal{E}_\alpha \times \nu = 0 & \text{on } \partial\Omega. \end{cases} \tag{6} $$

Of course, knowing $u_g$, while $\mathcal{E}_\alpha$ is in accordance with (6), we determine the electric field:

$$ E_\alpha := \mathcal{E}_\alpha + u_g. \tag{7} $$

Let us set:

$$ \begin{gathered} \Psi = H_0^1(\Omega), \\ \mathcal{H} = \{u \in H(\mathrm{curl}; \Omega); u \times \nu = 0 \text{ on } \partial\Omega\}. \end{gathered} $$

The vector field $\mathcal{E}_\alpha$ is sought in the space $\mathcal{H}$. Let us consider $q \in \Psi$ such that

$$ (\varepsilon_\alpha \mathrm{grad} q, \mathrm{grad} \psi)_{(L^2(\Omega))^3} = (\varepsilon_\alpha u_g, \mathrm{grad} \psi)_{(L^2(\Omega))^3} \quad \forall \psi \in \Psi, \tag{8} $$

and introduce the new unknown

$$ U_\alpha := \mathcal{E}_\alpha + \mathrm{grad} q, \tag{9} $$

as well as the space

$$ \mathcal{V} = \{v \in H(\mathrm{curl}; \Omega); (\varepsilon_\alpha v, \mathrm{grad} \psi)_{(L^2(\Omega))^3} = 0 \quad \forall \psi \in \Psi, v \times \nu = 0 \text{ on } \partial\Omega\}. $$

Here, the notation $(\cdot, \cdot)_{(L^2(\Omega))^3}$ denotes the usual Hermitian product of $(L^2(\Omega))^3$ and the associated norm will be represented by $\|\cdot\|_{(L^2(\Omega))^3}$. The space $\mathcal{V}$ is endowed with the norm equivalent to $\|\cdot\|_{H(\mathrm{curl}; \Omega)}$ and generated (see e.g. [7]) by the mapping $u \in \mathcal{V} \mapsto \|\mathrm{curl} u\|_{(L^2(\Omega))^3}$ according to hypotheses on $\Omega$ and to the definition of the real-valued parameter $\varepsilon_\alpha$. Following [14], we introduce the weak formulation defined below for $g \in TH^{-\frac{1}{2}}(\mathrm{div}; \partial\Omega)$ and therefore for $u_g$ taken as in (5).

Find $U_\alpha \in \mathcal{V}$ such that:

$$ (\frac{1}{\mu_\alpha} \mathrm{curl} U_\alpha, \mathrm{curl} v)_{(L^2(\Omega))^3} = -(\frac{1}{\mu_\alpha} \mathrm{curl} u_g, \mathrm{curl} v)_{(L^2(\Omega))^3} \quad \forall v \in \mathcal{V}. \tag{10} $$

The following statement, regarding (10) and established in [14], leads to the determination of a unique vector field $\mathcal{E}_\alpha$ in accordance with (6), due to the definition of the scalar potential $q$ in (9). It ensures then the existence and the uniqueness of $E_\alpha$ subject to (3), by taking into account (7).

**Theorem 2.1** For any $g \in TH^{-\frac{1}{2}}(\mathrm{div}; \partial\Omega)$, and therefore any $u_g$ defined as in (5), the formulation (10) has one and only one solution $U_\alpha \in \mathcal{V}$. Furthermore, there exists a constant $C > 0$ independent of $\alpha$ such that: $\|U_\alpha\|_{H(\mathrm{curl}; \Omega)} \le C \|g\|_{TH^{-\frac{1}{2}}(\mathrm{div}; \partial\Omega)}$.

Let us now derive a perturbation formula due to the presence of the imperfections in the domain $\Omega$, and directly based on (3)-(4).
## 2.2 Perturbation formula

In contrast with [14], we are here interested in a study of perturbations resulting directly from the limit model in electric field. Let us first notice the following preliminary remark.

**Remark 2.1** The vector fields $E_\alpha$ and $E_0$, subject to (3) and (4) respectively, are in particular such that:

$$ \mathrm{curl}(\frac{1}{\mu_{\alpha}} \mathrm{curl} E_{\alpha}) - \mathrm{grad}(\mathrm{div}(\varepsilon_{\alpha} E_{\alpha}) = 0 \text{ in } \Omega, \qquad (11) $$

$$ \mathrm{curl}(\frac{1}{\mu_0} \mathrm{curl} E_0) - \mathrm{grad}(\mathrm{div}(\varepsilon_0 E_0) = 0 \text{ in } \Omega. \qquad (12) $$

A description of the boundary perturbation in the electric field due to the presence of imperfections is summarized as below.

**Proposition 2.1** Let (1) be satisfied and $g \in TH^{-\frac{1}{2}}(\mathrm{div}; \partial\Omega)$. Let us denote by $w$ any smooth vector-valued function such that:

$$ \mathrm{curl}(\frac{1}{\mu_0} \mathrm{curl} w) - \mathrm{grad}(\mathrm{div}(\varepsilon_0 w)) = 0 \text{ in } W, \qquad (13) $$

where $W$ is an open neighborhood of $\Omega$. Let us consider the vector field $E_\alpha$, where $0 < \alpha < \alpha_0$, with the constant $\alpha_0$ depending on $\{B_j\}_{1\le j\le m}$, $\Omega$, $\{\mu_j\}_{0\le j\le m}$, $\{\varepsilon_j\}_{0\le j\le m}$ and $d_0$, but otherwise independent of $g$, $w$, and of the points $z_j$, $1 \le j \le m$. The vector field $E_\alpha$ satisfies:

$$ \int_{\partial\Omega} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma - \int_{\partial\Omega} \mathrm{curl} w \times \nu \cdot (\nu \times (E_{\alpha} \times \nu)) \, d\sigma - \mu_0 \int_{\partial\Omega} (\varepsilon_0 E_{\alpha}) \cdot \nu \, \mathrm{div}(w) \, d\sigma = \\ \int_{I_{\alpha}} (\frac{\mu_0}{\mu_{\alpha}} - 1) \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w \, dx - \mu_0 \int_{I_{\alpha}} (\varepsilon_{\alpha} - \varepsilon_0) E_{\alpha} \cdot \mathrm{grad}(\mathrm{div}(w)) \, dx. \qquad (14) $$

*Proof.* From (11), it follows that

$$ \int_{\Omega} \frac{1}{\mu_{\alpha}} \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w \, dx - \int_{\partial\Omega} \frac{1}{\mu_{\alpha}} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma + \int_{\Omega} \mathrm{div}(\varepsilon_{\alpha} E_{\alpha}) \, \mathrm{div}(w) \, dx = 0, $$

since $\int_{\partial\Omega} \mathrm{div}(\varepsilon_{\alpha} E_{\alpha}) w \cdot \nu \, d\sigma = 0$. Then,

$$ \int_{\Omega} \frac{1}{\mu_0} \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w \, dx + \int_{I_{\alpha}} (\frac{1}{\mu_{\alpha}} - \frac{1}{\mu_0}) \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w \, dx \\ - \int_{\partial\Omega} \frac{1}{\mu_{\alpha}} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma + \int_{\Omega} \mathrm{div}(\varepsilon_{\alpha} E_{\alpha}) \, \mathrm{div}(w) \, dx = 0, $$

i.e.,

$$ \int_{\Omega} E_{\alpha} \cdot \mathrm{curl}(\frac{1}{\mu_0} \mathrm{curl} w) \, dx + \int_{I_{\alpha}} (\frac{1}{\mu_{\alpha}} - \frac{1}{\mu_0}) \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w \, dx + \int_{\Omega} \mathrm{div}(\varepsilon_{\alpha} E_{\alpha}) \, \mathrm{div}(w) \, dx \\ - \int_{\partial\Omega} \frac{1}{\mu_0} E_{\alpha} \times \nu \cdot \mathrm{curl} w \, d\sigma - \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma = 0. $$
Thus,

$$ 
\begin{aligned}
& \int_{\Omega} E_{\alpha} \cdot \mathrm{curl}(\frac{1}{\mu_0} \mathrm{curl} w) dx + \int_{I_{\alpha}} (\frac{1}{\mu_{\alpha}} - \frac{1}{\mu_0}) \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w dx - \int_{\Omega} E_{\alpha} \cdot \mathrm{grad}(\mathrm{div}(\varepsilon_0 w)) dx \\
& \quad - \int_{I_{\alpha}} (\varepsilon_{\alpha} - \varepsilon_0) E_{\alpha} \cdot \mathrm{grad}(\mathrm{div} w) dx - \int_{\partial\Omega} \frac{1}{\mu_0} E_{\alpha} \times \nu \cdot \mathrm{curl} w d\sigma \\
& \quad - \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} E_{\alpha} \times \nu \cdot w d\sigma + \int_{\partial\Omega} (\varepsilon_0 E_{\alpha}) \cdot \nu \mathrm{div} w d\sigma = 0.
\end{aligned}
$$

By making now use of (13), it follows finally that

$$ 
\begin{aligned}
& - \int_{\partial\Omega} \frac{1}{\mu_0} E_{\alpha} \times \nu \cdot \mathrm{curl} w d\sigma - \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} E_{\alpha} \times \nu \cdot w d\sigma + \int_{\partial\Omega} (\varepsilon_0 E_{\alpha}) \cdot \nu \mathrm{div} w d\sigma \\
& \quad + \int_{I_{\alpha}} (\frac{1}{\mu_{\alpha}} - \frac{1}{\mu_0}) \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w dx - \int_{I_{\alpha}} (\varepsilon_{\alpha} - \varepsilon_0) E_{\alpha} \cdot \mathrm{grad}(\mathrm{div} w) dx = 0,
\end{aligned}
$$

and we obtain the relation

$$ 
\begin{aligned}
& \int_{\partial\Omega} \mathrm{curl} w \times \nu \cdot E_{\alpha} d\sigma - \int_{\partial\Omega} \mathrm{curl} E_{\alpha} \times \nu \cdot w d\sigma + \mu_0 \int_{\partial\Omega} (\varepsilon_0 E_{\alpha}) \cdot \nu \mathrm{div} w d\sigma \\
& \quad + \int_{I_{\alpha}} (\frac{\mu_0}{\mu_{\alpha}} - 1) \mathrm{curl} E_{\alpha} \cdot \mathrm{curl} w dx - \mu_0 \int_{I_{\alpha}} (\varepsilon_{\alpha} - \varepsilon_0) E_{\alpha} \cdot \mathrm{grad}(\mathrm{div} w) dx = 0,
\end{aligned}
$$

which leads then to (14). □

<figure><img src="image_2.png" /></figure>

The description summarized in Proposition 2.1 allows us to introduce the following perturbation formula.

**Proposition 2.2** Let (1) be satisfied and $g \in TH^{-\frac{1}{2}}(\mathrm{div}; \partial\Omega)$. Let us denote by $w$ any smooth vector-valued function satisfying (13). Let us consider the vector field $E_{\alpha}$, where $0 < \alpha < \alpha_0$, with the constant $\alpha_0$ depending on $\{B_j\}_{1 \le j \le m}$, $\Omega$, $\{\mu_j\}_{0 \le j \le m}$, $\{\varepsilon_j\}_{0 \le j \le m}$ and $d_0$, but otherwise independent of $g$, $w$, and of the points $z_j$, $1 \le j \le m$. The vector fields $E_{\alpha}$ and $E_0$ satisfy:

$$ 
\begin{aligned}
& \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} E_{\alpha} \times \nu \cdot w d\sigma - \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} w \times \nu \cdot (\nu \times (E_{\alpha} \times \nu)) d\sigma - \int_{\partial\Omega} (\varepsilon_0 E_{\alpha}) \cdot \nu \mathrm{div}(w) d\sigma \\
& = \alpha^3 \sum_{j=1}^{m} (\frac{1}{\mu_j} - \frac{1}{\mu_0}) \left[ M^j(\frac{\mu_0}{\mu_j}) \mathrm{curl} E_0(z_j) \right] \cdot \mathrm{curl} w(z_j) \quad (15) \\
& \quad + \alpha^3 \sum_{j=1}^{m} \varepsilon_0(\frac{\varepsilon_0}{\varepsilon_j} - 1) \left[ M^j(\frac{\varepsilon_0}{\varepsilon_j}) E_0(z_j) \right] \cdot \mathrm{grad}(\mathrm{div}(w))(z_j) + o(\alpha^3).
\end{aligned}
$$

Before concluding this section with some remarks on (15), let us present the main lines of the proof of Proposition 2.2 at first. In fact, the formula (15) is obtained from (14) and by making use of an asymptotic matching technique of expansions of $E_{\alpha}$ similar to the one developed in [4] in the context of the time-harmonic electromagnetic fields. Let us consider (14), where $E_{\alpha}$ is subject to (3), whereas the smooth vector-valued function $w$ satisfies (13). Also, the vector field $E_0$ is subject to (4).
Replacing curl $w$ and grad(div($w$)) by curl $w(z_j)$ and grad(div($w$))(z_j) in z_j + αB_j, with 1 ≤ j ≤ m fixed, we formally obtain from (14),

$$
\begin{align}
& \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} E_\alpha \times \nu \cdot w \, d\sigma - \int_{\partial\Omega} \frac{1}{\mu_0} \mathrm{curl} w \times \nu \cdot (\nu \times (E_\alpha \times \nu)) \, d\sigma - \int_{\partial\Omega} (\varepsilon_0 E_\alpha) \cdot \nu \, \mathrm{div}(w) \, d\sigma \nonumber \\
& = \sum_{j=1}^{m} \left(\frac{1}{\mu_j} - \frac{1}{\mu_0}\right) \left(\int_{z_j + \alpha B_j} \mathrm{curl} E_\alpha \, dx\right) \cdot \mathrm{curl} w(z_j) \tag{16} \\
& - \sum_{j=1}^{m} (\varepsilon_j - \varepsilon_0) \left(\int_{z_j + \alpha B_j} E_\alpha \, dx\right) \cdot \mathrm{grad}(\mathrm{div}(w))(z_j) + o(\alpha^3). \nonumber
\end{align}
$$

Let us consider now, for $1 \le j \le m$ fixed, expressions of the fields $E_\alpha$ and $\mathrm{curl} E_\alpha$, inside $z_j + \alpha B_j$ based on a matching asymptotic expansions technique. Let us denote by $y = \frac{x-z_j}{\alpha}$ the local variable. For the field $E_\alpha(x)$, we expect that it will differ appreciably from $E_0(x)$ for $x$ near $z_j$, but that it will differ little from $E_0(x)$ for $x$ far from $z_j$. This field is here represented by two different expansions, namely an inner expansion for $x$ near $z_j$ and an outer expansion for $x$ far from $z_j$. The outer expansion is expressed as,

$$
E_{\alpha}(x) = E_{0}(x) + \alpha^{\tau_{1}} E_{1}(x) + \alpha^{\tau_{2}} E_{2}(x) + \dots, \quad \text{for } |x - z_{j}| \gg \mathcal{O}(\alpha),
$$

where $0 < \tau_1 < \tau_2 < \dots$, and $E_1, E_2, \dots$, are to be found. The inner expansion is written as

$$
E_{\alpha}(z_j + \alpha y) = e_{\alpha}(y) = e_0(y) + \alpha e_1(y) + \alpha^2 e_2(y) + \dots, \quad \text{for } |y| = \mathcal{O}(1),
$$

where $e_0, e_1, \dots$, are to be found. The two expansions must be asymptotically equal in some “overlap” domain within which the stretched variable $y$ is large and $x-z_j$ is small. The matching constraint in this domain is

$$
E_0(x) + \alpha^{\tau_1} E_1(x) + \alpha^{\tau_2} E_2(x) + \dots \sim e_0(y) + \alpha e_1(y) + \alpha^2 e_2(y) + \dots
$$

The first matching condition, obtained from the terms of order $\alpha^0$, indicates that

$$
e_0(y) \rightarrow E_0(z_j) \quad \text{as} \quad |y| \rightarrow +\infty.
$$

By substituting the inner expansion into the two first relations of (3) and formally equating coefficients
of $\alpha^{-2}$ and $\alpha^{-1}$, it follows respectively that

$$
\mathrm{curl}_y \left( \frac{1}{\mu(y)} \mathrm{curl}_y e_0(y) \right) = 0 \quad \forall y \in \mathbb{R}^3, \quad \mathrm{div}_y (\varepsilon(y)e_0(y)) = 0 \quad \forall y \in \mathbb{R}^3,
$$

where

$$
\mu(y) = \begin{cases} \mu_0 & \text{if } y \in \mathbb{R}^3 \setminus \overline{B_j}, \\ \mu_j & \text{if } y \in B_j, \end{cases} \quad \varepsilon(y) = \begin{cases} \varepsilon_0 & \text{if } y \in \mathbb{R}^3 \setminus \overline{B_j}, \\ \varepsilon_j & \text{if } y \in B_j. \end{cases}
$$

The field $e_0$ is then subject to

$$
\left\{
\begin{array}{ll}
\mathrm{curl}_y \left( \frac{1}{\mu(y)} \mathrm{curl}_y e_0(y) \right) = 0 & \quad \forall y \in \mathbb{R}^3, \\
\\
\mathrm{div}_y (\varepsilon(y)e_0(y)) = 0 & \quad \forall y \in \mathbb{R}^3, \\
\\
e_0(y) \rightarrow E_0(z_j) & \quad \text{as } |y| \rightarrow +\infty.
\end{array}
\right.
\qquad (17)
$$
Let us represent by $y_1$, $y_2$ and $y_3$ the coordinates of $y \in \mathbb{R}^3$. With the notation used in [4], let us consider now the scalar potential $v_{1i}^{\frac{1}{\varepsilon}}$ (see Page 765 of [4]) such that

$$ \left\{ \begin{array}{ll} \operatorname{div}_y (\varepsilon(y) \operatorname{grad}_y v_{1i}^{\frac{1}{\varepsilon}}(y)) = 0 & \forall y \in \mathbb{R}^3, \\ v_{1i}^{\frac{1}{\varepsilon}}(y) \to y_i & \text{as } |y| \to +\infty. \end{array} \right. \qquad (18) $$

By denoting by $\{u_i\}_{1\le i\le 3}$ the canonical basis of $\mathbb{R}^3$, the asymptotic behavior of $\operatorname{grad}_y v_{1i}^{\frac{1}{\varepsilon}}$ is such that, $\operatorname{grad}_y v_{1i}^{\frac{1}{\varepsilon}}(y) = u_i + \frac{c}{|y|^3}y + o(\frac{1}{|y|^2})$, where the constant $c$ depends on $\varepsilon_0$, $\varepsilon_j$, and $B_j$. By considering, for $y \in \mathbb{R}^3$,

$$ e_0(y) = \sum_{i=1}^{3} (E_0(z_j) \cdot u_i) \operatorname{grad}_y v_{1i}^{\frac{1}{\varepsilon}}(y), \qquad (19) $$

we determine in fact $e_0$ subject to (17) since this field satisfies the relations $\operatorname{curl}_y e_0 = 0$, $\operatorname{div}_y (\varepsilon e_0) = 0$, and is such that $e_0(y) \to E_0(z_j)$ as $|y| \to +\infty$.

In a similar way, by setting $\tilde{h}_\alpha(x) = \operatorname{curl} E_\alpha(x)$, for $x \in \Omega$, it follows that the first order term, $\tilde{h}_0$, in the inner asymptotic expansion of this field is subject to

$$ \left\{ \begin{array}{ll} \operatorname{curl}_y \left(\frac{1}{\mu(y)} \tilde{h}_0(y)\right) = 0 & \forall y \in \mathbb{R}^3, \\ \operatorname{div}_y (\tilde{h}_0(y)) = 0 & \forall y \in \mathbb{R}^3, \\ \tilde{h}_0(y) \to \operatorname{curl} E_0(z_j) & \text{as } |y| \to +\infty. \end{array} \right. $$

Then the vector field $s_0$, defined in $\mathbb{R}^3$ as $s_0 := \frac{1}{\mu}\tilde{h}_0$, is such that

$$ \left\{ \begin{array}{ll} \operatorname{curl}_y (s_0(y)) = 0 & \forall y \in \mathbb{R}^3, \\ \operatorname{div}_y (\mu(y)s_0(y)) = 0 & \forall y \in \mathbb{R}^3, \\ s_0(y) \to \frac{1}{\mu_0} \operatorname{curl} E_0(z_j) & \text{as } |y| \to +\infty. \end{array} \right. $$

With again the notation used in [4], it follows that, for $y \in \mathbb{R}^3$,

$$ \frac{1}{\mu(y)} \tilde{h}_0(y) = s_0(y) = \sum_{i=1}^{3} \left( \frac{1}{\mu_0} \operatorname{curl} E_0(z_j) \cdot u_i \right) \operatorname{grad}_y v_{1i}^{\frac{1}{\mu}}(y), \qquad (20) $$

where the scalar potential $v_{1i}^{\frac{1}{\mu}}$ satisfies

$$ \left\{ \begin{array}{ll} \operatorname{div}_y (\mu(y) \operatorname{grad}_y v_{1i}^{\frac{1}{\mu}}(y)) = 0 & \forall y \in \mathbb{R}^3, \\ v_{1i}^{\frac{1}{\mu}}(y) \to y_i & \text{as } |y| \to +\infty. \end{array} \right. \qquad (21) $$

By using the asymptotic expression of $E_\alpha$ in accordance with (19) and the one of $\operatorname{curl} E_\alpha$ combined with (20), it results from (16) that

$$ \begin{aligned} & \int_{\partial\Omega} \frac{1}{\mu_0} \operatorname{curl} E_\alpha \times \nu \cdot w \, d\sigma - \int_{\partial\Omega} \frac{1}{\mu_0} \operatorname{curl} w \times \nu \cdot (\nu \times (E_\alpha \times \nu)) \, d\sigma - \int_{\partial\Omega} (\varepsilon_0 E_\alpha) \cdot \nu \operatorname{div}(w) \, d\sigma \\ &= \sum_{j=1}^{m} \left(\frac{1}{\mu_j} - \frac{1}{\mu_0}\right) \alpha^3 \left( \int_{B_j} \mu(y) \left[ \sum_{i=1}^{3} \left(\frac{1}{\mu_0} \operatorname{curl} E_0(z_j) \cdot u_i\right) \operatorname{grad}_y v_{1i}^{\frac{1}{\mu}}(y) \right] dy \right) \cdot \operatorname{curl} w(z_j) \quad (22) \\ & - \sum_{j=1}^{m} (\varepsilon_j - \varepsilon_0) \alpha^3 \left( \int_{B_j} \left[ \sum_{i=1}^{3} (E_0(z_j) \cdot u_i) \operatorname{grad}_y v_{1i}^{\frac{1}{\varepsilon}}(y) \right] dy \right) \cdot \operatorname{grad}(\operatorname{div}(w))(z_j) + o(\alpha^3). \end{aligned} $$
The notation in (2) allows us to mention that

$$ 
\left\{ 
\begin{aligned}
& \int_{B_j} \left[ \sum_{i=1}^{3} (E_0(z_j) \cdot u_i) \operatorname{grad} v_{1i}^{\frac{1}{\varepsilon}}(y) \right] dy = \frac{\varepsilon_0}{\varepsilon_j} M^j(\frac{\varepsilon_0}{\varepsilon_j}) E_0(z_j), \\
& \int_{B_j} \mu(y) \left[ \sum_{i=1}^{3} (\frac{1}{\mu_0} \operatorname{curl} E_0(z_j) \cdot u_i) \operatorname{grad} v_{1i}^{\frac{1}{\mu}}(y) \right] dy = M^j(\frac{\mu_0}{\mu_j}) \operatorname{curl} E_0(z_j).
\end{aligned}
\right. 
$$

By using in (22) these relations, it derives then (15).

Let us now make some remarks related to (15).

**Remark 2.2** In the particular context where the field w subject to (13) is furthermore divergence-free, we get from the perturbation formula (15), by neglecting the asymptotically small remainder term, the approximation formula of boundary measurements that was presented in [14], in combination with the model (3)-(4). The numerical investigations of [14] established that this approximation formula was particularly suitable for localizations of imperfections that are not purely electric.

**Remark 2.3** The formula (15) contains both the explicit information on the parameter $\mu_\alpha$ and the parameter $\varepsilon_\alpha$, contrary to the approximation formula presented in [14]. We can thus expect here, not only localizations of imperfections that are purely magnetic or of electromagnetic type with real-valued permittivity, as indicates Remark 2.2, but also the localization of imperfections that are purely electric. The context where the field w subject to (13) is furthermore such that

$$ \operatorname{grad}(\operatorname{div}(w)) \neq 0 \text{ in } W, \qquad (23) $$

will be hence considered.

**Remark 2.4** The relations (18) and (21) indicate (15) as being a formula which reflects a superposition of contrasts simultaneously induced by a conductivity problem in $\mu$ and a conductivity problem in $\varepsilon$. The polarization tensors intervening in (15) are the same as in the formula known from [5] in the frequency context for perturbations of the electromagnetic fields.

By combining the perturbation formula (15) with particular inversion algorithms, we shall present procedures (see Section 4) for effective localizations of imperfections.

# 3 Numerical Approximations

Let us assume for simplifying the presentation, here and in the following sections, that each inhomogeneity present in $\Omega$ is a polyhedron, as is $\Omega$. In simulations, we will make use of the formula (15) and consider the numerical approximation of $E_\alpha$ performed as in [14]. Namely, the conforming mesh of the domain $\bar{\Omega}$ consists of a collection $\mathcal{T}_\alpha$ made up of tetrahedra, obtained from a usual process of triangulation, in such a way that each inhomogeneity is entirely covered by a distinct subset of tetrahedra of
$\mathcal{T}_\alpha$. This collection is also such that there exists a constant $c > 0$ such that $\sup_{K \in \mathcal{T}_\alpha} \frac{h_K}{\partial K} \le c$, where $h_K$ denotes the diameter of the tetrahedron $K$ and $\varrho_K$ is the diameter of the largest sphere included in $K$. The mesh size $h$ of the domain, $h = \sup_{K \in \mathcal{T}_\alpha} h_K$, depends in particular on $\alpha$ in the sense that we need to have a (global) mesh as fine as its part associated with the smallest inhomogeneity.

The discrete formulation associated with (10) is then introduced with in particular the help of the edge elements (see Nédélec [16]) of the first order. Representing by $K$ a tetrahedron of $\mathcal{T}_\alpha$, let us consider

$$R^1(K) = \{u : K \to \mathbb{C}^3; \exists a, b \in \mathbb{C}^3, u(x) = a + b \times x, x = (x_1, x_2, x_3)^T \in K\},$$

where the superscript “T” denotes, here and in the next sections, the transpose. Let us set

$$\mathcal{R}_h = \{u_h \in H(\text{curl}; \Omega); u_h|_K \in R^1(K) \ \forall K \in \mathcal{T}_\alpha\},$$

and associate with $\mathcal{H}$ the discrete space

$$\mathcal{H}_h = \{u_h \in \mathcal{R}_h; u_h \times \nu = 0 \text{ on } \partial\Omega\},$$

also endowed, as in the case of $\mathcal{H}$, with $(\cdot, \cdot)_{H(\text{curl}; \Omega)}$. The expression of any element of $\mathcal{R}_h$ in each tetrahedron $K \in \mathcal{T}_\alpha$ can be written similarly as was done in [13] with $\mathbb{R}^3$-valued fields, for a practical implementation.

Let us next associate with $\Psi$ the discrete space

$$\Psi_h = \{\psi_h \in H^1(\Omega); \text{Re}(\psi_h)|_K, \text{Im}(\psi_h)|_K \in P_1(K) \ \forall K \in \mathcal{T}_\alpha, \psi_h = 0 \text{ on } \partial\Omega\},$$

where $P_1(K)$ is the space of polynomials of degree less than or equal to 1 defined in $K$. Finally, we consider

$$\mathcal{V}_h = \{v_h \in \mathcal{R}_h; (\varepsilon_\alpha v_h, \text{grad } \psi_h)_{(L^2(\Omega))^3} = 0 \ \forall v_h \in \Psi_h, v_h \times \nu = 0 \text{ on } \partial\Omega\},$$

the discrete space associated with $\mathcal{V}$, and endowed with the norm of $\mathcal{V}$. The discrete formulation associated with (10) consists of finding $U_h \in \mathcal{V}_h$ such that:

$$(\frac{1}{\mu_{\alpha}} \text{curl } U_h, \text{curl } v_h)_{(L^2(\Omega))^3} = -(\frac{1}{\mu_{\alpha}} \text{curl } u_g, \text{curl } v_h)_{(L^2(\Omega))^3} \ \forall v_h \in \mathcal{V}_h.$$

As observed in [14], although well-posed according in particular to properties (see e.g. Theorem 1. in [17]) of the used finite elements, to the assumptions on $\mu_\alpha$, $\varepsilon_\alpha$, as well as on the geometry of $\Omega$, it appears more suitable for a practical implementation to rewrite this formulation as follows.

Find $U_h \in \mathcal{H}_h$ such that:

$$\left\{ \begin{array}{l} (\frac{1}{\mu_{\alpha}} \text{curl } U_h, \text{curl } v_h)_{(L^2(\Omega))^3} = -(\frac{1}{\mu_{\alpha}} \text{curl } u_g, \text{curl } v_h)_{(L^2(\Omega))^3} \quad \forall v_h \in \mathcal{H}_h, \\ (\varepsilon_{\alpha} U_h, \text{grad } \psi_h)_{(L^2(\Omega))^3} = 0 \quad \forall \psi_h \in \Psi_h. \end{array} \right. \qquad (24)$$

Since $U_h$ represents an approximation of the vector field $U_\alpha$ satisfying (10), we introduce, following (9), the vector field $\mathcal{E}_h$ as being the discrete field associated with $\mathcal{E}_\alpha$ subject to (6), $U_h := \mathcal{E}_h + \text{grad } q$.
Next, by inserting this expression of $U_h$ into (24), and taking into account the definition of the scalar potential $q$ from (8), we now deal with a discrete formulation associated with (6) and well-posed following [14].

Find $\mathcal{E}_h \in \mathcal{H}_h$ such that:

$$ \left\{ \begin{array}{l} \left( \frac{1}{\mu_{\alpha}} \operatorname{curl} \mathcal{E}_h, \operatorname{curl} v_h \right)_{(L^2(\Omega))^3} = - \left( \frac{1}{\mu_{\alpha}} \operatorname{curl} u_g, \operatorname{curl} v_h \right)_{(L^2(\Omega))^3} \quad \forall v_h \in \mathcal{H}_h, \\ \left( \varepsilon_{\alpha} \mathcal{E}_h, \operatorname{grad} \psi_h \right)_{(L^2(\Omega))^3} = - \left( \varepsilon_{\alpha} u_g, \operatorname{grad} \psi_h \right)_{(L^2(\Omega))^3} \quad \forall \psi_h \in \Psi_h. \end{array} \right. \quad (25) $$

In numerical simulations, we will be concerned with the formulation (25) allowing us to approximate the electric field $E_\alpha$, through (7), by considering the vector field

$$ E_{\alpha}^{h} := \mathcal{E}_{h} + u_{g}. \quad (26) $$

In the numerical approximation of $E_\alpha$ and therefore in the computation of the discrete field $\mathcal{E}_h$, the scalar potential $q$ does not intervene.

For the evaluation of the integral terms of the formulation (25), we use a numerical integration method of order 2. A rectangular linear system results from this formulation and we are concerned in computations with the determination of a least squares solution obtained as in [14] with the help of a preconditioning technique based on an Incomplete Modified Gram-Schmidt (IMGS) factorization.

# 4 Localization Procedures

We describe here two localization procedures of imperfections, that result from the combination of the formula (15) with a Current Projection method or an Inverse Fourier method.

## 4.1 Procedure based on a Current Projection method

In this part, we present a localization procedure which can be used only in the context where the domain contains a single imperfection. Our aim in this context is to determine the center of the imperfection. Let us denote by $p = (p_1, p_2, p_3)^T$ this center, by $M$ the “rescaled” polarization tensor $(\frac{\mu_0}{\mu_1} - 1)M^1(\frac{\mu_0}{\mu_1})$, and by $N$ the “rescaled” polarization tensor $(\frac{\varepsilon_0}{\varepsilon_1} - 1)M^1(\frac{\varepsilon_0}{\varepsilon_1})$ of this imperfection. By neglecting the asymptotically small remainder term in (15), it follows that:

$$ \Gamma := \int_{\partial \Omega} \operatorname{curl} E_{\alpha} \times \nu \cdot w \, d\sigma - \int_{\partial \Omega} \operatorname{curl} w \times \nu \cdot (\nu \times g) \, d\sigma - \int_{\partial \Omega} \mu_0(\varepsilon_0 E_{\alpha}) \cdot \nu \operatorname{div}(w) \, d\sigma \quad (27) \\ \approx \alpha^3 (M \operatorname{curl} E_0(p)) \cdot \operatorname{curl} w(p) + \alpha^3 \mu_0 \varepsilon_0 (N E_0(p)) \cdot \operatorname{grad}(\operatorname{div}(w))(p), $$

where $g = E_\alpha \times \nu$, a datum also prescribed in (4), and $w$ is any smooth vector field satisfying (13).

Let us apply different currents for $g$ that correspond to the following background vector potentials:

$$ E_0^{(1)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ 0 \\ x_2 \end{pmatrix}, \quad E_0^{(2)}(x_1, x_2, x_3) = \begin{pmatrix} x_3 \\ 0 \\ 0 \end{pmatrix}, \quad E_0^{(3)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ x_1 \\ 0 \end{pmatrix}. $$
For each current $g^{(i)} = E_0^{(i)} \times \nu$, $1 \le i \le 3$, we use $g := g^{(i)}$ in (5) and compute through (25)-(26) the corresponding discrete electric field denoted here by $E_{\alpha,(i)}^h$, in association with $E_\alpha$. Next, we consider the test vector fields

$$w^{(1)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ 0 \\ x_2 \end{pmatrix}, \quad w^{(2)}(x_1, x_2, x_3) = \begin{pmatrix} x_3 \\ 0 \\ 0 \end{pmatrix},$$

$$w^{(3)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ x_1 \\ 0 \end{pmatrix}, \quad w^{(4)}(x_1, x_2, x_3) = \begin{pmatrix} x_2 x_3 \\ -x_1 x_3 \\ 0 \end{pmatrix},$$

in order to evaluate from the left-hand side of (27) the terms $\Gamma^{(j,i)}$, $1 \le j \le 4$,

$$\Gamma^{(j,i)} := \int_{\partial\Omega} \mathrm{curl} E_{\alpha,(i)}^h \times \nu \cdot w^{(j)} d\sigma - \int_{\partial\Omega} \mathrm{curl} w^{(j)} \times \nu \cdot (\nu \times g^{(i)}) d\sigma - \int_{\partial\Omega} \mu_0(\epsilon_0 E_{\alpha,(i)}^h) \cdot \nu \mathrm{div}(w^{(j)}) d\sigma, \\ \text{which take then the following form,}$$

$$\Gamma^{(j,i)} = \int_{\partial\Omega} \mathrm{curl} E_{\alpha,(i)}^h \times \nu \cdot w^{(j)} d\sigma - \int_{\partial\Omega} \mathrm{curl} w^{(j)} \times \nu \cdot (\nu \times g^{(i)}) d\sigma,$$

since each of these test vector fields is divergence-free. Each term $\Gamma^{(j,i)}$, $1 \le i \le 3$, $1 \le j \le 4$, is called the *measurement*, and its evaluation is achieved with the aid of a numerical integration method of order 2. From (27), it results that:

$$\Gamma^{(j,i)} \approx \alpha^3 M_{ji}, \quad \text{for } 1 \le i, j \le 3, \qquad (28)$$

after setting $M = (M_{pq})_{1 \le p, q \le 3}$. Once an approximation of the rescaled tensor $\alpha^3 M$ is thus obtained, we consider from (27) the following “linear system”,

$$\left\{ \begin{array}{l} \Gamma^{(4,1)} \approx \alpha^3 M_{11}p_1 + \alpha^3 M_{21}p_2 - 2\alpha^3 M_{31}p_3, \\ \Gamma^{(4,2)} \approx \alpha^3 M_{12}p_1 + \alpha^3 M_{22}p_2 - 2\alpha^3 M_{32}p_3, \\ \Gamma^{(4,3)} \approx \alpha^3 M_{13}p_1 + \alpha^3 M_{23}p_2 - 2\alpha^3 M_{33}p_3, \end{array} \right. \qquad (29)$$

which allows us to determine in a unique way the unknown $(p_1, p_2, p_3)^T$. This will always be possible since the matrix of the “system” (29) is invertible as its determinant is equal to $-2\det(\alpha^3 M)$ and $\alpha^3 M$ is invertible; of course, the polarization tensor $M^1(\mu_0/\mu_1)$ being positive definite.

The operated choices for the boundary currents and the test vector fields, as well as the resulting boundary measurements, are the same as those appearing in [14]. In conformity with Remark 2.2, we are thus concerned above, as in [14], with the same treatment then devoted in particular to the context where the single imperfection is purely magnetic.

Let us now consider the configuration where the single imperfection is purely electric. We apply different currents for $g$ that correspond to the following background vector potentials:

$$E_0^{(1)}(x_1, x_2, x_3) = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}, \quad E_0^{(2)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}, \quad E_0^{(3)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix},$$
and use, following (13) and (23), of course for $\mu_0\epsilon_0 = 1$, the test vector fields

$$w^{(1)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ x_1 x_2 \\ 0 \end{pmatrix}, \quad w^{(2)}(x_1, x_2, x_3) = \begin{pmatrix} 0 \\ 0 \\ x_2 x_3 \end{pmatrix},$$

$$w^{(3)}(x_1, x_2, x_3) = \begin{pmatrix} x_1 x_3 \\ 0 \\ 0 \end{pmatrix}, \quad w^{(4)}(x_1, x_2, x_3) = \begin{pmatrix} x_1 x_2 x_3 \\ x_1 x_2 x_3 \\ x_1 x_2 x_3 \end{pmatrix}.$$

Also for each current $g^{(i)} = E_0^{(i)} \times \nu$, $1 \le i \le 3$, we take $g := g^{(i)}$ in (5) and compute from (25)-(26) the corresponding discrete electric field $E_{\alpha,(i)}^h$, in association with $E_\alpha$. Next, we evaluate from the left-hand side of (27) the terms $\Gamma_{(j,i)}$, $1 \le j \le 4$,

$$\Gamma_{(j,i)} := \int_{\partial\Omega} \mathrm{curl} E_{\alpha,(i)}^h \times \nu \cdot w^{(j)} d\sigma - \int_{\partial\Omega} \mathrm{curl} w^{(j)} \times \nu \cdot (\nu \times g^{(i)}) d\sigma - \int_{\partial\Omega} \mu_0(\epsilon_0 E_{\alpha,(i)}^h) \cdot \nu \mathrm{div}(w^{(j)}) d\sigma,$$

with the aid of a numerical integration method of order 2. It derives from (27) that:

$$\Gamma_{(j,i)} \approx \alpha^3 \mu_0 \epsilon_0 N_{ji}, \quad \text{for } 1 \le i, j \le 3, \qquad (30)$$

after setting $N = (N_{pq})_{1 \le p,q \le 3}$. Once an approximation of the rescaled tensor $\alpha^3\mu_0\epsilon_0 N$ is thus determined, we consider from (27) the following “linear system”,

$$\left\{ \begin{array}{l} \Gamma_{(4,1)} \approx \alpha^3 \mu_0 \epsilon_0 N_{11}(p_2+p_3) + \alpha^3 \mu_0 \epsilon_0 N_{21}(p_1+p_3) + \alpha^3 \mu_0 \epsilon_0 N_{31}(p_1+p_2), \\ \Gamma_{(4,2)} \approx \alpha^3 \mu_0 \epsilon_0 N_{12}(p_2+p_3) + \alpha^3 \mu_0 \epsilon_0 N_{22}(p_1+p_3) + \alpha^3 \mu_0 \epsilon_0 N_{32}(p_1+p_2), \\ \Gamma_{(4,3)} \approx \alpha^3 \mu_0 \epsilon_0 N_{13}(p_2+p_3) + \alpha^3 \mu_0 \epsilon_0 N_{23}(p_1+p_3) + \alpha^3 \mu_0 \epsilon_0 N_{33}(p_1+p_2), \end{array} \right. \qquad (31)$$

written in a matrix form as

$$\begin{pmatrix} \Gamma_{(4,1)} \\ \Gamma_{(4,2)} \\ \Gamma_{(4,3)} \end{pmatrix} \approx \begin{pmatrix} \alpha^3 \mu_0 \epsilon_0 N_{11} & \alpha^3 \mu_0 \epsilon_0 N_{21} & \alpha^3 \mu_0 \epsilon_0 N_{31} \\ \alpha^3 \mu_0 \epsilon_0 N_{12} & \alpha^3 \mu_0 \epsilon_0 N_{22} & \alpha^3 \mu_0 \epsilon_0 N_{32} \\ \alpha^3 \mu_0 \epsilon_0 N_{13} & \alpha^3 \mu_0 \epsilon_0 N_{23} & \alpha^3 \mu_0 \epsilon_0 N_{33} \end{pmatrix} A \begin{pmatrix} p_1 \\ p_2 \\ p_3 \end{pmatrix},$$

where $A = \begin{pmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{pmatrix}$. The determination in a unique way of the unknown $(p_1, p_2, p_3)^T$ can then be achieved. The matrix of the “system” (31) is invertible as its determinant is equal to $2\det(\alpha^3\mu_0\epsilon_0 N)$ and $\alpha^3\mu_0\epsilon_0 N$ is invertible; the polarization tensor $M^1(\frac{\epsilon_0}{\epsilon_1})$ being positive definite.

Depending on the context where the imperfection is purely electric or not, we determine both an approximation of the rescaled tensor $\alpha^3\mu_0\epsilon_0 N$, or $\alpha^3 M$, and an approximation of the center of the imperfection, by considering therefore three background vector potentials and also four test vector fields.

Let us notice that an approximation of the order of magnitude of the diameter of the imperfection can be obtained, following (30) or (28), from one of the terms $\Gamma_{(i,i)}$, or $\Gamma^{(i,i)}$, $1 \le i \le 3$,
when in addition one of the diagonal coefficients of the tensor $\mu_0\varepsilon_0(\frac{\varepsilon_0}{\varepsilon_1} - 1)M^1(\frac{\varepsilon_0}{\varepsilon_1})$, or respectively $(\frac{\mu_0}{\mu_1} - 1)M^1(\frac{\mu_0}{\mu_1})$, is known.

In the configuration where the single imperfection is electromagnetic (with real permittivity), an approximation of its center can be determined from the measurements in (30)-(31) and an approximation of the order of magnitude of its diameter (from (30)) can be obtained when one of the diagonal coefficients of the tensor $\mu_0\varepsilon_0(\frac{\varepsilon_0}{\varepsilon_1} - 1)M^1(\frac{\varepsilon_0}{\varepsilon_1})$ is in addition known. Of course, as already established in [14], the measurements in (28)-(29) and the knowledge of one of the diagonal coefficients of the tensor $(\frac{\mu_0}{\mu_1} - 1)M^1(\frac{\mu_0}{\mu_1})$ lead also to the determination of similar approximations. Let us already mention that the numerical investigations that will be reported hereafter (see Subsection 5.2) deal uniquely with the context of measurements in (30)-(31) when interested in the present procedure.

## 4.2 Procedure based on an Inverse Fourier method

We are concerned in this part with a procedure devoted to the localization of a finite number $m$ ($m \ge 1$) of imperfections. This procedure makes use of the formula (15) and of the technique of Calderón [11] that reduces a localization problem to the calculation of an inverse Fourier transform.

Let us first reconsider the formula (15) as follows

$$ \Gamma := \int_{\partial\Omega} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma - \int_{\partial\Omega} \mathrm{curl} w \times \nu \cdot (\nu \times g) \, d\sigma - \int_{\partial\Omega} \mu_0(\varepsilon_0 E_{\alpha}) \cdot \nu \, \mathrm{div}(w) \, d\sigma \\ \approx \alpha^3 \sum_{j=1}^{m} \left(\frac{\mu_0}{\mu_j} - 1\right) \left[M^j\left(\frac{\mu_0}{\mu_j}\right) \mathrm{curl} E_0(z_j)\right] \cdot \mathrm{curl} w(z_j) \tag{32} \\ + \alpha^3 \sum_{j=1}^{m} \mu_0 \varepsilon_0 \left(\frac{\varepsilon_0}{\varepsilon_j} - 1\right) \left[M^j\left(\frac{\varepsilon_0}{\varepsilon_j}\right) E_0(z_j)\right] \cdot \mathrm{grad}(\mathrm{div}(w))(z_j), $$

where $g = E_{\alpha} \times \nu$.

For an arbitrary $\eta \in \mathbb{R}^3$, let us define $\beta$ and $\zeta$ in $\mathbb{R}^3$ such that:

$$ \left\{ \begin{array}{ll} \|\beta\|^2 = 1, & \beta \cdot \eta = 0, \\ \|\zeta\|^2 = 1, & \zeta \cdot \eta = \zeta \cdot \beta = 0, \end{array} \right. $$

with $\|\cdot\|$ denoting the usual norm associated with the Hermitian product on $\mathbb{C}^3$. Let $p = \eta + \gamma\beta$, where $\gamma$ is a complex number such that $\gamma = i\|\eta\|$, i.e. $p \cdot p = 0$. We also set $q = \eta - \gamma\beta$.

Let us consider the boundary current for which the corresponding background potential is given by

$$ E_0(x) = e^{i\mathbf{p} \cdot \mathbf{x}} \zeta, $$

namely $g(x) = (e^{i\mathbf{p} \cdot \mathbf{x}}\zeta) \times \nu(x)$, and use also, in accordance with (13), the test vector field

$$ w(x) = e^{i\mathbf{q} \cdot \mathbf{x}} \zeta. $$
It follows from (32) with these considerations that

$$
\begin{align*}
\Gamma &= \int_{\partial \Omega_m} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma - \int_{\partial \Omega} \mathrm{curl} w \times \nu \cdot (\nu \times g) \, d\sigma \\
&\approx \alpha^3 \sum_{j=1}^{m} \left( \frac{\mu_0}{\mu_j} - 1 \right) \left[ M^j \left( \frac{\mu_0}{\mu_j} \right) \mathrm{curl} E_0(z_j) \right] \cdot \mathrm{curl} w(z_j),
\end{align*}
$$

since the test vector field is divergence-free. By viewing the measurement as a function of η, we notice
that:

$$
\Gamma(\eta) \approx \alpha^3 \sum_{j=1}^{m} \left( -\left(\frac{\mu_0}{\mu_j} - 1\right) \left[ M^j \left(\frac{\mu_0}{\mu_j}\right) ((\eta + \gamma\beta) \times \zeta) \right] \cdot ((\eta - \gamma\beta) \times \zeta) \right) e^{i2\eta \cdot z_j}. \quad (33)
$$

If all the imperfections are balls — the tensors $M^j(c)$ being accordingly of the form $m^j(c)I_3$, with $I_3$ the $3 \times 3$ identity matrix, $m^j(c)$ a scalar depending on $c$ (cf. e.g. [8]), it follows then that

$$
\Gamma(\eta) \approx \alpha^3 \sum_{j=1}^{m} \left[ -\left(\frac{\mu_0}{\mu_j} - 1\right) m^j \left(\frac{\mu_0}{\mu_j}\right) (2\|\eta\|^2) \right] e^{i2\eta \cdot z_j}. \quad (34)
$$

Since ||η||² is a polynomial in the coordinates ηᵢ, 1 ≤ i ≤ 3, of η, we notice therefore from the right-
hand side of (34) the Fourier transform of a linear combination of derivatives of order less than or
equal to 2 of delta functions centered at the points -2zⱼ, 1 ≤ j ≤ m. Otherwise, the inverse Fourier
transform of Γ(η) is such that

$$
\check{\Gamma}(x) \approx \alpha^3 \sum_{j=1}^{m} L_j(\delta_{-2z_j})(x),
$$

where $L_j$ is a second order differential operator with constant coefficients depending on $m^j(\frac{\mu_0}{\mu_j})$. A numerical Fourier inversion of a sample of measurements should efficiently pin down the $z_j$'s.

When some of the imperfections are not balls, the expression in the right-hand side of (33) is the
Fourier transform of an operator of a more complicated kind acting on delta functions centered at the
points -2zⱼ, 1 ≤ j ≤ m.

We are in fact led to certain developments of [14]; the operated choices for the boundary current
and the test vector field, as well as the resulting boundary measurement, being the same as those
appearing in [14] for this configuration of imperfections involving magnetic contrasts.

In view of the configuration of purely electric imperfections, we consider the boundary current for
which the corresponding background potential is given by

$$
E_0(x) = e^{ip \cdot x} p,
$$

namely $g(x) = (e^{ip \cdot x} p) \times \nu(x)$, and use the test vector field

$$
w(x) = e^{i\mathbf{q} \cdot \mathbf{x}} \eta,
$$

under the condition

$$
\mu_0 \epsilon_0 = 1, \tag{35}
$$
since $w$ must satisfy (13). With these considerations, where in particular $\mathrm{curl} E_0$ vanishes, it follows
from (32) that

$$
\begin{align*}
\Gamma &= \int_{\partial\Omega} \mathrm{curl} E_{\alpha} \times \nu \cdot w \, d\sigma - \int_{\partial\Omega} \mathrm{curl} w \times \nu \cdot (\nu \times g) \, d\sigma - \int_{\partial\Omega} \mu_0(\varepsilon_0 E_{\alpha}) \cdot \nu \, \mathrm{div}(w) \, d\sigma \\
&\approx \alpha^3 \sum_{j=1}^{m} \mu_0 \varepsilon_0 \left(\frac{\varepsilon_0}{\varepsilon_j} - 1\right) \left[M^j\left(\frac{\varepsilon_0}{\varepsilon_j}\right) E_0(z_j)\right] \cdot \mathrm{grad}(\mathrm{div}(w))(z_j),
\end{align*}
$$

and hence,

$$
\Gamma \approx \alpha^3 \sum_{j=1}^{m} \mu_0 \varepsilon_0 \left(\frac{\varepsilon_0}{\varepsilon_j} - 1\right) \left[M^j \left(\frac{\varepsilon_0}{\varepsilon_j}\right) \left(e^{i\mathbf{p} \cdot z_j} p\right)\right] \cdot \left(- (q \cdot \eta) e^{i\mathbf{q} \cdot z_j} q\right).
$$

Let us view the measurement as a function of $\eta$:

$$
\Gamma(\eta) \approx \alpha^3 \sum_{j=1}^{m} \left( \mu_0 \varepsilon_0 \left( \frac{\varepsilon_0}{\varepsilon_j} - 1 \right) \left[ M^j \left( \frac{\varepsilon_0}{\varepsilon_j} \right) (\eta + \gamma\beta) \right] \cdot \left( - (\eta \cdot \eta)(\eta - \gamma\beta) \right) \right) e^{i2\eta \cdot z_j}. \quad (36)
$$

If all the imperfections are balls, it follows that

$$
\Gamma(\eta) \approx \alpha^3 \sum_{j=1}^{m} \left[ \mu_0 \varepsilon_0 \left(\frac{\varepsilon_0}{\varepsilon_j} - 1\right) m^j \left(\frac{\varepsilon_0}{\varepsilon_j}\right) (-2\|\eta\|^4) \right] e^{i2\eta \cdot z_j}, \quad (37)
$$

where $m^j(\frac{\varepsilon_0}{\varepsilon_j})$ is a scalar depending on $\frac{\varepsilon_0}{\varepsilon_j}$. The expression in the right-hand side of (37) is then
the Fourier transform of a linear combination of derivatives of order less than or equal to 4 of delta
functions centered at the points $-2z_j$, $1 \le j \le m$. Otherwise, the inverse Fourier transform of $\Gamma(\eta)$ is
such that

$$
\check{\Gamma}(x) \approx \alpha^3 \sum_{j=1}^{m} \mathcal{L}_j(\delta_{-2z_j})(x),
$$

where $\mathcal{L}_j$ is a fourth order differential operator with constant coefficients depending on $m^j(\frac{\varepsilon_0}{\varepsilon_j})$.

When some of the imperfections are not balls, the expression in the right-hand side of (36), where appears the term of the type $\|\eta\|^2 \mathcal{Q}(\eta) e^{i2\eta \cdot z_j}$ with $\mathcal{Q}(\eta)$ a polynomial function in $\eta$ of degree less than or equal to 2, corresponds to the Fourier transform of an operator of a more complicated kind acting on delta functions centered at the points $-2z_j$, $1 \le j \le m$.

The inversion process, then suitable both for configurations of imperfections that are purely electric or not, consists, from (36) or (33), of sampling $\Gamma(\eta)$ at some discrete set of points and then evaluating the discrete inverse Fourier transform of the corresponding sample. Namely, for each point $\eta$ (of such a set), we consider, depending on the context, $g(x) = (e^{i(\eta+\gamma\beta)\cdot x}\zeta) \times \nu(x)$ or $g(x) = (e^{i(\eta+\gamma\beta)\cdot x}(\eta + \gamma\beta)) \times \nu(x)$ as the boundary current in (5) and compute through (25)-(26) the corresponding discrete electric field, here denoted by $E_{\alpha}^{h}$. We then evaluate numerically, with the aid of an integration method of order 2, the measurement $\Gamma(\eta)$, by using of course $w(x) = e^{i(\eta-\gamma\beta)\cdot x}\zeta$ or $w(x) = e^{i(\eta-\gamma\beta)\cdot x}\eta$ as the test field in $\int_{\partial\Omega}\mathrm{curl}\,E_{\alpha}\times\nu\cdot w\,d\sigma - \int_{\partial\Omega}\mathrm{curl}\,w\times\nu\cdot(\nu\times g)\,d\sigma$ or $\int_{\partial\Omega}\mathrm{curl}\,E_{\alpha}\times\nu\cdot w\,d\sigma - \int_{\partial\Omega}\mathrm{curl}\,w\times\nu\cdot(\nu\times g)\,d\sigma - \int_{\partial\Omega}\mu_0(\varepsilon_0 E_{\alpha})\cdot\nu\,\mathrm{div}(w)\,d\sigma$ and by replacing $E_{\alpha}$ by $E_{\alpha}^{h}$ in this representation of the left-hand side of (33) or respectively of (36).
A possible way to choose a step size for sampling with respect to $\eta$ in the numerical simulations is the one introduced in [18]. It is then assumed that all the centers $z_j = (z_j^1, z_j^2, z_j^3)^T$ of the imperfections ($1 \le j \le m$) lie in a domain $[-K, K]^3$, where the bound $K$ is known, and that we are in possession of the sequence of data, measurements for $\eta = (\eta_1, \eta_2, \eta_3)^T \in [-\eta_{\max}, \eta_{\max}]^3$, on a regular grid made up of $n^3$ points with $\rho = \frac{2\eta_{\max}}{n}$ as the step size. In this process, centers $z_j, 1 \le j \le m$, are localized, after a rescaling (by $-\frac{1}{2}$), from the sequence of the modules of the terms that derive from the application of the discrete inverse Fourier transform to the sequence of data. This is achieved with at best (theoretically) a resolution of order $\frac{\pi}{2\eta_{\max}}$, and the link between $\eta_{\max}, K$ and $n$ is such that $\rho \approx \frac{1}{K}$ (see also [8]).

The present procedure aims then at localizing imperfections that are purely electric, under the background constraint (35). In view of the localization of imperfections that are purely magnetic, the procedure operates in exactly the same way as from [14] by making then use of the measurements in (33). It is also devoted to the configuration of imperfections that are electromagnetic (with real permittivity) from the measurements in (36) and in accordance with (35); of course, as already established in [14], the localization of such imperfections can also be achieved from the measurements in (33) without such a background constraint. Let us already mention that the numerical investigations that will be reported hereafter (see Subsection 5.3) deal uniquely with the context of measurements in (36), when concerned with the present procedure.

**Remark 4.1** The constraint (35) is not in fact restrictive since it remains possible to formulate the localization problem with the following normalized parameters,

$$ \mu_{\alpha}^{*} = \frac{\mu_{\alpha}}{\mu_{0}} = \begin{cases} \mu_{0}^{*} & \text{in } \Omega \setminus \overline{I_{\alpha}}, \\ \mu_{j}^{*} & \text{in } z_{j} + \alpha B_{j}, \end{cases} \quad \varepsilon_{\alpha}^{*} = \frac{\varepsilon_{\alpha}}{\varepsilon_{0}} = \begin{cases} \varepsilon_{0}^{*} & \text{in } \Omega \setminus \overline{I_{\alpha}}, \\ \varepsilon_{j}^{*} & \text{in } z_{j} + \alpha B_{j}, \end{cases} $$

where of course $\mu_j^* = \frac{\mu_j}{\mu_0}, \mu_0^* = 1, \varepsilon_j^* = \frac{\varepsilon_j}{\varepsilon_0}, 1 \le j \le m, \varepsilon_0^* = 1$, and then $\mu_0^* \varepsilon_0^* = 1$.

# 5 Numerical Simulations

We present in this part the numerical results obtained from the localization procedures previously introduced.

## 5.1 Computational configurations

Two settings for the polyhedral domain $\Omega$, having here the diameter and the shape of the unit ball, are distinctly considered: the case where $\Omega$ contains a single imperfection and when it contains multiple imperfections. In the first setting, the single imperfection is a polyhedron having the shape of a ball of center $p = (p_1, p_2, p_3)^T \in \Omega$ and of radius $\alpha$, the discretization of $\Omega$ being represented by

* $\mathcal{T}_{\alpha}^1$ when $p = (0.23, -0.31, 0.15)^T$, with $\alpha = 0.2$;

* $\mathcal{T}_{\alpha}^2$ when $p = (0.23, -0.31, 0.15)^T$, with $\alpha = 0.17$.
In the second setting, each imperfection is a polyhedron having the shape of a ball or of an ellipsoid. We represent now the discretization of $\bar{\Omega}$ by

* $\mathcal{T}_{\alpha}^{3}$ when $\Omega$ contains two ball-like shaped imperfections of centers $(0.23, -0.31, 0.15)^{T}$, $(-0.17, 0.43, -0.11)^{T}$, and of the same 'radius' $\alpha = 0.2$;

* $\mathcal{T}_{\alpha}^{4}$ when $\Omega$ contains three ball-like shaped imperfections of centers $(0.23, -0.31, 0.15)^{T}$, $(-0.17, -0.43, -0.11)^{T}$, $(-0.5, 0.25, 0.1)^{T}$, with respective 'radii' 0.18, 0.16 and 0.17. In this case, we denote by $\alpha$ the maximal radius and by $\alpha_{\min}$ the minimal radius: $\alpha = 0.18, \alpha_{\min} = 0.16$;

* $\mathcal{T}_{\alpha}^{5}$ when $\Omega$ contains three imperfections, one of which has the shape of a ball of radius 0.16 and of center $(0.23, -0.31, 0.15)^{T}$. The second one has the shape of an ellipsoid centered at $(-0.17, -0.43, -0.11)^{T}$ with semi-axes of lengths 0.16, 0.16, 0.18 in the directions $Ox, Oy, Oz$ respectively. The last imperfection is also ellipsoid-shaped, but centered at $(-0.5, 0.25, 0.1)^{T}$ with the 'semi-axes' (on $Oxy$) rotated about $Oz$ by an angle of $\frac{\pi}{4}$. The lengths of the 'semi-axes' of this imperfection are 0.16, 0.17 and 0.19. Now, $\alpha$ is the maximal value of the semi-axes lengths and the 'radius' of the first imperfection, $\alpha_{\min}$ is the minimal value of these quantities: $\alpha = 0.19, \alpha_{\min} = 0.16$.

Each one of these discretizations is associated with a conforming mesh made up of tetrahedra that takes implicitly into account the geometry of each imperfection. The resulting mesh size $h$ is here systematically smaller than the lowest of the 'radii' or 'semi-axes' lengths of the imperfections: $h < \alpha_{\min}$. In the following table, we give some characteristics of the mesh in each one of these configurations that were also used in [14], in the localization context of imperfections not purely electric.

<table><thead><tr><th></th><th>NK</th><th>NIE</th><th>NIV</th><th>nf</th><th>ne</th><th>h</th></tr></thead><tbody><tr><td>T<sup>1</sup><sub>α</sub></td><td>45101</td><td>49906</td><td>6643</td><td>3678</td><td>5517</td><td>0.17725</td></tr><tr><td>T<sup>2</sup><sub>α</sub></td><td>54368</td><td>60753</td><td>8215</td><td>3662</td><td>5493</td><td>0.15717</td></tr><tr><td>T<sup>3</sup><sub>α</sub></td><td>55847</td><td>62386</td><td>8425</td><td>3774</td><td>5661</td><td>0.15718</td></tr><tr><td>T<sup>4</sup><sub>α</sub></td><td>64765</td><td>72662</td><td>9872</td><td>3952</td><td>5928</td><td>0.14810</td></tr><tr><td>T<sup>5</sup><sub>α</sub></td><td>74093</td><td>83334</td><td>11363</td><td>4246</td><td>6369</td><td>0.14534</td></tr></tbody></table>

We have denoted by NK, NIE, NIV the number of tetrahedra, internal edges and internal vertices respectively. Also, nf, ne are respectively the number of boundary faces and boundary edges.

## 5.2 From the procedure based on the Current Projection method

By using the procedure described in Subsection 4.1, we present here numerical results related to the single imperfection setting. Since this procedure operates in particular, with (28)-(29), in exactly the same way as in [14] when, in particular, the imperfection is purely magnetic, we uniquely consider here
the purely electric and electromagnetic (with real permittivity) contexts in association with (30)-(31).
Hereafter, $\mu_0 = \varepsilon_0 = 1$.

While the measurements in (30)-(31) are known, namely for $\varepsilon_1 \neq \varepsilon_0$, $\varepsilon_1 > 0$, the reconstruction of the center of the imperfection will always be possible. For the reconstruction of the 'radius' $\alpha$, the measurements in (30) are needed, and our simulations require of course the numerical evaluation of one of the diagonal coefficients of the polarization tensor $M^1(\frac{\varepsilon_0}{\varepsilon_1})$ associated with the imperfection (such an evaluation is done as indicated in [8]).

The representations in Figs. 1-2, after considering in simulations the configurations $\mathcal{T}_\alpha^1, \mathcal{T}_\alpha^2$, concern the relative error on the center $p$ of the imperfection, $\frac{|p - p_h|_{\mathbb{R}^3}}{|p|_{\mathbb{R}^3}}$ (where $|\cdot|_{\mathbb{R}^3}$ denotes the infinity norm on $\mathbb{R}^3$), with $p_h$ the center of the localized imperfection. It is observed with respect to values used for $\varepsilon_1$ ($\varepsilon_1 > \varepsilon_0$) that this relative error presents an asymptotic behavior and in particular a deterioration for large values of $\varepsilon_1$. The relative error on the center obtained from $\mathcal{T}_\alpha^2$ is slightly smaller than the one resulting from $\mathcal{T}_\alpha^1$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Semi-log representation of the relative error on the center (percentage) with respect to some values of $\varepsilon_1$, when $\mu_1 = 1$, from $\mathcal{T}_\alpha^1$ (—) and $\mathcal{T}_\alpha^2$ (---).</figcaption></figure>

As indicates Fig. 2 in association with $\mathcal{T}_\alpha^1$, for different considerations of $\mu_1$, this relative error keeps the same behavior with respect to $\varepsilon_1$.

Fig. 3 presents the cross-sections at $x = p_1 = 0.23$, $y = p_2 = -0.31$ and $z = p_3 = 0.15$, of the localized imperfection resulting from $\mathcal{T}_\alpha^1$ when $\mu_1 = 1, \varepsilon_1 = 5$.

The experiments associated with Figs. 4-5 are also related to $\mathcal{T}_\alpha^1$ but deal with the contexts of higher electric contrast and electromagnetic contrast.

The results represented in Figs. 6-7 are obtained from the configuration $\mathcal{T}_\alpha^2$.

The experiment associated with Fig. 6 considers a purely electric contrast whereas the one associated with Fig. 7 is performed with an electromagnetic contrast.

The evaluation (based on $\mathcal{T}_\alpha^1, \mathcal{T}_\alpha^2$) with good numerical accuracy of the measurements required in (30)-(31) has contributed to obtain such concrete results.

By considering other values of $\mu_1$, and of $\varepsilon_1$ (where certain are associated with Fig. 2), results similar to those represented in Figs. 3-7 arise also from simulations with $\mathcal{T}_\alpha^1, \mathcal{T}_\alpha^2$.
<figure><img src="image_4.png" /><figcaption>Figure 2: Semi-log representation of the relative error on the center (percentage) with respect to some values of $\varepsilon_1$, for $\mu_1 = 1$ (—) and $\mu_1 = 5$ (---), from $\mathcal{T}_\alpha^1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Respective cross-sections at $x = p_1$, $y = p_2$ and $z = p_3$, from $\mathcal{T}_\alpha^1$ and with $\mu_1 = 1$, $\varepsilon_1 = 5$. Superposition of the original imperfection (—) whose center is marked by “+”, and of the localized imperfection (---) with its center marked by “x”.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Same as Fig. 3 except with $\mu_1 = 1$, $\varepsilon_1 = 250$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Same as Fig. 3 except with $\mu_1 = 3, \varepsilon_1 = 5$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Respective cross-sections at $x = p_1$, $y = p_2$ and $z = p_3$, from $\mathcal{T}_\alpha^2$ and with $\mu_1 = 1, \varepsilon_1 = 3$. Superposition of the original imperfection (—) whose center is marked by “+”, and of the localized imperfection (---) with its center marked by “x”.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Same as Fig. 6 except with $\mu_1 = 3, \varepsilon_1 = 5$.</figcaption></figure>
It follows that the procedure described in Subsection 4.1 leads, with (30)-(31), to effective localizations in purely electric and electromagnetic (with real permittivity) configurations.

## 5.3 From the procedure based on an Inverse Fourier method

We present here localization results obtained in the single imperfection context ($m = 1$) as well as from the multiple imperfections setting ($m > 1$), by making use of the procedure described in Subsection 4.2. Since this procedure operates, with (33), in the same way as in [14] when in particular the imperfections are purely magnetic, we uniquely consider here the purely electric and electromagnetic (with real permittivity) configurations in association with (36).

Recalling that the outcome of the procedure is the sequence of modules of the terms that derive from the application of the inverse Fourier transform reported in Subsection 4.2, the presentation of results will consist of representing, after a rescaling by $-\frac{1}{2}$, contour-plots based on this sequence, additionally enriched by a usual linear interpolation process. In order to overcome numerical instabilities that could be induced by the consideration of an arbitrarily large value of $\eta_{max}$ in the stage that consists of building the data of the procedure ($n^3$ numerical measurements), a cutoff process is incorporated (as in [14]). Hence, a threshold $\eta_*$ (independent of the centers and shapes of imperfections as well as of $\mu_\alpha$, $\varepsilon_\alpha$) is considered such that for $\eta$ with $\|\eta\| > \|(\eta_*, \eta_*, \eta_*)^T\|$, the associated datum is set equal to 0. All our numerical experiments are then described with respect to $\eta_{max}$, $n$ and $\eta_*$, in addition to the physical parameters $\mu_\alpha$, $\varepsilon_\alpha$. Furthermore, all the computations are performed in double precision arithmetic in each experiment (as it was also the case in Subsection 5.2). Systematically $\mu_0 = \varepsilon_0 = 1$, and suitable values for $\eta_*$ will result from computations.

At first, we present the results regarding the single imperfection context and compare some of these results with those of Subsection 5.2. Figs. 8-12 show the results obtained from the configurations $T_\alpha^1$ and $T_\alpha^2$. In the simulations, $\eta_{max} = 10$ and it has been considered $\rho = 2$ as the step size for sampling, i.e., $n = 10$. The expected order of resolution is then $\frac{\pi}{2\eta_{max}} \approx 0.157$. The mentioned cutoff process is required since this value of $\eta_{max}$ appears numerically large as observed from simulations.

As indicated Figs. 8-9, the localization from $T_\alpha^1$ is successfully achieved in purely electric and electromagnetic settings.

<figure><img src="image_10.png" /><figcaption>Figure 8: Contour-plot views respectively from the x-direction, the y-direction and the z-direction, by using $T_\alpha^1$. Here, $\mu_1 = 1$, $\varepsilon_1 = 5$, $\eta_{max} = 10$, $n = 10$ and $\eta_* = 3$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Same as for Fig. 8 except with $\mu_1 = 10, \varepsilon_1 = 5$.</figcaption></figure>

We notice from the experiments associated with Figs. 10-12 that the localization of a smaller imperfection is also achieved with good numerical accuracy, according to the fixed order of resolution.

<figure><img src="image_12.png" /><figcaption>Figure 10: Contour-plot views respectively from the x-direction, the y-direction and the z-direction, by using $T_\alpha^2$. Here, $\mu_1 = 1, \varepsilon_1 = 3, \eta_{max} = 10, n = 10$ and $\eta_\star = 3$.</figcaption></figure>

The experiment associated with Fig. 10 considers a purely electric contrast whereas the ones related to Figs. 11-12 are performed with electromagnetic contrasts.

<figure><img src="image_13.png" /><figcaption>Figure 11: Same as for Fig. 10 except with $\mu_1 = 10, \varepsilon_1 = 5$.</figcaption></figure>

By comparing the results presented here (see Figs. 8, 10 and 12) with those of the previous subsection (see Figs. 3, 6 and 7), we notice in this single imperfection context that the present localization procedure is less efficient than the one based on the Current Projection method. Namely, the required number of numerical measurements is very large and the obtained results are less accurate.
<figure><img src="image_14.png" /><figcaption>Figure 12: Same as for Fig. 10 except with $\mu_1 = 3, \varepsilon_1 = 5$.</figcaption></figure>

We keep however great interest in the present localization procedure, since devoted to the multiple imperfections context and hence potentially suitable for extended applications. We are thus also concerned with an inspection of the noise impact on the procedure.

When the numerical measurements are corrupted with random additive noise generated from the Gaussian normal distribution, we summarize through Figs. 13-15 some of the results that are obtained, for beginning, in the single imperfection context. For low noise levels, the localization is stable. The results start to become erroneous from incorporations of more than 20% of noise, as it can be observed for instance from Fig. 15.

<figure><img src="image_15.png" /><figcaption>Figure 13: Same as for Fig. 10 and incorporating 15% of random noise.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Same as for Fig. 10 and incorporating 20% of random noise.</figcaption></figure>

Let us consider now the multiple imperfections context. We use successively $\rho = 2, \frac{5}{4}$ and expect
<figure><img src="image_17.png" /><figcaption>Figure 15: Same as for Fig. 10 and incorporating 25% of random noise.</figcaption></figure>

as before the same order of resolution. The suited values for $\eta_*$ shall result again from simulations. In the presentation of results, when the z-direction, for instance, will be concerned, the software used for the present postprocessing will display, besides contour-plot obtained on the plane *Oxy*, horizontal and vertical segments for which the intersections correspond to centers of the original imperfections viewed on *Oxy*.

Figs. 16-18 show results obtained from the configuration $\mathcal{T}_{\alpha}^{3}$ and by fixing $\rho = 2$. The same number of measurements as in the single imperfection context is considered. The experiment associated with Fig. 16 concerns the localization of purely electric imperfections.

<figure><img src="image_18.png" /><figcaption>Figure 16: Contour-plot views respectively from the x-direction, the y-direction and the z-direction, by using $\mathcal{T}_{\alpha}^{3}$. Here, $\mu_{1} = \mu_{2} = 1$, $\varepsilon_{1} = \varepsilon_{2} = 5$, $\eta_{\max} = 10$, $n = 10$ and $\eta_{*} = 3$.</figcaption></figure>

Results similar to those associated with Fig. 16 are obtained from simulations by taking now $\mu_{1} = \mu_{2} = 1$, $\varepsilon_{1} = \varepsilon_{2} = 10$.

Fig. 17 presents the results obtained in the context of two imperfections one of which is purely electric and the other is electromagnetic.

When from $\mathcal{T}_{\alpha}^{3}$ the numerical measurements are corrupted with random additive noise, as indicates also the experiment associated with Fig. 18 and concerning electromagnetic imperfections, the localization is quite stable for low noise levels.

Regarding the experiments performed from $\mathcal{T}_{\alpha}^{4}$ and $\mathcal{T}_{\alpha}^{5}$, we use for $\rho$ a unique value, $\rho = \frac{5}{4}$, by taking into account that these configurations are associated with the same physical region of interest in the procedure. A larger value of $n$ than previously is considered though expecting the same order
<figure><img src="image_19.png" /><figcaption>Figure 17: Same as for Fig. 16 except with $\mu_1 = 1$, $\varepsilon_1 = 3$, $\mu_2 = 10$, $\varepsilon_2 = 5$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18: Same as for Fig. 16 except with $\mu_1 = 10$, $\varepsilon_1 = 5$, $\mu_2 = 5$, $\varepsilon_2 = 10$, and incorporating 5% of random noise.</figcaption></figure>

of resolution.

<figure><img src="image_21.png" /><figcaption>Figure 19: Contour-plot views respectively from the x-direction, the y-direction and the z-direction, by using $T_\alpha^4$. Here, $\mu_j = 1$, $\varepsilon_j = 5$ ($1 \le j \le 3$), $\eta_{\max} = 10$, $n = 16$ and $\eta_\star = 2.5$.</figcaption></figure>

Each of the experiments associated with Figs. 19-21 concerns the localization of imperfections that are purely electric or else electromagnetic. Similar results arise from simulations by using the parameters values of the experiment associated with Fig. 19 and by taking now $\mu_j = 5$, $\varepsilon_j = 5$ ($1 \le j \le 3$), as well as those of the experiment related to Fig. 21 but with $\mu_j = 10$, $\varepsilon_j = 5$ ($1 \le j \le 3$).

The results represented in Figs. 22-27 are related to the localization in the presence of purely electric and electromagnetic imperfections, or in the unique presence of purely electric imperfections.
<figure><img src="image_22.png" /><figcaption>Figure 20: Contour-plot views respectively from the x-direction, the y-direction and the z-direction, by using $T_\alpha^5$. Here, $\mu_1 = 1$, $\varepsilon_1 = 3$, $\mu_2 = 1$, $\varepsilon_2 = 5$, $\mu_3 = 1$, $\varepsilon_3 = 10$, $\eta_{\max} = 10$, $n = 16$ and $\eta_\star = 3$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 21: Same as for Fig. 20 except with $\mu_j = 5$, $\varepsilon_j = 10$ ($1 \le j \le 3$).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 22: Same as for Fig. 19 except with $\mu_1 = 5$, $\varepsilon_1 = 3$, $\mu_2 = 1$, $\varepsilon_2 = 5$, $\mu_3 = 5$, $\varepsilon_3 = 5$, and $\eta_\star = 3$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 23: Same as for Fig. 20 except with $\mu_1 = 1$, $\varepsilon_1 = 3$, $\mu_2 = 1$, $\varepsilon_2 = 5$, $\mu_3 = 10$, $\varepsilon_3 = 5$.</figcaption></figure>
When the numerical measurements are corrupted with random additive noise, we summarize through Figs. 24-27 some of the results based on $T_α^5$. For low noise levels, the localization is quite stable. The results start to become erroneous from incorporations of more than 15% of noise, as it can be observed for instance from Figs. 26 and 27.

<figure><img src="image_26.png" /><figcaption>Figure 24: Same as for Fig. 20 and incorporating 10% of random noise.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 25: Same as for Fig. 20 and incorporating 15% of random noise.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 26: Same as for Fig. 20 and incorporating 18% of random noise.</figcaption></figure>

As a first part of conclusions regarding the procedure used in this subsection, we notice, according to the fixed order of resolution, that the localization is also achieved with good numerical accuracy in the multiple imperfections context.

The localization accuracy can be improved by considering very large values for $\eta_{max}$. Such a consideration amplifies however the CPU time of the evaluation stage of numerical measurements.
<figure><img src="image_29.png" /><figcaption>Figure 27: Same as for Fig. 20 and incorporating 20% of random noise.</figcaption></figure>

Typically, even for a value of $\eta_{max}$ which is not very large, as in the previous experiments, we are concerned with a number of measurements which remains large and, even with the cutoff process of the Fourier domain that is applied, the evaluation stage of these is relatively costly (despite an efficient numerical computation associated with (25)). When compared in the single imperfection context with the procedure of Subsection 4.1, it results that the evaluation stage of measurements required by the present procedure is very expensive. We have been thus led to consider an order of resolution which is not very “fine” but appears pertinent, in order to perform numerical experiments with reasonable CPU time.

# 6 Conclusions

A framework in static regime has been presented for the localization in a three-dimensional bounded domain of purely electric imperfections, of purely magnetic imperfections and of electromagnetic imperfections (with real electric permittivity). This framework is based on a perturbation formula that results from a limit model in electric field, in combination with suited inversion algorithms. Two localization procedures using boundary measurements as data are thus described and numerically investigated; namely, the procedure defined from a Current Projection algorithm for the single imperfection context and the procedure defined from an Inverse Fourier method for the multiple imperfections context. Each of these procedures operating in particular in the same way as in [14] for configurations of purely magnetic imperfections while the measurements in (28)-(29) or (33) are required, the numerical experiments have therefore been essentially devoted here to purely electric settings and to electromagnetic (with real permittivity) settings by making use of the measurements in (30)-(31) or (36).

We notice that the procedure based on a Current Projection method is the most robust in the single imperfection case, whereas the second procedure seems well-suited to the context of large numbers of imperfections. Experiments with such numbers of imperfections were not performed since, according to configuration hypotheses, the imperfections should then be much smaller than those considered here, and the size of the rectangular system resulting from (25) would be in that case too large, prohibiting thus flexible numerical simulations because of the excessive memory storage and CPU time (evaluation
of measurements). However, with the help for instance of a reduced meshes process [15] in simulations
for achieving the evaluation stage of measurements, experiments based on the presented approach and
involving then much smaller imperfections shall be considered in a forthcoming study.

**Acknowledgments.** This work was partly supported by the IMMELEXT Project of the University of Henri Poincaré — Nancy-I and of the Regional Council of Lorraine in France.

References

[1] H. Ammari, E. Iakovleva, D. Lesselier and G. Perrusson, MUSIC-type electromagnetic imaging of a collection of small three-dimensional inclusions, *SIAM J. Sci. Comput.*, **29** (2007), 674-709.

[2] H. Ammari and H. Kang, Reconstruction of Small Inhomogeneities from Boundary Measurements, Lecture Notes in Mathematics, v. 1846, Springer-Verlag, Berlin, 2004.

[3] H. Ammari and H. Kang, Polarization and Moment Tensors with Application to Inverse Problems and Effective Medium Theory, Applied Mathematical Sciences, v. 162, Springer Science, New York, 2007.

[4] H. Ammari and A. Khelifi, Electromagnetic scattering by small dielectric inhomogeneities, *J. Math. Pures Appl.*, **82** (2003), 749-842.

[5] H. Ammari, S. Moskow and M.S. Vogelius, Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter, *ESAIM Control Optim. Calc. Var.*, **9** (2003), 49-66.

[6] H. Ammari, M. Vogelius and D. Volkov, Asymptotic formulas for perturbations in the electro-magnetic fields due to the presence of imperfections of small diameter II. The full Maxwell equations, *J. Math. Pures Appl.*, **80** (2001), 769-814.

[7] H. Ammari and D. Volkov, Asymptotic formulas for perturbations in the eigenfrequencies of the full Maxwell equations due to the presence of imperfections of small diameter, *Asymptot. Anal.*, **30** (2002), 331-350.

[8] M. Asch and S.M. Mefire, Numerical localization of electromagnetic imperfections from a perturbation formula in three dimensions, *J. Comput. Math.*, **26** (2008), 149-195.

[9] M. Asch and S.M. Mefire, Using reduced meshes for simulations of the localization of small electromagnetic inhomogeneities in a 3D bounded domain, *Int. J. Numer. Anal. Model.*, **6** (2009), 50-88.

[10] M. Brühl and M. Hanke, Numerical implementation of two noniterative methods for locating inclusions by impedance tomography, *Inverse Probl.*, **16** (2000), 1029-1042.

[11] A.P. Calderón, On an inverse boundary value problem, *Seminar on Numerical Analysis and its Applications to Continuum Physics, Soc. Brasileira de Matemática, Rio de Janeiro*, (1980), 65-73.
[12] D.J. Cedio-Fengya, S. Moskow and M.S. Vogelius, Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction, *Inverse Problems*, **14** (1998), 553-595.

[13] J. Laminie and S.M. Mefire, Three-dimensional computation of a magnetic field by mixed finite elements and boundary elements, *Appl. Numer. Math.*, **35** (2000), 221-244.

[14] S.M. Mefire, Three-dimensional numerical localization of imperfections based on a limit model in electric field and a limit perturbation model, *J. Comput. Math.*, **27** (2009), 495-524.

[15] S.M. Mefire, Associating a limit perturbation model in 3D with reduced meshes for simulations of the localization of certain electromagnetic inhomogeneities, *Int. J. Numer. Anal. Model.*, **7** (2010), 195-216.

[16] J.-C. Nédélec, Mixed finite elements in $R^3$, *Numer. Math.*, **35** (1980), 315-341.

[17] J.-C. Nédélec, Elements Finis mixtes incompressibles pour l'équation de Stokes dans $R^3$, *Numer. Math.*, **39** (1982), 97-112.

[18] D. Volkov, Numerical methods for locating small dielectric inhomogeneities, *Wave Motion*, **38** (2003), 189-206.