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
using BMesh, LMesh, TMeshes, LFEM
using LFilter

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

    # Load the problem
    m = Simply_supported2D(100,100,:solid2D)
    
    # Filter
    filter = Filter(m,3*1.0/100)

    # Exponent (SIMP)
    p = 3.0

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

    # Main loop
    for j=1:100

        # Filter 
        ρ = x2proj(x,filter)

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

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

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

        # 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

    # Filter 
    ρ = x2proj(x,filter)

    # Write solution
    Gmsh_init(output,m)

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

    # As requested by Dr. Olavo, export the displacements
    U,_ = Solve_linear(m;x=ρ,p=p)
    Gmsh_nodal_vector(m,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: m = Simply_supported2D(100,100,:solid2D)

3. Create the filter with radius 0.03: filter = Filter(m,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)