Home

Welcome to the LFEM wiki!

How to solve a simple Topology Optimization problem using LFEM

The main idea of the L* packages is to make it easy to implement a topology optimization program. To show a basic example, lets solve the traditional

#
# Min C(x)
# S.t
#      V(x)≤V̄
#
#      0<xmin ≤ x ≤1

where C is the compliance and V the volume. The limit volume V̄ is set as a fraction of the total volume. This problem can be solved by using an Optimal Criteria (OC) as explained by Bensdsoe and Kikuchi in their seminal paper of 1988.

The main function is the OC

#
# Min C(x)
# S.t
#      V(x)≤V̄
#
#      0<xmin ≤ x ≤1
#
function OC(m::Mesh,vf::Float64,x::Vector,dC::Vector,dV::Vector)
    
    
    # Initial interval
    λ1 = 0.0
    λ2 = 1E5
    
    # Moving limits
    δ = 0.1
    
    # Tolerance (volume constraint)
    tol = 1E-6
    
    # Relaxation
    η = 0.5
    
    # Minimal value
    x_min = 1E-3
    
    # Full volume
    full_volume = Volume(m,ones(m.bmesh.ne))
    
    # Limit volume
    Vlimit = vf*full_volume
    
    
    # Working array
    x_estim = copy(x)
   
    # External loop - OC
    for k=1:1000  
       
       # Middle point 
       λ = (λ1 + λ2)/2
          
       # Loop along the elements
       for j=1:m.bmesh.ne 
              
           # Amplification factor 
           β = -min(0.0,dC[j])/(λ*dV[j]) 
        
           # New value for this variable
           x_e = x[j]*(β^η)
            
           # Moving limits
           x_dir = min(x[j] + δ, 1.0) 
           x_esq = max(x[j] - δ, x_min) 
           
           # 
           #-------0[----|--x--|----]1---------  
           #           x-δ    x+δ 
           #
            
           # Check moving limits 
           x_estim[j] = max( min( x_e, x_dir), x_esq)
            
            
       end #j
         
       # Current volume
       V = Volume(m,x_estim)
 
       # Test constraint
       if abs(λ2 - λ1)<=tol
            break
       end
        
       # Adjust λ to match the volume constraint
       if (V > Vlimit)
             λ1 = λ
       else
             λ2 = λ
       end
 
        
    end # k
 
    # Return the x_estim
    return x_estim
    
end

Next, let's define the functions and derivatives needed to solve the optimization problem

#
# Volume
#
function Volume(m::Mesh,x::Vector)
   
    volume = 0.0
    
    # Loop for all elements
    for j=1:m.bmesh.ne
        
        # Area
        geo = m.geo_ele[j]
        A = m.geometries[geo].A
    
        # Length
        L = Length(m.bmesh,j)  

        # Add to the total volume
        volume = volume + L*A*x[j]
        
    end #j
    
    return volume
    
end

function dVolume(m::Mesh,x::Vector)
    
    # Output 
    D = zeros(m.bmesh.ne)

    # Loop for all elements
    for j=1:m.bmesh.ne
    
        # Area
        geo = m.geo_ele[j]
        A = m.geometries[geo].A
        
        # Length
        L = Length(m.bmesh,j)  

        # Local derivative
        D[j] = L*A
        
    end #j
    
    return D
    
end

and

function dCompliance(m::Mesh,U::Vector,x::Vector,p::Float64)
    
    D = zeros(m.bmesh.ne)

    for ele=1:m.bmesh.ne

        # Dofs
        dofs = DOFs(m.bmesh,ele)
   
        # Element displacement in global reference
        u = U[dofs] 

        # Element stiffness in global reference
        if contains(string(m.bmesh.etype),"truss")
           Ke = Keg_truss(m,ele)
        else
            error("bla")   
        end

        D[ele]  =  -p*x[ele]^(p-1)*dot(u,Ke,u)

    end
   
    return D

end

Main program

using LinearAlgebra
using BMesh, LMesh, TMeshes, LFEM

function main(vf=0.5,output="topo.sol")

    # Load the problem
    m = Simply_supported2D(6,6)
    
    # Exponent (SIMP)
    p = 3.0

    # Initial point
    x = vf*ones(m.bmesh.ne)

    # Main loop
    for j=1:100

        # Equilibrium
        U,_ = Solve_linear(m;x=x,p=p)

        # Sensitivities
        dV = dVolume(m,x)
        dC = dCompliance(m,U,x,p)

        # OC
        xn = OC(m,vf,x,dC,dV)

        # Stop criteria
        if norm(xn.-x).<1E-6
            println("Done")
            x .= xn
            break
        end

        # Roll over Beethoven
        x .= xn 

    end

    # Write solution
    Gmsh_init(output,m)

    # Write topology
    Gmsh_element_scalar(m,x,output,"Topology")

    # As requested by Dr. Olavo, export the displacements
    U,_ = Solve_linear(m;x=x,p=p)
    Gmsh_nodal_vector(m,U,output,"Displacement")

end