multiscale.jl

################################################################################
# A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, OCTOBER 1999
# MODIFIED FOR 3D MULTISCALE DESIGN VIA SURROGATE MODEL, LLNL, JULY 2018
#
# This work was produced under the auspices of the U.S. Department of Energy by
# Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
#
# This work was prepared as an account of work sponsored by an agency of the
# United States Government. Neither the United States Government nor Lawrence
# Livermore National Security, LLC, nor any of their employees makes any warranty,
# expressed or implied, or assumes any legal liability or responsibility for the
# accuracy, completeness, or usefulness of any information, apparatus, product, or
# process disclosed, or represents that its use would not infringe privately owned
# rights. Reference herein to any specific commercial product, process, or service
# by trade name, trademark, manufacturer, or otherwise does not necessarily
# constitute or imply its endorsement, recommendation, or favoring by the United
# States Government or Lawrence Livermore National Security, LLC. The views and
# opinions of authors expressed herein do not necessarily state or reflect those
# of the United States Government or Lawrence Livermore National Security, LLC,
# and shall not be used for advertising or product endorsement purposes.
#
# LLNL-CODE-757968
################################################################################

function multiscale(nelx, nely, nelz, volfrac, rmin, truss, Es, vs, minVF, maxVF, maxit)
    anim=Animation()

    # INITIALIZE
    x=zeros(nelx,nely,nelz)
    x[1:nelx, 1:nely, 1:nelz] .= volfrac;
    Gs = Es / (2*(1+vs));
    loop = 0; change = 1.0;
    nnx = nelx+1; nny = nely+1; nnz = nelz+1;
    # colormap(gray); caxis([0.0, 1.0]);
    # maxit=1
    # START ITERATION
    while change > 0.01 && loop < maxit
        loop = loop + 1;
        xold = x;
        # FE-ANALYSIS
        U = FE(nelx, nely, nelz, nnx, nny, nnz, x, truss, Es, vs, Gs);
        # display(U)
        # OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
        c = 0.0; 
        dc = zeros(nelx, nely, nelz);
        for elz = 1:nelz; 
            for ely = 1:nely 
                for elx = 1:nelx
                    KE   = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, false);
                    DKE  = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, true);
                    dofs = get_elem_dofs(nnx, nny, nnz, elx, ely, elz);
                    Ue = U[dofs];
                    c = c + Ue'*KE*Ue;
                    dc[elx,ely,elz] = -Ue'*DKE*Ue;
                end
            end
        end
        # FILTERING OF SENSITIVITIES
        dc   = check(nelx, nely, nelz, rmin, x, dc);
        # DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
        x    = OC(nelx, nely, nelz, x, volfrac, dc, minVF, maxVF);
        # PRINT RESULTS
        change = maximum((abs.(x-xold)));
        println(" It:$loop Obj:$c Vol:$(mean(x)) ch:$change ")
        
        if mod(loop,5)==0
            display(topplot3d(x,nelx, nely, nelz))
            topplot3d(x,nelx, nely, nelz)
            frame(anim)
        end

        # disp([" It.: " sprintf("#4i",loop) " Obj.: " sprintf("#10.4f",c)  " Vol.: " sprintf("#6.3f",sum(sum(sum(x)))/(nelx*nely*nelz)) " ch.: " sprintf("#6.3f",change )])
        # PLOT DENSITIES
        # viz3d(nelx, nely, nelz, x, volfrac, nelx==1);
        # SAVE PARAMETER VALUES (ELEMENT DENSITIES AND ROD DIAMETERS)
        # xOut = x[:]
        # # reshape(x,[],1);             
        # # save("-ascii","elVolFrac.txt", "xOut");
        # dOut = get_d(truss,x)[:]
        # # reshape(get_d(truss,x),[],1); 
        # # save("-ascii","elRodDiam.txt","dOut");
    end
    return x,anim
end

######### OPTIMALITY CRITERIA UPDATE ###########################################
function OC(nelx, nely, nelz, x, volfrac, dc, minVF, maxVF)
    l1 = 0; l2 = 100000; move = 0.2;
    xnew=zeros(size(x))
    while (l2-l1 > 1e-4)
        lmid = 0.5*(l2 + l1);
        xnew = max.(minVF, max.(x .- move, min.(maxVF, min.(x .+move,x.*sqrt.(abs.(0.0 .-dc./lmid))))));
        # xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.((0.0.-dc)./dv./lmid)))))

        if sum(sum(sum(xnew))) - volfrac*nelx*nely*nelz > 0;
            l1 = lmid;
        else
            l2 = lmid;
        end
    end
    return xnew
end

######### MESH-INDEPENDENCY FILTER #############################################
function check(nelx, nely, nelz, rmin, x, dc)
    dcn=zeros(size(dc));
    for elz = 1:nelz; 
        for ely = 1:nely; 
            for elx = 1:nelx
                sum = 0.0;
                for k = max(elz-round(rmin),1):min(elz+round(rmin),nelz)
                    for j = max(ely-round(rmin),1):min(ely+round(rmin),nely)
                        for i = max(elx-round(rmin),1):min(elx+round(rmin),nelx)
                            fac = rmin - sqrt((elx-i)^2+(ely-j)^2+(elz-k)^2);
                            sum = sum + max(0,fac);
                            dcn[elx,ely,elz] = dcn[elx,ely,elz] + max(0,fac)*x[Int(i),Int(j),Int(k)]*dc[Int(i),Int(j),Int(k)];
                        end
                    end
                end
                dcn[elx,ely,elz] = dcn[elx,ely,elz] / (x[elx,ely,elz]*sum);
            end; 
        end; 
    end
    return dcn
end

######### FE-ANALYSIS ##########################################################
function FE(nelx, nely, nelz, nnx, nny, nnz, x, truss, Es, vs, Gs)
    K = spzeros(3*nnx*nny*nnz, 3*nnx*nny*nnz);
    F = spzeros(3*nnx*nny*nnz);
    U = spzeros(3*nnx*nny*nnz);
    for elz = 1:nelz; 
        for ely = 1:nely; 
            for elx = 1:nelx
                KE   = get_KE(truss, x, Es, vs, Gs, elx, ely, elz, false);
                dofs = get_elem_dofs(nnx, nny, nnz, elx, ely, elz);
                # display(size(KE))
                # display(size(dofs))
                # display(size(K[dofs,dofs]))
                K[dofs,dofs] .= K[dofs,dofs] .+ KE;
            end; 
        end; 
    end
    # DEFINE LOADS AND SUPPORTS(HALF MBB-BEAM)
    coords = zeros(nnx*nny*nnz,3);
    n = 0;
    for k = 1:nnz; 
        for j = 1:nny; 
            for i = 1:nnx
                n = n+1; 
                coords[n,1] = i-1; 
                coords[n,2] = j-1; 
                coords[n,3] = k-1;
            end
        end 
    end
    midplane_nodes = findall(x->x==0, coords[:,2]);
    loaded_nodes   = intersect(findall(x->x==nelz, coords[:,3]), findall(x->x==0, coords[:,2]));
    fixed_nodes    = intersect(findall(x->x==0, coords[:,3]), findall(x->x==nely, coords[:,2]));
    fixeddofs      = zeros(Int,size(midplane_nodes,1) + 2*size(fixed_nodes,1));
    for i = loaded_nodes'
        F[3*(i-1)+3] = -1.0/nnx; 
    end
    n = 1;
    for i = midplane_nodes'
        for j=[2]
            fixeddofs[n] = 3*(i-1)+j; 
            n =n+1; 
        end 
    end
    for i = fixed_nodes' 
        for j=[1,3] 
            fixeddofs[n] = 3*(i-1)+j; 
            n =n+1; 
        end
    end
    alldofs   = 1:3*nnx*nny*nnz;
    freedofs  = setdiff(alldofs,fixeddofs);
    # display(K)

    # SOLVING
    U[freedofs]  .= K[freedofs,freedofs] \ Array(F[freedofs]);
    U[fixeddofs] .= 0;
    return U
end

######### ELEMENT AND NODE NUMBERING IN 3D MESH ################################
function get_num(nx, ny, nz, i, j, k)
    num = (nx*ny)*(k-1) + nx*(j-1) + i;
    return num
end

######### GLOBAL DOFS FOR A GIVEN ELEMENT ######################################
function  get_elem_dofs(nnx, nny, nnz, elx, ely, elz)
    n = get_num(nnx, nny, nnz, elx, ely, elz);
    N = [n; n+1; n+nnx+1; n+nnx; n+nnx*nny; n+nnx*nny+1; n+nnx*nny+nnx+1; n+nnx*nny+nnx];
    dofs = zeros(Int,24); 
    for j = 1:8; 
        for i = 1:3; 
            dofs[3*(j-1)+i] = Int(3*(N[j]-1)+i); 
        end
    end
    return dofs
end

######### INTEGRATE ELASTICITY TENSOR CE TO GET KE #############################
function get_KE(truss, x, Es, vs, Gs, i, j, k, deriv)
    KE = zeros(24,24);
    CE = get_CE(truss, x, Es, vs, Gs, i, j, k, deriv);
    
    for l = 1:8
        r = (sqrt(3)/3) * (-1 + 2*(sum([2,3,6,7].==l)>0)); rp = (1+r); rm = (1-r);
        s = (sqrt(3)/3) * (-1 + 2*(sum([3,4,7,8].==l)>0)); sp = (1+s); sm = (1-s);
        t = (sqrt(3)/3) * (-1 + 2*(sum([5,6,7,8].==l)>0)); tp = (1+t); tm = (1-t);
        DN = [ -sm*tm -rm*tm -rm*sm;
                sm*tm -rp*tm -rp*sm;
                sp*tm  rp*tm -rp*sp;
                -sp*tm  rm*tm -rm*sp;
                -sm*tp -rm*tp  rm*sm;
                sm*tp -rp*tp  rp*sm;
                sp*tp  rp*tp  rp*sp;
                -sp*tp rm*tp  rm*sp] ./ 8;
        B = DN * 2*Matrix(1.0I, 3, 3); 
        G = kron(B', Matrix(1.0I, 3, 3)); 
        KE = KE + G' * CE * G / 8;
    end
    # display(KE)
    return KE
end

######### DEFINE ELASTICITY TENSOR FOR DIFFERENT TRUSSES #######################
function get_CE(truss, x, Es, vs, Gs, i, j, k, D)
    p = x[i, j, k];
    if truss== "iso"       
        TM = iso_moduli
    elseif truss== "octet"  
        TM = octet_moduli
    elseif truss== "orc"  
        TM = orc_moduli
    elseif truss== "bound"  
        TM = bound_moduli
    else                          
        TM = simp_moduli
    end

    E, v, G = TM(p,Es,vs,Gs,false);  
    if D
        DE, Dv, DG = TM(p,Es,vs,Gs,true); 
    end
    if D == 0
        C1111 = E * (1.0 - v) / (1.0 - v - 2*v^2);
        C1122 = (E * v) / (1.0 - v - 2*v^2);
        C1212 = G;
    else # return the deriviatives instead
        C1111 = ((DE*(1-v)-E*Dv)*(1-v-2*v^2)-E*(1-v)*(-Dv-4*v*Dv)) / (1-v-2*v^2)^2;
        C1122 = ((DE*v+E*Dv)*(1-v-2*v^2)-E*v*(-Dv-4*v*Dv)) / (1-v-2*v^2)^2;
        C1212 = DG;
    end
    CE = [C1111   0     0     0   C1122   0     0     0   C1122;
            0   C1212   0   C1212   0     0     0     0     0  ;
            0     0   C1212   0     0     0   C1212   0     0  ;
            0   C1212   0   C1212   0     0     0     0     0  ;
        C1122   0     0     0   C1111   0     0     0   C1122;
            0     0     0     0     0   C1212   0   C1212   0  ;
            0     0   C1212   0     0     0   C1212   0     0  ;
            0     0     0     0     0   C1212   0   C1212   0  ;
        C1122   0     0     0   C1122   0     0     0   C1111];
    return CE
end

######### TRUSS-SPECIFIC MECHANICS MODELS ######################################
function iso_moduli(p, Es, vs, Gs, deriv)
    E = Es * (( 2.05292e-01 - 3.30265e-02*vs) * (p^(1-deriv)) * (1+0*deriv) + 
            ( 8.12145e-02 + 2.72431e-01*vs) * (p^(2-deriv)) * (1+1*deriv) + 
            ( 6.49737e-01 - 2.42374e-01*vs) * (p^(3-deriv)) * (1+2*deriv));
    v =      ( 2.47760e-01 + 1.69804e-02*vs) * (1-deriv) +  
        (-1.59293e-01 + 7.38598e-01*vs) * (p^(1-deriv)) * (1+0*deriv) + 
        (-1.86279e-01 - 4.83229e-01*vs) * (p^(2-deriv)) * (1+1*deriv) + 
        ( 9.77457e-02 + 7.26595e-01*vs) * (p^(3-deriv)) * (1+2*deriv);
    G = Gs * (( 1.63200e-01 + 1.27910e-01*vs) * (p^(1-deriv)) * (1+0*deriv) + 
            ( 6.00810e-03 + 4.13331e-01*vs) * (p^(2-deriv)) * (1+1*deriv) + 
                ( 7.22847e-01 - 3.56032e-01*vs) * (p^(3-deriv)) * (1+2*deriv));
    return E,v,G
end

function octet_moduli(p, Es, vs, Gs, deriv)
    E = Es * (( 1.36265e-01 - 1.22204e-02*vs) * (p^(1-deriv)) * (1+0*deriv) + 
        ( 8.57991e-02 + 6.63677e-02*vs) * (p^(2-deriv)) * (1+1*deriv) + 
            ( 7.39887e-01 - 6.26129e-02*vs) * (p^(3-deriv)) * (1+2*deriv));
    v = ( 3.29529e-01 + 1.86038e-02*vs) * (1-deriv) +  
        (-1.42155e-01 + 4.57806e-01*vs) * (p^(1-deriv)) * (1+0*deriv) +  
        (-3.29837e-01 + 5.59823e-02*vs) * (p^(2-deriv)) * (1+1*deriv) +  
        ( 1.41233e-01 + 4.72695e-01*vs) * (p^(3-deriv)) * (1+2*deriv);
    G = Gs * (( 2.17676e-01 + 7.22515e-02*vs) * (p^(1-deriv)) * (1+0*deriv) + 
        (-7.63847e-02 + 1.31601e+00*vs) * (p^(2-deriv)) * (1+1*deriv) + 
            ( 9.11800e-01 - 1.55261e+00*vs) * (p^(3-deriv)) * (1+2*deriv));
    return E,v,G
end

function orc_moduli(p, Es, vs, Gs, deriv)
    E = Es * (( 1.34332e-01 - 7.06384e-02*vs) * (p^(1-deriv)) * (1+0*deriv) + 
        ( 2.59957e-01 + 8.51515e-01*vs) * (p^(2-deriv)) * (1+1*deriv) + 
            ( 6.53902e-01 - 7.29803e-01*vs) * (p^(3-deriv)) * (1+2*deriv));
    v =      ( 3.38525e-01 + 7.04361e-03*vs) * (1-deriv) + 
            (-4.25721e-01 + 4.14882e-01*vs) * (p^(1-deriv)) * (1+0*deriv) + 
            (-7.68215e-02 + 5.58948e-01*vs) * (p^(2-deriv)) * (1+1*deriv) +  
            ( 1.64073e-01 + 3.98374e-02*vs) * (p^(3-deriv)) * (1+2*deriv);
    G = Gs * (( 1.96762e-01 + 1.66705e-01*vs) * (p^(1-deriv)) * (1+0*deriv) + 
            ( 1.30938e-01 + 1.72565e-01*vs) * (p^(2-deriv)) * (1+1*deriv) + 
            ( 6.45455e-01 - 2.87424e-01*vs) * (p^(3-deriv)) * (1+2*deriv));
    return E,v,G
end

function bound_moduli(p, Es, vs, Gs, deriv)
    Ks = 1.0 / (3*(1-2*vs));
    K = Ks + (1-p) / ( -1.0/Ks + p/(Ks + (4.0*Gs)/3.0) );
    G = Gs + (1-p) / ( -1.0/Gs + (2.0*p*(Ks+2.0*Gs)) / (5.0*Gs*(Ks+(4.0*Gs)/3.0)) );
    E = 9*K*G/(3*K+G);
    v = (3*K-2*G) / (2*(3*K+G));
    if deriv
        DK = (p - 1)/(((4*Gs)/3 + Ks)*(p/((4*Gs)/3 + Ks) - 1/Ks)^2) - 
            1/(p/((4*Gs)/3 + Ks) - 1/Ks);
        DG = 1/(1/Gs - (2*p*(2*Gs + Ks))/(5*Gs*((4*Gs)/3 + Ks))) + 
            (2*(2*Gs + Ks)*(p - 1))/(5*Gs*((4*Gs)/3 + Ks)*(1/Gs - 
            (2*p*(2*Gs + Ks))/(5*Gs*((4*Gs)/3 + Ks)))^2);
        DE = ( 9*(3*K+G)*(DK*G+K*DG) - 9*K*G*(3*DK+DG) ) / (3*K+G)^2;
        Dv = ( 2*(3*K+G)*(3*DK-2*DG) - 2*(3*K-2*G)*(3*DK+DG) ) / (2*(3*K+G))^2;
        G = DG;
        E = DE;
        v = Dv;
    end
    return E,v,G
end

function simp_moduli(p, Es, vs, Gs, deriv)
    E = Es * p^(3-deriv) * (1+2*deriv);
    v = vs * (1-deriv);
    G = Gs * p^(3-deriv) * (1+2*deriv);
    return E,v,G
end

######### TRUSS-SPECIFIC ROD DIAMETERS #########################################
function get_d(truss, p)
    if truss=="iso"
        d = 2.04920e-02 + 1.05076e+00*p - 1.59468e+00*(p.^2) + 1.09799e+00*(p.^3);
    elseif truss=="octet"
        d = 1.64505e-02 + 9.23773e-01*p - 1.61345e+00*(p.^2) + 1.23729e+00*(p.^3);
    elseif truss=="orc"
        d = 2.32950e-02 + 1.31602e+00*p - 2.28842e+00*(p.^2) + 1.90225e+00*(p.^3);
    else
        d = -1*ones(size(p));
    end
    return d
end
    
function topplot3d(xPhys,nelx, nely, nelz)
    ix=[]
    iy=[]
    iz=[]
    for j in 1:nely
        for i in 1:nelx
            for k in 1:nelz
                if (xPhys[i,j,k]> 0.5)
                    append!(ix,i)
                    append!(iy,j)
                    append!(iz,k)
                end
            end
        end
    end
    # r = 4.0
    # lim = FRect3D((-4,-4,-4*r),(8,8,8*r))
    return scatter(ix,iy,iz,color="black",label="",markersize =4, aspect_ratio=:equal,markerstrokealpha = 0.2,markeralpha = 0.6,markershape = :square,camera = (30, 60))#,markershape = :square
end

######### 3D VISUALIZATION #####################################################
# function viz3d(nelx, nely, nelz, x, volfrac, is2D)
#     y = zeros(nelx+2, nely+2, nelz+2); y(2:nelx+1, 2:nely+1, 2:nelz+1) = x;
#     if is2D; 
#         T=0; 
#         A=90; 
#         E=0; 
#     else; 
#         T=volfrac; 
#         A=142.5; 
#         E=30; 
#     end;
#     nf = nelx*nely*(nelz+1) + nelx*(nely+1)*nelz + (nelx+1)*nely*nelz; n = 0;
#     X = zeros(4,nf); Y = zeros(4,nf); Z = zeros(4,nf); C = zeros(1,nf);
#     for k = 1:nelz+1; 
#         for j = 1:nely+1; 
#             for i = 1:nelx+1;
#                 I = i-1; J = j-1; K = k-1; L = i+1; M = j+1; N = k+1;
#                 cz = max(y(L,M,k:N)); cy = max(y(L,j:M,N)); cx = max(y(i:L,M,N));
#                 dz = min(y(L,M,k:N)); dy = min(y(L,j:M,N)); dx = min(y(i:L,M,N));
#                 if cz > T && dz < T+is2D; n = n+1; C(1,n) = 1-cz; 
#                     X(:,n) = [I,i,i,I]'; Y(:,n) = [J,J,j,j]'; Z(:,n) = [K,K,K,K]';
#                 end
#                 if cy > T && dy < T+is2D; n = n+1; C(1,n) = 1-cy;
#                     X(:,n) = [I,i,i,I]'; Y(:,n) = [J,J,J,J]'; Z(:,n) = [K,K,k,k]';
#                 end
#                 if cx > T && dx < T+is2D; n = n+1; C(1,n) = 1-cx;
#                     X(:,n) = [I,I,I,I]'; Y(:,n) = [J,j,j,J]'; Z(:,n) = [K,K,k,k]';
#                 end
#             end; 
#         end; 
#     end
#     patch(X(:,1:n), Y(:,1:n), Z(:,1:n), C(1,1:n), 'EdgeColor', 'none');
#     view(A,E); 
#     axis equal; 
#     axis tight; 
#     axis off;
#      pause(1e-3);  
# end