Topology Optimisation using Level set Method

ZeDDKiing99 · July 3, 2025, 5:21pm

I am working based on the works of A.Laurain on “A Level set based structural opitmization code using FEniCS”.
I am trying to convert the code from legacy FEniCS to FEniCSx and try to simulate the same problems (A cantilever beam with point load at the free end). upon completion, I am facing an issue where my Function value which is supposed to decrease, keeps increasing in the optimisation process. I believe the issue is with how I defined the subdomain where the material exists (not completely sure).
Can anyone please go through the code snippet and tell me if I am making any mistake in the FEniCSx syntax.

def solve_pde(V, dx, ds, eps_er, bcd, mu, lmbda, Load, mesh):
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    S1 = 2.0*mu*inner(ufl.sym(grad(u)),ufl.sym(grad(v))) + lmbda*ufl.div(u)*ufl.div(v)
    a = form(S1*eps_er*dx(0) + S1*dx(1))
    A = assemble_matrix(a, bcs=[bcd])
    A.assemble()
    L = form(inner(as_vector([0.0, -1.0]), v) * ds(2))
    b = assemble_vector(L)
    # with b.localForm() as b_local:
    #     b_local.setValues(dofs, -1.0, addv=True)
    apply_lifting(b, [a], [bcd])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode = PETSc.ScatterMode.REVERSE)
    # set_bc(b, [bcd])
    bcd.set(b.array_w)
    U = Function(V)
    A.setOption(PETSc.Mat.Option.SPD, True)
    solver.setFromOptions()
    solver.setOperators(A)
    solver.solve(b, U.x.petsc_vec)
    return U

def shape_der(Vvec, u_vec , eps_er, mu, lmbda, dx, solver, Lag):
    xi = ufl.TestFunction(Vvec)
    rv = 0.0
    for u in u_vec:
        eu,Du,Dxi = ufl.sym(grad(u)),grad(u),grad(xi)
        S1 = 2*mu*(2*inner((Du.T)*eu,Dxi) -inner(eu,eu)*ufl.div(xi))\
         + lmbda*(2*inner( Du.T, Dxi )*ufl.div(u) - ufl.div(u)*ufl.div(u)*ufl.div(xi))
        form_to_assemble = form(-(eps_er*S1*dx(0) + S1*dx(1) + Lag*ufl.div(xi)*dx(1)))
        rv_contribution = assemble_vector(form_to_assemble)
        rv += rv_contribution
    th = Function(Vvec)
    rv.assemble()
    solver.solve(rv, th.x.petsc_vec)
    return th

def Omega(x):
    result = (x[0] >= 0.0) & (x[0] <= lx) & (x[1] >= 0.0) & (x[1] <= ly)
    phi_values = phi.x.array
    # Find dofs where phi < 0
    mask = phi_values < 0

    return result & mask 

main_loop():
    omega_cells = locate_entities(mesh, mesh.topology.dim, Omega)
    cell_tags = np.zeros(num_cells, dtype=dtype)
    cell_tags[omega_cells] = 1
    domains = meshtags(mesh, mesh.topology.dim, np.arange(num_cells), cell_tags)
    dx = Measure('dx',domain=mesh, subdomain_data=domains)

dokken · July 3, 2025, 6:06pm

The code is not complete, i.e. reproducible. Please modify to ensure that one can run it and reproduce your issue.

ZeDDKiing99 · July 4, 2025, 11:44am

from mpi4py import MPI 
from numpy import pi

try:
    from petsc4py import PETSc
    import dolfinx 
    if not dolfinx.has_petsc:
        print("This demo requires DOLFINX")
        exit(0)
except ModuleNotFoundError:
    print("This demo requires petsc4py.")
    exit(0)

import numpy as np, sys, os 
import ufl 
import dolfinx as la 
from dolfinx.mesh import locate_entities_boundary, meshtags, locate_entities
from dolfinx.fem import (Expression, Function, FunctionSpace, dirichletbc, form, locate_dofs_topological
                         , Constant, assemble_scalar, set_bc, functionspace, locate_dofs_geometrical)
from dolfinx.io import XDMFFile
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector
import gmsh
from ufl import grad, inner, Measure, as_vector, dot
from basix.ufl import element
from matplotlib import cm,pyplot as pp 
pp.switch_backend('Agg')
from dolfinx import log
from dolfinx.log import LogLevel
log.set_log_level(LogLevel.ERROR)

dtype = PETSc.ScalarType #Type:ignore

#------------------------Functions-----------------------------#
def solve_pde(V, dx, ds, eps_er, bcd, mu, lmbda, Load, mesh):
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    S1 = 2.0*mu*inner(ufl.sym(grad(u)),ufl.sym(grad(v))) + lmbda*ufl.div(u)*ufl.div(v)
    a = form(S1*eps_er*dx(0) + S1*dx(1))
    A = assemble_matrix(a, bcs=[bcd])
    A.assemble()
    L = form(inner(as_vector([0.0, -1.0]), v) * ds(2))
    b = assemble_vector(L)
    # with b.localForm() as b_local:
    #     b_local.setValues(dofs, -1.0, addv=True)
    apply_lifting(b, [a], [bcd])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode = PETSc.ScatterMode.REVERSE)
    # set_bc(b, [bcd])
    bcd.set(b.array_w)
    U = Function(V)
    A.setOption(PETSc.Mat.Option.SPD, True)
    solver.setFromOptions()
    solver.setOperators(A)
    solver.solve(b, U.x.petsc_vec)
    return U

def shape_der(Vvec, u_vec , eps_er, mu, lmbda, dx, solver, Lag):
    xi = ufl.TestFunction(Vvec)
    rv = 0.0
    for u in u_vec:
        eu,Du,Dxi = ufl.sym(grad(u)),grad(u),grad(xi)
        S1 = 2*mu*(2*inner((Du.T)*eu,Dxi) -inner(eu,eu)*ufl.div(xi))\
         + lmbda*(2*inner( Du.T, Dxi )*ufl.div(u) - ufl.div(u)*ufl.div(u)*ufl.div(xi))
        form_to_assemble = form(-(eps_er*S1*dx(0) + S1*dx(1) + Lag*ufl.div(xi)*dx(1)))
        rv_contribution = assemble_vector(form_to_assemble)
        rv += rv_contribution
    th = Function(Vvec)
    rv.assemble()
    solver.solve(rv, th.x.petsc_vec)
    return th

def hj(v,psi,lx,ly,Nx,Ny,beta):
    for k in range(10):
        Dym = Ny*np.repeat(np.diff(psi,axis=0),[2]+[1]*(Ny-1),axis=0)/ly 
        Dyp = Ny*np.repeat(np.diff(psi,axis=0),[1]*(Ny-1)+[2],axis=0)/ly
        Dxm = Nx*np.repeat(np.diff(psi),[2]+[1]*(Nx-1),axis=1)/lx 
        Dxp = Nx*np.repeat(np.diff(psi),[1]*(Nx-1)+[2],axis=1)/lx          
        g = 0.5*( v[0]*(Dxp + Dxm) + v[1]*(Dyp + Dym)) \
          - 0.5*(np.abs(v[0])*(Dxp - Dxm) + np.abs(v[1])*(Dyp - Dym)) 
        maxv = np.max(abs(v[0]) + abs(v[1]))
        dt  = beta*lx / (Nx*maxv)
        psi = psi - dt*g
    return psi 

def reinit(lx,ly,Nx,Ny,psi):
    Dxs = Nx*(np.repeat(np.diff(psi),[2]+[1]*(Nx-1),axis=1) \
          +np.repeat(np.diff(psi),[1]*(Nx-1)+[2],axis=1))/(2*lx) 
    Dys = Ny*(np.repeat(np.diff(psi,axis=0),[2]+[1]*(Ny-1),axis=0)\
          +np.repeat(np.diff(psi,axis=0),[1]*(Ny-1)+[2],axis=0))/(2*ly)                  
    signum = psi / np.power(psi**2 + ((lx/Nx)**2)*(Dxs**2+Dys**2),0.5)
    for k in range(0,2):      
        Dym = Ny*np.repeat(np.diff(psi,axis=0),[2]+[1]*(Ny-1),axis=0)/ly 
        Dyp = Ny*np.repeat(np.diff(psi,axis=0),[1]*(Ny-1)+[2],axis=0)/ly
        Dxm = Nx*np.repeat(np.diff(psi),[2]+[1]*(Nx-1),axis=1)/lx 
        Dxp = Nx*np.repeat(np.diff(psi),[1]*(Nx-1)+[2],axis=1)/lx               
        Kp = np.sqrt((np.maximum(Dxm,0))**2 + (np.minimum(Dxp,0))**2 \
           + (np.maximum(Dym,0))**2 + (np.minimum(Dyp,0))**2)
        Km = np.sqrt((np.minimum(Dxm,0))**2 + (np.maximum(Dxp,0))**2 \
           + (np.minimum(Dym,0))**2 + (np.maximum(Dyp,0))**2)           
        g  = np.maximum(signum,0)*Kp + np.minimum(signum,0)*Km
        psi  = psi - (0.5*lx/Nx)*(g - signum)
    return psi 

def comp_lsf(px,py,phi,phi_mat,dofsV_max):

    phi_array = phi.x.array
    for dof in range(dofsV_max):
        if np.rint(px[dof]) % 2 == 0:
            cx, cy = np.int_(np.rint([px[dof]/2, py[dof]/2]))
            phi_array[dof] = phi_mat[cy, cx]
        else:
            cx, cy = np.int_(np.floor([px[dof]/2, py[dof]/2]))
            phi_array[dof] = 0.25 * (phi_mat[cy, cx] + phi_mat[cy+1, cx] +
                                     phi_mat[cy, cx+1] + phi_mat[cy+1, cx+1])
    # Assign the modified array back to phi.x.array
    phi.x.array[:] = phi_array
    return phi

# def build_nullspace(V: FunctionSpace):
#     bs = V.dofmap.index_map_bs # no. of scalar components per dof (2 or 3)
#     length0 = V.dofmap.index_map.size_local # Scalar dof owned by MPI process
#     basis = [la.vector(V.dofmap.index_map, bs= bs, dtype = dtype) for i in range(6)] # creates 6 empty basis vectors
#     b = [b.array for b in basis]

#     dofs = [V.sub(i).dofmap.list.flatten() for i in range(3)]

#     for i in range(3):
#         b[i][dofs[i]] = 1.0 # sets 3 translational basis vectors to 1

#     x = V.tabulate_dof_coordinates()
#     dofs_block = V.dofmap.list.flatten()
#     x0, x1, x2 = x[dofs_block, 0], x[dofs_block,1], x[dofs_block, 2]
#     b[3][dofs[0]] = -x1
#     b[3][dofs[1]] = x0
#     b[4][dofs[0]] = x2
#     b[4][dofs[2]] = -x0
#     b[5][dofs[2]] = x1
#     b[5][dofs[1]] = -x2

#     la.orthonormalize(basis)

#     basis_petsc = [PETSc.Vec().createWithArray(x[: bs*length0], bsize = 3, comm = V.mesh.comm)  # type: ignore
#                    for x in b]
    
#     return PETSc.NullSpace().create(vectors=basis_petsc)

Lag,Nx,Ny,lx,ly = 40.0, 302, 151, 2.0, 1.0
XX,YY = np.meshgrid(np.linspace(0.0,lx,Nx+1),np.linspace(0.0,ly,Ny+1))  
mesh = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0.0, 0.0]), np.array([lx, ly])],
    [Nx, Ny],
    cell_type=dolfinx.mesh.CellType.triangle,
    diagonal=dolfinx.mesh.DiagonalType.crossed
)
def all_boundary(x):
    return np.full(x.shape[1], True, dtype=bool)
def left_boundary(x):
    return np.isclose(x[0],0.0)
def load_boundary(x):
    tol = 0.3
    return (np.isclose(x[0], lx)) & (np.abs(x[1] - ly/2) < tol)

all_boundary_facet = locate_entities_boundary(mesh, dim=mesh.topology.dim-1, marker=all_boundary)
boundary_facet = locate_entities_boundary(mesh, dim=mesh.topology.dim-1, marker=left_boundary)
load_facet = locate_entities_boundary(mesh, dim=mesh.topology.dim-1, marker=load_boundary)
all_boundary_tag = np.zeros(len(all_boundary_facet), dtype=dtype)
left_mask = np.isin(all_boundary_facet, boundary_facet)
load_mask = np.isin(all_boundary_facet, load_facet)
all_boundary_tag[left_mask] = 1
all_boundary_tag[load_mask] = 2
boundary_tag = meshtags(
    mesh, mesh.topology.dim-1, all_boundary_facet, all_boundary_tag
)
ds = Measure('ds', subdomain_data=boundary_tag)

Vvec = functionspace(mesh, element("Lagrange", mesh.topology.cell_name(), 1, shape=(mesh.geometry.dim,)))
bcd = dirichletbc(
    np.zeros(2, dtype=dtype), 
    locate_dofs_topological(Vvec, entity_dim=mesh.topology.dim-1,entities=boundary_facet), V=Vvec
)

Load = np.array([[lx,ly/2]])
Load_mag = as_vector([0.0, -1.0])

phi_mat =  -np.cos(8.0/lx*pi*XX) * np.cos(4.0*pi*YY) - 0.4\
  
#-------------------------INIT END----------------------------------#

eps_er, E, nu = 0.001, 1.0, 0.3
# mu, lmbda = Constant(E/(2*(1+nu))), Constant(E*nu/((1+nu)*(1-2*nu)))
mu, lmbda = (E/(2*(1+nu))),(E*nu/((1+nu)*(1-2*nu)))

beta0_init, ls, ls_max, gamma, gamma2 = 0.5, 0, 3, 0.8, 0.8
beta0 = beta0_init
beta = beta0
ItMax, It, stop = int(1.5*Nx), 0, False 
J = np.zeros(ItMax)

#---------------------------assigning values end------------------#

V = functionspace(mesh, ("Lagrange", 1))
VolUnit = Function(V)
VolUnit.interpolate(lambda x: np.ones(x.shape[1]))
U = [0]*len(Load)

gdim = mesh.geometry.dim
dofsV = V.tabulate_dof_coordinates()
dofsVvec = Vvec.tabulate_dof_coordinates()
px, py = (dofsV[:,0]/lx)*2*Nx, (dofsV[:,1]/ly)*2*Ny
pxvec, pyvec = (dofsVvec[:,0]/lx)*2*Nx, (dofsVvec[:,1]/ly)*2*Ny
pxvec, pyvec = np.repeat(pxvec,2) , np.repeat(pyvec,2)
dofsV_max, dofsVvec_max = ((Nx+1)*(Ny+1) + Nx*Ny) * np.array([1,2])

#-----------------------Setting Coordinates-----------------------#
from scipy.interpolate import griddata
phi = Function(V)
phi = comp_lsf(px, py, phi, phi_mat, dofsV_max)
# def Omega(x):
#     return (x[0] >= 0.0) & (x[0] <= lx) & (x[1] >= 0.0) & (x[1] <= ly) & phi.x.array(x) < 0.0
def Omega(x):
    result = (x[0] >= 0.0) & (x[0] <= lx) & (x[1] >= 0.0) & (x[1] <= ly)
    phi_values = phi.x.array
    # Find dofs where phi < 0
    mask = phi_values < 0

    return result & mask  # Indices of dofs in Omega

num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
dX = Measure('dx')
n = ufl.FacetNormal(mesh)

theta, xi = ufl.TrialFunction(Vvec), ufl.TestFunction(Vvec)
av = form((inner(grad(theta),grad(xi)) +0.1*inner(theta,xi))*dX\
         + 1.0e4*(inner(dot(theta,n),dot(xi,n)) * (ds(0)+ds(1)+ds(2))))
Av = assemble_matrix(av)
Av.assemble()

opts = PETSc.Options()
opts["ksp_type"] = "cg" 
opts["ksp_rtol"] = 1.0e-8 
opts["pc_type"] = "gamg"

solver = PETSc.KSP().create(mesh.comm)
solver.setFromOptions()
solver.setOperators(Av)

#--------------------------------MAIN LOOP-----------------------------#

while It < ItMax and stop == False:
    omega_cells = locate_entities(mesh, mesh.topology.dim, Omega)
    cell_tags = np.zeros(num_cells, dtype=dtype)
    cell_tags[omega_cells] = 1
    domains = meshtags(mesh, mesh.topology.dim, np.arange(num_cells), cell_tags)
    dx = Measure('dx',domain=mesh, subdomain_data=domains)

    for k in range(len(Load)):
        U[k] = solve_pde(Vvec, dx, ds, eps_er, bcd, mu, lmbda, Load[k], mesh)     #CHECK SOLVE_PDE
    compliance = 0
    for u in U:
        eU = ufl.sym(grad(u))
        S1 = 2.0 * mu * inner(eU, eU) + lmbda * ufl.tr(eU)**2
        integrand = eps_er * S1 * dx(0) + S1 * dx(1)
        L = form(integrand)
        compliance += assemble_scalar(L)

    vol = form(VolUnit * dx(1))
    Vol = assemble_scalar(vol)
    J[It] = compliance + Lag * Vol

    #--------------------LINE SEARCH-------------------------#
    if It>0 and J[It] > J[It-1] and ls < ls_max:
        ls += 1
        beta *= gamma 
        phi_mat, phi = phi_mat_old, phi_old
        phi_mat = hj(th_mat, phi_mat, lx, ly, Nx, Ny, beta)
        phi= comp_lsf(px, py, phi, phi_mat, dofsV_max)
        print("Line search iteration: %s" % ls)
    else:
        print('************ ITERATION NUMBER %s' % It)                   
        print('Function value        : %.2f' % J[It])
        print('Compliance            : %.5f' % compliance)
        print('Volume fraction       : %.2f' % (Vol/(lx*ly)))

        if ls == ls_max:
            beta0 = max(beta0 * gamma2, 0.1*beta0_init)  
        if ls == 0:
            beta0 = min(beta0 / gamma2, 1)
        
        ls, beta, It = 0, beta0, It+1

        th = shape_der(Vvec, U, eps_er, mu, lmbda, dx, solver, Lag)
        th_array = th.x.array
        th_mat = [np.zeros((Ny+1,Nx+1)), np.zeros((Ny+1,Nx+1))]
        for dof in range(0, dofsVvec_max, 2):
            if np.rint(pxvec[dof]) %2 == .0:
                cx, cy = np.int_(np.rint([pxvec[dof]/2,pyvec[dof]/2]))
                th_mat[0][cy,cx] = th_array[dof]
                th_mat[1][cy,cx] = th_array[dof+1]
        phi_old, phi_mat_old = phi, phi_mat
        phi_mat = hj(th_mat, phi_mat, lx,ly,Nx,Ny, beta)
        if np.mod(It,5) == 0:
            phi_mat = reinit(lx,ly,Nx,Ny,phi_mat)
        phi = comp_lsf(px,py,phi,phi_mat,dofsV_max)

        #--------------------stopping criteria-----------------#
        if It>20 and max(abs(J[It-5:It]-J[It-1]))<2.0*J[It-1]/Nx**2: 
            stop = True 
        #-----------------plot geometry-----------------------#
        if np.mod(It,10)==0 or It==1 or It==ItMax or stop==True:   
            import os
            result = np.where((phi_mat >= -10) & (phi_mat <= 0), 1, 0)
            pp.close()     
            pp.contourf(result,[-0.5, 0.5, 1.5],\
                cmap=cm.get_cmap('bone'))
            # pp.axes().set_aspect('equal','box')
            pp.show()
            if os.path.exists('it_' + str(It) + '.pdf'):
                os.remove('it_' + str(It) + '.pdf')
            pp.savefig('./it_' + str(It) + '.pdf', bbox_inches='tight')

Topic		Replies	Views
Level set based domain definition in FEniCS (3D) mesh	2	449	September 5, 2023
Subdomains with different materials, solver not converging dolfinx	3	147	May 21, 2024
Problem when visualising modes with different materials General	1	236	November 13, 2023
Interpolation issue for topology optimisation dolfinx	2	427	June 8, 2022
Using different time steps in different subdomains variational formulation	2	503	August 25, 2022

Topology Optimisation using Level set Method

Related topics