Adaptive Reconstruction for Electrical Impedance Tomography with a Piecewise Constant Conductivity

TL;DR

This paper introduces an adaptive numerical method for electrical impedance tomography that reconstructs piecewise constant conductivities using Tikhonov regularization, a posteriori error estimators, and mesh refinement, with proven convergence.

Contribution

It develops a novel adaptive algorithm combining residual-based error estimation with a Modica-Mortola penalty for improved conductivity reconstruction.

Findings

The adaptive method converges to a solution of the continuous optimality system.

Numerical examples demonstrate the effectiveness and convergence of the proposed algorithm.

The approach effectively reconstructs piecewise constant conductivities from boundary measurements.

Abstract

In this work we propose and analyze a numerical method for electrical impedance tomography of recovering a piecewise constant conductivity from boundary voltage measurements. It is based on standard Tikhonov regularization with a Modica-Mortola penalty functional and adaptive mesh refinement using suitable a posteriori error estimators of residual type that involve the state, adjoint and variational inequality in the necessary optimality condition and a separate marking strategy. We prove the convergence of the adaptive algorithm in the following sense: the sequence of discrete solutions contains a subsequence convergent to a solution of the continuous necessary optimality system. Several numerical examples are presented to illustrate the convergence behavior of the algorithm.

Equations

$$\text{mbox}SOLVE\to \text{mbox}ESTIMATE\to \text{mbox}MARK\to \text{mbox}REFINE.$$

$${\mathcal{F}}_{\epsilon }(\sigma )=\epsilon {\int }_{\Omega }|\nabla \sigma {|}^2\mathrm{d}x+\frac{1}{\epsilon }{\int }_{\Omega }W(\sigma )\mathrm{d}x,$$

$$W(s)=(s-c_0{)}^2(s-c_1{)}^2,$$

$$\underset{\sigma \in \tilde{\mathcal{A}}}{\inf }\left\{{\mathcal{J}}_{\epsilon }(\sigma )=\frac{1}{2}\Vert U(\sigma )-U^{\delta }{\Vert }^2+\frac{\tilde{\alpha }}{2}{\mathcal{F}}_{\epsilon }(\sigma )\right\},$$

$$\left\{\begin{aligned} \begin{array}[]{ll}-\nabla\cdot(\sigma\nabla u)=0&\mbox{ in }\Omega,\\[4.30554pt] u+z_{l}\sigma\frac{\partial u}{\partial n}=U_{l}&\mbox{ on }e_{l},l=1,2,\ldots,L,\\[4.30554pt] \int_{e_{l}}\sigma\frac{\partial u}{\partial n}\mathrm{d}s=I_{l}&\mbox{ for }l=1,2,\ldots,L,\\[4.30554pt] \sigma\frac{\partial u}{\partial n}=0&\mbox{ on }\partial\Omega\backslash\cup_{l=1}^{L}e_{l}.\end{array}\end{aligned}\right.$$

$$\mathcal{A}:=\{\sigma \in L^{\infty }(\Omega ):\sigma =c_0+(c_1-c_0){\chi }_{{\Omega }_1}\},$$

$$\Vert (u,U){\Vert }_{\mathbb{H}}^2=\Vert u{\Vert }_{H^1(\Omega )}^2+\Vert U{\Vert }^2.$$

$$\Vert (u,U){\Vert }_{\mathbb{H},\ast }^2=\Vert \nabla u{\Vert }_{L^2(\Omega )}^2+\sum _{l=1}^L\Vert u-U_l{\Vert }_{L^2(e_l)}^2.$$

$$a(\sigma ,(u,U),(v,V))=⟨I,V⟩\forall (v,V)\in \mathbb{H},$$

$$a(\sigma ,(u,U),(v,V))=(\sigma \nabla u,\nabla v)+\sum _{l=1}^L{z}_l^{-1}(u-U_l,v-V_l{)}_{L^2(e_l)},$$

$$\mathrm{BV}(\Omega )=\left\{v\in L^1(\Omega ):|v{|}_{\mathrm{TV}(\Omega )}<\infty \right\},$$

$$|v{|}_{\mathrm{TV}(\Omega )}:=\sup \left\{{\int }_{\Omega }v\nabla \cdot \mathbf{\varphi }\mathrm{d}x:\mathbf{\varphi }\in (C_c^1(\Omega ){)}^d,\Vert \mathbf{\varphi }(x)\Vert \le 1\right\}.$$

$$\underset{\sigma \in \mathcal{A}}{\min }\left\{\mathcal{J}(\sigma )=\frac{1}{2}\Vert U(\sigma )-U^{\delta }{\Vert }^2+\alpha |\sigma {|}_{\mathrm{TV}(\Omega )}\right\},$$

$$\mathcal{F}_{\varepsilon}(z):=\left\{\begin{array}[]{ll}\varepsilon\|\nabla z\|^{2}_{L^{2}(\Omega)}+\frac{1}{\varepsilon}\int_{\Omega}W(z)\,{\rm d}x&\mbox{if}~{}z\in H^{1}(\Omega),\\ +\infty&\mbox{otherwise},\end{array}\right.$$

$$\mathcal{F}(z):=\left\{\begin{array}[]{ll}2c_{W}|z|_{\rm TV(\Omega)}&\mbox{if}~{}z\in\mathrm{BV}(\Omega)\cap\mathcal{A},\\ +\infty&\mbox{otherwise}.\end{array}\right.$$

$$\underset{\sigma \in \tilde{\mathcal{A}}}{\inf }\left\{{\mathcal{J}}_{\epsilon }(\sigma )=\frac{1}{2}\Vert U(\sigma )-U^{\delta }{\Vert }^2+\frac{\tilde{\alpha }}{2}{\mathcal{F}}_{\epsilon }(\sigma )\right\},$$

$$\tilde{\mathcal{A}}:=\left\{\sigma \in H^1(\Omega ):c_0\le \sigma (x)\le c_1\text{mbox}a.e.x\in \Omega \right\}.$$

$$\underset{n\to \infty }{\lim }\Vert \left(u({\sigma }_n)-u({\sigma }^{\ast }),U({\sigma }_n)-U({\sigma }^{\ast })\right){\Vert }_{\mathbb{H}}=0.$$

$${\sigma }_n\rightharpoonup {\sigma }^{\ast }\text{mbox}weaklyin{H}^1(\Omega ),{\sigma }_n\to {\sigma }^{\ast }\text{mbox}in{L}^1(\Omega ),{\sigma }_n\to {\sigma }^{\ast }\text{mbox}a.e.in\Omega .$$

$${\mathcal{J}}_{\epsilon }({\sigma }^{\ast })\le \underset{n\to \infty }{\mathrm{liminf}}{\mathcal{J}}_{\epsilon }({\sigma }_n)\le \underset{n\to \infty }{\lim }{\mathcal{J}}_{\epsilon }({\sigma }_n)=m_{\epsilon }.$$

$$a(\sigma ,(p,P),(v,V))=⟨U(\sigma )-U^{\delta },V⟩\forall (v,V)\in \mathbb{H}.$$

$${\mathcal{J}}_{\epsilon }^{\prime }(\sigma )[\mu ]=\tilde{\alpha }\left[\epsilon (\nabla \sigma ,\nabla \mu )+\frac{1}{2\epsilon }(W^{\prime }(\sigma ),\mu )\right]-(\mu \nabla u(\sigma ),\nabla p(\sigma )).$$

$$\{\begin{array}{rl}& a({\sigma }^{\ast },(u^{\ast },U^{\ast }),(v,V))=⟨I,V⟩\forall (v,V)\in \mathbb{H},\\ & a({\sigma }^{\ast },(p^{\ast },P^{\ast }),(v,V))=⟨U^{\ast }-U^{\delta },V⟩\forall (v,V)\in \mathbb{H},\\ & \tilde{\alpha }\epsilon (\nabla {\sigma }^{\ast },\nabla (\mu -{\sigma }^{\ast }))+\frac{\tilde{\alpha }}{2\epsilon }(W^{\prime }({\sigma }^{\ast }),\mu -{\sigma }^{\ast })-((\mu -{\sigma }^{\ast })\nabla u^{\ast },\nabla p^{\ast })\ge 0\forall \mu \in \tilde{\mathcal{A}},\end{array}$$

$$V_{\mathcal{T}}=\left\{v\in C(\overline{\Omega }):v{|}_T\in P_1(T)\forall T\in \mathcal{T}\right\},$$

$${\tilde{\mathcal{A}}}_{\mathcal{T}}:=V_{\mathcal{T}}\cap \tilde{\mathcal{A}}.$$

$$a({\sigma }_{\mathcal{T}},(u_{\mathcal{T}},U_{\mathcal{T}}),(v_{\mathcal{T}},V))=⟨I,V⟩\forall (v_{\mathcal{T}},V)\in {\mathbb{H}}_{\mathcal{T}}.$$

$$J_{\epsilon ,\mathcal{T}}({\sigma }_{\mathcal{T}})=\frac{1}{2}\Vert U_{\mathcal{T}}({\sigma }_{\mathcal{T}})-U^{\delta }{\Vert }^2+\frac{\tilde{\alpha }}{2}{\mathcal{F}}_{\epsilon }({\sigma }_{\mathcal{T}}).$$

$$\{\begin{array}{rl}& a({\sigma }_{\mathcal{T}}^{\ast },(u_{\mathcal{T}}^{\ast },U_{\mathcal{T}}^{\ast }),(v,V))=⟨I,V⟩\forall (v,V)\in {\mathbb{H}}_{\mathcal{T}},\\ & a({\sigma }_{\mathcal{T}}^{\ast },(p_{\mathcal{T}}^{\ast },P_{\mathcal{T}}^{\ast }),(v,V))=⟨U_{\mathcal{T}}^{\ast }-U^{\delta },V⟩\forall (v,V)\in {\mathbb{H}}_{\mathcal{T}},\\ & \tilde{\alpha }\epsilon (\nabla {\sigma }_{\mathcal{T}}^{\ast },\nabla (\mu -{\sigma }_{\mathcal{T}}^{\ast }))+\frac{\tilde{\alpha }}{2\epsilon }(W^{\prime }({\sigma }_{\mathcal{T}}^{\ast }),\mu -{\sigma }_{\mathcal{T}}^{\ast })-((\mu -{\sigma }_{\mathcal{T}}^{\ast })\nabla u_{\mathcal{T}}^{\ast },\nabla p_{\mathcal{T}}^{\ast })\ge 0\forall \mu \in {\tilde{\mathcal{A}}}_{\mathcal{T}}.\end{array}$$

$$\Vert (u_{\mathcal{T}}^{\ast },U_{\mathcal{T}}^{\ast }){\Vert }_{\mathbb{H}}+\Vert (p_{\mathcal{T}}^{\ast },P_{\mathcal{T}}^{\ast }){\Vert }_{\mathbb{H}}\le c(\Vert I\Vert +\Vert U^{\delta }\Vert ),$$

$$R_{T,1}({\sigma }_{\mathcal{T}}^{\ast },u_{\mathcal{T}}^{\ast })$$

Full text

Keywords: electrical impedance tomography, piecewise constant conductivity, Modica-Mortola functional, a posteriori error estimator, adaptive finite element method, convergence analysis

1 Introduction

Electrical impedance tomography (EIT) aims at recovering the electrical conductivity distribution of an object from voltage measurements on the boundary. It has attracted much interest in medical imaging, geophysical prospecting, nondestructive evaluation and pneumatic oil pipeline conveying etc. A large number of reconstruction algorithms have been proposed; see, e.g., [37, 36, 26, 35, 31, 30, 24, 15, 21, 40, 39, 18, 3, 34, 33, 57, 28, 54, 51] for a rather incomplete list. One prominent idea underpinning many imaging algorithms is regularization, especially Tikhonov regularization [29]. In practice, they are customarily implemented using the continuous piecewise linear finite element method (FEM) on quasi-uniform meshes, due to its flexibility in handling spatially variable coefficients and general domain geometry. The convergence analysis of finite element approximations was carried out in [22, 47, 27].

In several practical applications, the physical process is accurately described by the complete electrode model (CEM) [14, 50]. It employs nonstandard boundary conditions to capture characteristics of the experiment. In particular, around the electrode edges, the boundary condition changes from the Neumann to Robin type, which, according to classical elliptic regularity theory [23], induces weak solution singularity around the electrode edges; see, e.g., [45] for an early study. Further, the low-regularity of the sought-for conductivity distribution, especially that enforced by a nonsmooth penalty, e.g., total variation, can also induce weak interior singularities of the solution. Thus, a (quasi-)uniform triangulation of the domain can be inefficient for resolving these singularities, and the discretization errors around electrode edges and internal interfaces can potentially compromise the reconstruction accuracy. These observations motivate the use of an adaptive strategy to achieve the desired accuracy in order to enhance the overall computational efficiency.

For direct problems, the mathematical theory of AFEM, including a posteriori error estimation, convergence and computational complexity, has advanced greatly [1, 44, 53, 12]. A common adaptive FEM (AFEM) consists of the following successive loops:

[TABLE]

The module ESTIMATE employs the given problem data and computed solutions to provide computable quantities on the local errors, and distinguishes different adaptive algorithms.

In this work, we develop an adaptive EIT reconstruction algorithm with a piecewise constant conductivity. In practice, the piecewise constant nature is commonly enforced by a total variation penalty. However, it is challenging for AFEM treatment (see e.g., [5] for image denoising). Thus, we take an indirect approach based on a Modica-Mortola type functional:

[TABLE]

where the constant $\varepsilon>0$ is small and $W(s):\mathbb{R}\rightarrow\mathbb{R}$ is the double-well potential, i.e.,

[TABLE]

with $c_{0},c_{1}>0$ being two known values that the conductivity $\sigma$ can take. The functional $\mathcal{F}_{\varepsilon}$ $\Gamma$-converges to the total variation semi-norm [42, 43, 2]. The corresponding regularized least-squares formulation reads

[TABLE]

where $\tilde{\alpha}>0$ is a regularization parameter; see Section 2 for further details. In this work, we propose a posteriori error estimators and an adaptive reconstruction algorithm of the form (1.1) for (1.3) based on a separate marking using three error indicators in the module MARK; see Algorithm 3.1. Further, we give a convergence analysis of the algorithm, in the sense that the sequence of state, adjoint and conductivity generated by the adaptive algorithm contains a convergent subsequence to a solution of the necessary optimality condition. The technical proof consists of two steps: Step 1 shows the subsequential convergence to a solution of a limiting problem, and Step 2 proves that the solution of the limiting problem satisfies the necessary optimality condition. The main technical challenges in the convergence analysis include the nonlinearity of the forward model, the nonconvexity of the double well potential and properly treating the variational inequality. The latter two are overcome by pointwise convergence of discrete minimizers and Lebesgue’s dominated convergence theorem, and AFEM analysis techniques for elliptic obstacle problems, respectively. The adaptive algorithm and its convergence analysis are the main contributions of this work.

Last, we situate this work in the existing literature. In recent years, several adaptive techniques, including AFEM, have been applied to the numerical resolution of inverse problems. In a series of works [6, 7, 8, 9], Beilina et al studied the AFEM in a dual weighted residual framework for parameter identification. Feng et al [20] proposed a residual-based estimator for the state, adjoint and control by assuming convexity of the cost functional and high regularity on the control. Li et al [38] derived a posteriori error estimators for recovering the flux and proved their reliability; see [55] for a plain convergence analysis. Clason et al [17] studied functional a posteriori estimators for convex regularized formulations. Recently, Jin et al [32] proposed a first AFEM for Tikhonov functional for EIT with an $H^{1}(\Omega)$ penalty, and also provided a convergence analysis. This work extends the approach in [32] to the case of a piecewise constant conductivity. There are a number of major differences between this work and [32]. First, the $H^{1}(\Omega)$ penalty in [32] facilitates deriving the a posteriori estimator on the conductivity $\sigma$, by completing the squares and suitable approximation argument, which is not directly available for the Mordica-Mortola functional $\mathcal{F}_{\varepsilon}$. Second, we develop a sharper error indicator associated with the crucial variational inequality than that in [32], by a novel constraint preserving interpolation operator [13]; see the proof of Theorem 5.5, which represents the main technical novelty of this work. Third, Algorithm 3.1 employs a separate marking for the estimators, instead of a collective marking in [32], which automatically takes care of different scalings of the estimators.

The rest of this paper is organized as follows. In Section 2, we introduce the complete electrode model, and the regularized least-squares formulation. In Section 3, we give the AFEM algorithm. In Section 4, we present extensive numerical results to illustrate its convergence and efficiency. In Section 5, we present the lengthy technical convergence analysis. Throughout, $\langle\cdot,\cdot\rangle$ and $(\cdot,\cdot)$ denote the inner product on the Euclidean space and $(L^{2}(\Omega))^{d}$, respectively, by $\|\cdot\|$ the Euclidean norm, and occasionally abuse $\langle\cdot,\cdot\rangle$ for the duality pairing between the Hilbert space $\mathbb{H}$ and its dual space. The superscript $\rm t$ denotes the transpose of a vector. The notation $c$ denotes a generic constant, which may differ at each occurrence, but it is always independent of the mesh size and other quantities of interest.

2 Regularized approach

This part describes the regularized approach for recovering piecewise constant conductivities.

2.1 Complete electrode model (CEM)

Let $\Omega$ be an open bounded domain in $\mathbb{R}^{d}$ $(d=2,3)$ with a polyhedral boundary $\partial\Omega$. We denote the set of electrodes by $\{e_{l}\}_{l=1}^{L}$, which are line segments/planar surfaces on $\partial\Omega$ and satisfy $\bar{e}_{i}\cap\bar{e}_{k}=\emptyset$ if $i\neq k$. The applied current on electrode $e_{l}$ is denoted by $I_{l}$, and the vector $I=(I_{1},\ldots,I_{L})^{\mathrm{t}}\in\mathbb{R}^{L}$ satisfies $\sum_{l=1}^{L}I_{l}=0$, i.e., $I\in\mathbb{R}_{\diamond}^{L}:=\{V\in\mathbb{R}^{L}:\sum_{l=1}^{L}V_{l}=0\}$. The electrode voltage $U=(U_{1},\ldots,U_{L})^{\mathrm{t}}$ is normalized, i.e., $U\in\mathbb{R}_{\diamond}^{L}$. Then the CEM reads: given the conductivity $\sigma$, positive contact impedances $\{z_{l}\}_{l=1}^{L}$ and input current $I\in\mathbb{R}_{\diamond}^{L}$, find $(u,U)\in H^{1}(\Omega)\otimes\mathbb{R}_{\diamond}^{L}$ such that [14, 50]

[TABLE]

The inverse problem is to recover the conductivity $\sigma$ from a noisy version $U^{\delta}$ of the electrode voltage $U(\sigma^{\dagger})$ (for the exact conductivity $\sigma^{\dagger}$) corresponding to one or multiple input currents.

Below the conductivity $\sigma$ is assumed to be piecewise constant, i.e., in the admissible set

[TABLE]

where the constants $c_{1}>c_{0}>0$ are known, $\overline{\Omega}_{1}\subset\Omega$ is an unknown open set with a Lipschitz boundary and $\chi_{\Omega_{1}}$ denotes its characteristic function. We denote by $\mathbb{H}$ the space $H^{1}(\Omega)\otimes\mathbb{R}_{\diamond}^{L}$ with its norm given by

[TABLE]

A convenient equivalent norm on the space $\mathbb{H}$ is given below.

Lemma 2.1**.**

On the space $\mathbb{H}$, the norm $\|\cdot\|_{\mathbb{H}}$ is equivalent to the norm $\|\cdot\|_{\mathbb{H},*}$ defined by

[TABLE]

The weak formulation of the model (2.1) reads [50]: find $(u,U)\in\mathbb{H}$ such that

[TABLE]

where the trilinear form $a(\sigma,(u,U),(v,V))$ on $\mathcal{A}\times\mathbb{H}\times\mathbb{H}$ is defined by

[TABLE]

where $(\cdot,\cdot)_{L^{2}(e_{l})}$ denotes the $L^{2}(e_{l})$ inner product. For any $\sigma\in\mathcal{A}$, $\{z_{l}\}_{l=1}^{L}$ and $I\in\Sigma_{\diamond}^{L}$, the existence and uniqueness of a solution $(u,U)\in\mathbb{H}$ to (2.2) follows from Lemma 2.1 and Lax-Milgram theorem.

2.2 Regularized reconstruction

For numerical reconstruction with a piecewise constant conductivity, the total variation (TV) penalty is popular. The conductivity $\sigma$ is assumed to be in the space $\mathrm{BV}(\Omega)$ of bounded variation [4, 19], i.e.,

[TABLE]

equipped with the norm $\|v\|_{\mathrm{BV}(\Omega)}=\|v\|_{L^{1}(\Omega)}+|v|_{\mathrm{TV}(\Omega)}$, where

[TABLE]

Below we discuss only one dataset, since the case of multiple datasets is similar. Then Tikhonov regularization leads to the following minimization problem:

[TABLE]

The scalar $\alpha>0$ is known as a regularization parameter. It has at least one minimizer [46, 22].

Since $\sigma$ is piecewise constant, by Lebesgue decomposition theorem [4], the TV term $|\sigma|_{\rm TV(\Omega)}$ in (2.3) reduces to $\int_{S_{\sigma}}|[\sigma]|\,{\rm d}\mathcal{H}^{d-1}$, where $S_{\sigma}$ is the jump set, $[\sigma]=\sigma^{+}-\sigma^{-}$ denotes the jump across $S_{\sigma}$ and $\mathcal{H}^{d-1}$ refers to the $(d-1)$-dimensional Hausdorff measure. The numerical approximation of (2.3) requires simultaneously treating two sets of different Hausdorff dimensions (i.e., $\Omega$ and $S_{\sigma}$), which is very challenging. Thus, we replace the TV term $|\sigma|_{\rm TV(\Omega)}$ in (2.3) by a Modica–Mortola type functional [43]

[TABLE]

where $\varepsilon$ is a small positive constant controlling the width of the transition interface, and $W:\mathbb{R}\rightarrow\mathbb{R}$ is the double-well potential given in (1.2). The functional $\mathcal{F}_{\varepsilon}$ was first proposed to model phase transition of two immiscible fluids in [11]. It is connected with the TV semi-norm as follows [43, 42, 2]; see [10] for an introduction to $\Gamma$-convergence.

Theorem 2.1**.**

With $c_{W}=\int_{c_{0}}^{c_{1}}\sqrt{W(s)}\,{\rm d}s$, let

[TABLE]

Then $\mathcal{F}_{\varepsilon}$ $\Gamma$-converges to $\mathcal{F}$ in $L^{1}(\Omega)$ as $\varepsilon\rightarrow 0^{+}$. Let $\{\varepsilon_{n}\}_{n\geq 1}$ and $\{v_{n}\}_{n\geq 1}$ be given sequences such that $\varepsilon_{n}\rightarrow 0^{+}$ and $\{\mathcal{F}_{\varepsilon_{n}}(v_{n})\}_{n\geq 1}$ is bounded. Then $v_{n}$ is precompact in $L^{1}(\Omega)$.

The proposed EIT reconstruction method reads

[TABLE]

where $\widetilde{\alpha}=\alpha/c_{W}$, and the admissible set $\widetilde{\mathcal{A}}$ is defined as

[TABLE]

Now we recall a useful continuity result of the forward map [22, Lemma 2.2], which gives the continuity of the fidelity term in the functional $\mathcal{J}_{\varepsilon}$. See also [31, 18] for related continuity results.

Lemma 2.2**.**

Let $\{\sigma_{n}\}_{n\geq 1}\subset\widetilde{\mathcal{A}}$ satisfy $\sigma_{n}\rightarrow\sigma^{\ast}$ in $L^{1}(\Omega)$. Then

[TABLE]

Lemma 2.2 implies that $\{\mathcal{J}_{\varepsilon}\}_{\varepsilon>0}$ are continuous perturbations of $\mathcal{J}$ in $L^{1}(\Omega)$. Then the stability of $\Gamma$-convergence [2, Proposition 1(ii)] [10, Remark 1.4] and Theorem 2.1 indicate that $\mathcal{J}_{\varepsilon}$ $\Gamma$-converges to $\mathcal{J}$ with respect to $L^{1}(\Omega)$, and $\mathcal{J}_{\varepsilon}$ can (approximately) recover piecewise constant conductivities. Next we show the existence of a minimizer to $\mathcal{J}_{\varepsilon}$.

Theorem 2.2**.**

For each $\varepsilon>0$, there exists at least one minimizer to problem (2.4).

Proof.

Since $\mathcal{J}_{\varepsilon}$ is nonnegative, there exists a minimizing sequence $\{\sigma_{n}\}_{n\geq 1}\subset\widetilde{\mathcal{A}}$ such that $\mathcal{J}_{\varepsilon}(\sigma_{n})\rightarrow m_{\varepsilon}:=\inf_{\sigma\in\widetilde{\mathcal{A}}}\mathcal{J}_{\varepsilon}(\sigma)$. Thus, $\sup_{n}\|\nabla\sigma_{n}\|_{L^{2}(\Omega)}<\infty$, which, along with $c_{0}\leq\sigma_{n}\leq c_{1}$, yields $\|\sigma_{n}\|_{H^{1}(\Omega)}\leq c$. Since $\widetilde{\mathcal{A}}$ is closed and convex, there exist a subsequence, relabeled as $\{\sigma_{n}\}_{n\geq 1}$, and some $\sigma^{\ast}\in\widetilde{\mathcal{A}}$ such that

[TABLE]

Since $W(s)\in C^{2}[c_{0},c_{1}]$, $\{W(\sigma_{n})\}_{n\geq 1}$ is uniformly bounded in $L^{\infty}(\Omega)$ and converges to $W(\sigma^{\ast})$ almost everywhere in $\Omega$. By Lebesgue’s dominated convergence theorem [19, p. 28, Theorem 1.19], $\int_{\Omega}W(\sigma_{n})\,{\rm d}x\rightarrow\int_{\Omega}W(\sigma^{\ast})\,{\rm d}x.$ By Lemma 2.2 and the weak lower semi-continuity of the $H^{1}(\Omega)$-seminorm, we obtain

[TABLE]

Thus $\sigma^{\ast}$ is a global minimizer of the functional $\mathcal{J}_{\varepsilon}$. ∎

To obtain the necessary optimality system of (2.4), we use the standard adjoint technique. The adjoint problem for (2.2) reads: find $(p,P)\in\mathbb{H}$ such that

[TABLE]

By straightforward computation, the Gâteaux derivative $\mathcal{J}_{\varepsilon}^{\prime}(\sigma)[\mu]$ of the functional $\mathcal{J}_{\varepsilon}$ at $\sigma\in\widetilde{\mathcal{A}}$ in the direction $\mu\in H^{1}(\Omega)$ is given by

[TABLE]

Then the minimizer $\sigma^{*}$ to problem (2.4) and the respective state $(u^{*},U^{*})$ and adjoint $(p^{*},P^{*})$ satisfy the following necessary optimality system:

[TABLE]

where the variational inequality at the last line is due to the box constraint in the admissible set $\widetilde{\mathcal{A}}$. The optimality system (2.8) forms the basis of the adaptive algorithm and its convergence analysis.

3 Adaptive algorithm

Now we develop an adaptive FEM for problem (2.4). Let $\mathcal{T}_{0}$ be a shape regular triangulation of $\overline{\Omega}$ into simplicial elements, each intersecting at most one electrode surface $e_{l}$, and $\mathbb{T}$ be the set of all possible conforming triangulations of $\overline{\Omega}$ obtained from $\mathcal{T}_{0}$ by the successive use of bisection. Then $\mathbb{T}$ is uniformly shape regular, i.e., the shape-regularity of any mesh $\mathcal{T}\in\mathbb{T}$ is bounded by a constant depending only on $\mathcal{T}_{0}$ [44, 52]. Over any $\mathcal{T}\in\mathbb{T}$, we define a continuous piecewise linear finite element space

[TABLE]

where $P_{1}(T)$ consists of all linear functions on $T$. The space $V_{\mathcal{T}}$ is used for approximating the state $u$ and adjoint $p$, and the discrete admissible set $\widetilde{\mathcal{A}}_{\mathcal{T}}$ for the conductivity is given by

[TABLE]

Given $\sigma_{\mathcal{T}}\in\widetilde{\mathcal{A}}_{\mathcal{T}}$, the discrete analogue of problem (2.2) is to find $(u_{\mathcal{T}},U_{\mathcal{T}})\in\mathbb{H}_{\mathcal{T}}\equiv V_{\mathcal{T}}\otimes\mathbb{R}^{L}_{\diamond}$ such that

[TABLE]

Then we approximate problem (2.4) by minimizing the following functional over $\widetilde{\mathcal{A}}_{\mathcal{T}}$:

[TABLE]

Then, similar to Theorem 2.2, there exists at least one minimizer $\sigma_{\mathcal{T}}^{*}$ to (3.2), and the minimizer $\sigma_{\mathcal{T}}^{\ast}$ and the related state $(u^{*}_{\mathcal{T}},U^{*}_{\mathcal{T}})\in\mathbb{H}_{\mathcal{T}}$ and adjoint $(p^{*}_{\mathcal{T}},P^{*}_{\mathcal{T}})\in\mathbb{H}_{\mathcal{T}}$ satisfy

[TABLE]

Further, $(u_{\mathcal{T}}^{\ast},U_{\mathcal{T}}^{\ast})$ and $(p_{\mathcal{T}}^{\ast},P_{\mathcal{T}}^{\ast})$ depend continuously on the problem data, i.e.,

[TABLE]

where the constant $c$ can be made independent of $\alpha$ and $\varepsilon$.

To describe the error estimators, we first recall some useful notation. The collection of all faces (respectively all interior faces) in $\mathcal{T}\in\mathbb{T}$ is denoted by $\mathcal{F}_{\mathcal{T}}$ (respectively $\mathcal{F}_{\mathcal{T}}^{i}$) and its restriction to the electrode $\bar{e}_{l}$ and $\partial\Omega\backslash\cup_{l=1}^{L}e_{l}$ by $\mathcal{F}_{\mathcal{T}}^{l}$ and $\mathcal{F}_{\mathcal{T}}^{c}$, respectively. A face/edge $F$ has a fixed normal unit vector $\boldsymbol{n}_{F}$ in $\overline{\Omega}$ with $\boldsymbol{n}_{F}=\boldsymbol{n}$ on $\partial\Omega$. The diameter of any $T\in\mathcal{T}$ and $F\in\mathcal{F}_{\mathcal{T}}$ is denoted by $h_{T}:=|T|^{1/d}$ and $h_{F}:=|F|^{1/(d-1)}$, respectively. For the solution $(\sigma^{\ast}_{\mathcal{T}},(u_{\mathcal{T}}^{\ast},U_{\mathcal{T}}^{\ast}),(p_{\mathcal{T}}^{\ast},P_{\mathcal{T}}^{\ast}))$ to problem (3.3), we define two element residuals for each element $T\in\mathcal{T}$ and two face residuals for each face $F\in\mathcal{F}_{\mathcal{T}}$ by

[TABLE]

where $[\cdot]$ denotes the jump across interior face $F$. Then for any element $T\in\mathcal{T}$, we define the following three error estimators

[TABLE]

with $q=d/(d-1)$. The estimator $\eta_{\mathcal{T},1}(\sigma_{\mathcal{T}}^{\ast},u_{\mathcal{T}}^{\ast},U_{\mathcal{T}}^{\ast},T)$ is identical with the standard residual error indicator for the direct problem: find $(\tilde{u},\tilde{U})\in\mathbb{H}$ such that

[TABLE]

It differs from the direct problem in (2.8) by replacing the conductivity $\sigma^{*}$ with $\sigma^{*}_{\mathcal{T}}$ instead, and is a perturbation of the latter case. The perturbation is vanishingly small in the event of the conjectured (subsequential) convergence $\sigma^{*}_{\mathcal{T}}\rightarrow\sigma^{*}$. The estimator $\eta_{\mathcal{T},2}(\sigma_{\mathcal{T}}^{\ast},p_{\mathcal{T}}^{\ast},P_{\mathcal{T}}^{\ast},T)$ admits a similar interpretation. These two estimators are essentially identical with that for the $H^{1}(\Omega)$ penalty in [32], and we refer to [32, Section 3.3] for a detailed heuristic derivation. The estimator $\eta_{\mathcal{T},3}(\sigma_{\mathcal{T}}^{\ast},u_{\mathcal{T}}^{\ast},p_{\mathcal{T}}^{\ast},T)$ is related to the variational inequality in the necessary optimality condition (2.8), and roughly provides a quantitative measure how well it is satisfied. The estimator (including the exponent $q$) is motivated by the convergence analysis; see the proof of Theorem 5.5 and Remark 5.2 below. It represents the main new ingredient for problem (2.4), and differs from that for the $H^{1}(\Omega)$ penalty in [32].

Remark 3.1**.**

The estimator $\eta_{k,3}$ improves that in [32], i.e.,

[TABLE]

in terms of the exponents on $h_{T}$ and $h_{F}$. This improvement is achieved by a novel constraint preserving interpolation operator defined in (5.13) below.

Now we can formulate an adaptive algorithm for (2.4); see Algorithm 3.1. Below we indicate the dependence on the mesh $\mathcal{T}_{k}$ by the subscript $k$, e.g., $\mathcal{J}_{\varepsilon,k}$ for $\mathcal{J}_{\varepsilon,\mathcal{T}_{k}}$.

The MARK module selects a collection of elements in the mesh $\mathcal{T}_{k}$. The condition (3.5) covers several commonly used marking strategies, e.g., maximum, equidistribution, modified equidistribution, and Dörfler’s strategy [49, pp. 962]. Compared with a collective marking in AFEM in [32], Algorithm 3.1 employs a separate marking to select more elements for refinement in each loop, which leads to fewer iterations of the adaptive process. The error estimators may also be used for coarsening, which is relevant if the recovered inclusions change dramatically during the iteration. However, the convergence analysis below does not carry over to coarsening, and it will not be further explored.

Last, we give the main theoretical result: for each fixed $\varepsilon>0$, the sequence of discrete solutions $\{\sigma_{k}^{\ast},(u_{k}^{\ast},U_{k}^{\ast}),(p_{k}^{\ast},P_{k}^{\ast})\}_{k\geq 0}$ generated by Algorithm 3.1 contains a subsequence converging in $H^{1}(\Omega)\times\mathbb{H}\times\mathbb{H}$ to a solution of system (2.8). The proof is lengthy and technical, and thus deferred to Section 5.

Theorem 3.1**.**

The sequence of discrete solutions $\{\sigma_{k}^{\ast},(u_{k}^{\ast},U_{k}^{\ast}),(p_{k}^{\ast},P_{k}^{\ast})\}_{k\geq 0}$ by Algorithm 3.1 contains a subsequence $\{\sigma_{k_{j}}^{\ast},(u_{k_{j}}^{\ast},U_{k_{j}}^{\ast}),(p_{k_{j}}^{\ast},P_{k_{j}}^{\ast})\}_{j\geq 0}$ convergent to a solution $(\sigma^{\ast},(u^{\ast},U^{\ast}),(p^{\ast},P^{\ast}))$ of system (2.8):

[TABLE]

4 Numerical experiments and discussions

Now we present numerical results to illustrate Algorithm 3.1 on a square domain $\Omega=(-1,1)^{2}$. There are sixteen electrodes $\{e_{l}\}_{l=1}^{L}$ (with $L=16$) evenly distributed along $\partial\Omega$, each of length $1/4$. The contact impedances $\{z_{l}\}_{l=1}^{L}$ are all set to unit. We take ten sinusoidal input currents, and for each voltage $U(\sigma^{\dagger})\in\mathbb{R}^{L}_{\diamond}$, generate the noisy data $U^{\delta}$ by

[TABLE]

where $\epsilon$ is the (relative) noise level, and $\{\xi_{l}\}_{l=1}^{L}$ follow the standard normal distribution. Note that $\epsilon=\text{1e-2}$ refers to a relatively high noise level for EIT. The exact data $U(\sigma^{\dagger})$ is computed using a much finer uniform mesh, to avoid the most obvious form of “inverse crime”.

In the experiments, we fix $K$ (the number of refinements) at $15$, $q$ (exponent in $\eta_{k,3}^{q}$) at $2$, and $\varepsilon$ (the functional $\mathcal{F}_{\varepsilon}$) at 1e-2. The marking strategy (3.5) in the module MARK selects a minimal refinement set $\mathcal{M}_{k}{:=\cup_{i=1}^{3}\mathcal{M}_{k}^{i}}\subseteq\mathcal{T}_{k}$ such that

[TABLE]

with a threshold $\theta=0.7$. The refinement is performed with one popular refinement strategy, i.e., newest vertex bisection [41]. Specifically, it connects the midpoint $x_{T}$, as a newest vertex, of a reference edge $F$ of an element $T\in\mathcal{T}_{k}$ to the opposite node of $F$, and employs two edges opposite to the midpoint $x_{T}$ as reference edges of the two newly created triangles in $\mathcal{T}_{k+1}$. Problem (3.1)-(3.2) is solved by a Newton type method; see Appendix A for the detail. The conductivity on $\mathcal{T}_{0}$ is initialized to $\sigma_{0}=c_{0}$, and then for $k=1,2,\ldots$, $\sigma_{k-1}^{*}$ (defined on $\mathcal{T}_{k-1}$) is interpolated to $\mathcal{T}_{k}$ to warm start the optimization. The regularization parameter $\widetilde{\alpha}$ in (2.4) is determined in a trial-and-error manner. All computations are performed using MATLAB 2018a on a personal laptop with 8.00 GB RAM and 2.5 GHz CPU.

The first set of examples are concerned with two inclusions.

Example 1**.**

The background conductivity $\sigma_{0}(x)=1$.

• $\rm(i)$

The true conductivity $\sigma^{\dagger}$ is given by $\sigma_{0}(x)+\chi_{B_{1}}(x)+\chi_{B_{2}}(x)$, with $B_{1}$ and $B_{2}$ denote two circles centered at $(0,0.5)$ and $(0,-0.5)$, respectively, both with a radius 0.3.

• $\rm(ii)$

The true conductivity $\sigma^{\dagger}$ is given by $\sigma_{0}(x)+1+1.2e^{-\frac{25(x_{1}^{2}+(x_{2}-0.5)^{2})}{2}}+1.2e^{-\frac{25(x_{1}^{2}+(x_{2}+0.5)^{2})}{2}}$, i.e., two Gaussian bumps centered at $(0,0.5)$ and $(0,-0.5)$.

• $\rm(iii)$

The true conductivity $\sigma^{\dagger}$ is given by $\sigma_{0}(x)+5\chi_{B_{1}}(x)+5\chi_{B_{2}}(x)$, with $B_{1}$ and $B_{2}$ denote two circles centered at $(0,0.5)$ and $(0,-0.5)$, respectively, both with a radius 0.3.

The numerical results for Example 1(i) with $\epsilon=\text{1e-3}$ and $\epsilon=\text{1e-2}$ are shown in Figs. 1–5, where d.o.f. denotes the degree of freedom of the mesh. It is observed from Fig. 1 that with both uniform and adaptive refinements, the final recoveries have comparable accuracy and capture well the inclusion locations.

Next we examine the adaptive refinement process more closely. In Figs. 2 and 3, we show the meshes $\mathcal{T}_{k}$ during the iteration and the corresponding recoveries $\sigma_{k}$ for Example 1(i) at two noise levels $\epsilon=\text{1e-3}$ and $\epsilon=\text{1e-2}$, respectively. On the coarse mesh $\mathcal{T}_{0}$, the recovery has very large errors and can only identify one component and thus fails to correctly identify the number of inclusions, due to the severe under-resolution of both state and conductivity. Nonetheless, Algorithm 3.1 can correctly recover the two components with reasonable accuracy after several adaptive loops, and accordingly, the support of the recovery is gradually refined with its accuracy improving steadily. In particular, the inclusion locations stabilize after several loops, and thus coarsening of the mesh seems unnecessary. Throughout, the refinement occurs mainly in the regions around the electrode edges and internal interface, which is clearly observed for both noise levels. This is attributed to the separable marking strategy, which allows detecting different sources of singularities simultaneously. In Fig. 4, we display the evolution of the error indicators for Example 1(i) with $\epsilon=\text{1e-3}$. The estimators play different roles: $\eta_{k,1}^{2}$ and $\eta_{k,2}^{2}$ indicate the electrode edges during first iterations and then also internal interface, whereas throughout $\eta_{k,3}^{2}$ concentrates on the internal interface. Thus, $\eta_{k,1}^{2}$ and $\eta_{k,2}^{2}$ are most effective for resolving the state and adjoint, whereas $\eta_{k,3}^{2}$ is effective for detecting internal jumps of the conductivity. The magnitude of $\eta_{k,2}^{2}$ is much smaller than $\eta_{k,1}^{2}$, since the boundary data $U^{\delta}-U(\sigma_{k})$ for the adjoint is much smaller than the input current $I$ for the state. Thus, a simple collective marking strategy (i.e., $\eta_{k}^{2}=\eta_{k,1}^{2}+\eta_{k,2}^{2}+\eta_{k,3}^{2}$) may miss the correct singularity, due to their drastically different scalings. In contrast, the separate marking in (3.5) can take care of the scaling automatically.

In Fig. 5, we plot the $L^{2}(\Omega)$ and $L^{1}(\Omega)$ errors of the recoveries versus d.o.f. $N$, where the recovery on the corresponding finest mesh is taken as the reference (since the recoveries by the adaptive and uniform refinements are slightly different; see Fig. 1). Due to the discontinuity of the sought-for conductivity, the $L^{1}(\Omega)$ norm is especially suitable for measuring the convergence. The convergence of the algorithm is clearly observed for both adaptive and uniform refinements. Further, with a fixed d.o.f., AFEM gives more accurate results than the uniform one in both error metrics. These observations show the computational efficiency of the adaptive algorithm.

Examples 1(ii) and (iii) are variations of Example 1(i), and the results are presented in Figs. 6–9. The proposed approach assumes a piecewise constant conductivity with known lower and upper bounds. Example 1(ii) does not fulfill the assumption, since the true conductivity $\sigma^{\dagger}$ is not piecewise constant. Thus the algorithm can only produce a piecewise constant approximation to the exact one. Nonetheless, the inclusion support is reasonably identified. When the noise level $\epsilon$ increases from 1e-3 to 1e-2, the reconstruction accuracy deteriorates significantly; see Fig. 6. Example 1(iii) involves high contrast inclusions, which are well known to be numerically more challenging. This is clearly observed in Fig. 8, where the recovery accuracy is inferior, especially for the noise level $\epsilon=\text{1e-2}$. However, the adaptive refinement procedure works well similarly as the preceding examples: the refinement occurs mainly around the electrode edges and inclusion interface; see Figs. 7 and 9 for the details.

Now we consider one more challenging example with four inclusions.

Example 2**.**

The true conductivity $\sigma^{\dagger}$ is given by $\sigma_{0}(x)+\sum_{i=1}^{4}\chi_{B_{i}}(x)$, with the circles $B_{i}$ centered at $(0.6,\pm 0.6)$ and $(-0.6,\pm 0.6)$, respectively, all with a radius 0.2, and the background conductivity $\sigma_{0}(x)=1$.

The numerical results for Example 2 are given in Figs. 10–12. The results are in excellent agreement with the observations from Example 1: The algorithm converges steadily as the adaptive iteration proceeds, and with a low noise level, it can accurately recover all four inclusions, showing clearly the efficiency of the adaptive approach. The refinement is mainly around the electrode edge and interval interface.

5 Proof of Theorem 3.1

The lengthy and technical proof is divided into two steps: Step 1 shows the convergence to an auxiliary minimization problem over a limiting admissible set in Section 5.1, and Step 2 shows that the solution of the auxiliary problem satisfies the necessary optimality system (2.8) in section 5.2. The overall proof strategy is similar to [32], and hence we omit relevant arguments.

5.1 Auxiliary convergence

Since the two sequences $\{\mathbb{H}_{k}\}_{k\geq 0}$ and $\{\widetilde{\mathcal{A}}_{k}\}_{k\geq 0}$ generated by Algorithm 3.1 are nested, we may define

[TABLE]

Clearly $\mathbb{H}_{\infty}$ is a closed subspace of $\mathbb{H}$. For the set $\widetilde{\mathcal{A}}_{\infty}$, we have the following result [32, Lemma 4.1].

Lemma 5.1**.**

$\widetilde{\mathcal{A}}_{\infty}$* is a closed convex subset of $\widetilde{\mathcal{A}}$.*

Over the limiting set $\widetilde{\mathcal{A}}_{\infty}$, we define an auxiliary limiting minimization problem:

[TABLE]

where $(u_{\infty},U_{\infty})\in\mathbb{H}_{\infty}$ satisfies

[TABLE]

By Lemma 2.1 and Lax-Milgram theorem, problem (5.2) is well-posed for any fixed $\sigma_{\infty}\in\mathcal{A}_{\infty}$. The next result gives the existence of a minimizer to (5.1)–(5.2).

Theorem 5.1**.**

There exists at least one minimizer to problem (5.1)–(5.2).

Proof.

Let $\{\sigma_{k}^{\ast},(u_{k}^{\ast},U_{k}^{\ast})\}_{k\geq 0}$ be the sequence of discrete solutions given by Algorithm 3.1. Since $c_{1}\in\widetilde{\mathcal{A}}_{k}$ for all $k$, by (3.4), $\mathcal{J}_{\varepsilon,k}(\sigma_{k}^{\ast})\leq\mathcal{J}_{\varepsilon,k}(c_{1})\leq c$, and thus $\{\sigma_{k}^{\ast}\}_{k\geq 0}$ is uniformly bounded in $H^{1}(\Omega)$. By Lemma 5.1 and Sobolev embedding, there exist a subsequence, denoted by $\{\sigma_{k_{j}}^{\ast}\}_{j\geq 0}$, and some $\sigma^{\ast}\in\widetilde{\mathcal{A}}_{\infty}$ such that

[TABLE]

Next we introduce a discrete analogue of problem (5.2) with $\sigma_{\infty}=\sigma^{\ast}$: find $(u_{k_{j}},U_{k_{j}})\in\mathbb{H}_{k_{j}}$ such that

[TABLE]

By Lemma 2.1, Cea’s lemma and the construction of the space $\mathbb{H}_{\infty}$, the solution $(u^{\ast}_{\infty},U_{\infty}^{\ast})\in\mathbb{H}_{\infty}$ of (5.2) with $\sigma_{\infty}=\sigma^{\ast}$ satisfies

[TABLE]

Taking the test function $(v,V)=(u_{k_{j}}-u_{k_{j}}^{\ast},U_{k_{j}}-U_{k_{j}}^{\ast})\in\mathbb{H}_{k_{j}}$ in the first line of (3.3) and (5.4) and then applying the Cauchy-Schwarz inequality lead to

[TABLE]

In view of (5.5), pointwise convergence in (5.3) and Lebesgue’s dominated convergence theorem,

[TABLE]

This and Lemma 2.1 imply $\|(u_{k_{j}}-u_{k_{j}}^{\ast},U_{k_{j}}-U_{k_{j}}^{\ast})\|_{\mathbb{H}}\rightarrow 0.$ Then, (5.5) and the triangle inequality imply

[TABLE]

Meanwhile, repeating the argument of Theorem 2.2 gives

[TABLE]

next we apply a density argument. For any $\sigma_{\infty}\in\widetilde{\mathcal{A}}_{\infty}$, by the construction of the space $\mathbb{H}_{\infty}$, there exists a sequence $\{\sigma_{k}\}_{k\geq 0}\subset\bigcup_{k\geq 0}\widetilde{\mathcal{A}}_{k}$ such that $\sigma_{k}\rightarrow\sigma_{\infty}$ in $H^{1}(\Omega)$. Repeating the preceding argument gives $\|U(\sigma_{k})-U^{\delta}\|^{2}\rightarrow\|U(\sigma_{\infty})-U^{\delta}\|^{2}$ and $\int_{\Omega}W(\sigma_{k})\,{\rm d}x\rightarrow\int_{\Omega}W(\sigma_{\infty})\,{\rm d}x$. Now (5.6), the weak lower semicontinuity of the $H^{1}(\Omega)$-norm, (5.7) and the minimizing property of $\sigma_{k}^{\ast}$ to $\mathcal{J}_{\varepsilon,k}$ over the set $\widetilde{\mathcal{A}}_{k}$ imply

[TABLE]

Since $\sigma^{\ast}\in\widetilde{\mathcal{A}}_{\infty}$, $\sigma_{\infty}^{\ast}:=\sigma^{\ast}$ is a minimizer of $\mathcal{J}_{\varepsilon,\infty}$ over $\widetilde{\mathcal{A}}_{\infty}$. ∎

Further, we have the following auxiliary convergence.

Theorem 5.2**.**

The sequence of discrete solutions $\{\sigma_{k}^{\ast},(u^{\ast}_{k},U_{k}^{\ast})\}_{k\geq 0}$ to problem (3.2) contains a subsequence $\{\sigma_{k_{j}}^{\ast},(u_{k_{j}}^{\ast},U_{k_{j}}^{\ast})\}_{j\geq 0}$ convergent to a minimizer $(\sigma_{\infty}^{\ast},(u_{\infty}^{\ast},U_{\infty}^{\ast}))$ to problem (5.1)–(5.2):

[TABLE]

Proof.

The convergence of $(u_{k_{j}}^{*},U_{k_{j}}^{*})$ was already proved in Theorem 5.1. Taking $\sigma_{\infty}=\sigma_{\infty}^{\ast}$ in (5.1) gives $\lim_{j\rightarrow\infty}\mathcal{J}_{\varepsilon,k_{j}}(\sigma_{k_{j}}^{\ast})=\mathcal{J}_{\varepsilon,\infty}(\sigma^{\ast}_{\infty})$. By (5.6) and (5.7), we have $\|\nabla\sigma_{k_{j}}^{\ast}\|^{2}_{L^{2}(\Omega)}\rightarrow\|\nabla\sigma_{\infty}^{\ast}\|^{2}_{L^{2}(\Omega)}$. Thus, the sequence $\{\sigma_{k_{j}}^{\ast}\}_{j\geq 0}$ converges to $\sigma_{\infty}^{\ast}$ in $H^{1}(\Omega)$. ∎

Next we consider the convergence of the sequence $\{(p_{k}^{\ast},P_{k}^{\ast})\}_{k\geq 0}$. With a minimizer $(\sigma_{\infty}^{\ast},(u^{\ast}_{\infty},U^{\ast}_{\infty}))$ to problem (5.1), we define a limiting adjoint problem: find $(p^{\ast}_{\infty},P^{\ast}_{\infty})\in\mathbb{H}_{\infty}$ such that

[TABLE]

By Lemma 2.1 and Lax-Milgram theorem, (5.9) is uniquely solvable. We have the following convergence result for $(p_{\infty}^{\ast},P_{\infty}^{\ast})$. The proof is identical with [32, Theorem 4.5], and hence omitted.

Theorem 5.3**.**

Under the condition of Theorem 5.2, the subsequence of adjoint solutions $\{(p_{k_{j}}^{\ast},P_{k_{j}}^{\ast})\}_{j\geq 0}$ generated by Algorithm 3.1 converges to the solution $(p_{\infty}^{\ast},P_{\infty}^{\ast})$ of problem (5.9):

[TABLE]

5.2 Proof of Theorem 3.1

Theorem 3.1 follows directly by combining Theorems 5.2-5.3 in Section 5.1 and Theorems 5.4-5.5 below. The proof in this part relies on the marking condition (3.5). First, we show that the limit $(\sigma^{\ast}_{\infty},(u^{\ast}_{\infty},U^{\ast}_{\infty}),(p^{\ast}_{\infty},P^{\ast}_{\infty}))$ solves the variational equations in (2.8).

Theorem 5.4**.**

The solutions $(\sigma_{\infty}^{*},u_{\infty}^{*},U_{\infty}^{*})$ and $(p_{\infty}^{*},P_{\infty}^{*})$ to problems (5.1)-(5.2) and (5.9) satisfy

[TABLE]

Proof.

The proof is identical with [32, Lemma 4.8], using Theorems 5.2-5.3, and hence we only give a brief sketch. By [32, Lemma 3.5], for each $T\in\mathcal{T}_{k}$ with its face $F$ (intersecting with $e_{l}$), there hold

[TABLE]

where the notation $D_{T}$ is defined below. Then by the marking condition (3.5), [32, Lemma 4.6] implies that for each convergent subsequence $\{\sigma^{\ast}_{k_{j}},(u^{\ast}_{k_{j}},U^{\ast}_{k_{j}}),(p^{\ast}_{k_{j}},P^{\ast}_{k_{j}})\}_{j\geq 0}$ from Theorems 5.2 and 5.3, there hold

[TABLE]

Last, by [32, Lemma 4.7] and Theorems 5.2-5.3, the argument of [32, Lemma 4.8] completes the proof. ∎

Remark 5.1**.**

The argument of Theorem 5.4 dates back to [49], and the main tools include the Galerkin orthogonality of the residual operator, the Lagrange and the Scott-Zhang interpolation operators [16, 48], the marking condition (3.5) and a density argument. Further, the error estimators $\eta_{k,1}(\sigma_{k}^{*},u_{k}^{*},U_{k}^{*})$ and $\eta_{k,2}(\sigma_{k}^{*},p_{k}^{*},P_{k}^{*})$ emerge in the proof and are then employed in the module ESTIMATE of Algorithm 3.1.

Next we prove that the limit $(\sigma^{\ast}_{\infty},(u^{\ast}_{\infty},U^{\ast}_{\infty}),(p^{\ast}_{\infty},P^{\ast}_{\infty}))$ satisfies the variational inequality in (2.8). The proof relies crucially on a constraint preserving interpolation operator. We denote by $D_{T}$ the union of elements in $\mathcal{T}$ with a non-empty intersection with an element $T\in\mathcal{T}$, and by $\omega_{F}$ the union of elements in $\mathcal{T}$ sharing a common face/edge with $F\in\mathcal{F}_{\mathcal{T}}$. Let

[TABLE]

The set $\mathcal{T}_{k}^{+}$ consists of all elements not refined after the $k$-th iteration, and all elements in $\mathcal{T}_{k}^{0}$ are refined at least once after the $k$-th iteration. Clearly, $\mathcal{T}_{l}^{+}\subset\mathcal{T}_{k}^{+}$ for $l<k$. We also define a mesh-size function $h_{k}:\overline{\Omega}\rightarrow\mathbb{R}^{+}$ almost everywhere

[TABLE]

where ${T}^{i}$ denotes the interior of an element $T\in\mathcal{T}_{k}$, and ${F}^{i}$ the relative interior of an edge $F\in\mathcal{F}_{k}$. It has the following property [49, Corollary 3.3]:

[TABLE]

The next result gives the limiting behaviour of the maximal error indicator $\eta_{k,3}$.

Lemma 5.2**.**

Let $\{(\sigma^{\ast}_{k},(u^{\ast}_{k},U^{\ast}_{k}),(p^{\ast}_{k},P^{\ast}_{k}))\}_{k\geq 0}$ be the sequence of discrete solutions generated by Algorithm 3.1. Then for each convergent subsequence $\{\sigma^{\ast}_{k_{j}},(u^{\ast}_{k_{j}},U^{\ast}_{k_{j}}),(p^{\ast}_{k_{j}},P^{\ast}_{k_{j}})\}_{j\geq 0}$, there holds

[TABLE]

Proof.

The inverse estimate and scaled trace theorem imply that for each $T\in\mathcal{T}_{k}$ (with its face $F$)

[TABLE]

With the choice $q=d/(d-1)$, combining these two estimates gives

[TABLE]

where $c$ depends on $\widetilde{\alpha}$ and $\varepsilon$ in $\mathcal{F}_{\varepsilon}$. Next, for the subsequence $\{\sigma^{\ast}_{k_{j}},(u^{\ast}_{k_{j}},U^{\ast}_{k_{j}}),(p^{\ast}_{k_{j}},P^{\ast}_{k_{j}})\}_{j\geq 0}$, let $\widetilde{T}_{j}^{3}\in\mathcal{M}_{k_{j}}^{3}$ be the element with the largest error indicator $\eta_{k_{j},3}(\sigma_{k_{j}}^{*},u_{k_{j}}^{*},p_{k_{j}}^{*},T)$. Since $D_{\widetilde{T}_{j}^{i}}\subset\Omega_{k_{j}}^{0}$, (5.10) implies

[TABLE]

By (5.11), Cauchy-Schwarz inequality and triangle inequality, there holds

[TABLE]

By Theorems 5.2 and 5.3, Lebesgue’s dominated convergence theorem, the choice $q=d/(d-1)\leq 2$ and Hölder inequality, we obtain $\|W^{\prime}(\sigma_{k_{j}}^{\ast})-W^{\prime}(\sigma_{\infty}^{\ast})\|^{q}_{L^{q}(\Omega)}\rightarrow 0$ and $\|\nabla(\sigma_{k_{j}}^{\ast}-\sigma_{\infty}^{\ast})\|^{q}_{L^{q}(\Omega)}\rightarrow 0$. Then the absolute continuity of the norm $\|\cdot\|_{L^{q}(\Omega)}$ with respect to Lebesgue measure and (5.12) complete the proof. ∎

Due to a lack of Galerkin orthogonality for variational inequalities, we employ a local $L^{r}$-stable interpolation operator of Clément/Chen-Nochetto type. Let $\mathcal{N}_{k}$ be the set of all interior nodes of $\mathcal{T}_{k}$ and $\{\phi_{x}\}_{x\in\mathcal{N}_{k}}$ be the nodal basis functions in $V_{k}$. For each $x\in\mathcal{N}_{k}$, the support of $\phi_{x}$ is denoted by $\omega_{x}$, i.e., the union of all elements in $\mathcal{T}_{k}$ with a non-empty intersection with $x$. Then we define $\Pi_{k}:L^{1}(\Omega)\rightarrow V_{k}$ by

[TABLE]

Clearly, $\Pi_{k}v\in\widetilde{\mathcal{A}}_{k}$ if $c_{0}\leq v\leq c_{1}$ a.e. $x\in\Omega$. The definition is adapted from [13] (for elliptic obstacle problems) by replacing the maximal ball $\Delta_{x}\subset\omega_{x}$ centered at an interior node $x$ by $\omega_{x}$. $\Pi_{k}v$ satisfies following properties; see Appendix B for a proof.

Lemma 5.3**.**

For any $v\in W^{1,r}(\Omega)$, there hold for all $r\in[1,+\infty]$, any $T\in\mathcal{T}_{k}$ and any $F\subset\partial T$,

[TABLE]

Last we show that the limit $(\sigma^{\ast}_{\infty},(u^{\ast}_{\infty},U^{\ast}_{\infty}),(p^{\ast}_{\infty},P^{\ast}_{\infty}))$ satisfies the variational inequality in (2.8).

Theorem 5.5**.**

The solutions $(\sigma_{\infty}^{*},u_{\infty}^{*},U_{\infty}^{*})$ and $(p_{\infty}^{*},P_{\infty}^{*})$ to problems (5.1)-(5.2) and (5.9) satisfy

[TABLE]

Proof.

The proof is lengthy, and we break it into five steps.

Step i. Derive a preliminary variational inequality. We relabel the subsequence $\{\sigma_{k_{j}}^{\ast},(u_{k_{j}}^{\ast},U_{k_{j}}^{\ast}),(p_{k_{j}}^{\ast},P_{k_{j}}^{\ast})\}_{j\geq 0}$ in Theorems 5.2 and 5.3 as $\{\sigma_{k}^{\ast},(u_{k}^{\ast},U_{k}^{\ast}),(p_{k}^{\ast},P_{k}^{\ast})\}_{k\geq 0}$. Let $I_{k}$ be the Lagrange interpolation operator on $V_{k}$, and let $\alpha^{\prime}=\widetilde{\alpha}\varepsilon$ and $\alpha^{\prime\prime}=\frac{\widetilde{\alpha}}{2\varepsilon}$. For any $\mu\in\widetilde{\mathcal{A}}\cap C^{\infty}(\overline{\Omega})$, $I_{k}\mu\in\widetilde{\mathcal{A}}_{k}$ and let $\nu=\mu-I_{k}\mu$. Direct computation gives

[TABLE]

where the last inequality is due to the variational inequality in (3.3) with $\mu_{k}=I_{k}\mu$.

Step ii. Bound the ${\rm I}$. By elementwise integration by parts, Hölder inequality, the definition of the estimator $\eta_{k,3}$ and Lemma 5.3 with $r=q^{\prime}$ (with $q^{\prime}$ being the conjugate exponent of $q$),

[TABLE]

Thus, for any $k>l$, by (discrete) Hölder’s inequality and the finite overlapping property of the patches $D_{T}$, due to uniform shape regularity of the meshes $\mathcal{T}_{k}\in\mathbb{T}$, there holds

[TABLE]

Since $W^{\prime}(s)\in C^{1}[c_{0},c_{1}]$, by the pointwise convergence of $\{\sigma_{k}^{\ast}\}_{k\geq 0}$ in Theorem 5.2 and Lebesgue’s dominated convergence theorem, we deduce

[TABLE]

Since $q=d/{(d-1)}\leq 2$, the sequence $\{W^{\prime}(\sigma_{k}^{\ast})\}_{k\geq 0}$ is uniformly bounded in $L^{q}(\Omega)$. By Theorems 5.2 and 5.3, the sequences $\{\sigma_{k}^{\ast}\}_{k\geq 0}$, $\{u_{k}^{\ast}\}_{k\geq 0}$ and $\{p_{k}^{\ast}\}_{k\geq 0}$ are uniformly bounded in $H^{1}(\Omega)$. Thus, (5.11) and (5.10), and Hölder inequality give

[TABLE]

Then by the error estimate of $I_{k}$ [16],

[TABLE]

By (5.10), $c\|h_{l}\chi_{\Omega^{0}_{l}}\|_{L^{\infty}(\Omega)}\|\mu\|_{W^{2,q^{\prime}}(\Omega)}\rightarrow 0$ as $l\rightarrow\infty$. Since $\mathcal{T}_{l}^{+}\subset\mathcal{T}_{k}$ for $k>l$, (3.5) implies

[TABLE]

By Lemma 5.2, for any small $\varepsilon>0$, we can choose $k_{1}>l_{1}$ for some large fixed $l_{1}$ such that whenever $k>k_{1}$,

[TABLE]

Consequently,

[TABLE]

Step iii. Bound the term ${\rm II}$. For the term $\rm II$, elementwise integration and Hölder inequality yield

[TABLE]

By the scaled trace theorem, local inverse estimate, $L^{q^{\prime}}$-stability of $\Pi_{k}$ in Lemma 5.3, local quasi-uniformity and interpolation error estimate for $I_{k}$ [16], we deduce that for $k>l$

[TABLE]

Since \big{(}\sum_{T\in\mathcal{T}_{k}\setminus\mathcal{T}_{l}^{+}}\eta^{q}_{k,3}(\sigma_{k}^{\ast},u_{k}^{\ast},p_{k}^{\ast},T)\big{)}^{1/q}\leq c, cf. (5.16), there holds

[TABLE]

Now by repeating the argument for the term $\rm I$, we obtain

[TABLE]

Step iv. Take limit in preliminary variational inequality. Using (5.15) and the $H^{1}(\Omega)$-convergence of $\{\sigma_{k}^{\ast}\}_{k\geq 0}$ in Theorem 5.2, we have for each $\mu\in\widetilde{\mathcal{A}}\cap C^{\infty}(\overline{\Omega})$

[TABLE]

Further, the uniform boundedness on $\{u_{k}^{\ast}\}_{k\geq 0}$ in $H^{1}(\Omega)$ and the convergence of $\{p^{\ast}_{k}\}_{k\geq 0}$ to $p_{\infty}^{\ast}$ in $H^{1}(\Omega)$ in Theorem 5.3 yield

[TABLE]

This and Theorem 5.2 imply

[TABLE]

In the splitting

[TABLE]

the arguments for (5.20) directly yields

[TABLE]

The boundedness on $\{u_{k}^{\ast}\}_{k\geq 0}$ in $H^{1}(\Omega)$, pointwise convergence of $\{\sigma_{k}^{\ast}\}_{k\geq 0}$ of Theorem 5.2 and Lebesgue’s dominated convergence theorem imply

[TABLE]

Hence, there holds

[TABLE]

Now by passing both sides of (5.14) to the limit and combining the estimates (5.17)-(5.21), we obtain

[TABLE]

Step v. Density argument. By the density of $C^{\infty}(\overline{\Omega})$ in $H^{1}(\Omega)$ and the construction via a standard mollifier [19], for any $\mu\in\widetilde{\mathcal{A}}$ there exists a sequence $\{\mu_{n}\}\subset\widetilde{\mathcal{A}}\cap C^{\infty}(\overline{\Omega})$ such that $\|\mu_{n}-\mu\|_{H^{1}(\Omega)}\rightarrow 0$ as $n\rightarrow\infty$. Thus, $(\nabla\sigma^{\ast}_{\infty},\nabla\mu_{n})\rightarrow(\nabla\sigma^{\ast}_{\infty},\nabla\mu)$, $(W^{\prime}(\sigma_{\infty}^{\ast}),\mu_{n})\rightarrow(W^{\prime}(\sigma_{\infty}^{\ast}),\mu)$, and $(\mu_{n}\nabla u^{\ast}_{\infty},\nabla p_{\infty}^{\ast})\rightarrow(\mu\nabla u^{\ast}_{\infty},\nabla p_{\infty}^{\ast})$, after possibly passing to a subsequence. The desired result follows from the preceding two estimates. ∎

Remark 5.2**.**

The computable quantity $\eta_{k,3}(\sigma_{k}^{*},u_{k}^{*},p_{k}^{*},T)$ emerges naturally from the proof, i.e., the upper bounds on ${\rm I}$ and ${\rm II}$, which motivates its use as the a posteriori error estimator in Algorithm 3.1.

Acknowledgements

The authors are grateful to an anonymous referee and the boarder member for the constructive comments, which have significantly improved the presentation of the paper. The work of Y. Xu was partially supported by National Natural Science Foundation of China (11201307), Ministry of Education of China through Specialized Research Fund for the Doctoral Program of Higher Education (20123127120001) and Natural Science Foundation of Shanghai (17ZR1420800).

Appendix A The solution of the variational inequality

Now we describe an iterative method for minimizing the energy functional

[TABLE]

Let $p(z)=(z-c_{0})(z-c_{1})$. Then one linearized approximation $p_{L}(z,z_{k})$ reads (with $\delta z=z-z_{k}$)

[TABLE]

Upon substituting the approximation $p_{L}(z,z_{k})$ for $p(z)$ and linearizing the forward map $U(\sigma)$, we obtain the following surrogate energy functional (with $\delta\sigma=\sigma-\sigma_{k}$ being the increment and $\delta U=U^{\delta}-U(\sigma_{k})$)

[TABLE]

The treatment of the double well potential term $\int_{\Omega}W(\sigma){\rm d}x$ is in the spirit of the classical majorization-minimization algorithm in the following sense (see [56] for a detailed derivation)

[TABLE]

This algorithm is known to have excellent numerical stability. Upon ignoring the box constraint on the conductivity $\sigma$, problem (A.1) is to find $\delta\sigma\in H^{1}(\Omega)$ such that

[TABLE]

This equation can be solved by an iterative method for the update $\delta\sigma$ (with the box constraint treated by a projection step). Note that $U^{\prime}(\sigma_{k})$ and $U^{\prime}(\sigma_{k})^{*}$ can be implemented in matrix-free manner using the standard adjoint technique. In our experiment, we employ the conjugate gradient method to solve the resulting linear systems, preconditioned by the sparse matrix corresponding to $\tilde{\alpha}\varepsilon(\nabla\delta\sigma,\nabla\phi)+\frac{\tilde{\alpha}}{\varepsilon}(p^{\prime}(\sigma_{k})^{2}\delta\sigma,\phi)$.

Appendix B Proof of Lemma 5.3

The proof follows that in [13, 25]. By Hölder inequality and $h_{T}^{d}\leq|\omega_{x}|$ for each node $x\in T$,

[TABLE]

The desired $L^{r}$-stability follows from the estimate $\|\phi_{x}\|_{L^{r}(T)}\leq ch_{T}^{d/r}$, by the local quasi-uniformity of the mesh. In view of the definition (5.13), $\Pi_{k}\zeta=\zeta$ for any $\zeta\in\mathbb{R}$. By local inverse estimate, the $L^{r}$-stability of $\Pi_{k}$, standard interpolation error estimate [16] and local quasi-uniformity,

[TABLE]

Similarly,

[TABLE]

By the scaled trace theorem, for any $F\subset\partial T$, there holds

[TABLE]

Then (B.1) and (B.2) complete the proof of the lemma.