Inverse problems and invisibility cloaking for FEM models and resistor networks

In this paper we consider inverse problems for resistor networks and for models obtained via the Finite Element Method (FEM) for the conductivity equation. These correspond to discrete versions of the inverse conductivity problem of Calder\'on. We characterize FEM models corresponding to a given triangulation of the domain that are equivalent to certain resistor networks, and apply the results to study nonuniqueness of the discrete inverse problem. It turns out that the degree of nonuniqueness for the discrete problem is larger than the one for the partial differential equation. We also study invisibility cloaking for FEM models, and show how an arbitrary body can be surrounded with a layer so that the cloaked body has the same boundary measurements as a given background medium.

1. Introduction 1.1. FEM and resistor networks. In this paper we consider inverse problems for FEM models of the (anisotropic) conductivity equation, and for resistor networks. Let Ω ⊂ R 2 be the interior of a simple polygon, let σ ∈ L ∞ (Ω, R 2×2 ) be a symmetric positive definite matrix function, and consider the Dirichlet problem for the conductivity equation If one is given a triangulation of Ω, then a Finite Element Method (FEM) model, based on piecewise affine basis functions, for this Dirichlet problem is a matrix equation where u = (u 1 · · · u N ) t ∈ R N corresponds to values of the FEM solution at the vertices of the triangulation, f = (f 1 · · · f B ) t corresponds to the boundary value f at the boundary vertices, and the matrix (L L ) depends on σ and the triangulation and has rows whose elements add up to zero.
On the other hand, if one considers the triangulation of Ω as a graph, and if to each edge of the graph one assigns a resistor with a given conductivity, then the Dirichlet problem for this resistor network is to find a solution u = (u 1 · · · u N ) t , corresponding to voltages at each vertex, of the matrix equation where f = (f 1 · · · f B ) t corresponds to voltages at the boundary vertices. Also here, the rows of the matrix ( L L ), which depends on the resistors, add up to zero. Thus, formally, the FEM model and resistor networks lead to the same kind of matrix equation. This has also been pointed out in [59], see also [1]. In this paper we will show that there is a sort of equivalence between the two. Given an anisotropic conductivity in a triangulated domain Ω, one can assign resistors (with possibly negative conductivity) to each edge in the triangulation, in such a way that the equations for the FEM model and resistor network are the same (and therefore have the same solution). Conversely, given a resistor network corresponding to a triangulated domain Ω ⊂ R 2 , one can find an anisotropic conductivity (not necessarily positive definite) which is constant in each triangle, such that the equations for the resistor network and the FEM model again coincide. Moreover, in the two dimensional case we characterize all such conductivities for the FEM model that correspond to the given resistor network. In the three dimensional case the number of degrees of freedom in the FEM and resistor networks are analyzed in [1].
In this paper, we actually do not consider the inverse problem of recovering a conductivity from boundary measurements in the discrete models. Rather, we will use the equivalence of FEM models and resistor networks to describe a new type of nonuniqueness in the inverse problem for the FEM model. Our results indicate that the degree of nonuniqueness in the FEM model is strictly larger than that in the continuous case, and might be different from the case of the related Finite Volume Method (FVM) which will be discussed below. We will also describe invisibility cloaking constructions in the discrete case. These results highlight a certain difference when compared with known results for the continuum inverse problem, i.e. the inverse problem for the conductivity equation ∇ · σ∇u = 0, which is a partial differential equation (PDE). To discuss the difference, let us first recall some results related to the PDE problem.
1.2. Inverse problem for the conductivity equation. The inverse conductivity problem, also known as Calderón's inverse problem, is the question whether an unknown conductivity distribution inside a domain can be determined from voltage and current measurements made on its boundary. The measurements correspond to the knowledge of the Dirichlet-to-Neumann map Λ σ (the continuous voltage-to-current map) associated with a matrix conductivity σ = (σ jk ) 2 j,k=1 . This map takes the Dirichlet boundary values of the solution of the conductivity equation to the corresponding Neumann boundary values, In the classical theory of the problem, the conductivity σ is measurable and bounded uniformly from above and below by positive constants. We say that the conductivity σ is isotropic if it is a scalar function times the identity matrix, and otherwise anisotropic. We first discuss known results for isotropic conductivities. The inverse conductivity problem was proposed by Calderón [12] in 1980. In dimensions n ≥ 3, Sylvester and Uhlmann [74] proved unique identifiability of conductivities which are C 2 smooth, and Haberman and Tataru [40] proved uniqueness for C 1 conductivities. For earlier results for nonsmooth conductivities, see [33,69]. In two dimensions, the first global uniqueness result is due to Nachman [65] for conductivities with two derivatives by the ∂ technique. Astala and Päivärinta [7] established uniqueness for L ∞ conductivities which are bounded from above and below by positive constants. See [75] for further references, and for numerical implementations see [71,63,64,47].
In the case of anisotropic, i.e. matrix-valued, conductivities that are uniformly bounded from above and below, there is a known obstruction to uniqueness in the Calderón problem [50]. If Φ : Ω → Ω is a diffeomorphism with Φ| ∂Ω = Id, it follows that where Φ * σ is the pushforward of σ under Φ (DΦ is the Jacobian matrix), .
In two dimensions [73,65,53,4,34] it is known that the Dirichlet-to-Neumann map Λ σ determines a regular conductivity tensor σ up to such a diffeomorphism. This means that one can obtain an image of the interior of Ω in deformed coordinates. Thus, the inverse problem is not uniquely solvable, but the nonuniqueness of the problem can be characterized. The same problem in dimensions n ≥ 3 is open, see [55,53,52,24,22,23] for partial results. The inverse conductivity problem in the two dimensional case has also been studied with degenerate, that is, non-regular, conductivities. Such conductivities appear in physical models where the medium varies continuously from a perfect conductor to a perfect insulator. As an example, we may consider the case where the conductivity goes to zero or to infinity near ∂D where D ⊂ Ω is a smooth open set. The paper [5] characterizes the kinds of degeneracies for which one has nonuniqueness in the inverse problem (the limit of visibility). The same paper also determines the cases where it is even possible to coat of an arbitrary object so that it appears like a homogeneous body in all static boundary measurements (the limit of invisibility cloaking). Thus, the degree of nonuniqueness in the inverse conductivity equation is quite well understood in two dimensions. Surprisingly, in the discrete problems that we study in this paper a new type of nonuniqueness appears.
1.3. Equivalence of FEM models and resistor networks. We will now state the first main theorems of this paper concerning the equivalence of FEM models and resistor networks. To keep this introduction brief we will not describe the discrete models in full detail, but rather refer to Sections 2 and 3 for precise definitions. For the FEM model, let Ω ⊂ R 2 be the interior of a simple polygon and let T be a triangulation of Ω. We consider conductivity matrices from the class P C(T ) = {σ ∈ L ∞ (Ω, R 2×2 ) ; σ is a symmetric constant matrix in the interior of each triangle}, and also from the class P C + (T ) which consists of those matrices in P C(T ) that are strictly positive definite in each triangle. Below, we define a FEM model to be the pair (T , σ) where T is a triangulation of a simple polygon and σ ∈ P C(T ). The FEM model for the conductivity equation leads to an equation of the form (1.1), where the matrix on the left hand side will be denoted by L(σ). If L(σ) is invertible, we say that 0 is not a Dirichlet eigenvalue of the FEM model (if σ ∈ P C + (T ) this is always true). In this case, there are natural boundary measurements for the FEM model that are given by the Dirichlet-to-Neumann map where B is the number of boundary vertices.
We next consider resistor networks. By definition, a resistor network is a graph (V, E) together with a distinguished set V B ⊂ V a map γ : E → R.
Here V = {x 1 , . . . , x N } is some finite set, called the set of vertices, E ⊂ {{x, y} ; x, y ∈ V } is the set of edges, V B is the set of boundary vertices, and γ(e) is thought of as the conductivity of a resistor that occupies the edge e. The Dirichlet problem for the resistor network is a matrix equation of the form (1.2), where the matrix on the left is denoted by L(γ). If L(γ) is invertible we say that 0 is not a Dirichlet eigenvalue of the resistor network, and one has boundary measurements for the resistor network given by a Dirichlet-to-Neumann map It turns out that there is a canonical way of going from a given FEM model to a resistor network having the same boundary measurements. Theorem 1.1. Let Ω ⊂ R 2 be the interior of a simple polygon, let T be a triangulation of Ω having vertices V , boundary vertices V B , and edges E, and let σ ∈ P C + (T ). There is a unique function γ such that the resistor Even if any FEM model with given triangulation naturally gives rise to a resistor network, there are many ways of going from the resistor network back to the FEM model. This may be described by going from the resistor network to a "double resistor network" where the resistor corresponding to each interior edge has been split in two parts. The splitting is determined by a real valued function α defined on the set of interior edges. It turns out that there is an essentially unique way of going from this double resistor network back to a FEM model with conductivity matrix σ α that is constant on each triangle of the original triangulation and has the same Dirichlet-to-Neumann map as the original FEM model. The detailed construction of σ α is given in Section 4. We emphasize that the conductivities arising in this way are not always positive definite, but we will also give conditions that guarantee positive definiteness (see Lemma 4.3).
be a resistor network that corresponds to a triangulation T , and assume that 0 is not a Dirichlet eigenvalue. For any real valued function α defined on the set of interior edges, there is a symmetric matrix function σ α ∈ P C(T ) such that L(σ α ) = L(γ) and Λ FEM σα = Λ γ . Moreover, any σ ∈ P C(T ) that satisfies L(σ) = L(γ) and Λ FEM σ = Λ γ is of the form σ = σ α for some such function α.
We mention here that inverse problems for resistor networks have been extensively studied using graph theoretical methods [9,10,19]. In [9,10] the authors modeled the Finite Volume Method with isotropic conductivities by using resistor networks. A closely related discretized version of the inverse conductivity problem, combined with regularization and an analysis related to homogenization theory, has been studied in [20]. Our work differs from these results in that we study the nonuniqueness which can arise in the FEM model for anisotropic conductivities.
The problems studied in [9,10,19] focused on two different types of problems. First, one can model a body with a fixed resistor network, e.g. a square lattice corresponding to Euclidean coordinates or a circular lattice corresponding to polar coordinates [15,16,17,18,35] and ask whether one can reconstruct the resistors from boundary measurements. Second, one can consider arbitrary graphs with a subset of the nodes declared as the boundary nodes and use the fact that the resistor network can be modified by applying a so-called ∆ − Y transformation (or the triangle-star transformation) in the graph without changing the boundary measurements. In this transformation the graph structure, in particular the number of the nodes, changes. In these studies one considers e.g. the planarity of the graph and the characterization of equivalence classes of the graphs that appear the same in boundary measurements [15,16,17].
Our paper falls closer to the first category of problems as we consider triangulations of a given domain which correspond to a fixed abstract graph. The results of the paper are also related to inverse problems for quantum graphs, see e.g. [3,51].
1.4. Nonuniqueness for FEM models. To consider nonuniqueness in the inverse problem for FEM models, let us start from the familiar nonuniqueness related to the diffeomorphism invariance of the conductivity equation ∇ · σ∇u = 0. As pointed out above, the Dirichlet-to-Neumann maps for σ and a pushforward conductivity Φ * σ coincide if Φ : Ω → Ω is a diffeomorphism fixing the boundary. If one approximates the conductivity equation by the FEM model and considers the corresponding resistor network, then the resistor network is an abstract graph which may be embedded in R 2 in several ways. These different embeddings (which may for instance move the interior vertices around) can be thought to correspond to the diffeomorphisms Φ, and the use of resistor networks "factors out" the diffeomorphism. Thus the graph structure associated with the resistor network can be considered as the analogue of the underlying Riemannian manifold structure associated with the anisotropic conductivity, cf. [52,53,55]. Since the diffeomorphism invariance is inherent in the continuous problem, similar type of nonuniqueness occurs when using FVM to study anisotropic conductivities.
We will next formalize the analogue of diffeomorphism nonuniqueness. Two triangulations T and T of Ω are said to be isomorphic via a map Φ if Φ : Ω → Ω is an orientation preserving homeomorphism that gives a 1-1 correspondence between triangles of T and T and is affine on each triangle. If σ ∈ P C + (T ) and if Φ is an isomorphism of triangulations T and T , we define the pushforward of σ via Φ by where y is in the interior of some triangle of T (see Fig. 1).
Theorem 1.3. Let T and T be two triangulations of Ω ⊂ R 2 that are isomorphic via some Φ : Ω → Ω with Φ| ∂Ω = Id. If σ ∈ P C + (T ), then we have Thus diffeomorphism invariance in the FEM model looks basically the same as in the continuum model. However, there is another type of nonuniqueness in the discrete model which is not present in the continuous case. As mentioned above, any FEM model with given triangulation naturally gives rise to a resistor network, but there are many ways of going from the resistor network back to the FEM model. This involves going from the resistor network to a "double resistor network" where the resistor corresponding to each interior edge has been split in two parts. The different splittings will be parametrized by real valued functions β on the set of interior edges, and these splittings account for the other type of nonuniqueness.
To make this precise, we consider the directed graphs associated with triangulations. Let T be a triangulation of Ω and let σ ∈ P C + (T ). The nondirected graph associated with T is the pair (V, E) where V = {x 1 , . . . , x N } is the set of all vertices of triangles in T and E is the set of non-directed edges. The directed graph associated with T is the pair (V, D) where V is the set of vertices as above, and D is the set of all ordered pairs (x j , x k ) such that the line segment from x j to x k is a positively oriented edge of some triangle in T . (The triangles of T and their boundaries have an orientation induced by the standard orientation of R 2 .) We write E B for the set of non-directed boundary edges (those edges that belong to only one triangle), and E I for the set of non-directed interior edges (those edges that belong to exactly two triangles). We assume that the elements of V are labeled so that the set of boundary vertices (those vertices that lie on ∂Ω) is Let β : E I → R be a real valued function, and let σ be the matrix function in Ω satisfying the following conditions: (i) σ is a constant matrix in the interior of each triangle in T , (ii) if e = {x j , x k } ∈ E I and j < k, and if T and T are the two triangles of T having both x j and x k as vertices such that ∂T has the same orientation as the directed edge d = (x j , x k ) (so then ∂T has the opposite orientation from d), then , and if T is the unique triangle which has both x j and x k as vertices, then It will turn out (see Section 4) that there is indeed a unique matrix function σ in Ω that satisfies the above conditions, and we write S β (σ) = σ. Clearly S β (σ) ∈ P C(T ), and we also have S 0 (σ) = σ. Note that here the conductivity S β (σ) may not be positive definite, even if σ ∈ P C + (T ), but it is always such that 0 is not a Dirichlet eigenvalue of the FEM problem and there is a well defined Dirichlet-to-Neumann map. The next result makes precise the second type of nonuniqueness that is present in the discrete model. Theorem 1.4. Let T be a triangulation of Ω and let σ ∈ P C + (T ). For any real valued function β on E I , there exists a unique conductivity matrix σ = S β (σ) ∈ P C(T ) that satisfies (i)-(iii). One also has . We note that in the deformations considered in Theorem 1.4, the graph structure associated with the conductivities σ and S β (σ) is the same. Finally, it is instructive to state an "infinitesimal" version of the above two results.
If t is sufficiently small, let T t be the triangulation of Ω that is isomorphic to T via Φ t . Let also β be a map E I → R. Then In this result, one may think that Φ t moves around the vertices of the triangulation slightly, and t → S tβ (σ) deforms the conductivity slightly according to the double resistor network picture (see Fig. 2). We emphasize that if σ is positive definite and t > 0 is small enough, then also the matrix (Φ t ) * S tβ (σ) is positive definite. Roughly speaking, this shows that the degree of the nonuniqueness in inverse problems for discrete FEM models, even when the graph structure associated with the triangulation is fixed, is strictly larger than than the one for the conductivity equation considered as a partial differential equation.
It remains an open problem to see whether the types of nonuniqueness described in Theorems 1.3 and 1.4 are the only ones admitted by a FEM discretization. In particular, the analysis in this paper, based on FEM models, does not allow to take full advantage of the powerful theory of inverse problems on graphs studied in [15,16,17,18,19,21,35]. We note that these techniques give a characterization for the equivalence class of the resistor networks that correspond to the given boundary measurements. However, this characterization is based on the ∆ − Y transformations that may change triangular graphs to non-triangular ones. Due to this, at this moment we do not have an effective characterization of the triangular resistor networks (corresponding to the FEM models studied in this paper) that correspond to the given boundary measurements.
A closely related technique to FEM is the Finite Volume Method (FVM) which was used in [9,10] to study the reconstruction of isotropic conductivities in Euclidean and polar coordinates. As we have indicated above, the goal of this paper is to study the nonuniqueness issues and cloaking phenomenon for anisotropic models in FEM. Both FVM and FEM have their own merits and are used in different applications. In FVM the flux entering a given element is identical to that leaving the element and thus FVM is a conservative method for flux. On the other hand, FEM is compatible to the weak formulation of the elliptic equation making it easy to analyze the convergence of the approximation when discretization becomes finer. The FEM approach is also easily adapted to irregularly shaped meshes and electrode models. For these reasons, FEM has been used extensively in the study of the inverse conductivity problem using statistical inference and regularization theory, see e.g. [38,39,44,45] and references therein.
Unlike FEM, FVM is based on balancing the fluxes through interfaces of the elements. Thus it is an interesting question if inverse problems based on the FVM discretization of the anisotropic conductivity equation have a smaller degree of nonuniqueness than inverse problems based on the FEM discretization. It seems that the nonuniqueness related to the diffeomorphisms (see Fig. 1) appears also for the FVM models, but it is not clear whether there is an analogy of the nonuniqueness related to the splitting of the resistors (see Fig. 2) for the FVM models, or not. However, this question is outside the context of this work.
1.5. Discrete nonuniqueness and cloaking. It is a well known fact that the problem of determining an anisotropic conductivity is unique only up to pull-backs by diffeomorphisms. In the discrete setting, this corresponds to different embeddings of a planar graph into R 2 . When considering a FEM model and its corresponding resistor networks there is an additional nonuniqueness. Namely, the piecewise constant conductivity on each triangle can be distributed differently along the adjacent edges. In the last section of the paper we will also discuss discrete cloaking layers where any resistor network Ω * can be embedded in a cloaking resistor network without being detected and the corresponding FEM models.
Structure of paper. Section 1 is the introduction and contains the statements of the main results. In Section 2 we describe the FEM model for the conductivity equation and also prove Theorem 1.3. Resistor networks are discussed in Section 3. The correspondence between FEM models and resistor networks is explained in Section 4, where also Theorems 1.1, 1.2, 1.4 and 1.5 are proved. The final section describes invisibility cloaking constructions for FEM models and resistor networks. . . , T R } where each T j is a closed triangle in R 2 , and two triangles T j and T k for j = k are either disjoint or their intersection is one vertex or one full edge of each triangle. The domain of T is the set and we have Ω = R j=1 T j . Two triangulations T = {T 1 , . . . , T R } and T = {T 1 , . . . , T R } with domains Ω and Ω are said to be isomorphic if R = R and there is an orientation preserving homeomorphism Φ : Ω → Ω that is affine on each T j and maps each T j onto some T j (the map Φ is called an isomorphism).
We denote the set of all vertices of triangles T ∈ T by V = {x 1 , . . . , x N }. We assume that the points are labeled so that Given a triangulation T of a domain Ω, we will consider conductivity matrices from the class in the interior of each triangle}.
We also write P C + (T ) for the subspace of matrices that are positive definite on each triangle.
Below, we say that v : Ω → R is a piecewise affine function (associated with the triangulation T ) if v is a continuous function on Ω and the restriction v j | T k is an affine function for each triangle T k ∈ T . The set of piecewise affine functions is denoted by A(Ω) = A T (Ω). Moreover, we denote by v k , k = 1, 2, . . . , N , the specific basis of the space A(Ω) for which v k has value 1 at the vertex x k and vanishes at all other vertices x j , j = k. Also, write for the space of piecewise affine functions (associated with the triangulation T ) on the boundary.
Next we will formulate Dirichlet problem for the FEM model. To motivate it, consider the Dirichlet problem for the anisotropic conductivity equation where Ω ⊂ R 2 is a bounded domain, and σ ∈ L ∞ (Ω; R 2×2 ) is a symmetric positive definite matrix. Weak solutions are functions in Let Ω be the interior of a simple polygon, and fix a triangulation of Ω with vertices {x 1 , . . . , x N } of which {x 1 , . . . , x B } lie on ∂Ω. Let A(Ω) be the set of continuous functions Ω → R 2 which are affine on each triangle. Recall that v k (x) is the unique function in A(Ω) which is 1 at x k and 0 at all the other vertices and Note that A(Ω) ⊂ H 1 (Ω) and the gradient of a function in A(Ω) is constant in the interior of each triangle. If f ∈ A(∂Ω), the FEM approximation to (2.1) asks to find a function u ∈ A(Ω) such that We refer to [14] for more information on FEM approximations.

FEM model as a matrix equation.
Let f ∈ A(∂Ω) be a given piecewise affine function on the boundary and let u ∈ A(Ω) be the solution of equation (2.2). Writing f ( We write L = L(σ) = (l jk (σ)) N j,k=1 for the matrix above. Next we record some well-known properties of the matrix L.
where L is symmetric. The rows of L L add up to zero: If σ is positive definite (in the sense that σ(x)ξ · ξ ≥ c|ξ| 2 for a.e. x ∈ Ω and for all ξ ∈ R 2 where c > 0), then L is positive definite (thus invertible) and l jj > 0 for j ≥ B + 1.
Proof. Clearly L has the above form, and if j ≥ B + 1 then Thus also L is positive definite and invertible, and for j ≥ B + 1 we have It will be useful to consider the FEM approximation for conductivities which may not be positive definite. If σ ∈ L ∞ (Ω; R 2×2 ) is a symmetric matrix, we say that 0 is not an eigenvalue of the FEM problem if the matrix L = L(σ) is invertible. In this case, for any f = (f 1 · · · f B ) t the FEM problem (2.3) has a unique solution u = L −1 f , and one has a well defined Dirichlet-to-Neumann map Λ FEM σ defined below.
For the relation to resistor networks below, one might ask for the condition l jk ≤ 0 for j = k (this would ensure that no resistor has negative conductivity). This is not true in general, but in the following case the non-negativity of elements of the matrix L can be guaranteed.
Lemma 2.2. If σ(x) is isotropic, i.e. a scalar function multiplied with the identity matrix, and if the triangulation of Ω is such that no triangle has an angle > π/2, then l jk ≤ 0 for j = k.
Proof. It is enough to consider l jk for j ≥ B + 1 and k = j. One has l jk = 0 unless x j and x k are neighboring vertices, and in that case there are exactly two triangles T 1 and T 2 which have both x j and x k as vertices. Since no angle is > π/2, elementary geometric considerations imply that in these two triangles ∇v j · ∇v k ≤ 0. Thus also σ∇v j · ∇v k ≤ 0 in these triangles.

2.3.
Dirichlet-to-Neumann map for FEM model. We will next define a Dirichlet-to-Neumann map Λ FEM σ related to the FEM model. Assume that 0 is not a Dirichlet eigenvalue for the FEM model, so (2.2) has a unique solution with all f ∈ A(∂Ω). We define the Dirichlet-to-Neumann map Λ FEM σ : A(∂Ω) → A(∂Ω) to be the unique map given by the pairing where u ∈ A(Ω) is the solution of (2.2), and w is any function in A(Ω) with w| ∂Ω = g. Here we have written Writing u(x) = N j=1 u j v j (x), u j = u(x j ) ∈ R, we obtain an explicit formula for the Dirichlet-to-Neumann map: 2.4. Diffeomorphism invariance. If σ ∈ P C(T ) and if Φ is an isomorphism of two triangulations T and T , we define the pushforward of σ via Φ by where y is in the interior of some triangle of T . Since Φ is affine on each triangle, it follows that Φ * σ ∈ P C(T ). Moreover, if σ ∈ P C + (T ) then Φ * σ ∈ P C + (T ). and therefore (writing ∇ as a column vector) The result follows.
Proof of Theorem 1.3. This is just a special case of the previous lemma where σ ∈ P C + (T ).

Resistor networks
To properly state the equivalence of the FEM model and resistor networks, we need some concepts from graph theory (see [21]). By a graph we mean a set (V, E) where V (vertices) is a finite set, and E (edges) is a set of unordered pairs of elements in V . A path from x ∈ V to y ∈ V is a sequence of vertices (x 0 , x 1 , . . . , x k ) where x 0 = x, x k = y, and {x j , x j+1 } ∈ E for 0 ≤ j ≤ k − 1. We only consider graphs which are connected (any two vertices can be joined by a path) and simple (no self-loops, i.e. E has no elements of the form {x, x}).

Definition.
A resistor network is a graph (V, E) together with a distinguished set V B ⊂ V and a function γ : E → R.
We will label elements of V so that V B = {x 1 , . . . , x B } (the set of boundary vertices) and V V B = {x B+1 , . . . , x N } (the set of interior vertices). If e = {x j , x k } ∈ E we write γ jk := γ(e) (the conductivity of a resistor located on the edge e), with the convention that γ jk = 0 if {x j , x k } / ∈ E. Given a function f : V B → R, we would like to find a function u : V → R satisfying Here N (x) = {y ∈ V ; {x, y} ∈ E} is the set of neighboring vertices of x. If the graph is such that (V, E) can be drawn as a straight-line graph in the plane, by Kirchhoff's law the function u corresponds to a voltage potential in the resistor network induced by the voltage f applied at the boundary vertices.
Writing f (x j ) = f j and u(x j ) = u j , the problem (3.1) corresponds to the matrix equation where u = (u 1 · · · u N ) t and f = (f 1 · · · f B ) t , and L = ( l jk ) N j,k=1 with It is clear that the matrix L has the form where L is symmetric, and for j ≥ B + 1 one has N k=1 l jk = 0.
Further, if γ(e) > 0 for all e ∈ E then l jj > 0 for all j.  Proof. Since the conductivities are positive, L is diagonally dominant: This implies that L is positive semidefinite, as a symmetric diagonally dominant matrix with nonnegative diagonal entries (this is a consequence of Gershgorin's circle theorem). For invertibility of L we note that if there is strict inequality in If L is invertible, we say that 0 is not an eigenvalue of (3.1). In this case it is possible to define the Dirichlet-to-Neumann map of the resistor network as the operator acting on functions f : V B → R by where u is the unique solution of (3.1).

Correspondence of FEM models and resistor networks
We wish to relate the FEM model to resistor networks. To do this we need to consider resistor networks which can be drawn as the triangulation of a simple polygon in R 2 . A (planar straight-line) drawing of a graph (V, E) is a pair of maps ϕ : V → R 2 and ψ : E × [0, 1] → R 2 , where ϕ is injective and ψ satisfies the following conditions: (1) if e = {x, y} ∈ E then ψ e = ψ(e, · ) is a line segment with end points ϕ(x) and ϕ(y), (2) ψ e ([0, 1]) does not intersect ϕ(V ) except at the end points, (3) for any e, e ∈ E with e = e the open line segments ψ e ((0, 1)) and ψ e ((0, 1)) do not intersect.
The components of R 2 ψ(E × [0, 1]) are called the faces of the graph (there is exactly one unbounded face). The number of faces is independent of the embedding by Euler's formula. A graph which admits a drawing is called planar. Planarity can be checked directly from the abstract graph (for instance Kuratowski's condition). Also, there is no loss of generality in using straight-line embeddings by Fary's theorem. See [21] for these facts.

Definition. A triangular resistor network is a resistor network (V, E) such that
(1) there is a drawing of (V, E) for which every face (except the unbounded one) is a triangle, (2) the interior of the complement of the unbounded face, which we denote by Ω (the domain), is such that Ω is connected, Remark. There is a notion of triangulation in graph theory which only uses the abstract graph, but this seems to be different from what we need here.
If σ is an anisotropic conductivity in a triangulated set Ω, we denote by L(σ) the matrix in (2.3). Also, if (V, V B , E, γ) is a resistor network we denote by L(γ) the matrix in (3.2). It is easy to see that anisotropic conductivities give rise to resistor networks in this way. Let Ω ⊂ R 2 be the interior of a simple polygon, and let σ ∈ L ∞ (Ω; R 2×2 ) be a symmetric matrix function such that 0 is not an eigenvalue of (2.2). If a triangulation of Ω with vertices V , boundary vertices V B = V ∩ ∂Ω, and edges E is given, then there is a unique function γ such that the triangular resistor network R = (V, V B , E, γ) satisfies L(γ) = L(σ) and Λ γ = Λ FEM σ . Proof. With the notation in Section 2, we choose Note that Ω σ∇v j · ∇v k dx = 0 if x k / ∈ N (x j ). Using the fact that N k=1 Ω σ∇v j · ∇v k dx = Ω σ∇v j · ∇(1) dx = 0, we have l jj = x k ∈N (x j ) γ jk = Ω σ|∇v j | 2 dx and thus L(γ) = L(σ). For the equality of Dirichlet-to-Neumann maps we also choose This implies x k ∈N (x j ) γ jk = Ω σ|∇v j | 2 dx as above, so Λ γ = Λ FEM σ .
Proof of Theorem 1.1. This is a special case of the previous theorem when σ ∈ P C + (T ).
Remark. If σ is a positive definite matrix then L(γ) is positive definite and the diagonal entries l jj are positive, but apparently some of the conductivities γ jk may be negative unless additional restrictions are imposed.
In the converse direction there are many conductivities σ corresponding to a given resistor network. However, if one restricts to conductivities which are constant inside each triangle, it is possible to characterize all such σ in terms of directed resistor networks. To describe the characterization, we let (V, E) be a resistor network and consider a splitting into two parts of all the resistors associated with interior edges of E. Informally, this is done as follows: (1) Let E I = {e ∈ E ; at least one endpoint of e is an interior vertex}. We now describe the construction of piecewise constant conductivities which correspond to a given triangular resistor network (V, E, V B , γ) with domain Ω ⊂ R 2 . Let D I be some choice of directed edges corresponding to E I as above, and let α be a function D I → R. We wish to find a symmetric matrix function σ = σ α satisfying the following conditions: (c) if x j and x k are neighboring boundary vertices, and if T is the unique triangle which has both x j and x k as vertices, then The next result shows that such a conductivity σ = σ α always exists and it is uniquely determined by the above conditions. Proof. Let α be a function D I → R. We define α also for boundary edges as follows: if e = {x j , x k } ∈ E where both x j and x k are boundary vertices, let T be the unique triangle having both x j and x k as vertices, let e be the directed edge corresponding to e such that ∂T has the same orientation as e , and define α(e ) = 1 2 γ jk . We first show the existence of σ α satisfying (a)-(c). Choose any triangle T , and let x j , x k , x l be its vertices. If e jk is the directed edge corresponding to e jk = {x j , x k }, define if ∂T has the same orientation as e jk , 1 2 γ jk − α(e jk ) otherwise. (4.1) The numbers µ kl and µ lj are defined analogously. Consider the three equations Since σ and each ∇v m are constant in T , these are equivalent with σ∇v l · ∇v j = −µ lj /|T |.
Writing σ = a b b c in T for some real constants a, b, c, (4.3) becomes a linear system of three equations in three unknowns (the constants a, b, c). The corresponding homogeneous system is σ∇v j · ∇v k = 0, σ∇v k · ∇v l = 0, σ∇v l · ∇v j = 0. This only has the zero solution σ = 0, since for instance ∇v k and ∇v l are linearly independent so for any w ∈ R 2 one obtains σ∇v j · w = σ∇v j · (w 1 ∇v k + w 2 ∇v l ) = 0.
Consequently σ∇v j = 0, and analogously σ∇v k = σ∇v l = 0 which implies σ = 0. This shows that the original system (4.2) has a unique solution, and so we obtain a constant matrix σ in T satisfying (4.2).
We denote by σ = σ α the matrix in Ω which is defined by the above procedure in each triangle. Clearly (a) is satisfied, and also (b) and (c) are valid by the definition of σ α . Also, if σ is another matrix satisfying (a)-(c) then σ − σ satisfies (a)-(c) with all right hand sides replaced by zero, and localizing to each triangle and using the uniqueness of solutions to 3 × 3 systems as above implies that σ = σ.
Note that one has L(σ) = L(γ) and Λ FEM This statement follows by writing out the definitions of the matrices and Dirichlet-to-Neumann maps and using that N m=1 Ω σ∇v j · ∇v m dx = 0. To prove (4.4), if at least one of x j , x k is an interior vertex, there are two triangles T and T having both x j and x k as vertices. We choose T to be the triangle such that ∂T has the same orientation as the directed edge e corresponding to {x j , x k }. One has by the construction of σ α . On the other hand, if both x j and x k are boundary vertices, there is a unique triangle T having both these points as vertices, and again by construction Thus L(σ) = L(γ) and Λ FEM σ = Λ γ . Finally we prove that for any symmetric matrix σ which is constant in each open triangle and satisfies L(σ) = L(γ) and Λ FEM σ = Λ γ , one has σ = σ α for some α : D I → R. This is easy at this point. The assumption implies that (4.4) holds. Let e ∈ D I , and let e = {x j , x k } be the corresponding undirected edge where x j is an interior vertex. There are two triangles T and T having both x j and x k as vertices, chosen so that ∂T has the same orientation as e . Note that by (4.4) T σ∇v j · ∇v k dx + T σ∇v j · ∇v k dx = −γ jk .
We define It follows that with this definition of α, σ satisfies conditions (a) and (b) above. If x j and x k are boundary vertices, there is a unique triangle T having both these points as vertices and therefore by (4.4) Thus σ also satisfies (c). Since σ is uniquely determined by (a)-(c), it follows that σ = σ α .
Proof of Theorem 1.2. Theorem 1.2 follows from the previous theorem after fixing some choice of directed interior edges D I corresponding to the set E I of non-directed interior edges, and after identifying the function α : Remark. An isomorphism of two directed graphs (V, D) and (V , D ) is a bijective map φ : V → V for which (x, y) ∈ D if and only if (φ(x), φ(y)) ∈ D . Note that an isomorphism Φ of two triangulations T and T induces an isomorphism of the associated directed graphs, and conversely two triangulations T and T whose associated graphs are isomorphic via some φ are isomorphic via a map Φ that is uniquely determined by ϕ. Indeed, the map Φ is uniquely determined by its values at the vertex points.
The next result gives a condition that characterizes when the matrix σ α above is positive semidefinite. Lemma 4.3. The conductivity matrix σ constructed in Theorem 1.2 is positive semidefinite if and only if on each triangle T with vertices x j , x k , and x l , the conductivities µ j,k , µ l,j and µ k,l that appear in the proof of the theorem and are associated with the directed edges of the triangle are chosen to satisfy Proof. Let T be a triangle with vertices which we assume without loss of generality are x 1 , x 2 , and x 3 . Assume first the special case that T is a a equilateral triangle, and choose coordinates so that where l > 0. Then |T | = l 2 / √ 3, and in T we have Solving the equation (4.3) for σ = a b b c we get that With these entries, the trace and determinant of σ satisfy tr(σ) = 4 √ 3 (µ 1,2 + µ 2,3 + µ 3,1 ), det (σ) = 4(µ 1,2 µ 2,3 + µ 1,2 µ 3,1 + µ 2,3 µ 3,1 ).
Since σ is real and symmetric, we have that σ is positive semidefinite if and only if the trace and determinant are ≥ 0. This is equivalent with the conditions in the statement. Now for a general triangle T , find an affine diffeomorphism F such that T := F (T ) is an equilateral triangle. Then make a change of variables using To end this section, we show how the double resistor network construction can be used to prove Theorems 1.4 and 1.5.
Proof of Theorem 1.4. Assume that we are given a triangulation T of a domain Ω ⊂ R 2 that has associated graph (V, E) with interior edges E I and boundary edges E B , and that we are also given a conductivity matrix σ ∈ P C + (T ) and a function β : E I → R.
Recall from Section 1 that we are looking for a conductivity σ = S β (σ) that satisfies the following three conditions: (i) σ is a constant matrix in the interior of each triangle in T , (ii) if e = {x j , x k } ∈ E I and j < k, and if T and T are the two triangles of T having both x j and x k as vertices such that ∂T has the same orientation as the directed edge d = (x j , x k ) (so then ∂T has the opposite orientation from d), then , and if T is the unique triangle which has both x j and x k as vertices, then We define a set of directed edges corresponding to E I by We will then choose σ = σ α as in Theorem 4.2 for some function α : D I → R that will depend on β.
Let us begin by constructing the resistor network (V, E, V B , γ) that is equivalent with the FEM model as in Theorem 4.1 and satisfies L(γ) = L(σ) and Λ γ = Λ FEM σ . This was obtained by taking (See the proof of Theorem 4.1.) Let now e = {x j , x k } ∈ E I where j < k, and let T and T be the two triangles of T having both x j and x k as vertices such that ∂T has the same orientation as the directed edge e = (x j , x k ) (so ∂T has the opposite orientation from e ). Looking at the supports of ∇v j and ∇v k , we have We make the choice By using the expression for −γ jk above, we have The above construction defines a function α : D I → R, and Theorem 4.2 guarantees that there is a conductivity σ α ∈ P C(T ) that satisfies conditions S tβ (σ) . If Φ t is as in the statement of the theorem and if t is sufficiently small, then there is a triangulation T t of Ω that is isomorphic to T via Φ t . Also (Φ t ) * S tβ (σ) is in P C + (T ) for t small. Theorem 1.3 then implies that This proves the result.

Cloaking for discrete and continuous models
By invisibility cloaking one means the possibility, both theoretical and practical, of shielding a region or object from detection via electromagnetic fields. The invisibility cloaking examples, see [25,26,32,62,61,56,68], have been the starting point for transformation optics, a very topical subject in recent studies in mathematics, physics, and material science.
To consider invisibility cloaking results for discrete models, let us recall the corresponding results for the anisotropic conductivity equation (1.3) in Ω ⊂ R 2 , where the conductivity σ = [σ jk (x)] 2 j,k=1 is a measurable function which values are symmetric, positive definite matrixes. We say that a conductivity σ is regular if there are c 1 , c 2 > 0 such that the matrix σ(x) satisfies If conductivity is not regular, it is said to be degenerate.
Let us first recall the counterexamples for the uniqueness of the inverse problem that have a close connection to the invisibility cloaking. In 2003, before the appearance of practical possibilities for cloaking, it was shown in [31,32] that passive objects can be coated with a layer of material with a degenerate conductivity which makes the object undetectable by the electrostatic boundary measurements. These constructions were based on the blowup maps and gave counterexamples for the uniqueness of inverse conductivity problem in the three and higher dimensional cases. In two dimensional case, the mathematical theory of the cloaking examples for conductivity equation have been studied in [48,49,54,66].
The interest in cloaking was raised in particular in 2006 when it was realized that practical cloaking constructions are possible using so-called metamaterials which allow fairly arbitrary specification of electromagnetic material parameters. The construction of Leonhardt [56] was based on conformal mapping on a non-trivial Riemannian surface. At the same time, Pendry et al [68] proposed a cloaking construction for Maxwell's equations using a blow-up map and the idea was demonstrated in laboratory experiments [70]. Cloaking constructions for the conductivity equation based on the mathematical models in [49,31] have also been recently verified experimentally in [76].
Let Σ = Σ(Ω) be the class of measurable matrix valued functions σ : Ω → M , where M is the set of symmetric positive semi-definite matrices. Instead of defining the Dirichlet-to-Neumann operator which may not be well defined for these conductivities, we consider the corresponding quadratic forms.
Definition. Let h ∈ H 1/2 (∂Ω). The Dirichlet-to-Neumann quadratic form corresponding to the conductivity σ ∈ Σ(Ω) is given by and the infimum is taken over real valued u ∈ L 1 (Ω) such that ∇u ∈ L 1 (Ω) 3 and u| ∂Ω = h. In the case where Q σ [h] < ∞ and A σ [u] reaches its minimum at some u, we say that u is a W 1,1 (Ω) solution of the conductivity problem.
In the case when σ is smooth, bounded from below and above by positive constants, we can define a Dirichlet-to-Neumann map (1.4). For such regular conductivities Q σ [h] is the quadratic form corresponding the Dirichlet-to-Neumann map, that is, where ds is the length measure on ∂Ω. Physically, Q σ [h] corresponds to the power needed to keep voltage h at the boundary. Let us next consider various counterexamples for the solvability of inverse conductivity problem with degenerate conductivities. 5.1. Invisibility cloaking for the continuous model. In this subsection we recall the invisibility cloaking example based on the coordinate transformation by a blow-up map [32]. Let us consider the general background conductivity σ. We aim to coat an arbitrary body with a layer of exotic material so that the coated body appears in measurements the same as the background conductivity σ. Usually one is interested in the case when the background conductivity σ is equal to the constant γ = 1. However, we consider here a more general case and assume that σ is a L ∞ -smooth conductivity in a domain Ω ⊂ R 2 that contains the closed disc D(2). Here, D(ρ) is an open 2-dimensional disc of radius ρ and center zero and D(ρ) is its closure. Also, assume that σ(z) ≥ c 0 I, c 0 > 0. Consider the homeomorphism given by Using the map (5.4), we define a singular conductivity where η(z) = [η jk (x)] is any symmetric measurable matrix satisfying c 1 I ≤ η(z) ≤ c 2 I with c 1 , c 2 > 0. The conductivity F * σ is called the cloaking conductivity obtained from the transformation map F and background conductivity σ and η(z) is the conductivity of the cloaked (i.e. hidden) object. In dimensions n ≥ 3 it shown in 2003 in [31,32] that the Dirichlet-to-Neumann map corresponding to H 1 (Ω) solutions for the conductivity (5.5) coincide with the Dirichlet-to-Neumann map for σ when σ = 1. In 2008, the analogous result for Dirichlet-to-Neumann maps was proven in the twodimensional case in [49]. Using the above formalism of Dirichlet-to-Neumann forms it was shown in [5] that the quadratic form Q σ coincides with Q σ . This means that when σ is a regular conductivity in the domain Ω containing the disc D(2), we can modify the conductivity in D(2) to a conductivity σ so that D(1) contains the object with conductivity η and the conductivity in D(2)\D(1) is such that for the obtained conductivity σ in Ω all measurements on ∂Ω coincide with those corresponding to σ, that is, the object η is hidden in the conducting medium and one can not detect from the boundary even the fact that something is hidden. Due to the presence of degenerate coefficients in the partial differential equation, numerical simulation of cloaking is a challenging (and interesting) problem, see e.g. [42,58,67] and references therein. Next we will consider cloaking starting directly from the FEM model for the background space and construct cloaking models using only the FEM models and the corresponding resistor networks. We study in this paper only the static (i.e. zero frequency) cloaking and leave open the problem how the corresponding constructions could be done for the Helmholtz equation or time-harmonic Maxwell's equations.

5.2.
Invisibility cloaking for the discrete model. In this subsection we construct resistor networks and their equivalent FEM models with cloaking properties. Consider first the resistor network with graph given by x 3 x 2 x 5 x 4 x 6 We call this graph consisting of four equilateral triangles the graph of type B (B refers to background) and say that (x 1 x 3 x 5 ) is the basic triangle of the graph. We consider the graph B as a resistor network by assigning resistivity 1 on all edges of the graph.
Also, we consider the second graph, shown below, that is an extension of the above graph of type B.
x 1 x 8 x 2 x 3 x 9 x 7 x 5 x 4 x 6 x 10 This graph is formed by choosing points x 10 , x 11 , and x 12 to be on the line segments connecting the points x 2 , x 6 , and x 4 to the centroid A of the triangle (x 2 , x 6 , x 4 ) so that the points x 10 , x 11 , and x 12 have the equal distances to A.
Note that in this graph no triangle has an angle > π/2. We consider also this graph as a resistor network by assigning the resistivity by the following rules. On regular lines the resistivity is 1, on dashed lines the resistivity is zero, on bold lines the resistivity is 2, and on the inner triangle the dotted edges have arbitrary strictly positive resistivity. Observe that the edge x 5 x 12 has the resistivity 1. We call this graph consisting of 16 triangles the graph of type C (C refers to cloak) and say that (x 1 x 3 x 5 ) is the basic triangle corresponding to the center triangle (x 10 x 11 x 12 ).
In the following, we consider one fixed graph of type B and another graph of type C that is an extension of the graph B, and refer to these graphs and the corresponding resistor networks and as the graph B and C and the resistor network B and C.
By using the construction given in the proofs of Theorems 1.2 and 4.2 with α j,k = γ j,k /2 we can construct a piecewise constant positive semi-definite conductivity σ C ∈ P C(T C ) with triangularization T C given by the graph C. Note that the conductivity σ C in everywhere non-zero. This gives a FEM model that corresponds to the resistor network C. Note that the constant conductivity σ B = 1 with triangularization T B given by the graph B gives a FEM model that corresponds to the resistor network B.
Proposition 5.1 means that in the FEM model corresponding to σ C the center triangle is cloaked. Below we use this conductivity to hide more general bodies.

5.3.
Cloaking of a more complicated body. Let us now consider the four triangles T 1 , T 2 , T 3 , T 4 that form the graph B with T 1 being the basic triangle and let triangle Ω ⊂ R 2 be the union of these triangles.
Similarly to the above let T B be the triangulation {T 1 , T 2 , T 3 , T 4 } of Ω and σ B be the conductivity that is constant 1 in all these triangles.
Let us also consider the graph C with the same boundary vertexes as the graph B and let T 0 denote the center triangle of the graph C. Recall that above we assigned arbitrary positive resistivities on the edges of the center Figure 3. In the figure a general body Ω * is coated with a cloaking structure. The red triangles are a triangulation of Ω * . This triangulation is extended with blue edges to form a triangulation of the triangle T 0 that is the center triangle in the resistor network corresponding to the graph C. Transforming the resistor network to a FEM model, we obtain a conductivity where the piecewise constant conductivity σ * in the domain Ω * is surrounded with a positive-semidefinite conductivity that cloaks the presence of the conductivity in Ω * . triangle of the graph C and constructed a FEM model with conductivity σ C ∈ P C(T C ).
Let Ω * ⊂ R 2 be a convex polygonal domain having some triangulation T * and assume that Ω * ⊂ T int 0 where T 0 is the center triangle of the graph C. Let σ * ∈ P C + (T * ) be an arbitrary piecewise constant conductivity in Ω * .
Let T 0 denote the triangulation of T 0 that contains all triangles in T * and all triangles in T 0 \Ω * that are obtained by connecting all boundary vertexes of Ω * to the vertexes of T 0 with such line segments that intersect Ω * only at one end point of the line segment. In Figure 3 the red triangles correspond to the triangles in T * and the blue line segments and the boundary edges of Ω * define the other triangles in T 0 .
The pair (Ω * , σ * ) corresponds to a body we want to hide from observations. Let us define in T 0 the conductivity σ 0 ∈ P C + (T 0 ) so that σ 0 coincides with σ * in Ω * and is an arbitrary positive definite piecewise constant conductivity in T 0 \ Ω * . Let T be a triangulation of Ω that is the union of T C \ {T 0 } and T 0 . Finally, let σ ∈ P C( T ) by the conductivity that coincides with σ 0 in T 0 and with the cloaking conductivity σ C in Ω \ T 0 . The FEM model ( T , σ) corresponds to a resistor network obtained by starting from a resistor network of type C, with some positive resistivities on the edges of the central triangle T 0 , and adding some additional vertexes and edges inside the triangle T 0 . In the resistor network corresponding the FEM model ( T , σ) the solution of the equation (3.1) is constant in all vertexes that are inside or on the boundary of the triangle T 0 . Thus the exactly same arguments that were used to prove Prop. 5.1 give the following result: Theorem 5.2. The conductivity σ ∈ P C( T ) is positive semidefinite and satisfies σ| Ω * = σ * . Moreover, for the FEM models ( T , σ) and (T B , σ B ) all boundary measurements coincide in the sense that (5.8) Summarizing the above construction, we can construct a piecewise constant conductivity σ ∈ P C( T ) for which the Dirichlet-to-Neumann maps for the FEM models ( T , σ) and (T B , σ B ) coincide and ( T , σ) contains the body (Ω * , σ * ). This example can be considered as discrete invisibility cloaking of the body (Ω * , σ * ). Note that in this example the "cloaking layer" around the body (Ω * , σ * ) is quite thick, but analogously to the PDE cloaking examples, one can consider a piecewise affine diffeomorphism Φ : Ω → Ω 1 to some suitable domain Ω 1 , such that Φ is identity in Ω * . Using such diffeomorphisms one can obtain other cloaking examples where the domain Ω 1 is a smaller neighborhood of Ω * . However, we do not consider any such particular examples in this paper.
Note that by Lemma 4.3 the conductivity matrices in σ on some triangles will be only positive semi-definite but not strictly positive definite. Analogues to the continuous PDE-based cloaking in [30,49] we will remedy this by constructing an approximate cloaking resistor network which will allow us to obtain a positive definite equivalent FEM model. Indeed, by setting in the above constructions all the dashed lines in graphs of type C to have conductivity > 0, we obtain a positive definite conductivity σ( ) ∈ P C + ( T ). Hence, approximate cloaking for FEM models is possible with positive definite conductivities.