Using a BESO like algorithm

The manuscript Zhi Hao Zuo, Yi Min Xie, A simple and compact Python code for complex 3D topology optimization, 2015, proposes a very simple algorithm to account for discrete values of relative density. The complete computer code using LFEM and LFilter is

'''julia

Min C(x)

S.t

V(x)≤V̄

x = xmin OR x = 1

using LinearAlgebra, ProgressMeter using BMesh, LMesh, TMeshes, LFEM using LFilter

vf is the volume fraction

function main(vf=0.3,output="topo_beso.sol")

Load the problem
nx = 40*2
ny = 25*2
Lx = 8.0
Ly = 5.0
mesh = Cantilever_beam_bottom2D(nx,ny,:solid2D)

# It just makes sense to filter if using elasticity 
Get_eclass(mesh)==:solid || throw("main:: this example is for :solid2D or :solid3D elements")

# Element size
te = max(Lx/nx,Ly/ny)

# Filter
filter = Filter(mesh,3*te)

# Material parametrization and derivative
kparam(xe::Float64,p=3.0)=xe^p
dkparam(xe::Float64,p=3.0)=p*xe^(p-1)

# BESO parameters

# Evolution rate
er = 0.01

# Number of iterations
niter = 150

# Number of elements
ne = Get_ne(mesh)

# Initial point
ρ = ones(ne)

# Volume of each element
V = Volumes(mesh,ones(ne))

# Total volume
volume_full = sum(V)

# Target volume
Vast = vf*volume_full

# Sensitivity index in the last iteration
ESED_F_ANT = zeros(ne)

# Main loop
 for j=1:niter

    # Actual volume
    volume_actual = sum(ρ.*V)

    # Volume for this iteration
    if volume_actual > Vast
       vol = max(Vast, volume_actual*(1.0 - er))
    elseif volume_actual < Vast
       vol = min(Vast,volume_actual*(1 + er))
    else
        vol = Vast
    end
    
    # Equilibrium
    U,F,_ = Solve_linear(mesh,ρ,kparam)

    # Compliance sensitivities
    dC = Sensitivity(mesh,U,ρ,dkparam)

    # ESED
    SN = dC ./ V
    
    # Filter the ESEDs
    ESED_F =  LFilter.x2ρ(SN,filter)

    # Mean value using the last iteration
    ESED_F_media = (ESED_F .+ ESED_F_ANT)./2

    # Store the value for the next iteration
    ESED_F_ANT .= ESED_F_media 

    # Update the relative densities
    ρ = BESO(ρ, ESED_F_media, vol, V)
    

end #j 

# Write solution
Gmsh_init(output,mesh)

# Write topology
Gmsh_element_scalar(mesh,ρ,output,"Topology")

end

Volume

function Volumes(mesh::Mesh,x::Vector)

# Number of elements
ne = Get_ne(mesh)

# Output vector
V = zeros(ne)

# Loop for all elements
for ele in mesh
    
    # Volume of this element
    vol = Volume_element(mesh,ele)

    # Add to the total volume
    V[ele] =  vol*x[ele]
    
end #ele

return V

end

Sensitivities - Compliance

function Sensitivity(mesh::Mesh,U::Vector,x::Vector,dkparam::Function)

# Number of elements
ne = Get_ne(mesh)

# Output vector
D = zeros(ne)

# Loop over elements in the mesh
for ele in mesh

    # Dofs
    dofs = DOFs(mesh,ele)

    # Element displacement 
    ul = U[dofs] 

    # Element stiffness 
    Ke = Local_K(mesh,ele)

    D[ele]  =  dkparam(x[ele])*dot(ul,Ke,ul)

end

return D

end

BESO as in Zhi Hao Zuo, Yi Min Xie, A simple and compact Python code for complex 3D topology optimization, 2015

function BESO(x, D, Vlim, V, tol=1E-6, xmin=1E-3, xmax=1.0)

# Limit values
D_min = minimum(D)
D_max = maximum(D)

# Make a copy
xn = copy(x)

# Main loop (bisection)
while (D_max - D_min)/D_max > tol

    # Threshold value
    corte = (D_min + D_max)/2 

    # Initialize the volume
    volume = 0.0

    # For each variable
    for i in LinearIndices(x)

        # Round down or up, relative to the threshold
        if D[i] < corte

            xn[i] = xmin  

        elseif D[i] > corte

            xn[i] = xmax 

        end #if

        # Add the volume
        volume += xn[i]*V[i]
    
    end #i

    # Bisection
    if volume > Vlim
        D_min = corte
    else
        D_max = corte
    end

end # while

# Return the new distribution of design variables
return xn

end

'''