Discrete Total Variation of the Normal Vector Field as Shape Prior with Applications in Geometric Inverse Problems

An analogue of the total variation prior for the normal vector field along the boundary of piecewise flat shapes in 3D is introduced. A major class of examples are triangulated surfaces as they occur for instance in finite element computations. The analysis of the functional is based on a differential geometric setting in which the unit normal vector is viewed as an element of the two-dimensional sphere manifold. It is found to agree with the discrete total mean curvature known in discrete differential geometry. A split Bregman iteration is proposed for the solution of discretized shape optimization problems, in which the total variation of the normal appears as a regularizer. Unlike most other priors, such as surface area, the new functional allows for piecewise flat shapes. As two applications, a mesh denoising and a geometric inverse problem of inclusion detection type involving a partial differential equation are considered. Numerical experiments confirm that polyhedral shapes can be identified quite accurately.


Introduction
The total variation (TV) functional is popular as a regularizer in imaging and inverse problems; see for instance Rudin, Osher, Fatemi, 1992;Chan, Golub, Mulet, 1999;Bachmayr, Burger, 2009;Langer, 2017 andVogel, 2002, Chapter 8. It is most commonly applied to functions with values in R or R n . In the companion paper Bergmann, Herrmann, et al., 2019, we introduced the total variation of the normal vector field n along smooth surfaces Γ ⊂ R 3 : |n| T V (Γ) := Γ |(D Γ n) ξ 1 | 2 g + |(D Γ n) ξ 2 | 2 g 1/2 ds. (1.1) In contrast to the setting of real-or vector-valued functions, the normal vector field is manifold-valued with values in the sphere S 2 = {v ∈ R 3 : |v| 2 = 1}. In (1.1), D Γ n denotes the derivative (push-forward) of n, and {ξ 1 (s), ξ 2 (s)} is an arbitrary orthonormal basis (w.r.t. the Euclidean inner product in the embedding Γ ⊂ R 3 ) of the tangent spaces T s Γ along Γ. Finally, | · | g denotes the norm induced by a Riemannian metric on S 2 . It was shown in Bergmann, Herrmann, et al., 2019 that (1.1) can be alternatively expressed as where k 1 and k 2 are the principal curvatures of the surface. In this paper, we discuss a discrete variant of (1.1) tailored to piecewise flat surfaces Γ h , where (1.1) does not apply. In contrast with the smooth setting, the total variation of the piecewise constant normal vector field n is concentrated in jumps across edges between flat facets. We therefore propose the following discrete total variation of the normal, ( 1.2) Here E denotes an edge of length |E| between facets, and d(n + E , n − E ) is the geodesic distance between the two neighboring normal vectors.
We investigate (1.2) in Section 2. It turns out to coincide with the discrete total mean curvature known in discrete differential geometry. Subsequently, we discuss the utility of this functional as a prior in shape optimization problems cast in the form Minimize (u(Ω h ), w.r.t. the vertex positions of the discrete shape Ω h with boundary Γ h . (1.3) Here u(Ω h ) denotes the solution of the problem specific partial differential equation (PDE), which depends on the unknown domain Ω h . Moreover, represents a loss function, such as a least squares function. In particular, (1.3) includes geometric inverse problems, where one seeks to recover a shape Ω h ⊂ R 3 representing, e.g., the location of a source or inclusion inside a given, larger domain, or the geometry of an inclusion or a scatterer. Numerical experiments confirm that |n| DT V (Γ h ) , as a shape prior, can help to identify polyhedral shapes.
Similarly as for the case of smooth surfaces discussed in Bergmann, Herrmann, et al., 2019, solving discrete shape optimization problems (1.3) is challenging due to the non-trivial dependency of n on the vertex positions of the discrete surface Γ h , as well as the non-smoothness of |n| DT V (Γ h ) . We therefore propose in Section 3 a version of the split Bregman method proposed in Goldstein, Osher, 2009, an algorithm from the alternating direction method of multipliers (ADMM) class in which the jumps in the normal vector are treated as a separate variable. The particularity here is that the normal vector has values in S 2 and thus the jump, termed d, is represented by a logarithmic map in the appropriate tangent space. An outstanding feature of the proposed splitting is that the two subproblems, the minimization w.r.t. the vertex coordinates representing the discrete surface and w.r.t. d, are directly amenable to numerical algorithms.
Although many optimization algorithms have been recently generalized to Riemannian manifolds, see, e.g., Bačák, 2014;, the Riemannian split Bregman method for manifolds proposed in this and the companion paper Bergmann, Herrmann, et al., 2019 is new to the best of our knowledge. Its detailed investigation will be postponed to future work. For a general overview of optimization on manifolds, we refer the reader to Absil, Mahony, Sepulchre, 2008. We anticipate that our method can be applied to other non-smooth problems involving manifold-valued total variation functionals as well. Examples falling into this class have been introduced for instance in Lellmann et al., 2013;Bergmann, Tenbrinck, 2018. An alternative splitting scheme, the so-called half-quadratic minimization, was introduced by Bergmann, Chan, et al., 2016.
The structure of the paper is as follows. In the following section we provide an analysis of the discrete total variation of the normal (1.2) and its properties. We also compare it to geometric functionals appearing elsewhere in the literature. In particular, we provide a numerical comparison between (1.2) and surface regularization for a mesh denoising problem. Section 3 is devoted to the formulation of an ADMM method which generalizes the split Bregman algorithm to the manifoldvalued problem (1.3). In section 4, we describe an inclusion detection problem of type (1.3), motivated by geophysical applications. We also provide implementation details in the finite element framework FEniCS. Corresponding numerical results are presented in Section 5.

Discrete Total Variation of the Normal
From this section onwards we assume that Γ h ⊂ R 3 is a piecewise flat, compact, orientable surface without boundary, which consists of a finite number of flat facets with straight sided edges between facets. Consequently, Γ h can be thought of as a mesh consisting of polyhedral cells with a consistently oriented outer unit normal. We also assume this mesh to be geometrically conforming, i.e., there are no hanging nodes. A frequent situation is that Γ h is the boundary mesh of a geometrically conforming volume mesh with polyhedral cells, representing a volume domain Ω h ⊂ R 3 . In our numerical example in Section 5, we will utilize a volume mesh consisting of tetrahedra, whose surface mesh consists of triangles; see  occurs in bands of width 2ε around the edges. Such an approximation can be constructed, for instance, by a level-set representation of Γ h by means of a signed distance function Φ. Then a family of smooth approximations Γ ε can be obtained as zero level sets of mollifications Φ ϕ ε for sufficiently small ε. Here ϕ ε is the standard Friedrichs mollifier in 3D and denotes convolution. A construction of this type is used, for instance, in Gómez, Hernández, López, 2005;Bonito, Demlow, Nochetto, 2019. An alternative to this procedure is the so-called Steiner smoothing, where Γ ε is taken to be the boundary of the Minkowski sum of Ω h with the ball B ε (0) ⊂ R 3 ; see for instance Sullivan, 2008, Section 4.4.
Proof. Let us denote the vertices in Γ h by V and its edges by E. Since mollification is local, the normal vector is constant in the interior of each facet minus its collar, which is of order ε. Consequently, changes in the normal vector are confined to a neighborhood of the skeleton. We decompose this area into the disjoint unioṅ Here I E,ε are the transition regions around edge E where the normal vector is modified due to mollification, and B V,ε are the regions around vertex V . On I E,ε , we can arrange the basis ξ 1,2 to be aligned and orthogonal to E so that holds, which can be easily evaluated as an iterated integral. In each stripe in I E,ε perpendicular to E, n ε changes monotonically along the geodesic path between n + E and n − E , so that the integral along this stripe yields the constant d(n + E , n − E ). Since the length of I E,ε parallel to E is |E| up to terms of order ε, we obtain The contributions to |n ε | T V (Γε) from integration over B V,ε are of order ε since |(D Γε n ε ) ξ 1 | 2 g + |(D Γε n ε ) ξ 2 | 2 g 1/2 is of order ε −1 and the area of B V,ε is of order ε 2 . This yields the claim.

Comparison with Prior
Work for Discrete Surfaces. The functional (2.1) has been used previously in the literature. We mention that it fits into the framework of total variation of manifold-valued functions defined in Giaquinta, Mucci, 2007;Lellmann et al., 2013. Specifically in the context of discrete surfaces, we mention Sullivan, 2005 where the term H E := |E| Θ E appears as the total mean curvature of the edge E. Here Θ E is the exterior dihedral angle, which agrees with d(n + E , n − E ), see (2.2). Consequently, (2.1) can be written as E H E . Moreover, (2.1) appears as a regularizer in Wu et al., 2015 within a variational model for mesh denoising but the geodesic distances are approximated for the purpose of numerical solution. We also mention the recent Pellis et al., 2019 where (2.1) appears as a measure of visual smoothness of discrete surfaces. Particular emphasis is given to the impact of the mesh connectivity. In our study, the mesh connectivity will remain fixed and only triangular surface meshes are considered in the numerical experiments.
In addition, we are aware of Zhang et al., 2015;Zhong et al., 2018 was proposed in the context of variational mesh denoising. Notice that in contrast to (2.1), (2.4) utilizes the Euclidean as opposed to the geodesic distance between neighboring normals and is therefore an underestimator for (2.1).
Once again, we are not aware of any work in which (2.1) or its continuous counterpart (1.1) were used as a prior in shape optimization or geometric inverse problems involving partial differential equations.
2.2. Properties of the Discrete Total Variation of the Normal. In this section we investigate some properties of the discrete total variation of the normal. As can be seen directly from (2.1), a scaling in which Γ h is replaced by δΓ h for some δ > 0 yields This is the same behavior observed, e.g., for the total variation of scalar functions defined on two-dimensional domains. Consequently, when studying optimization problems involving (2.1), we need to take precautions to avoid that Γ h degenerates to a point. This can be achived either by imposing a constraint, e.g., on the surface area, or by considering tracking problems in which an additional loss term appears.

Simple Minimizers of the Discrete Total Variation of the Normal.
In this section, we investigate minimizers of |n| DT V (Γ h ) subject to an area constraint. More precisely, we consider the following problem. Given a triangulated surface mesh consisting of vertices V , edges E and facets F , find the mesh with the same connectivity, which To the best of our knowledge, a precise characterization of the minimizers of (2.5) is an open problem and the solution depends on the connectivity; compare the observations in Pellis et al., 2019, Section 4. That is, different triangulations of the same (initial) mesh, e.g., a cube, may yield different minimizers. We also refer the reader to Alexa, Wardetzky, 2011 for a related observation in discrete mean curvature flow.
We do have, however, the following partial result. For the proof, we exploit that (2.1) coincides with the discrete total mean curvature and utilize results from discrete differential geometry. The reader may wish to consult Meyer et al., 2003;Polthier, 2005;Wardetzky, 2006;Bobenko, Springborn, 2007;Crane et al., 2013.
Theorem 2.2. The icosahedron and the cube with crossed diagonals are stationary for (2.5) within the class of triangulated surfaces Γ h of constant area and identical connectivity.
Proof. Let us consider the Lagrangian associated with (2.5), Here x i ∈ R 3 denote the coordinates of vertex #i and N V is the total number of vertices of the triangular surface mesh. Notice that the normal vectors n ± E , edge lengths |E| and facet areas |F | depend on these coordinates. The gradient of (2.6) w.r.t. x i can be represented as (2.7) see for instance Crane et al., 2013. Here N (i) denotes the index set of vertices adjacent to vertex #i. For any j ∈ N (i), E ij denotes the edge between vertices #i and #j. Moreover, α ij and β ij are the angles as illustrated in Figure 2.3.
For the icosahedron with surface area A 0 , all edges have length |E ij | = A0 Moreover, since all facets are unilaterial triangles, α ij = β ij = π/3 holds. Finally, the exterior dihedral angles d(n + Eij , n − Eij ) are all equal to arccos( We remark that (2.1) and thus (2.7) is not differentiable when one or more of the angles d(n + Eij , n − Eij ) are zero. This is the case for the cube with crossed diagonals, see Figure 2.3. However, the right hand side in (2.7) still provides a generalized derivative of L in the sense of Clarke. In contrast to the icosahedron, the cube has two types of vertices. When x i is the center vertex of one of the lateral surfaces, then d(n + Eij , n − Eij ) = 0 and α ij = x i , independently of the value of the Lagrange multiplier µ. Now when x i is a vertex of "corner type", we need to distinguish two types of edges. Along the three edges leading to neighbors of the same type, we have a exterior dihedral angle of d(n + Eij , n − Eij ) = π/2, length |E ij | = (A 0 /6) 1/2 and α ij = β ij = π/2. Along the three remaining edges leading to surface centers, we have d(n + Eij , n − Eij ) = 0 and α ij = β ij = π/4. Thus for vertices of "corner type", it is straightforward to verify that 0 belongs to the generalized (partial) differential of L at (x 1 , . . . , holds, which is true for the obvious choice of µ. Numerical experiments indicate that the icosahedron as well as the cube are not only stationary points, but also local minimizers of (2.5). We can thus conclude that the discrete objective (2.1) exhibits different minimizers than its continuous counterpart (1.1) for smooth surfaces. In particular, (2.1) admits and promotes piecewise flat minimizers such as the cube. This is in accordance with observations made in Pellis et al., 2019, Section 3.2 that optimal meshes typically exhibit a number of zero dihedral angles. This property sets our functional apart from other functionals previously used as priors in shape optimization and geometric inverse problems. For instance, the popular surface area prior is well known to produce smooth shapes; see the numerical experiments in section 2.2.3 below.

Comparison of Discrete and Continuous Total Variation of the
Normal. In this section we compare the values of (1.1) and (2.1) for a sphere Γ, and a sequence of discretized spheres Γ h . For comparison, we choose Γ to have the same surface area as the cube in the previous section, i.e., we use r = 3/(2π) as the radius. It is easy to see that since the principal curvatures of a sphere Γ of radius r are k 1 = k 2 = 1/r, (1.1) becomes which amounts to 4 √ 3π ≈ 12.2799 for the sphere under consideration.
To compare this to the discrete total variation of the normal, we created a sequence of triangular meshes Γ h of this sphere with various resolutions using Gmsh and evaluated (2.1) numerically. The results are shown in Table 2.2. They reveal a  (2π), their values of (2.1) and the ratio between (2.1) and (1.1). The value of the latter is |n| T V (Γ) = 4 √ 3π ≈ 12.2799. N V , N E and N T denote the number of vertices, edges, and triangles of the respective mesh.
factor of approximately √ 2 between the discrete and continuous functionals for the sphere. To explain this discrepancy, recall that the principal curvatures of the sphere are k 1 = k 2 = 1/r. This implies that the derivative map D Γ n has rank two everywhere. Discretized surfaces behave fundamentally different in the following respect. Their curvature is concentrated on the edges, and one of the principal curvatures (the one in the direction along the edge) is always zero. So even for successively refined meshes, e.g., of the sphere, one is still measuring only one principal curvature at a time. We are thus led to the conjecture that the limit of (2.1) for sucessively refined meshes is the "anisotropic", yet still intrinsic measure Γ |k 1 | + |k 2 | ds, whose value for the sphere in Table 2.1 is 4 √ 6π ≈ 17.3664. The factor √ 2 can thus be attributed to the ratio between the 1 -and 2 -norms of the vector (1, 1) . This observation is in accordance with the findings in Pellis et al., 2019, Section 1.2.
One could consider an "isotropic" version of (2.1) in which the dihedral angles across all edges meeting at any given vertex are measured jointly. These alternatives will be considered elsewhere.

Discrete Total Variation
Compared to Surface Area Regularization. In this section we consider a specific instance of the general problem (1.3) and compare our discrete TV functional with the surface area regularizer. We begin with a triangular surface mesh Γ h of a box Ω = (−1, 1) × (−1.5, 1.5) × (−2, 2) and add normally distributed noise to the coordinate vector of each vertex in average normal direction of the adjacent triangles with zero mean and standard deviation σ = 0.2 times the average edge length. We denote the noisy vertex positions as x V and utilize a simple least-squares functional as our loss function and consider the following mesh denoising problem, w.r.t. the vertex positions x V of the discrete surface Γ h . (2.8) Here the sum runs over the vertices of Γ h . For comparison, we also consider a variant Minimize where we utilize the total surface area as prior.
A numerical approach to solve the non-smooth problem (2.8) will be discussed in section 3. By contrast, problem (2.9) is a fairly standard smooth discrete shape optimization problem and we solve it using a simple shape gradient descent scheme.
The details how to obtain the shape derivative and shape gradient are the same as described in section 4.2 for problem (2.8). Figure 2.4 shows the numerical solutions of (2.8) and (2.9) for various choices of the regularization parameters β and γ, respectively. The initial guess for both problems is a sphere with the same connectivity as Γ h . We can clearly see that our functional (2.1) achieves a very good reconstruction of the original shape for a proper choice of β. By contrast, the surface area regularization requires a relatively large choice of γ in order to reasonably reduce the noise, which in turn leads to a significant shrinkage of the surface and a rounding of the sharp features. Top row: original box and desired outcome of the noise reduction, noisy box with vertex coordinates x V used for the data fidelity term, and sphere with same connectivity used as the initial guess; middle row: results for total variation of the normal (2.8) with β = 10 −2 , 10 −3 , 10 −4 ; bottom row: results for surface area regularization (2.9) with γ = 0.02, 0.01, 0.005.

Discrete Split Bregman Iteration
In this section, we develop an optimization scheme to solve the non-smooth problem (1.2). To this end, we adapt the well-known split Bregman method to our setting. This leads to a discrete realization of the approach presented in Bergmann, Herrmann, et al., 2019, section 4. Recall that combining (1.2) with (1.3) results in the problem where E are the edges of the unknown part Γ h of the boundary ∂Ω h . We will consider a concrete example in Section 4.1.
Notice that the second term in the objective in (3.1) is non-differentiable whenever n + E = n − E occurs on at least one edge. Following the classical split Bregman approach, we introduce a splitting in which the variation of the normal vector becomes an independent variable. Since this variation is confined to edges, where the normal vector jumps (without loss of generality) from n + E to n − E , this new variable becomes Here log n + E n − E denotes the logarithmic map, which specifies the unique tangent vector at the point n + E such that the geodesic departing from n + E in that direction will reach n − E at unit time. The logarithmic map is well-defined whenever n + Together with the set of Lagrange multipliers b E ∈ T n + E S 2 , we define the Augmented Lagrangian pertaining to (3.1) and (3.2) as The vectors d and b are simply the collections of their entries d E , b E ∈ T n + E S 2 , three components per edge E. Hence, since the tangent space T n + E S 2 changes between shape updates, the respective quantities have to be parallely transported, which is a major difference to ADMM methods in Euclidean or Hilbert spaces.
We state the split Bregman iteration in Algorithm 3.1. Perform several gradient steps for

5:
Parallely transport the multiplier estimate b Set k := k + 1 9: end while We now address the individual steps of algorithm 3.1 in more detail, i.e., the successive minimization with respect to the unknown vertices of Ω h and d, followed by an explicit update for the multiplier b.
Step 4 is the minimization of (3.3) with respect to the unknown vertex positions of Ω h . To this end, we employ a gradient descent scheme, where we compute the sensitivities with respect to those node positions discretely, see Section 4.2 for more details. Following Goldstein, Osher, 2009, an approximate minimization suffices, and thus only a certain number of steepest descent steps are performed. After Ω (k) h has been updated to Ω for each edge E.

An EIT Model Problem and its Implementation in FEn-iCS
In this section we address some details concerning the implementation of Algorithm 3.1 in the finite element framework FEniCS (version 2018.2.dev0), Logg, Mardal, Wells, et al., 2012;Alnaes, Blechta, et al., 2015. For concreteness, we elaborate on a particular reduced loss function (u(Ω), Ω) where the state u(Ω) arises from a PDE modeling a geological electrical impedance tomography (EIT) problem with Robin-type far field boundary conditions. We introduce the problem under consideration first and discuss implementation details and derivative computations later on.

EIT Model Problem.
Electrical impedance tomography (EIT) problems are a prototypical class of inverse problems. Common to these problems is the task of reconstructing the internal conductivity inside a volume from boundary measurements of electric potentials or currents. These problems are both nonlinear and severely ill-posed and require appropriate regularization; see for instance Santosa, Vogelius, 1990;Cheney, Isaacson, Newell, 1999;Chung, Chan, Tai, 2005. Traditionally, EIT problems are modeled with Neumann (current) boundary conditions and the internal conductivity is an unknown function across the entire domain. In order to focus on the demonstration of the utility of (1.2) as a regularizer in geometric inverse problems, we consider a simplified situation in which we seek to Consequently, the unknowns are the vertex positions of the interface of the inclusion. As a perfect conductor shields its interior from the electric field, there is no necessity to mesh and simulate the interior of the inclusion. However, we mention that our methodology can be extended also to interface problems, non-perfect conductors and other geometric inverse problems.
The perfect conductor is modeled via a homogenous Neumann condition on the unknown interior boundary Γ 1 of the domain Ω. To overcome the non-uniqueness of the electric potential, we employ Robin boundary conditions on the exterior boundary Γ 2 . The use of homogeneous Robin boundary conditions to model the far field is well-established for geological EIT problems; see, e.g., Helfrich-Schkarbanenko, 2011. We use them here also for current injection.
The geometry of our model is shown in Figure 4.1, where Γ 1 is the unknown boundary of the perfect conductor and Γ 2 is a fixed boundary where currents are injected and measurements are taken. We assume that i = 1, . . . , r ∈ N experiments are conducted, each resulting in a measured electric potential z i ∈ CG 1 (Γ 2 ), the finite element space consisting of piecewise linear, globally continuous functions on the outer boundary Γ 2 . Experiment #i is conducted by applying the right hand side source f i ∈ DG 0 (Γ 2 ), which is the characteristic function of one of the colored regions shown in Figure 4.1. Here, DG 0 denotes the space of piecewise constant functions. We then seek to reconstruct the interface of the inclusion Γ 1 by solving the following regularized least-squares problem of type (1.3), with respect to the vertex positions of Γ 1 . Here u i ∈ CG 1 (Ω h ) is the computed electric field for source f i . Hence, the problem features r PDE constraints with identical operator but different right hand sides.
As detailed in Section 4.2, we compute the shape derivative of the least-squares objective and the PDE constraint separately from the shape derivative of the regularization term. To evaluate the former, we utilize a classical adjoint approach. To this end, we consider the Lagrangian for p i ∈ CG 1 (Ω h ). The differentiation w.r.t. u i leads to the following adjoint problem for p i : (4.3) The above adjoint PDE was implemented by hand. Since all forward and adjoint problems are governed by the same differential operator, we assemble the associated stiffness matrix once and solve the state and adjoint equations via an ILU-preconditioned conjugate gradient method.
Provided that u i and p i solve the respective state and adjoint equations, the directional derivative of (u(Ω h ), Ω h ) coincides with the partial directional derivative of F (u 1 , . . . , u r , p 1 , . . . , p r , Ω h ), both with respect to the vertex positions. In practice, we evaluate the latter using the coordinate derivative functionality of FEniCS as described in the following subsection.

Discrete Shape Derivative.
We now focus on computing the sensitivity of finite element functionals, when mesh vertices x of Ω h are moved in accordance to x ε = x + εV Ω h with V Ω h ∈ CG 3 1 (Ω h ). As discussed in Ham et al., 2018, a convenient way to compute this within the finite element world is by tapping into the transformation of the reference element to the physical one. Hence, we use the symbol d (u(Ω h ), Ω h )[V Ω h ] for this object and we obtain it using the coordinate derivative functionality, first introduced in FEniCS release 2018.2.dev0.
Our split Bregman scheme requires the shape derivative of (3.3), which is given by where originates from the splitting approach (3.3). Because our design variable is Γ 1 only, we introduce the extension P Ω h (V Γ1 ) of V Γ1 ∈ CG 3 1 (Γ 1 ) to the volume Ω h by padding with zeros. Furthermore, a reduction to boundary only sensitivities can also be motivated from considering shape derivatives in the continuous setting, see Bergmann, Herrmann, et al., 2019, Section 3.
is computed via the adjoint approach as explained above, In order to employ this AD functionality, (4.5) needs to be given as a UFL form, a domain specific language based on Python, which forms the native language of the FEniCS framework, see Alnaes, Logg, et al., 2014. Such a UFL representation is easy to achieve if all mathematical expressions are finite element functions. Notice that d and b in (4.5) are constant functions on the edges of the boundary mesh representing Γ 1 . We can thus represent them in the so called HDivTrace space of lowest order in FEniCS. From the directional derivatives (4.4), we pass to a shape gradient on the surface w.r.t. a scaled H 1 (Γ 1 ) scalar product by solving a variational problem. This problem involves the weak form of a Laplace-Beltrami operator with potential term and it finds W Γ1 ∈ CG 1 (Γ 1 ) 3 such that holds for all test functions V Γ1 ∈ CG 1 (Γ 1 ) 3 . The previous procedure provides us with a shape gradient W Γ1 on the surface Γ 1 alone. In order to propagate this information into the volume Ω h , we solve the following mesh deformation equation: find W Ω h ∈ CG 1 (Ω h ) 3 such that for all test functions V Ω h ∈ CG 1 (Ω h ) 3 with zero Dirichlet boundary conditions, where W Ω h is subject to the Dirichlet boundary condition W Ω h = W Γ1 on Γ 1 and W Ω h = 0 on Γ 2 . Subsequently, the vertices of the mesh are moved in the direction of W Ω h .

Intrinsic Formulation Using Co-Normal Vectors.
We recall that our functional of interest (1.2) is formulated in terms of the unit outer normal n of the oriented surface Γ 1 . This leads to the term (4.5) inside the augmented Lagrangian (3.3). In order to utilize the differentiation capability of FEniCS w.r.t. vertex coordinates, we need to represent (4.5) in terms of an integral. Since the edges are the interior facets of the surface mesh for Γ 1 , and d and b can be represented as constant on edges as explained above, (4.5) can indeed be written as an integral w.r.t. the interior facet measure dS on Γ 1 . Then, however, the outer normal vectors appearing in the term log n + domain boundary of Ω boundary of true inclusion Γ 1 boundary of [−0.4, 0.4] 3 initial guess for Γ 1 boundary of B 0.5 (0) ⊂ R 3 number of measurements r = 48 Robin coefficient in (4.1) α = 10 −5 split Bregman parameter λ = 10 −5 standard deviation of noise σ = 0 or σ = 0.34 · 10 −2 regularization parameter . . . for total variation regularization β = 10 −6 for surface area regularization γ = 5 · 10 −5 , 2 · 10 −5 shape step size 10 2

Numerical Results
In this section we present numerical results obtained with Algorithm 3.1 for the geological impedance tomography model problem described in the previous section. The data of the problem are given in Table 5.1 and the initial guess of the inclusion Γ 1 , as well as the true inclusion, are shown in Figure 5.1. The state u and adjoint state p were discretized using piecewise linear, globally continuous finite elements on a tetrahedral grid of Ω minus the volume enclosed by Γ 1 . The mesh has 4429 vertices and 19 384 tetrahedra. Regarding the shape optimization problem of Algorithm 3.1, we perform 10 gradient steps per split Bregman iteration combined with an Armijo linesearch with starting step size of 10 2 . Also, we stop the whole algorithm, i.e., the outer Bregman iteration, when the initial gradient of the above mentioned shape optimization problem has a norm below 10 −7 in the sense of (4.6).
In Figure 5.2, we show the results obtained in the noise-free setting (top row) and with noise (bottom row). In the latter case, normally distributed random noise is added with zero mean and standard deviation σ = 0.34 · 10 −2 per degree of freedom of z i on Γ 2 for each of the r = 48 simulations of the forward model (4.1). The amount of noise is considerable when put in relation to the average range of values for the simulated states, which is Due to mesh corruption, we have to remesh Ω h at some point in the cases with noise. Afterwards, we start again Algorithm 3.1 with the remeshed Ω h as new initial guess.
For comparison, we also provide results obtained for a related problem in Figure 5.2, using the popular surface area regularization with the same data otherwise. For the surface area regularization, β |n| T V (Γ1) is replaced by γ Γ1 ds = γ F |F |, where F are the facets of Γ 1 . Because the problem is smooth in this case, we apply a shape gradient scheme directly rather than a split Bregman scheme and terminate as soon as the norm of the gradient falls below 5 · 10 −8 . The regularization parameters β and γ are selected by hand in each case. Automatic parameter selection strategies can clearly be applied here as well, but this is out of the scope of the present paper. As is expected and well known, the use of surface area regularization leads to results in which the identified inclusion Γ 1 is smoothed out. This can be explained by the observation that the gradient based minimization of the surface area yields a mean curvature flow. By contrast, our novel prior (2.1) allows for piecewise flat shapes and thus the interface Γ 1 is closely reconstructed in the noise-free situation. Even in the presence of noise, the reconstruction is remarkably good. In particular, the flat lateral surfaces and sharp edges can be identified quite well.

Conclusions
In this paper we introduced a discrete analogue of the total variation prior for the normal vector field as shown in Bergmann, Herrmann, et al., 2019. While we are currently unable to characterize all minimizers of its discrete counterpart, we showed that the icosahedron and a cube with crossed diagonals are stationary under an area constraint. We conjecture that the full set of minimizers is much richer than this, in particular when the connectivity is included as design unknown. It has been Figure 5.2. Top row: setting without noise; left: total variation regularization, β = 10 −6 and 90 iterations; middle: surface area regularization with γ = 5 · 10 −5 and 1129 iterations; right: surface area regularization with γ = 2 · 10 −5 and 978 iterations. Bottom row: setting with noise; left: total variation regularization with β = 10 −6 and 173 iterations with remeshing after iteration 121; middle: surface area regularization with γ = 5 · 10 −5 and 1016 iterations with remeshing after iteration 539; right: surface area regularization with γ = 2 · 10 −5 and 987 iterations with remeshing after iteration 308.
argued in Pellis et al., 2019, Section 3.3 that minimal energy is achieved for meshes which are not triangular, but whose faces are approximately rectangular.
We proposed, described and implemented a split Bregman (ADMM) scheme for the numerical solution of shape optimization problems involving the discrete total variation of the normal. In contrast to a Euclidean ADMM as proposed for instance in Goldstein, Osher, 2009, the normal vector data belongs to the sphere S 2 . Therefore, the formulation of the ADMM method requires concepts from differential geometry. In particular, the discrete setting utilizes logarithmic maps and parallel transport of tangent vectors. An analysis of the ADMM scheme is beyond the scope of this paper and will be presented elsewhere.
We demonstrate the utility of the discrete total variation of the normal as a shape prior in a geometric inverse problem, in which we aim to detect a polyhedral inclusion. Unlike the popular surface area regularization, our prior allows for piecewise flat shapes.

A The Sphere as a Riemannian Manifold
In this section we provide some useful formulas for the sphere S 2 = {n ∈ R 3 : |n| 2 = 1} equipped with the Riemannian metric obtained from the pull back of the Euclidean metric from the ambient space R 3 . We are going to represent points n ∈ S 2 by vectors in R 3 . Moreover, we identify the tangent space at n with the twodimensional subspace T n S 2 = {ξ ∈ R 3 : ξ n = 0}. We utilize the Riemannian metric g(a, b) = a b in T n S 2 and the norm |a| g = (a a) 1/2 . The geodesic distance between any two n, n ∈ S 2 is given by d(n, n ) = arccos(n n ). (A.1) The geodesic curve γ( · ; n, ξ) : R → S 2 departing from n ∈ S 2 in the direction of ξ ∈ T n S 2 is given by γ(t; n, ξ) = cos t |ξ| g n + sin t |ξ| g ξ |ξ| g .
The exponential map is thus given by exp n ξ = γ(1; n, ξ) = cos |ξ| g n + sin |ξ| g ξ |ξ| g . (A. 3) The logarithmic map is the inverse of the exponential map w.r.t. to the tangent direction ξ. In other words, ξ = log n n holds if any only if ξ is the unique element in T n S 2 such that exp n ξ = n holds. The logarithmic map is well-defined whenever n = −n holds. In this case, we have log n n = d(n, n ) n − (n n ) n |n − (n n ) n| g . (A.4) Finally we require the concept of parallel transport of a tangent vector from one tangent space to another, along the unique shortest geodesic connecting the base points. Specifically, the parallel transport P n→n : T n S 2 → T n S 2 along the unique shortest geodesic γ( · ; n, log n n ) connecting n and n = −n is given by P n→n (ξ) = ξ − ξ (log n n ) d 2 (n, n ) (log n n + log n n) = ξ + cos(|v| g ) u − u − sin(|v| g ) n u ξ, (A.5) see for instance Hosseini, Uschmajew, 2017 andPersch, 2018, Section 2.3.1, repectively. Here we used the abbreviations v = log n n , |v| g = d(n, n ) and u = v |v|g . To see that both expressions in (A.5) coincide -after plugging in the definition of the geodesic distance (A.1)-it remains to show that − n n log n n |log n n | g + 1 − (n n ) 2 n = log n n |log n n| g .
which holds true since the norm of the logarithmic map is |log n n | g = |n − n n n| g = (n n ) − (n n ) = 1 − (n n ) = |log n n| g .
Hence multiplying with the denominator of the first term in (A.5) yields the equality with the second term, since using the definition of the logarithmic map we obtain (n n ) n − (n n ) 2 n − (1 − (n n ) 2 ) n = n − (n n ) n .