Skip to content
Snippets Groups Projects
multimaterial.jl 36.5 KiB
Newer Older
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020


#############################################multimaterial############################################################

#Based on: https://link.springer.com/article/10.1007/s00158-013-0999-1
function multitop(nx,ny,tol_out,tol_f,iter_max_in,iter_max_out,p,q,e,v,rf)
    anim=Animation()
    alpha = zeros(Float64,nx*ny,p)
    for i = 1:p
        alpha[:,i] .= v[i]
    end
    
    # MAKE FILTER
    H,Hs = make_filter(nx,ny,rf)
    obj=0
    change_out = 2*tol_out; iter_out = 0;
    while (iter_out < iter_max_out) && (change_out > tol_out)
        alpha_old = copy(alpha)
        for a = 1:p
            for b = a+1:p
                obj,alpha = bi_top(a,b,nx,ny,p,v,e,q,alpha,H,Hs,iter_max_in);
            end
        end
        iter_out = iter_out + 1;
        change_out = norm(alpha[:].-alpha_old[:],Inf);
        display("Iter:$iter_out Obj.:$obj change:$change_out");
        #  UPDATE FILTER
        if (change_out < tol_f) && (rf>3)
            tol_f = 0.99*tol_f; rf = 0.99*rf; 
            H,Hs = make_filter(nx,ny,rf);
        end
        #  SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
        if mod(iter_out,5)==0
            I = make_bitmap(p,nx,ny,alpha);
            display(RGB.(I[:,:,1],I[:,:,2],I[:,:,3]))
            I = permutedims(I, [3, 1, 2])
            img = colorview(RGB, I)
            heatmap(img,xaxis=nothing,yaxis=nothing,aspect_ratio=:equal)
            frame(anim)
        end
    end
    return anim,alpha
end

function bi_top(a,b,nx,ny,p,v,e,q,alpha_old,H,Hs,maxIter)
    alpha = copy(alpha_old)
    nu = 0.3;
    
    # PREPARE FINITE ELEMENT ANALYSIS
    A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
    A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
    B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
    B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
    KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
    
    nelx=nx
    nely=ny
    
    nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx)
    edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,nelx*nely,1)
    edofMat = repeat(edofVec,1,8).+repeat([0 1 2*nely.+[2 3 0 1] -2 -1],nelx*nely,1)
    iK = convert(Array{Int},reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1))
    jK = convert(Array{Int},reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1))
    
    # DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
    F = sparse([2],[1],[-1.0],2*(nely+1)*(nelx+1),1)
    U = zeros(2*(nely+1)*(nelx+1),1)
    fixeddofs = union(1:2:2*(nely+1),2*(nelx+1)*(nely+1))
    alldofs = 1:(2*(nely+1)*(nelx+1))
    freedofs = setdiff(alldofs,fixeddofs)
    
    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

        sK = reshape(KE[:]*E[:]',64*nelx*nely,1)
        K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
        @timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
    

        # Objective function and sensitivity analysis
        ce = sum((U[edofMat]*KE).*U[edofMat],dims=2)
        o = sum(sum(E.*ce))
        
        # FILTERING OF SENSITIVITIES 
        dc = 0.0 .-(q .*(e[a].-e[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.2;
        r = ones(nx*ny,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
            lmid = 0.5*(l2+l1)
            alpha_a = max.(l,min.(u,alpha[:,a].*sqrt.(-dc./lmid)));
            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

########

function multitop_compliant(nx,ny,tol_out,tol_f,iter_max_in,iter_max_out,p,q,e,v,rf,Load,Support,Spring,DOut)
    anim=Animation()
    alpha = zeros(Float64,nx*ny,p)
    for i = 1:p
        alpha[:,i] .= v[i]
    end
    
    # MAKE FILTER
    H,Hs = make_filter(nx,ny,rf)
    obj=0
    change_out = 2*tol_out; iter_out = 0;
    while (iter_out < iter_max_out) && (change_out > tol_out)
        alpha_old = copy(alpha)
        for a = 1:p
            for b = a+1:p
                obj,alpha = bi_top_compliant(a,b,nx,ny,p,v,e,q,alpha,H,Hs,iter_max_in);
            end
        end
        iter_out = iter_out + 1;
        change_out = norm(alpha[:].-alpha_old[:],Inf);
        display("Iter:$iter_out Obj.:$obj change:$change_out");
        #  UPDATE FILTER
        if (change_out < tol_f) && (rf>3)
            tol_f = 0.99*tol_f; rf = 0.99*rf; 
            H,Hs = make_filter(nx,ny,rf);
        end
        #  SCREEN OUT TEMPORAL TOPOLOGY EVERY 5 ITERATIONS
        if mod(iter_out,5)==0
            I = make_bitmap_compliant(p,nx,ny,alpha);
            display(RGB.(I[:,:,1],I[:,:,2],I[:,:,3]))
            I = permutedims(I, [3, 1, 2])
            img = colorview(RGB, I)
            heatmap(img,xaxis=nothing,yaxis=nothing,aspect_ratio=:equal)
            frame(anim)
        end
    end
    return anim,alpha
end

function bi_top_compliant(a,b,nx,ny,p,v,e,q,alpha_old,H,Hs,maxIter,Load,Support,Spring,DOut)
    alpha = copy(alpha_old)
    nu = 0.3;
    
    # PREPARE FINITE ELEMENT ANALYSIS
    A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
    A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
    B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
    B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
    KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
    
    nelx=nx
    nely=ny
    
    nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx)
    edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,nelx*nely,1)
    edofMat = repeat(edofVec,1,8).+repeat([0 1 2*nely.+[2 3 0 1] -2 -1],nelx*nely,1)
    iK = convert(Array{Int},reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1))
    jK = convert(Array{Int},reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1))

    

    # # DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
    # F = sparse([2],[1],[-1.0],2*(nely+1)*(nelx+1),1)
    # U = zeros(2*(nely+1)*(nelx+1),1)
    # fixeddofs = union(1:2:2*(nely+1),2*(nelx+1)*(nely+1))
    # alldofs = 1:(2*(nely+1)*(nelx+1))
    # freedofs = setdiff(alldofs,fixeddofs)

    # DEFINE LOADS AND SUPPORTS 
    F = sparse(2 .*Load[:,1] .-2 .+ Load[:,2],ones(size(Load,1)),Load[:,3],2*(nely+1)*(nelx+1),2)
    DofDOut = 2 * DOut[1] - 2 +DOut[2]; #only one
Loading
Loading full blame...