Shape and topology optimization for variational inequalities with pointwise boundary observation

TL;DR

This paper addresses shape and topology optimization problems involving variational inequalities with unilateral conditions and boundary observations, employing regularization and penalization techniques within the Hamiltonian framework.

Contribution

It introduces a novel approach combining regularization and penalization methods for variational inequalities in shape and topology optimization.

Findings

Effective numerical examples demonstrate the method's viability.

The approach successfully handles unilateral boundary conditions.

The techniques improve the stability and accuracy of the optimization process.

Abstract

We study optimal design problems involving variational inequalities with unilateral conditions in the domain and pointwise boundary observation. We use regularizing and penalization tehniques in the setting of the Hamiltonian approach to shape/topology optimization problems. Numerical examples are also included.

Equations

$$\Omega ={\Omega }_g=\mathrm{int}\{\mathbf{x}\in D:g(\mathbf{x})\le 0\},g\in \mathcal{F}.$$

$$g+\lambda h,\lambda \in \mathbb{R},g,h\in \mathcal{F},$$

$${\mathbf{x}}_0\in \partial {\Omega }_g,\forall g\in \mathcal{F}.$$

$$\underset{g\in \mathcal{F}}{\min }J(g,y_g),$$

$$A{y}_g=f\mathrm{in}{\Omega }_g,$$

$$B{y}_g=0\mathrm{on}\partial {\Omega }_g.$$

$$\nabla g(\mathbf{x})\ne 0,\forall g\in \mathcal{F},\forall \mathbf{x}\in D,withg(\mathbf{x})=0.$$

$$g>0\mathrm{on}\partial D,\forall g\in \mathcal{F},$$

$$z_1^{\prime }(t)$$

$$z_2^{\prime }(t)$$

$$\left(z_1(0),z_2(0)\right)$$

$$-\Delta y_g+\beta (y_g)$$

$$y_g$$

$$0\in \beta (0).$$

$$-\Delta y_{\epsilon }+{\beta }_{\epsilon }(y_{\epsilon })+\frac{1}{\epsilon }{H}_{\epsilon }(g)y_{\epsilon }$$

$$y_{\epsilon }$$

$${\beta }_{\epsilon }(y)={\int }_{-\infty }^{\infty }{\beta }^{\epsilon }(y-\epsilon s)\rho (s)ds,$$

$$\underset{g\in \mathcal{F}}{\min }\left\{J_1={\left|\frac{\partial y_{\epsilon }}{\partial \mathbf{n}}({\mathbf{x}}_0)-\alpha \right|}^2\right\},$$

$$\underset{g\in \mathcal{F}}{\min }\left\{J_2=\sum _{i=0}^{l-1}{\left|\frac{\partial y_{\epsilon }}{\partial \mathbf{n}}\left(\mathbf{z}\left(\frac{i{T}_g}{l}\right)\right)-{\alpha }_i\right|}^2\right\},$$

$$-\Delta u+{\beta }_{\epsilon }^{\prime }(y_{\epsilon })u+\frac{1}{\epsilon }{H}_{\epsilon }(g)u+\frac{1}{\epsilon }{H}_{\epsilon }^{\prime }(g)y_{\epsilon }h$$

$$u$$

$$w_1^{\prime }$$

$$w_2^{\prime }$$

$$w_1(0)$$

$$\underset{\lambda \to 0}{\lim }\frac{T_{\lambda }-T_g}{\lambda }=-\frac{w_2(T_g)}{z_2^{\prime }(T_g)}$$

$$-\Delta (y_{\epsilon }(g+\lambda h)-y_{\epsilon }(g))+{\beta }_{\epsilon }(y_{\epsilon }(g+\lambda h))-{\beta }_{\epsilon }(y_{\epsilon }(g))$$

$$+\frac{1}{\epsilon }{H}_{\epsilon }(g+\lambda h)(y_{\epsilon }(g+\lambda h)-y_{\epsilon }(g))$$

$$-\frac{1}{\epsilon }\left[H_{\epsilon }(g+\lambda h)-H_{\epsilon }(g)\right]y_{\epsilon }(g)=0.$$

$$L$$

$$\underset{\lambda \to 0}{\lim }\frac{1}{\lambda }\left[{\left|\frac{\partial y_{\epsilon }(g+\lambda h)}{\partial {\mathbf{n}}_{\lambda }}({\mathbf{x}}_0)-\alpha \right|}^2-{\left|\frac{\partial y_{\epsilon }(g+\lambda h)}{\partial \mathbf{n}}({\mathbf{x}}_0)-\alpha \right|}^2\right].$$

$$\underset{\lambda \to 0}{\lim }\frac{1}{\lambda }\left[\frac{\partial y_{\epsilon }(g+\lambda h)}{\partial {\mathbf{n}}_{\lambda }}({\mathbf{x}}_0)-\frac{\partial y_{\epsilon }(g+\lambda h)}{\partial \mathbf{n}}({\mathbf{x}}_0)\right]$$

$$M=\left[-\frac{\partial y_{\epsilon }(g)}{\partial \mathbf{n}}\frac{\partial h}{\partial \mathbf{n}}\frac{1}{|\nabla g|}+\frac{\nabla y_{\epsilon }(g)\cdot \nabla h}{|\nabla g|}\right]({\mathbf{x}}_0).$$

$$\underset{\lambda \to 0}{\lim }\frac{1}{\lambda }\left[\frac{1}{|\nabla (g+\lambda h)|}-\frac{1}{|\nabla (g)|}\right]({\mathbf{x}}_0)=-\frac{\partial h}{\partial \mathbf{n}}({\mathbf{x}}_0)\frac{1}{|\nabla g({\mathbf{x}}_0){|}^2}.$$

$$L+2M\left(\frac{\partial y_{\epsilon }(g)}{\partial \mathbf{n}}({\mathbf{x}}_0)-\alpha \right)=$$

$$2\left(\frac{\partial y_{\epsilon }(g)}{\partial \mathbf{n}}({\mathbf{x}}_0)-\alpha \right)\left[\frac{\partial u}{\partial \mathbf{n}}-\frac{\partial y_{\epsilon }(g)}{\partial \mathbf{n}}\frac{\partial h}{\partial \mathbf{n}}\frac{1}{|\nabla g|}+\frac{\nabla y_{\epsilon }(g)\cdot \nabla h}{|\nabla g|}\right]({\mathbf{x}}_0).$$

$$2\left(\frac{\partial y_{\epsilon }(g)}{\partial \mathbf{n}}({\mathbf{x}}_0)-\alpha \right)\frac{\partial u}{\partial \mathbf{n}}({\mathbf{x}}_0)=<p,\mu >=<p,\frac{1}{\epsilon }{H}_{\epsilon }^{\prime }(g)y_{\epsilon }h>,$$

Full text

Keywords: optimal design, pointwise observation, nonsmooth variational inequalities.

MSC: 35J87, 49M41

1 Introduction

We consider a fixed domain approach, based on level functions for the (unknown) geometry and on the penalization/regularization of the state system. The shape and topology optimization problem is associated to nonlinear second order elliptic equations with Dirichlet boundary conditions. The optimal design problem can be approximated, via penalization, by an elliptic optimal control problem (see [12]).

Our methodology is underpinned by functional variations [13], [10]. For linear and nonlinear state systems, it is known that both closing and creating new holes is possible. Moreover, we prove differentiability and we apply standard gradient algorithms in the optimization process, that are known to converge locally, [3]. For a detailed discussion of optimal control problems governed by nonlinear and nonsmooth state systems, see [2], [14].

In the setting of functional variations, one assumes that all the admissible domains $\Omega$ are contained in a given bounded domain $D\subset\mathbb{R}^{2}$. We denote by $\mathcal{O}$ their family. We assume that this is generated by a family $\mathcal{F}\subset\mathcal{C}(\overline{D})$ of admissible level functions, as open sublevel sets:

[TABLE]

We recall that functional variations have the form

[TABLE]

which enables the work in functional spaces instead of the very specialized geometric variations.

To generate domains, we add to (1) a condition that selects just the component satisfying ($\mathbf{x}_{0}\in D$ is fixed)

[TABLE]

Notice that, for any admissible $g\in\mathcal{F}$, we have $g(\mathbf{x}_{0})=0$.

A general shape and topology optimization problem has the form (for $\Omega_{g}\in\mathcal{O}$ defined as above):

[TABLE]

In (4) - (6), $A$ is an elliptic partial differential operator, $B$ is a boundary operator, defining the state system and $J$ is some cost functional. The operators may be linear or nonlinear and appropriate assumptions will be imposed as necessity appears.

Assume now that $\mathcal{F}\subset\mathcal{C}^{1}(\overline{D})$ and

[TABLE]

Assume, as well, that

[TABLE]

which ensures that $\partial\Omega_{g}\cap\partial D=\emptyset$. If, moreover, the open cone $\mathcal{F}\subset\mathcal{C}^{2}(D)$ is given by (7), (8), then the solution of the Hamiltonian system below is periodic and gives a global parameterization of $\partial\Omega_{g}$

[TABLE]

by [10], [6]. We may take $I_{g}=[0,T_{g}]$, where $T_{g}>0$ denotes the principal period in (9)-(11). Moreover, the dependence $g\in\mathcal{F}\to T_{g}$ has differentiability properties, due to [10] if the family of shape functions satisfies $\mathcal{F}\subset\mathcal{C}^{2}(D)$ and 7), (8).

In every $\Omega\in\mathcal{O}$, $\Omega=\Omega_{g}\subset D$, $g\in\mathcal{F}$, we consider the variational inequality:

[TABLE]

where $f\in L^{2}(D)$ is given and $\beta\subset\mathbb{R}\times\mathbb{R}$ is a maximal monotone graph satisfying

[TABLE]

The nonlinear boundary value problem (12), (13) includes as an example the obstacle problem with obstacle $\varphi=0$. In case the obstacle $\varphi:D\to\mathbb{R}$ is some other function (negative in $D$, to ensure compatibility with (13) for any $\Omega\in\mathcal{O}$), then $\beta(y_{g}-\varphi)$ should be written in (12).

It is known, [4], that (12), (13) has a unique solution $y_{g}\in H^{2}(\Omega_{g})\cap H_{0}^{1}(\Omega_{g})$ with $\beta(y_{g})\in L^{2}(\Omega_{g})$, if $\partial\Omega_{g}$ is in $\mathcal{C}^{1,1}$. Therefore the formulation (4), (12), (13) makes sense. The regularity of $\partial\Omega_{g}$ may be obtained from the implicit functions theorem since $\mathcal{F}\subset\mathcal{C}^{2}(D)$ and (7) is satisfied (we get even $\mathcal{C}^{2}$ regularity for $\partial\Omega_{g}$).

We penalize and regularize (12) as follows ($\varepsilon>0$ is “small”):

[TABLE]

Above, $\beta_{\varepsilon}$ combines the Yosida approximation $\beta^{\varepsilon}$ with Friedrichs mollifiers:

[TABLE]

with $\rho\in\mathcal{C}_{0}^{\infty}(R)$, $\mathrm{supp}(\rho)\subset[-1,1]$, $\rho(-s)=\rho(s)$, $\rho(s)\geq 0$ and $\int_{-\infty}^{\infty}\rho(s)ds=1$.

Concerning $H(\cdot)$, it denotes the Heaviside function on $\mathbb{R}$ and its maximal monotone extension in $\mathbb{R}\times\mathbb{R}$ as well (i.e., we take $H(0)=[0,1]$). The notation $H_{\varepsilon}$ is obtained from $H\subset\mathbb{R}\times\mathbb{R}$ again via Yosida approximation and Friedrichs mollifiers, as in (17). Clearly, (15), (16) has a unique solution $y_{\varepsilon}\in H^{2}(D)\cap H_{0}^{1}(D)$ if $\partial D$ is supposed to be of class $\mathcal{C}^{1,1}$.

This distributed penalization technique, using $H,H_{\varepsilon}$, is due to Natori and Kawarada [7], in the setting of free boundary problems, and employed in shape optimization, for instance, in [13], for linear elliptic equations.

We indicate an approximation result, following [11], under weak regularity assumptions (we do not require $g\in\mathcal{F}$, for the moment) on the domain $\Omega$ (class $\mathcal{C}$, i.e., continuous boundary or the segment property, see [1], Def. 4.2, p. 66, [15]). We recall here that $\Omega$ is of class $\mathcal{C},$ if there exists a locally finite open cover of $\partial\Omega$ and an associated sequence of vectors $\{y_{j}\}$, not null, such that for any $x\in\partial\Omega$ there is $j$ with $x+ty_{j}\in\Omega$ for $0<t<1$ and this remains valid on the associated neighborhood of $x$, in the open cover of $\partial\Omega$.

Theorem 1**.**

If $\Omega=\Omega_{g}$ is of class $\mathcal{C}$, then, on a subsequence, $y_{\varepsilon}|_{\Omega_{g}}\to y_{g}$ weakly in $H^{1}(\Omega_{g})$ and $\Delta y_{\varepsilon}\to\Delta y_{g}$, $\beta_{\varepsilon}(y_{\varepsilon})\to\beta(y_{g})$ weakly in $L^{2}(\Omega_{g})$. If $g\in\mathcal{F}$, then $\frac{\partial y_{\varepsilon}}{\partial\mathbf{n}}\to\frac{\partial y_{g}}{\partial\mathbf{n}}$ in $H^{-1/2}(\partial\Omega_{g})$ weakly.

In the next section we discuss a special choice of the cost functional (4), namely in the so called pointwise observation case, that describes realistic, but computationally difficult situations.

2 Pointwise observation

In this section, an important role in the justification of the following computations and arguments is played by the topological stability result from [18], for small functional variational perturbations. Somewhat paradoxical, this stability result allows topological changes in optimal design since the corresponding steps in the descent gradient algorithms are not small, in general.

A natural and maybe the easiest case of pointwise boundary observation (taking into acount (9)-(11)) is the following performance index:

[TABLE]

where $\alpha\in\mathbb{R}$ is given, $y_{\varepsilon}$ is the solution of (15), (16) and $\mathbf{x}_{0}\in\partial\Omega$ is the initial condition in (9)-(11). Other points may be taken into account as well, see (19). The normal derivative refers to the domain $\Omega_{g}$ and can be expressed as $\frac{\nabla y_{\varepsilon}\cdot\nabla g}{|\nabla g|}$ in any point of $D$. Notice that (18) makes sense due to

Remark 1**.**

If $\mathcal{F}\subset\mathcal{C}^{2}$ and $f\in H^{1}(D)$, then $y_{\varepsilon}\in H^{3}(D)$ since $\beta_{\varepsilon}(y_{\varepsilon})$ and $H_{\varepsilon}(y_{\varepsilon})$ are in $H^{2}(D)$ due to the regularity of $\beta_{\varepsilon},H_{\varepsilon}$ and we know that $y_{\varepsilon}\in H^{2}(D)$. Consequently, $\Delta y_{\varepsilon}\in H^{1}(D)$ and we get $y_{\varepsilon}\in H^{3}(D)$, by standard regularity results for the Laplace equation and $y_{\varepsilon}\in\mathcal{C}^{1}(\overline{D})$ by the Sobolev theorem. That is $\frac{\partial y_{\varepsilon}}{\partial\mathbf{n}}$ has a pointwise sense.

We have denoted by $T_{g}$ the main period of the solution of (9)-(11)) corresponding to $g\in\mathcal{F}$. Let $l$ be a given natural number and divide $[0,T_{g}]$ (and implicitly $\partial\Omega_{g}$) in $l$ parts starting from the initial condition (11) (that we may arbitrarily choose on $\partial\Omega_{g}$, in this case). The described process is like a discretization procedure. We define in this way a preassigned number of observation points on $\partial\Omega_{g}$ and another example of pointwise cost functional is

[TABLE]

where $\alpha_{i}\in\mathbb{R}$ are given. A variant of (19) is to choose the observation points just on some part of $\partial\Omega_{g}$ , by fixing the indices $l_{1}\leq i\leq l_{2}$, with $l_{1},l_{2}$ given natural numbers, in the corresponding sum (or by using another appropriate subset of indices).

Of course, one can define pointwise observation functionals, by using points in $\Omega$, not on $\partial\Omega$ (and under the additional condition that $g\in\mathcal{F}$ are negative in such points). But this type of observation is less realistic, in general (it is difficult to place sensors inside a body) and we shall consider just the minimization of tracking type functionals as in (18), (19).

In the sequel, we study the differentiability properties of (18), (for (19) the arguments are quite similar). We recall first the following two results from [9], [16], [10], [11]. We consider the perturbations of the parameters (controls) $g+\lambda h$ and the system below (20)-(23) characterizes the variations of the corresponding solutions of the state system (15), (16), (9)-(11). Furthermore, $\mathbf{z}_{\lambda}$ denotes the solution of (9)-(11) associated to $g+\lambda h$, $T_{\lambda}$ denotes the main period of $\mathbf{z}_{\lambda}$ and we shall also write $y_{\varepsilon}(g),y_{\varepsilon}(g+\lambda h)$ when we emphasize the dependence on the level functions describing the geometry.

Proposition 1**.**

*The system in variations corresponding to (15), (16),

(9)-(11) is:*

[TABLE]

[TABLE]

where $u=\lim_{\lambda\rightarrow 0}\frac{y_{\varepsilon}(g+\lambda h)-y_{\varepsilon}(g)}{\lambda}$, $\mathbf{w}=[w_{1},w_{2}]=\lim_{\lambda\rightarrow 0}\frac{\mathbf{z}_{\lambda}-\mathbf{z}}{\lambda}$ and the limits exist in $H^{2}(D)$, respectively $\mathcal{C}^{1}([0,T_{g}])^{2}$.

Proposition 2**.**

Under the assumptions (7), (8), if $T_{\lambda}$ denotes the main period of the solution of (9)-(11), we have:

[TABLE]

if $z_{2}^{\prime}(T_{g})\neq 0$.

Remark 2**.**

If $z_{1}^{\prime}(T_{g})\neq 0$, the limit is $-\frac{w_{1}(T_{g})}{z_{1}^{\prime}(T_{g})}$. In general, we denote by $\theta(g,r)$ this limit. Either the last condition in Proposition 2 or this condition here follow by (7). In case both conditions are valid, it is known that the two mentioned results coincide (this is a consequence of Prop.1 in [17]).

The directional derivatives of cost functionals can be computed as follows.

We consider the cost (18) and we notice (under the notations of Prop. 1 and by Rem.1) that the computations to obtain (20) can be slightly refined as follows ($\varepsilon$ is fixed here):

[TABLE]

Dividing (24) by $\lambda\neq 0$ and using Prop.1, we get $\Delta\frac{y_{\varepsilon}(g+\lambda h)-y_{\varepsilon}(g)}{\lambda}$ bounded in $H^{1}(D)$ due to the regularity properties of $\beta_{\varepsilon},H_{\varepsilon}$, that is $\frac{y_{\varepsilon}(g+\lambda h)-y_{\varepsilon}(g)}{\lambda}\rightarrow u$ strongly in $\mathcal{C}^{1}(D)$, due to the Sobolev theorem. Then, we can compute the derivative of the pointwise cost (18). First, we notice that:

[TABLE]

In (25), we have decomposed the first paranthesis in the product of the sum by the difference of the squared terms and we apply Prop.1 and the above argument, to pass to the limit. We continue with the other term:

[TABLE]

where $n_{\lambda}$ is the normal in $\mathbf{x}_{0}$ to $\partial\Omega_{\lambda}$, $\Omega_{\lambda}$ being the domain associated to $g+\lambda h$ via (1), (3). After the decomposition of the difference of the squares, the important part to be investigated is:

[TABLE]

where we apply that $g=0$ on $\partial\Omega$ and $g+\lambda h=0$ on $\partial\Omega_{\lambda}$ and the well known formula for the normal derivative. By some calculations, we get the limit:

[TABLE]

For the first term above, we have used that:

[TABLE]

Combining (25) - (28) with $M$, we infer:

Proposition 3**.**

Assume that $\mathcal{F}\subset\mathcal{C}^{2}$, $f\in H^{1}(D)$ and (7), (8) are valid. Then, the directional directional derivative of the cost (18) with respect to functional variations $g+\lambda h$, $g,\;h$ satisfying the above assumptions and $g(\mathbf{x}_{0})=0,\;h(\mathbf{x}_{0})=0$ is given by:

[TABLE]

A similar calculation can be performed for the pointwise cost functional (19), term by term.

Denote by $\mu\in L^{2}(D)$ the last term in (20). In fact, at least $H^{1}_{0}(D)$ regularity is valid for $\mu$, under our assumptions, and $u\in H^{3}(D)$ depends linearly and continuously on $\mu$ in these spaces. Then, $\frac{\partial u}{\partial\mathbf{n}}$ and the first term in (29) depend linearly and continuously on $\mu\in H^{1}_{0}(D)$. There is a unique $p\in H^{-1}(D)$ such that the linear continuous functional on the left side below, can be expressed as:

[TABLE]

where $<\cdot,\cdot>$ denotes the pairing in $H^{-1}(D)\times H^{1}_{0}(D)$. We define the adjoint equation by rewriting (30) in the form:

[TABLE]

where $u\in H^{3}(D)\cap H^{1}_{0}(D)$ has to be considered an arbitrary test function. This approach, called the transposition method, is due to J.L. Lions, see [8], Ch.2.5.4, where a related cost functional is considered. Moreover, one adds the null Dirichlet condition to (31), due to the employed class of test functions and a formal interpretation of equation (31) is:

[TABLE]

in $D$, where $\delta_{\mathbf{x}_{0}}$ is the Dirac distribution concentrated in $\mathbf{x}_{0}$ and $P(y_{\varepsilon})$ is generated starting from the linear form in $u$:

[TABLE]

We have proved:

Proposition 4**.**

The adjoint equation associated to the cost functional (18) and the state system (15), (16) is given by (31) and it has a unique transposition solution $p\in H^{-1}(D)$.

In the next section, we shall employ a direct approach based on Prop.3, in the finite dimensional optimization problem obtained via discretization. Due to the very weak regularity properties of the adjoint state discussed in Prop.4, the usual method to compute the gradient of the cost via the adjoint system, seems difficult to apply (but we indicate a regularized variant in one example).

3 Discretization and numerical tests

We assume $D\subset\mathbb{R}^{2}$ to be polygonal domain and $\mathcal{T}_{h}$ a mesh in $D$ of size $h>0$, each element $T$ of the mesh is a triangle.

We set the $\mathbb{P}_{1}$ Lagrangian finite element spaces

[TABLE]

Let $\varphi_{h}\in W_{h}$ be an approximation of the obstacle $\varphi$. We set

[TABLE]

The discretized state equation, associated to (15), (16), in variational formulation, is: find $y_{h}\in\mathcal{K}_{h}$ such that

[TABLE]

with $f\in L^{2}(D)$. Here $H_{\eta}\in\mathcal{C}^{1}(\mathbb{R})$, $\eta>0$, is a regularization of the Heaviside function, given by

[TABLE]

For the following example of maximal monotone graph $\beta\subset\mathbb{R}\times\mathbb{R}$,

[TABLE]

we introduce the approximation $\beta_{\eta,\varepsilon_{2}}\in\mathcal{C}^{1}(\mathbb{R})$, ($\varepsilon_{2}>0$) given by

[TABLE]

We take $\eta>\varepsilon$ and $\varepsilon_{2}>\varepsilon$ to avoid some numerical difficulties, see [11].

Each $g_{h}\in W_{h}$ can be written as $g_{h}=\sum_{i=1}^{n}G_{i}\phi_{i}$ where $(\phi_{i})_{1\leq i\leq n}$ is the basis of the hat functions and $n=\mathrm{dim}(W_{h})$. We introduce $\mathcal{J}_{1}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ defined by

[TABLE]

where $G=(G_{i})_{1\leq i\leq n}$, see(18).

We denote $I=\{1,\dots,n\}$ and let $I_{0}\subset I$ be the indices of the vertices $A_{i}$ of $\mathcal{T}_{h}$ such that $\|\mathbf{x}_{0}-A_{i}\|<hC$ with $C\geq 2$. As an example, for the triangulation and $\mathbf{x}_{0}$ of Test 1, we get $I_{0}$ of 13 elements for $C=2$ and 33 elements for $C=3$.

Another choice for $I_{0}$ is the set of indices of the vertices of the triangle $T\in\mathcal{T}_{h}$ containing $\mathbf{x}_{0}$. In this case, the parameter $C$ is not needed.

If $i\in I\setminus I_{0}$, then $\mathbf{x}_{0}$ doesn’t belong to the support of $\phi_{i}$ and

[TABLE]

Let $g_{h}^{0}\in W_{h}$ be a finite element approximation of the initial guess $g^{0}$ in the subsequent algorithm, such that $g_{h}^{0}(\mathbf{x}_{0})=0$ and $g_{h}^{0}=\sum_{i=1}^{n}G_{i}^{0}\phi_{i}$. For $g_{h}=\sum_{i\in I_{0}}G_{i}^{0}\phi_{i}+\sum_{i\in I\setminus I_{0}}R_{i}\phi_{i}$, $R_{i}\in\mathbb{R}$, we get

[TABLE]

In other words, we can choose freely $R_{i}$, for all $i\in I\setminus I_{0}$.

For $i\in I\setminus I_{0}$, we set the elliptic problem: find $u_{i}\in V_{h}$ such that

[TABLE]

For each $i\in I\setminus I_{0}$, the matrix of the linear system (40) is the same, just the right-hand side is different. We can solve effectively by factorization of the matrix, see [5], Chap. 4. We point out that $u_{i}\in V_{h}$, so the linear system has $\dim(V_{h})$ equations.

Corollary 1**.**

Let $g_{h}$ in $W_{h}$ and we assume that $g_{h}(\mathbf{x}_{0})=0$. Then

[TABLE]

Proof.

For given $g_{h}$, we set $y_{h}$ the solution of (35). The partial derivative $\partial_{i}\mathcal{J}_{1}(G)$ is in fact $J_{1}^{\prime}(g_{h})\phi_{i}$. The last one can be computed using (29), where $y_{\varepsilon}$, $u$, $h$ are replaced with $y_{h}$, $u_{i}$, $\phi_{i}$ respectively. Here $\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})$, $\frac{\partial u_{i}}{\partial\mathbf{n}}(\mathbf{x}_{0})$ $\frac{\partial\phi_{i}}{\partial\mathbf{n}}(\mathbf{x}_{0})$ mean

[TABLE]

respectively. ∎

The optimization problem to be solved is

[TABLE]

and we have for all $i\in I\setminus I_{0}$

[TABLE]

The steepest descent direction for (42) is given by

[TABLE]

and we can write the Algorithm:

Step 1. Let $g_{h}^{0}=\sum_{i=1}^{n}G_{i}^{0}\phi_{i}$ be the initial guess, such that $g_{h}^{0}(\mathbf{x}_{0})=0$. Set $R_{i}^{0}=G_{i}^{0}$ for $i\in I\setminus I_{0}$ and put $k=0$.

Step 2. For $i\in I\setminus I_{0}$, compute $u_{i}$ using (40), then compute $\nabla\mathbb{J}_{1}(R^{k})$ using (41), where $R^{k}=(R_{i}^{k})_{i\in I\setminus I_{0}}$.

Step 3. Find $\lambda_{k}>0$ by line search

[TABLE]

Step 4. Update

[TABLE]

Step 5. We stop the algorithm when $|\mathbb{J}_{1}(R^{k+1})-\mathbb{J}_{1}(R^{k})|<tol$. If not, go to Step 2.

We notice that the values at the vertices of index $i\in I_{0}$ of the parameterization function $g_{h}^{k}=\sum_{i\in I_{0}}G_{i}^{0}\phi_{i}+\sum_{i\in I\setminus I_{0}}R_{i}^{k}\phi_{i}$ are fixed, in order to satisfy the condition $g_{h}^{k}(\mathbf{x}_{0})=0$, for all $k\geq 0$. In fact, $g_{h}^{k}$ is fixed in a neighborhood of $\mathbf{x}_{0}$.

Test 1.

The domain is $D=]0,1[\times]0,1[$, the obstacle is $\varphi:D\rightarrow\mathbb{R}$, $\varphi=-0.5$, $f:D\rightarrow\mathbb{R}$, $f=-100$ and $\alpha=0$. We use a mesh for $D$ of 26870 vertices, 53138 triangles and size $h=1/150$. The following values have been used: penalization parameter $\varepsilon=0.0001$, parameters for smoothed Heaviside and $\beta$ are $\eta=0.05$, $\varepsilon_{2}=0.01$. The starting domain is a disk of radius $r_{0}=0.25$ and center $(0.5,0.5)$, while $\mathbf{x}_{0}=(0.25,0.5)$. For the stopping criterion, we use $tol=10^{-6}$. For the cases a), b), c) we use the descent direction $-\nabla\mathbb{J}_{1}$ and for the case d) we use the simplified (partial) direction $-y_{h}p_{h}$.

Case a). We use $\alpha=0$ and $C=2$; $I_{0}$ has 13 vertices.

The steepest descent direction algorithm stops after 5 iterations and the objective functions decreases from $36.82$ to $2.80\times 10^{-7}$, see Figure 1. The evolution of $\Omega$ is prezented in Figure 3. The boundary of the final $\Omega$ is plotted in Figure 2.

In Figure 4, we show the final state solution $y_{h}$. At the final iteration, we obtained

[TABLE]

then $\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})-\alpha=-0.000529=\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})$.

In Figure 4, we see that the obtained optimal geometry “generates” a saddle point for the optimal $y_{h}$ and $\mathbf{x}_{0}$ is close to this saddle point. That’s why the final cost has such a low value. Notice however that the norm of $\nabla y_{h}(\mathbf{x}_{0})\neq 0$, that is the tangential derivative of the optimal state in $\mathbf{x}_{0}$ is not zero and the Dirichlet boundary condition for $y_{h}$ is slightly violated due to the numerical errors (but the numerical minimization is correctly performed).

Case b). We use $\alpha=0$ and $I_{0}$ with 3 vertices.

If $I_{0}$ is the set of indices of the vertices of the triangle $T\in\mathcal{T}_{h}$ containing $\mathbf{x}_{0}$, the algorithm stops after 3 iterations, the values of the objective functions are: $36.82$ (initial), $3.40515$, $0.000325$, $7.05\times 10^{-7}$ (final).

We see in Figure 5 that the final $\partial\Omega$ has two connected boundaries (topological optimization), but one is very small. It is a triangular polygonal line containing $\mathbf{x}_{0}$. The triangle, having as boundary the polygonal line, represents a hole for the optimal $\Omega$. We point out that this triangle is not an element of the triangulation, otherwise we get $g_{h}=0$ in this triangle, which contradicts (7).

At the final iteration, we have

[TABLE]

then $\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})-\alpha=-0.000839=\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})$.

Case c). We use $\alpha=1$ and $C=2$; $I_{0}$ has 13 elements.

The algorithm stops after 3 iterations, the values of the objective functions are: $25.69$ (initial), $1.45939$, $0.000395$, $1.26\times 10^{-7}$ (final).

At the final iteration, we obtained

[TABLE]

then $\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})-\alpha=-0.000356$. Some images for the final solution are plotted in Figure 6. Although in Figure 6, we also see a saddle point for $y_{h}$, $\mathbf{x}_{0}$ is sufficiently far from it such that $\frac{\partial y_{h}}{\partial\mathbf{n}}(\mathbf{x}_{0})-1=-0.000356$.

We see that $\frac{\nabla g_{h}(\mathbf{x}_{0})}{|\nabla g_{h}(\mathbf{x}_{0})|}$ is the same for the three cases a), b), c). It is correct, because $g_{h}$ is not modified during the algorithm in the vertices $i\in I_{0}$. We have that $\mathbf{x}_{0}$ belongs to $\cup_{i\in I_{0}}\mathrm{supp}(\phi_{i})$. The set $I_{0}$ in the case b) is included in the set $I_{0}$ of cases a) and c). The triangulation, $\mathbf{x}_{0}$ and the initial $g_{h}$ are the same for the three cases.

The holdall domain $D$ is defined by the user. For rectangular shape, we can built a mesh composed only of right triangles or rectangular elements. Here, we use an unstructured triangulation (irregular grid) which can be employed also for non rectangular $D$.

Case d). We use $\alpha=0$, $C=2$ and $-y_{h}p_{h}$ as descent direction.

We can increase the regularity in the adjoint equation (32) by using $\zeta_{\varepsilon_{1}}(\mathbf{x}-\mathbf{x}_{0})$, a mollifier approximation of the Dirac functional $\delta_{\mathbf{x}_{0}}$. We write $\zeta_{\varepsilon_{1}}(\mathbf{x})=\frac{1}{\varepsilon_{1}^{2}}\zeta\left(\frac{\mathbf{x}}{\varepsilon_{1}}\right)$, where $\zeta\in\mathcal{C}^{\infty}(\mathbb{R}^{2})$ has support in the unitary disk centered at $(0,0)$ and $\int_{\mathbb{R}^{2}}\zeta(\mathbf{x})\,d\mathbf{x}=1$, $\varepsilon_{1}>0$.

An approximate version of the adjoint state equation (32) is: find $p\in H_{0}^{1}(D)$ such that

[TABLE]

Taking $v=u$ in (44) and using (20), we get that

[TABLE]

and the right-hand side is an approximation of L defined by (25). Since $H_{\varepsilon}^{\prime}\geq 0$, the left-hand side is negative, for $h=-y_{\varepsilon}p$ and can be used in the gradient-type algorithms, inspired by [10], although it has a “partial” character.

All the data are the same as in the case a), but here we use the discretization (45) and in the Algorithm, the simplified descent direction $-y_{h}p_{h}$ instead of $-\nabla\mathbb{J}_{1}$. That is, we use just a part of the gradient computed in the previous section (but the descent properties are preserved according to Fig. 7). The objective functions decreases from $36.82$ to $1.15\times 10^{-7}$, see Figure 7. The Algorithm stops after 4 iterations.

The finite element approximation of the adjoint state is: find $p_{h}\in V_{h}$ such that

[TABLE]

For the mollifier approximation of the Dirac function, we use $\varepsilon_{1}=0.05$.

We notice in Figure 8 that the final $\partial\Omega$ has three connected components of the boundary (topological optimization) and the smallest one is containing $\mathbf{x}_{0}$. That is, the final $\Omega$ has two holes.

Test 2.

We can use similar computations, when we have pointwise observations at a finite number of fixed points $\mathbf{x}_{0}^{j}$, $1\leq j\leq\ell$. We assume that $g(\mathbf{x}_{0}^{j})=0$, for all $j$. The objective function is the sum of the expression (18) computed at each point $\mathbf{x}_{0}^{j}$. The set $I_{0}$ is, in this case, the indices of the vertices $A_{i}$ of $\mathcal{T}_{h}$ such that there exists $1\leq j\leq\ell$, $\|\mathbf{x}_{0}^{j}-A_{i}\|<hC$ with $C\geq 2$. The descent direction is the sum of the terms in (43) computed at each point $\mathbf{x}_{0}^{j}$.

We have $\ell=3$, the observation points are:

[TABLE]

For $C=3$ we get $I_{0}$ of 96 elements. The others parameters are the same as in Test 1 and $\alpha=0$. The algorithm stops after 5 iterations, the objective functions decreases from $326.12$ to $1.64\times 10^{-5}$, see Figure 9.

In Figure 10 we plot the initial and final $\partial\Omega$ and in Figure 11 we show the final state solution $y_{h}$. Final $\Omega$ seems to be not connected, while the support of $y_{h}$ is connected. This phenomena appears since in the state variational inequality we use in the penalization term $H_{\eta}(g_{h})$ which is an approximation of the characteristic function of the final $\Omega$, that is the approximations yield important perturbations on the result without affecting its essence.