This page was generated from unit-7-optimization/01_Shape_Derivative_Levelset.ipynb.

7.1 Shape optimization via the shape derivative

In this tutorial we discuss the implementation of shape optimization algorithms to solve

$$\underset{\mathrm{\Omega }}{min}J(\mathrm{\Omega }):={\int }_{\mathrm{\Omega }}f(x)dx,f\in C^1({\mathbf{R}}^d).$$

The analytic solution to this problem is

$${\mathrm{\Omega }}^{\ast }:=\{x\in {\mathbf{R}}^d:f(x)\le 0\}.$$

This problem is solved by fixing an initial guess ${\mathrm{\Omega }}_0\subset {\mathbf{R}}^d$ and then considering the problem

$$\underset{X}{min}J((\text{id}+X)({\mathrm{\Omega }}_0)),$$

where $X:{\mathbf{R}}^d\to {\mathbf{R}}^d$ are (at least) Lipschitz vector fields. We approximate $X$ by a finite element function.

Initial geometry and shape function

We choose $f$ as

$$\begin{array}{r}\begin{array}{rl}f(x,y)=& (\sqrt{(x-a{)}^2+b{y}^2}-1)(\sqrt{(x+a{)}^2+b{y}^2}-1)\\ & (\sqrt{b{x}^2+(y-a{)}^2}-1)(\sqrt{b{x}^2+(y+a{)}^2}-1)-\epsilon .\end{array}\end{array}$$

where $\epsilon =0.001,a=4/5,b=2$. The corresponding zero level sets of this function are as follows. The green area indicates where $f$ is negative.

Let us set up the geometry.

# Generate the geometry
from ngsolve import *

# load Netgen/NGSolve and start the gui
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry

geo = SplineGeometry()
geo.AddCircle((0,0), r = 2.5)

ngmesh = geo.GenerateMesh(maxh = 0.08)
mesh = Mesh(ngmesh)
mesh.Curve(2)

<ngsolve.comp.Mesh at 0x7fb587c54e50>

#gfset.vec[:]=0
#gfset.Set((0.2*x,0))
#mesh.SetDeformation(gfset)
#scene = Draw(gfzero,mesh,"gfset", deformation = gfset)

# define the function f and its gradient
a =4.0/5.0
b = 2
f = CoefficientFunction((sqrt((x - a)**2 + b * y**2) - 1) \
                * (sqrt((x + a)**2 + b * y**2) - 1) \
                * (sqrt(b * x**2 + (y - a)**2) - 1) \
                * (sqrt(b * x**2 + (y + a)**2) - 1) - 0.001)

# gradient of f defined as vector valued coefficient function
grad_f = CoefficientFunction((f.Diff(x),f.Diff(y)))

# vector space for shape gradient
VEC = H1(mesh, order=1, dim=2)

# grid function for deformation field
gfset = GridFunction(VEC)
gfX = GridFunction(VEC)

Draw(gfset)

BaseWebGuiScene

Shape derivative

$$DJ(\mathrm{\Omega })(X)={\int }_{\mathrm{\Omega }}f\text{div}(X)+\mathrm{\nabla }f\cdot Xdx$$

# Test and trial functions
PHI, PSI = VEC.TnT()

# shape derivative
dJOmega = LinearForm(VEC)
dJOmega += (div(PSI)*f + InnerProduct(grad_f, PSI) )*dx

Bilinear form

$$(\phi ,\psi )\mapsto b_{\mathrm{\Omega }}(\phi ,\psi ):={\int }_{\mathrm{\Omega }}(\mathrm{\nabla }\phi +\mathrm{\nabla }{\phi }^{\mathrm{\top }}):\mathrm{\nabla }\psi +\phi \cdot \psi dx.$$

to compute the gradient $X:=\text{grad}J(\mathrm{\Omega })$ by

$$b_{\mathrm{\Omega }}(X,\psi )=DJ(\mathrm{\Omega })(\psi )\text{ for all }\psi \in H^1(\mathrm{\Omega })$$

# bilinear form for H1 shape gradient; aX represents b_\Omega
aX = BilinearForm(VEC)
aX += InnerProduct(grad(PHI) + grad(PHI).trans, grad(PSI)) * dx
aX += InnerProduct(PHI, PSI) * dx

The first optimisation step

Fix an initial domain ${\mathrm{\Omega }}_0$ and define

$${\mathcal{J}}_{{\mathrm{\Omega }}_0}(X):=J((\text{id}+X)({\mathrm{\Omega }}_0)).$$

Then we have the following relation between the derivative of ${\mathcal{J}}_{{\mathrm{\Omega }}_0}$ and the shape derivative of $J$:

$$D{\mathcal{J}}_{{\mathrm{\Omega }}_0}(X_n)(X)=DJ((\text{id}+X_n)({\mathrm{\Omega }}_0))(X\circ (\text{id}+X_n{)}^{-1}).$$

Here

• $(\text{id}+X_n)({\mathrm{\Omega }}_0)$ is current shape

Now we note that $\phi \mapsto \phi \circ (\text{id}+X_n{)}^{-1}$ maps the finite element space on $(\text{id}+X_n)({\mathrm{\Omega }}_0)$ to the finite element space on ${\mathrm{\Omega }}_0$. Therefore the following are equivalent:

• assembling $\phi \mapsto D{\mathcal{J}}_{{\mathrm{\Omega }}_0}(X_n)(\phi )$ on the fixed domain ${\mathrm{\Omega }}_0$

• assembling $\phi \mapsto DJ((\text{id}+X_n)({\mathrm{\Omega }}_0))(\phi )$ on the deformed domain $(\text{id}+X_n)({\mathrm{\Omega }}_0)$.

# deform current domain with gfset
mesh.SetDeformation(gfset)

# assemble on deformed domain
aX.Assemble()
dJOmega.Assemble()

mesh.UnsetDeformation()
# unset deformation

Now let’s make one optimization step with step size ${\alpha }_1>0$:

$${\mathrm{\Omega }}_1=(\text{id}-{\alpha }_1{X}_0)({\mathrm{\Omega }}_0).$$

# compute X_0
gfX.vec.data = aX.mat.Inverse(VEC.FreeDofs(), inverse="sparsecholesky") * dJOmega.vec

print("current cost ", Integrate(f*dx, mesh))
print("Gradient norm ", Norm(gfX.vec))

alpha = 20.0 / Norm(gfX.vec)

gfset.vec[:]=0
scene = Draw(gfset)
# input("A")
gfset.vec.data -= alpha * gfX.vec
mesh.SetDeformation(gfset)
#draw deformed shape
scene.Redraw()
# input("B")

print("cost after gradient step:", Integrate(f, mesh))
mesh.UnsetDeformation()
scene.Redraw()

current cost  63.58544060454411
Gradient norm  441.59290658661877

cost after gradient step: 5.023994232440709

Optimisation loop

Now we can set up an optimisation loop. In the second step we compute

$${\mathrm{\Omega }}_2=(\text{id}-{\alpha }_0{X}_0-{\alpha }_1{X}_1)({\mathrm{\Omega }}_0)$$

by the same procedure as above.

import time


iter_max = 50
gfset.vec[:] = 0
mesh.SetDeformation(gfset)
scene = Draw(gfset,mesh,"gfset")

converged = False

alpha =0.11#100.0 / Norm(gfX.vec)
# input("A")
for k in range(iter_max):
    mesh.SetDeformation(gfset)
    scene.Redraw()
    Jold = Integrate(f, mesh)
    print("cost at iteration ", k, ': ', Jold)

    # assemble bilinear form
    aX.Assemble()

    # assemble shape derivative
    dJOmega.Assemble()

    mesh.UnsetDeformation()

    gfX.vec.data = aX.mat.Inverse(VEC.FreeDofs(), inverse="sparsecholesky") * dJOmega.vec

    # step size control

    gfset_old = gfset.vec.CreateVector()
    gfset_old.data = gfset.vec

    Jnew = Jold + 1
    while Jnew > Jold:

        gfset.vec.data = gfset_old
        gfset.vec.data -= alpha * gfX.vec

        mesh.SetDeformation(gfset)

        Jnew = Integrate(f, mesh)

        mesh.UnsetDeformation()

        if Jnew > Jold:
            print("reducing step size")
            alpha = 0.9*alpha
        else:
            print("linesearch step accepted")
            alpha = alpha*1.5
            break

    print("step size: ", alpha, '\n')

    time.sleep(0.1)
    Jold = Jnew

cost at iteration  0 :  63.58544060454411
linesearch step accepted
step size:  0.165

cost at iteration  1 :  -0.7382545763570348
linesearch step accepted
step size:  0.2475

cost at iteration  2 :  -0.8176867576439005
linesearch step accepted
step size:  0.37124999999999997

cost at iteration  3 :  -0.9325774655646434
linesearch step accepted
step size:  0.556875

cost at iteration  4 :  -1.0701609035569957
linesearch step accepted
step size:  0.8353125

cost at iteration  5 :  -1.1521932913592035
linesearch step accepted
step size:  1.25296875

cost at iteration  6 :  -1.154748226275979
reducing step size
reducing step size
linesearch step accepted
step size:  1.52235703125

cost at iteration  7 :  -1.1549550076062747
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.4982276723046877

cost at iteration  8 :  -1.15504717135492
reducing step size
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.3270326873287925

cost at iteration  9 :  -1.1554533885077471
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.4511102435940346

cost at iteration  10 :  -1.155597635089538
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.4281101462330694

cost at iteration  11 :  -1.1556252662136504
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.4054746004152756

cost at iteration  12 :  -1.1558537564268163
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.3831978279986936

cost at iteration  13 :  -1.1559568694538345
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.3612741424249144

cost at iteration  14 :  -1.156152895177702
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.3396979472674797

cost at iteration  15 :  -1.1561632124479113
reducing step size
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.1866173613229614

cost at iteration  16 :  -1.1563069254537834
reducing step size
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0510285285313932

cost at iteration  17 :  -1.1563735494124863
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.034369726354171

cost at iteration  18 :  -1.1564005065055163
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0179749661914572

cost at iteration  19 :  -1.1564106607179754
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0018400629773228

cost at iteration  20 :  -1.1564510867511106
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9859608979791321

cost at iteration  21 :  -1.1565106488237795
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9703334177461629

cost at iteration  22 :  -1.1565853953839687
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9549536330748863

cost at iteration  23 :  -1.156665227033129
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9398176179906492

cost at iteration  24 :  -1.1567406736299617
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.027690565272775

cost at iteration  25 :  -1.1567415931350467
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0114016698132016

cost at iteration  26 :  -1.1567527690291632
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9953709533466624

cost at iteration  27 :  -1.1567746927196771
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9795943237361178

cost at iteration  28 :  -1.1568044739349084
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9640677537049005

cost at iteration  29 :  -1.1568390508947293
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9487872798086778

cost at iteration  30 :  -1.15687395189781
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9337490014237104

cost at iteration  31 :  -1.1569061440231116
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0210545330568275

cost at iteration  32 :  -1.1569107199788815
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0048708187078768

cost at iteration  33 :  -1.1569189782060194
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.988943616231357

cost at iteration  34 :  -1.156930421933171
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.97326885991409

cost at iteration  35 :  -1.1569442348717474
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9578425484844518

cost at iteration  36 :  -1.1569589035206531
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9426607440909731

cost at iteration  37 :  -1.1569732064699092
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0307995236634793

cost at iteration  38 :  -1.1569759485587745
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0144613512134133

cost at iteration  39 :  -1.156980232839626
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9983821387966807

cost at iteration  40 :  -1.1569860444980924
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9825577818967535

cost at iteration  41 :  -1.1569932002791836
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.96698424105369

cost at iteration  42 :  -1.1570010762794725
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.05739726759221

cost at iteration  43 :  -1.1570011715183581
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0406375209008738

cost at iteration  44 :  -1.1570020152391864
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0241434161945948

cost at iteration  45 :  -1.1570041949268655
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.0079107430479106

cost at iteration  46 :  -1.157007847452682
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9919353577706012

cost at iteration  47 :  -1.1570130270308265
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9762131823499371

cost at iteration  48 :  -1.1570192444424463
reducing step size
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  0.9607402034096908

cost at iteration  49 :  -1.1570260548839857
reducing step size
reducing step size
reducing step size
linesearch step accepted
step size:  1.050569412428497

Using SciPy optimize toolbox

We use the package scipy.optimize and wrap the shape functions around it.

We first setup the initial geometry as a circle of radius r = 2.5 (as before)

# The code in this cell is the same as in the example above.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry

import numpy as np

ngsglobals.msg_level = 1

geo = SplineGeometry()
geo.AddCircle((0,0), r = 2.5)

ngmesh = geo.GenerateMesh(maxh = 0.08)
mesh = Mesh(ngmesh)
mesh.Curve(2)

# define the function f

a =4.0/5.0
b = 2
f = CoefficientFunction((sqrt((x - a)**2 + b * y**2) - 1) \
                * (sqrt((x + a)**2 + b * y**2) - 1) \
                * (sqrt(b * x**2 + (y - a)**2) - 1) \
                * (sqrt(b * x**2 + (y + a)**2) - 1) - 0.001)

# Now we define the finite element space VEC in which we compute the shape gradient

# element order
order = 1

VEC = H1(mesh, order=order, dim=2)

# define test and trial functions
PHI = VEC.TrialFunction()
PSI = VEC.TestFunction()


# define grid function for deformation of mesh
gfset = GridFunction(VEC)
gfset.Set((0,0))

# only for new gui
#scene = Draw(gfset, mesh, "gfset")

#if scene:
#    scene.setDeformation(True)

# plot the mesh and visualise deformation
#Draw(gfset,mesh,"gfset")
#SetVisualization (deformation=True)

 Generate Mesh from spline geometry
 Curve elements, order = 2

Shape derivative

$$DJ(\mathrm{\Omega })(X)={\int }_{\mathrm{\Omega }}f\text{div}(X)+\mathrm{\nabla }f\cdot Xdx$$

fX = LinearForm(VEC)

# analytic shape derivative
fX += f*div(PSI)*dx + (f.Diff(x,PSI[0])) *dx + (f.Diff(y,PSI[1])) *dx

Bilinear form for shape gradient

Define the bilinear form

$$(\phi ,\psi )\mapsto b_{\mathrm{\Omega }}(\phi ,\psi ):={\int }_{\mathrm{\Omega }}(\mathrm{\nabla }\phi +\mathrm{\nabla }{\phi }^{\mathrm{\top }}):\mathrm{\nabla }\psi +\phi \cdot \psi dx$$

to compute the gradient $X=:\text{grad}J$ by

$$b_{\mathrm{\Omega }}(X,\psi )=DJ(\mathrm{\Omega })(\psi )\text{ for all }\psi \in H^1(\mathrm{\Omega })$$

# Cauchy-Riemann descent CR
CR = False

# bilinear form for gradient
aX = BilinearForm(VEC)
aX += InnerProduct(grad(PHI)+grad(PHI).trans, grad(PSI))*dx + InnerProduct(PHI,PSI)*dx

## Cauchy-Riemann regularisation

if CR == True:
    gamma = 100
    aX += gamma * (PHI.Deriv()[0,0] - PHI.Deriv()[1,1])*(PSI.Deriv()[0,0] - PSI.Deriv()[1,1]) *dx
    aX += gamma * (PHI.Deriv()[1,0] + PHI.Deriv()[0,1])*(PSI.Deriv()[1,0] + PSI.Deriv()[0,1]) *dx


aX.Assemble()
invaX = aX.mat.Inverse(VEC.FreeDofs(), inverse="sparsecholesky")

gfX = GridFunction(VEC)

Wrapping the shape function and its gradient

Now we define the shape function $J$ and its gradient grad(J) use the shape derivative $DJ$. The arguments of $J$ and grad(J) are vector in ${\mathbf{R}}^d$, where $d$ is the dimension of the finite element space.

def J(x_):

    gfset.Set((0,0))
    # x_ NumPy vector
    gfset.vec.FV().NumPy()[:] += x_

    mesh.SetDeformation(gfset)

    cost_c = Integrate (f, mesh)

    mesh.UnsetDeformation()

    Redraw(blocking=True)

    return cost_c

def gradJ(x_, euclid = False):

    gfset.Set((0,0))
    # x_ NumPy vector
    gfset.vec.FV().NumPy()[:] += x_

    mesh.SetDeformation(gfset)

    fX.Assemble()

    mesh.UnsetDeformation()

    if euclid == True:
        gfX.vec.data = fX.vec
    else:
        gfX.vec.data = invaX * fX.vec

    return gfX.vec.FV().NumPy().copy()

Gradient descent and Armijo rule

We use the functions $J$ and grad J to define a steepest descent method.We use scipy.optimize.line_search_armijo to compute the step size in each iteration. The Arjmijo rule reads

$$J((\text{Id}+{\alpha }_k{X}_k)({\mathrm{\Omega }}_k))\le J({\mathrm{\Omega }}_k)-c_1{\alpha }_k\Vert \mathrm{\nabla }J({\mathrm{\Omega }}_k){\Vert }^2,$$

where $c_1:=1e-4$ and ${\alpha }_k$ is the step size.

# import scipy linesearch method
from scipy.optimize.linesearch import line_search_armijo

def gradient_descent(x0, J_, gradJ_):

    xk_ = np.copy(x0)

    # maximal iteration
    it_max = 50
    # count number of function evals
    nfval_total = 0

    print('\n')
    for i in range(1, it_max):

        # Compute a step size using a line_search to satisfy the Wolf
        # compute shape gradient grad J
        grad_xk = gradJ_(xk_,euclid = False)

        # compute descent direction
        pk = -gradJ_(xk_)

        # eval cost function
        fold = J_(xk_)

        # perform armijo stepsize

        if CR == True:
            alpha0 = 0.15
        else:
            alpha0 = 0.11



        step, nfval, b = line_search_armijo(J_, xk_, pk = pk, gfk = grad_xk, old_fval = fold, c1=1e-4, alpha0 = alpha0)
        nfval_total += nfval

        # update the shape and print cost and gradient norm
        xk_ = xk_ - step * grad_xk
        print('Iteration ', i, '| Cost ', fold, '| grad norm', np.linalg.norm(grad_xk))

        mesh.SetDeformation(gfset)
        scene.Redraw()
        mesh.UnsetDeformation()

        if np.linalg.norm(gradJ_(xk_)) < 1e-4:

            #print('#######################################')
            print('\n'+'{:<20}  {:<12} '.format("##################################", ""))
            print('{:<20}  {:<12} '.format("### success - accuracy reached ###", ""))
            print('{:<20}  {:<12} '.format("##################################", ""))
            print('{:<20}  {:<12} '.format("gradient norm: ", np.linalg.norm(gradJ_(xk_))))
            print('{:<20}  {:<12} '.format("n evals f: ", nfval_total))
            print('{:<20}  {:<12} '.format("f val: ", fold) + '\n')
            break
        elif i == it_max-1:

            #print('######################################')
            print('\n'+'{:<20}  {:<12} '.format("#######################", ""))
            print('{:<20}  {:<12} '.format("### maxiter reached ###", ""))
            print('{:<20}  {:<12} '.format("#######################", ""))
            print('{:<20}  {:<12} '.format("gradient norm: ", np.linalg.norm(gradJ_(xk_))))
            print('{:<20}  {:<12} '.format("n evals f: ", nfval_total))
            print('{:<20}  {:<12} '.format("f val: ", fold) + '\n')

/tmp/ipykernel_88245/197104950.py:2: DeprecationWarning: Please use `line_search_armijo` from the `scipy.optimize` namespace, the `scipy.optimize.linesearch` namespace is deprecated.
  from scipy.optimize.linesearch import line_search_armijo

Now we are ready to run our first optimisation algorithm:

gfset.vec[:]=0
x0 = gfset.vec.FV().NumPy() # initial guess = 0
scene = Draw(gfset, mesh, "gfset")
gradient_descent(x0, J, gradJ)

Iteration  1 | Cost  63.58544060454411 | grad norm 441.59290658661877
Iteration  2 | Cost  -0.7382545763570351 | grad norm 4.795270097704366
Iteration  3 | Cost  -0.7857727667400566 | grad norm 4.810253357263991
Iteration  4 | Cost  -0.8340100892582087 | grad norm 4.778142295070978
Iteration  5 | Cost  -0.8819829476436056 | grad norm 4.690302282301047
Iteration  6 | Cost  -0.9285214898272736 | grad norm 4.5397982631673415
Iteration  7 | Cost  -0.9723556895806401 | grad norm 4.322976857862516
Iteration  8 | Cost  -1.0122531578877725 | grad norm 4.041051815268497
Iteration  9 | Cost  -1.0471883778786093 | grad norm 3.701224549853241
Iteration  10 | Cost  -1.076504443665547 | grad norm 3.3168858328915554
Iteration  11 | Cost  -1.1000162626592398 | grad norm 2.906202755156926
Iteration  12 | Cost  -1.1180150904153732 | grad norm 2.4897881641333446
Iteration  13 | Cost  -1.1311719155347586 | grad norm 2.0875804970132634
Iteration  14 | Cost  -1.1403746530625036 | grad norm 1.7159371235340435
Iteration  15 | Cost  -1.1465562098318822 | grad norm 1.3857186789182068
Iteration  16 | Cost  -1.1505622552842771 | grad norm 1.1024720780232573
Iteration  17 | Cost  -1.1530807674166712 | grad norm 0.8664390301856686
Iteration  18 | Cost  -1.154625445865094 | grad norm 0.6744733131571226
Iteration  19 | Cost  -1.1555550744962926 | grad norm 0.5214023587336313
Iteration  20 | Cost  -1.1561072737628182 | grad norm 0.40125724712233707
Iteration  21 | Cost  -1.1564329926583787 | grad norm 0.30812647847844515
Iteration  22 | Cost  -1.1566250729528902 | grad norm 0.2366460842080817
Iteration  23 | Cost  -1.156739221824035 | grad norm 0.18221681468365375
Iteration  24 | Cost  -1.1568082482769877 | grad norm 0.1410442451728767
Iteration  25 | Cost  -1.156851212982562 | grad norm 0.11008000257674355
Iteration  26 | Cost  -1.1568790835671514 | grad norm 0.08691942912250207
Iteration  27 | Cost  -1.156898134112397 | grad norm 0.06968543512602463
Iteration  28 | Cost  -1.1569119439709799 | grad norm 0.05692156240844273
Iteration  29 | Cost  -1.1569225568515866 | grad norm 0.04750014325185838
Iteration  30 | Cost  -1.156931146061588 | grad norm 0.040550313810789974
Iteration  31 | Cost  -1.1569383917279348 | grad norm 0.03540437114758971
Iteration  32 | Cost  -1.1569446953786267 | grad norm 0.03155781707996517
Iteration  33 | Cost  -1.1569502999928782 | grad norm 0.028637497719290292
Iteration  34 | Cost  -1.1569553580085303 | grad norm 0.02637399371759419
Iteration  35 | Cost  -1.1569599695029114 | grad norm 0.024577086183132767
Iteration  36 | Cost  -1.1569642036044872 | grad norm 0.023114746696247402
Iteration  37 | Cost  -1.1569681109470848 | grad norm 0.021896245450960167
Iteration  38 | Cost  -1.1569717306087155 | grad norm 0.02085939054432362
Iteration  39 | Cost  -1.1569750939270342 | grad norm 0.019961362418991883
Iteration  40 | Cost  -1.1569782269881668 | grad norm 0.019172877193450753
Iteration  41 | Cost  -1.1569811514952684 | grad norm 0.018471912284186168
Iteration  42 | Cost  -1.156983887154983 | grad norm 0.017845413165404594
Iteration  43 | Cost  -1.1569864496512077 | grad norm 0.017278988193844186
Iteration  44 | Cost  -1.1569888554425787 | grad norm 0.01676565064074166
Iteration  45 | Cost  -1.1569911182298955 | grad norm 0.016298795918093117
Iteration  46 | Cost  -1.1569932500627187 | grad norm 0.01587314690341795
Iteration  47 | Cost  -1.1569952619910004 | grad norm 0.015485165897309769
Iteration  48 | Cost  -1.1569971633976466 | grad norm 0.01512970090315433
Iteration  49 | Cost  -1.1569989636230704 | grad norm 0.014804371462622527

#######################
### maxiter reached ###
#######################
gradient norm:        0.01450646124756406
n evals f:            49
f val:                -1.1569989636230704

L-BFGS method

Now we use the L-BFGS method provided by SciPy. The BFGS method requires the shape function $J$ and its gradient grad J. We can also specify additional arguments in options, such as maximal iterations and the gradient tolerance.

In the BFGS method we replace

$$b_{\mathrm{\Omega }}(X,\phi )=\mathrm{\nabla }J(\mathrm{\Omega })(\phi )$$

by

$$H_{\mathrm{\Omega }}(X,\phi )=\mathrm{\nabla }J(\mathrm{\Omega })(\phi ),$$

where $H_{\mathrm{\Omega }}$ is an approximation of the second shape derivative at $\mathrm{\Omega }$. On the discrete level we solve

$$\begin{array}{r}\begin{array}{rl}\text{ solve }{B}_n{X}_k& =\mathrm{\nabla }J({\mathrm{\Omega }}_k)\\ & \\ \text{ update }{s}_k& :=-{\alpha }_k{p}_k\\ y_k& :=\mathrm{\nabla }{J}_h({\mathrm{\Omega }}_{k+1})-\mathrm{\nabla }{J}_h({\mathrm{\Omega }}_k)\\ B_{k+1}& :=B_k+\frac{y_k\otimes y_k}{y_k\cdot s_k}+\frac{B_k{s}_k\otimes B_k{s}_k}{B_k{y}_k\cdot B_k{s}_k}.\end{array}\end{array}$$

from scipy.optimize import minimize

x0 = gfset.vec.FV().NumPy()

# options for optimiser

options = {"maxiter":1000,
              "disp":True,
              "gtol":1e-10}

# we use quasi-Newton method L-BFGS
minimize(J, x0, method='L-BFGS-B', jac=gradJ, options=options)

RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N =         7760     M =           10
 This problem is unconstrained.

At X0         0 variables are exactly at the bounds

At iterate    0    f= -1.15700D+00    |proj g|=  9.63932D-04

At iterate    1    f= -1.15700D+00    |proj g|=  1.78419D-03

At iterate    2    f= -1.15703D+00    |proj g|=  1.15992D-03

At iterate    3    f= -1.15706D+00    |proj g|=  3.29684D-03
  ys=-1.811E-04  -gs= 2.687E-03 BFGS update SKIPPED

At iterate    4    f= -1.15707D+00    |proj g|=  2.17392D-03

At iterate    5    f= -1.15709D+00    |proj g|=  9.20828D-04

At iterate    6    f= -1.15709D+00    |proj g|=  1.22511D-03

At iterate    7    f= -1.15709D+00    |proj g|=  1.26433D-03

At iterate    8    f= -1.15709D+00    |proj g|=  1.34061D-03

 Bad direction in the line search;
   refresh the lbfgs memory and restart the iteration.

At iterate    9    f= -1.15710D+00    |proj g|=  1.52892D-03

At iterate   10    f= -1.15710D+00    |proj g|=  4.54092D-04

At iterate   11    f= -1.15710D+00    |proj g|=  4.52665D-04

At iterate   12    f= -1.15711D+00    |proj g|=  7.15400D-04

At iterate   13    f= -1.15711D+00    |proj g|=  7.15400D-04

           * * *

Tit   = total number of iterations
Tnf   = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip  = number of BFGS updates skipped
Nact  = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F     = final function value

           * * *

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
 7760     13     95      2     1     0   7.154D-04  -1.157D+00
  F =  -1.1571107484943068

CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH

 Warning:  more than 10 function and gradient
   evaluations in the last line search.  Termination
   may possibly be caused by a bad search direction.

      fun: -1.1571107484943068
 hess_inv: <7760x7760 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 2.53834124e-06, -6.12335884e-04,  6.18489141e-04, ...,
       -5.50879208e-06,  4.86565638e-06, -5.41181817e-06])
  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 94
      nit: 13
     njev: 94
   status: 0
  success: True
        x: array([-5.12537914e-04,  6.98505623e-01, -6.98455974e-01, ...,
        4.64527232e-03, -1.45348167e-03,  4.56399822e-03])

Improving mesh quality via Cauchy-Riemann equations

In the previous section we computed the shape gradient grad J:= X of $J$ at $\mathrm{\Omega }$ via

$${\int }_{\mathrm{\Omega }}\mathrm{\nabla }X:\mathrm{\nabla }\phi +X\cdot \phi dx=DJ(\mathrm{\Omega })(\phi )\mathrm{\forall }\phi \in H^1(\mathrm{\Omega }{)}^d.$$

This may lead to overly stretched or even degenerated triangles. One way to improve this is to modify the above equation by

$${\int }_{\mathrm{\Omega }}\mathrm{\nabla }X:\mathrm{\nabla }\phi +X\cdot \phi +\gamma BX\cdot B\phi dx=DJ(\mathrm{\Omega })(\phi )\mathrm{\forall }\phi \in H^1(\mathrm{\Omega }{)}^d,$$

where

$$\begin{array}{r}B:=\left(\begin{array}{cc}-{\partial }_x& {\partial }_y\\ {\partial }_y& {\partial }_x\end{array}\right)\end{array}$$

The two equations $BX=0$ are precisely the Cauchy-Riemann equations $-{\partial }_x{X}_1+{\partial }_y{X}_2=0$ and ${\partial }_y{X}_1+{\partial }_x{X}_2$. So the bigger $\gamma$ the more angle preserving is the gradient. So by adding the B term we enforce conformaty with strength $\gamma$.

This only means we need to add the $\alpha$ term to the above bilinear form aX:

alpha = 100

aX += alpha * (PHI.Deriv()[0,0] - PHI.Deriv()[1,1])*(PSI.Deriv()[0,0] - PSI.Deriv()[1,1]) *dx
aX += alpha * (PHI.Deriv()[1,0] + PHI.Deriv()[0,1])*(PSI.Deriv()[1,0] + PSI.Deriv()[0,1]) *dx