Adding filter and projection to our simple Topology Optimization program.

Recalling the previous example. If we intend to solve the same problem using elasticity (:solid2D or :solid3D elements), then we may face some common problems in topology optimization: checkerboard and mesh-dependency.

A common approach to solve both problems is to use spatial filter (with or without projection). This functionality is provided by package LFilter (in this repository). We only have to change the main function of our previous program

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

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

    # Load the problem
    mesh = Simply_supported2D(50,50,:solid2D)

    # It just makes sense to filter if not using truss
    Get_eclass(mesh)==:solid || throw("main:: this example is for :solid2D or :solid3D elements")
    
    # Filter
    filter = Filter(mesh,3*1.0/100)

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

    # Number of elements
    ne = Get_ne(mesh)

    # Initial point
    x = vf*ones(ne)

    # Main loop
    @showprogress "Optimizing..." for j=1:100

        # Filter 
        ρ = x2proj(x,filter)

        # Equilibrium
        U,_ = Solve_linear(mesh,ρ,kparam)

        # Sensitivities
        dV = dVolume(mesh,ρ)
        dC = dCompliance(mesh,U,ρ,dkparam)

        # Convert Sensitivities back to x
        dproj2dx!(x,dV,filter)
        dproj2dx!(x,dC,filter)
        

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

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

        # Roll over Beethoven
        x .= xn 

    end

    # Filter 
    ρ = x2proj(x,filter)

    # Write solution
    Gmsh_init(output,mesh)

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

    # As requested by Dr. Olavo, export the displacements
    U,_ = Solve_linear(mesh,ρ,kparam)
    Gmsh_nodal_vector(mesh,U,output,"Displacement")

end

The main modifications with respect to previous example are:

1. Load the new package: using LFilter

2. Use solid elements in our mesh: mesh = Simply_supported2D(100,100,:solid2D)

3. Create the filter with radius 0.03: filter = Filter(mesh,3*1.0/100)

4. Filter field x to obtain field ρ

5. Use field ρ to evaluate equilibrium and sensitivities

6. Fix sensitivities back to the original field: dproj2dx!(x,dV,filter) and dproj2dx!(x,dC,filter)

Putting all together

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

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

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

    # Load the problem
    mesh = Simply_supported2D(50,50,:solid2D)

    # It just makes sense to filter if not using truss
    Get_eclass(mesh)==:solid || throw("main:: this example is for :solid2D or :solid3D elements")
    
    # Filter
    filter = Filter(mesh,3*1.0/100)

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

    # Number of elements
    ne = Get_ne(mesh)

    # Initial point
    x = vf*ones(ne)

    # Main loop
    @showprogress "Optimizing..." for j=1:100

        # Filter 
        ρ = x2proj(x,filter)

        # Equilibrium
        U,_ = Solve_linear(mesh,ρ,kparam)

        # Sensitivities
        dV = dVolume(mesh,ρ)
        dC = dCompliance(mesh,U,ρ,dkparam)

        # Convert Sensitivities back to x
        dproj2dx!(x,dV,filter)
        dproj2dx!(x,dC,filter)
        

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

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

        # Roll over Beethoven
        x .= xn 

    end

    # Filter 
    ρ = x2proj(x,filter)

    # Write solution
    Gmsh_init(output,mesh)

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

    # As requested by Dr. Olavo, export the displacements
    U,_ = Solve_linear(mesh,ρ,kparam)
    Gmsh_nodal_vector(mesh,U,output,"Displacement")

end

function OC(mesh::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
    
    # Number of elements
    ne = Get_ne(mesh)

    # Full volume
    full_volume = Volume(mesh,ones(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 in mesh
              
           # 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(mesh,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

#
# Volume
#
function Volume(mesh::Mesh,x::Vector)
   
    volume = 0.0
    
    # Loop for all elements
    for ele in mesh
        
        # Volume of this element
        vol = Volume_element(mesh,ele)

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

function dVolume(mesh::Mesh,x::Vector)
    
    # Number of elements
    ne = Get_ne(mesh)

    # Output 
    D = zeros(ne)

    # Loop for all elements
    for ele in mesh
    
        # Local derivative
        D[ele] = Volume_element(mesh,ele)

        
    end #ele
    
    return D
    
end


function dCompliance(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 in global reference
        ug = U[dofs] 

        # If needed, rotate to local reference
        ul = To_local(ug,mesh,ele)
        
        # Element stiffness in Local reference
        Ke = Local_K(mesh,ele)

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

    end
   
    return D

end