multimaterial.jl

            I = make_bitmap_3d(p,nx,ny,nz,alpha)
            threshold=0.05
            display(topplotmulti3d(nx,ny,nz,I,threshold))
            frame(anim)
        end
    end
    return anim,alpha
end

function bi_top3D_compliantKE(a,b,nx,ny,nz,p,v,KEs,Es,q,alpha_old,H,Hs,maxIter,getProblem=inverter)
    alpha = copy(alpha_old)
    nu = 0.3;

    
    
    # PREPARE FINITE ELEMENT ANALYSIS
    # KE=lk_H8(nu)
    
    nelx=nx
    nely=ny
    nelz=nz
    nele = nelx*nely*nelz
    ndof = 3*(nelx+1)*(nely+1)*(nelz+1)

    nodenrs = reshape(1:(1+nelx)*(1+nely)*(1+nelz),1+nely,1+nelx,1+nelz)
    edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1].+1,nelx*nely*nelz,1)
    edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nely+1)*(nelx+1).+[0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1]],nelx*nely*nelz,1)
        
    iK = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
    jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
    
    # # USER-DEFINED LOAD DOFs
    # il = [0 nelx]; jl = [nely nely]; kl = [0 0];
    # loadnid = kl .*(nelx+1) .*(nely+1) .+il .*(nely .+1) .+(nely+1 .-jl);
    # loaddof = 3 .*loadnid[:] .- 2;
    # din = loaddof[1]; dout = loaddof[2];
    # F = sparse(loaddof,[1;2],[1;-1],ndof,2);
    
    # # USER-DEFINED SUPPORT FIXED DOFs
    # # Top face
    # m = Matlab.meshgrid(0:nelx,0:nelz)
    # iif=m[1]
    # kf=m[2]
    # fixednid = kf .*(nelx+1) .*(nely+1) .+iif .*(nely+1) .+1;
    # fixeddof_t = 3 .*fixednid[:] .-1;
    # # Side face
    # m = Matlab.meshgrid(0:nelx,0:nely)
    # iif=m[1]
    # jf=m[2]
    # fixednid = iif .*(nely+1) .+(nely+1 .-jf);
    # fixeddof_s = 3*fixednid[:];
    # # Two pins
    # iif = [0;0]; jf = [0;1]; kf = [nelz;nelz];
    # fixednid = kf .*(nelx+1) .*(nely .+1) .+iif .*(nely .+1) .+(nely+1 .-jf)
    # fixeddof_p = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2] # DOFs
    # # Fixed DOFs
    # fixeddof = union(fixeddof_t,fixeddof_s);
    # fixeddof = union(fixeddof, fixeddof_p);
    # fixeddof = sort(fixeddof);
    
    # # Displacement vector
    # U = zeros(ndof,2);
    # freedofs = setdiff(1:ndof,fixeddof)

    U,F,freedofs,din,dout=getProblem(nelx,nely,nelz)

    
    
    o=0;alpha_a=0
    # inner iteration
    for i =0: (maxIter-1)
        # FE-ANALYSIS
        # E = e[1]*alpha[:,1].^q
        # for phase = 2:p
        #     E = E + e[phase]*alpha[:,phase].^q;
        # end

        KEE = KEs[1][:]*alpha[:,1][:]'.^q
        for phase = 2:p
            KEE = KEE + KEs[phase][:]*alpha[:,phase][:]'.^q
        end

        # EE = mean(Es[1])*alpha[:,1].^q
        # for phase = 2:p
        #     EE = EE + mean(Es[phase])*alpha[:,phase].^q;
        # end

        

        sK = reshape(KEE[:],24*24*nelx*nely*nelz,1)
        

        # sK = reshape(KE[:]*E[:]',24*24*nelx*nely*nelz,1)
        K = sparse(vec(iK),vec(jK),vec(sK));
        K[din,din] = K[din,din] + 0.1;
        K[dout,dout] = K[dout,dout] + 0.1;
        K = (K+K')/2;

        # @timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
        @timed U[freedofs,:] = K[freedofs,freedofs] \ Array(F[freedofs,:])
        U1 = U[:,1]
        U2 = U[:,2]

        KE=reshape(sum(KEE,dims=2),24,24)
        # ce=zeros(nelx* nely* nelz)
        # for i=1:Int(nelx* nely* nelz)
        #     KE = KEs[1].*alpha[i,1]
        #     for phase = 2:p
        #         KE = KE + KEs[phase].*alpha[i,phase]
        #     end
        #     ce[i]=-sum((U1[edofMat][i,:]'*KE).*U2[edofMat][i,:]')
        # end

        ce=zeros(nelx* nely* nelz)
        for i=1:Int(nelx* nely* nelz)
            KE = KEs[1].*alpha[i,1]
            for phase = 2:p
                KE = KE + KEs[phase].*alpha[i,phase]
            end
            ce[i]=-sum((U1[edofMat][i,:]'*KE).*U2[edofMat][i,:]')
        end
    

        # Objective function and sensitivity analysis
        # ce = -sum((U1[edofMat]*KE).*U2[edofMat],dims=2)

        # o = sum(sum(E.*ce))
        o  = U[dout,1]
        
        # FILTERING OF SENSITIVITIES 
        # dc = 0.0 .-(q .*(e[a].-e[b]).*alpha[:,a].^(q-1)).*ce;
        dc = 0.0 .-(q .*mean(Es[a].-Es[b]).*alpha[:,a].^(q-1)).*ce;
        dc = H*(alpha[:,a].*dc)./Hs./max.(1e-3,alpha[:,a])
        dc = min.(dc,0.0)

        # UPDATE LOWER AND UPPER BOUNDS OF DESIGN VARIABLES
        move = 0.1;
        r = ones(nx*ny*nz,1)
  
        for k = 1:p
            #if (k ~= a) && (k ~= b)
            if (k != a) && (k != b)
                # display("k $k a $a b $b")
                r = r .- alpha[:,k]
            end
        end
        l = max.(0,alpha[:,a] .-move)
        u = min.(r,alpha[:,a] .+move)
        # OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES
        l1 = 0; l2 = 1e9
        # while (l2-l1)/(l1+l2) > 1e-3
        while (l2-l1)/(l1+l2) > 1e-4 && l2 > 1e-40
            lmid = 0.5*(l2+l1)
            # alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
            alpha_a = max.(l,min.(u,alpha[:,a].*((max.(1e-10,-dc./lmid)).^0.3 )));
            if sum(alpha_a) > nx*ny*v[a]
                l1 = lmid
            else
                l2 = lmid 
            end
        end
        alpha[:,a] = alpha_a
        alpha[:,b] = r .-alpha_a

    end
    return o,alpha
end