microstructure_topX.jl

    
        
    
    for j = 1:3
        ufixed[3:4,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[nelx;0];
        ufixed[7:8,j] = [e0[1,j] e0[3,j]/2 ; e0[3,j]/2 e0[2,j]]*[0;nely];
        ufixed[5:6,j] = ufixed[3:4,j] .+ ufixed[7:8,j];
    end
    wfixed = [repeat(ufixed[3:4,:],nely-1,1); repeat(ufixed[7:8,:],nelx-1,1)];
    

    
    ## INITIALIZE ITERATION
    qe = Array{Any,2}(undef, 3, 3);
    Q = zeros(3,3);
    dQ = Array{Any,2}(undef, 3, 3);
    x = volfrac.*ones(nely,nelx)
    
    for i = 1:nelx
        for j = 1:nely
            vall=3
            if optType=="poisson"
                vall=6
            end
            if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2) < min(nelx,nely)/vall
                x[j,i] = volfrac/2.0;
            end
        end
    end
    xPhys = copy(x);
    
    #display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
    loop=0
    
    function FE(xPhys)
        ## FE-ANALYSIS
        sK = reshape(KE[:]*(Emin .+ xPhys[:]'.^penal*(80 .- Emin)),64*nelx*nely,1);
        
        K = sparse(vec(iK),vec(jK),vec(sK)); 
        K = (K.+K')./2.0;
        Kr = vcat(hcat(K[d2,d2] , K[d2,d3]+K[d2,d4]),hcat((K[d3,d2]+K[d4,d2]),(K[d3,d3]+K[d4,d3]+K[d3,d4]+K[d4,d4])));
        U[d1,:] .= ufixed;
        U[[d2;d3],:] = Kr\(-[K[d2,d1]; K[d3,d1]+K[d4,d1]]*ufixed-[K[d2,d4]; K[d3,d4]+K[d4,d4]]*wfixed);
        U[d4,:] = U[d3,:]+wfixed;
        
        
        ## OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
        for i = 1:3
            for j = 1:3
                U1 = U[:,i]; U2 = U[:,j];
                qe[i,j] = reshape(sum((U1[edofMat]*KE).*U2[edofMat],dims=2),nely,nelx)./(nelx*nely);
                Q[i,j] = sum(sum((Emin .+ xPhys.^penal*(E0 .-Emin)).*qe[i,j]));
                dQ[i,j] = penal*(E0-Emin)*xPhys.^(penal-1).*qe[i,j];
            end
        end
    end
    
    # function minimize volume
    function F(x1::Vector, grad::Vector)
        
        xPhys=reshape(x1,nely,nelx)
        FE(xPhys)

        
        dv = ones(nely,nelx)
        if ft == 2
            dv[:] = H*(dv[:]./Hs)
        end

        grad[:] .=  dv[:]

        #return (sum(x1) - volfrac*nel)
        return (sum(x1))
    end

    function G33(x1::Vector, grad::Vector)

        #shear
        c=-Q[3,3];
        dc=-dQ[3,3];
            
          
        ## FILTERING/MODIFICATION OF SENSITIVITIES
        if ft == 1
            dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
        elseif ft == 2
            dc[:] = H*(dc[:]./Hs);
        end
        
        grad[:] .= dc[:]
        
        return c +CStar[3,3]
    end
    
    function G11(x1::Vector, grad::Vector)

        #shear
        c=-(Q[1,1]);
        dc=-(dQ[1,1]);
            
          
        ## FILTERING/MODIFICATION OF SENSITIVITIES
        if ft == 1
            dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
        elseif ft == 2
            dc[:] = H*(dc[:]./Hs);
        end
        
        grad[:] .= dc[:]
        
        return c +CStar[1,1]
    end
    
    function G12(x1::Vector, grad::Vector)

        #shear
        c=-(Q[1,2]);
        dc=-(dQ[1,2]);
            
          
        ## FILTERING/MODIFICATION OF SENSITIVITIES
        if ft == 1
            dc[:] = H*(x[:].*dc[:])./Hs./max.(1e-3,x[:]);
        elseif ft == 2
            dc[:] = H*(dc[:]./Hs);
        end
        
        grad[:] .= dc[:]
        
        return c +CStar[1,2]
    end

    


    opt = Opt(:LD_MMA, nel)
    opt.lower_bounds = fill(1e-9,nel)
    opt.upper_bounds = fill(1,nel)
    opt.xtol_rel = 1e-6
    opt.maxeval = maxEval

    opt.min_objective = F
    

    inequality_constraint!(opt, (x,gg) -> G33(x,gg), 1e-3)
    # inequality_constraint!(opt, (x,gg) -> G11(x,gg), 1e-3)
    inequality_constraint!(opt, (x,gg) -> G12(x,gg), 1e-6)

    display(@time (minf,minx,ret) = optimize(opt, xPhys[:]))
    numevals = opt.numevals # the number of function evaluations
    display("got $minf after $numevals iterations (returned $ret)")
    
    xPhys=reshape(minx,nely,nelx)
    display(heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing, foreground_color_subplot=colorant"white",fc=:grays))
    
    display(Q)
    
    return xPhys
end