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 Bensdso e 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