using LinearAlgebra
using Plots
import JSON
using SparseArrays
using Images
using StaticArrays, BenchmarkTools

Based on:https://www.researchgate.net/publication/229061561_Bone_remodeling_as_a_hybrid_cellular_automaton_optimization_process

nelx=25;
nely=50;
X=ones(nelx,nely)*0.85;

fxF= zeros(Int,2,3);
fxF[1,1]=2 .*((nely)/2*(nelx+1));
fxF[1,2]=-100;
fxF[2,1]=2 .*((nely)/2*(nelx+1)+1);
fxF[2,2]=-100;
fxF[:,3].=1;

fxU=zeros(Int,3,3);
fxU[1,1]=2 .*(nelx+1);
fxU[2,1]=2 .*(nelx+1)+1;
fxU[3,1]=2 .*(nelx+1)*(nely+1);
fxU[:,3].=1;

Ymod=20 ;
nu=0.3;
p=3;
#equilibrium value
tgS=ones(nelx,nely)*5.0;

rho=2e-3;
w=0.25;
U=2261;
M=nelx*nely*2.0*1e-3;
tgS=ones(nelx,nely).*((1-w)/w*U/M*rho/p);
tgS=ones(nelx,nely).*1.8084;

#tgS=ones(nelx,nely)*(w*rho/p);

s=0.02;
scond=2;
# 3 #strict surface condition
# 2 #soft surface condition
# proportional gain cp cy
ct=0.2;
cp=0.2/tgS[1];
cp=0.4;
cI=0.4/tgS[1];
cD=0.2/tgS[1];
cII=0.4;
k=[ct,cp,cI,cD,cII];#0.2/1.8084
#Extended Moore (N=24), MvonN (N=12), Moore (N=8) ,von Newmann (N=4), > 24 # Global (nelx*ney-1)
neighbors=8;
# bcond==2 %periodic bcond bcond==1 %fixed bcond
bcond=1;
itmax=120;
vtol=2.5e-3;
ctype= "e";# f' 'p' 'i' 'd' 'e' 'n' 'pi' 'pd' 'id' 'pid'
ttype= 2; #or 2??
T=10;#????
movflg= 0;
matflg=0;
filestr="";
xfilter=true; #effective relative mass
sfilter=true; #EFFECTIVE STIMULUS

function hcastruc1(X, fxF, fxU, Ymod, nu, p, tgS, s, scond, k, neighbors, bcond, itmax, vtol, ctype, ttype, T, xfilter, sfilter)
    states=[]
    anim=Animation()
    # INPUT PARAMETERS
    nely, nelx = size(X);

    # CORRECT DENSITY
    X[X .< 0.001] .= 0.001;
    X[X .> 1.000] .= 1.000;
    # START HCA ALGORITHM

    println("running hca...")

    # FE-ANALYSIS
    SED = fea(copy(X), fxF, fxU, Ymod, nu, p);
    
    
    # MECHANICAL STIMULUS
    S=zeros(size(X))
    if ttype==1
        S = SED;
    else
        S = SED./X;
    end
    
    # PLOT
    it = 0;
    SE = sum(sum(SED,dims=1));
    M = sum(sum(X,dims=1));
    
    println("t: $it, SE: $SE, M: $M")
    
#     hfig = imagesc(-X,[-1 0]); axis equal; axis tight; axis off;
#     str = [ " t: " sprintf("#4i",it) ...
#         " SE: " sprintf("#10.4f",SE) ...
#         " M: " sprintf("#6.3f",M)];
#     title(str); pause(1e-6); disp(str);
    
    
    efX=zeros(size(X))
    efS=zeros(size(X))
    
    # EFFECTIVE RELATIVE MASS
    if xfilter
        if bcond==1 #fixed bcond
            efX = doavgf(X,neighbors,0);
        elseif bcond==2 #periodic bcond bcond==1 #fixed bcond
            efX = doavgp(X,neighbors);
        end
    else
        efX .= X;
    end
    # EFFECTIVE STIMULUS
    if sfilter
        if bcond==1 #fixed bcond
            efS = doavgf(S,neighbors,0);
        elseif bcond==2 #periodic bcond
            efS = doavgp(S,neighbors);
        end
    else
        efS = copy(S);
    end
    # EFFECTIVE ERROR SIGNAL
    efE = zeros(size(efS));
    efE[efS .> tgS] .= efS[efS .> tgS] .- (1+s)*tgS[efS .> tgS];
    efE[efS .< tgS] .= efS[efS .< tgS] .- (1-s)*tgS[efS .< tgS];
    # FOR INTEGRAL CONTROL: CREATE HISTORY OF STATES FOR INTEGRAL CONTROL
    efEhis = zeros(size(efE,1),size(efE,2),T+2);
    efEsum = zeros(size(efE,1),size(efE,2));
    # FOR CONVERGENCE
    Mtol = Inf;
    # FOR CONVERGENCE PLOT
    ITplt = 0; SEplt = SE; Mplt = M;
    # FOR DERIVATIVE CONTROL
    efEold = efE; #default
    #efEold = zeros(size(efE)); #modified
    # LOOP
    Xnew = copy(X);
    Mold = 0;
    while (it<itmax && Mtol>vtol)
        it = it + 1;
        # SURFACE CONDITION
        xid = efE .!=0;
        sneighbors = 4;
        if scond==2 #soft surface condition
            if bcond==1 #fixed bcond
                sfX = doavgf(X,sneighbors,0); #average in N+1 neighbors
            elseif bcond==2 #periodic bcond
                sfX = doavgf(X,sneighbors,0); #average in N+1 neighbors
            end
            xid=(sfX.<1.000) .& (sfX.>0.001)
            
#             x=[1 2 3;5 3 1]
#             m=(x.<3) .& (x.>1)
#             yyy=zeros(Int,size(sfX));zzz=zeros(Int,size(sfX))
#             yyy[sfX.<1.000 ].=1;zzz[sfX.>0.001].=1
#             xid = yyy .& zzz
            #xid .= (sfX.<1.000 .& sfX.>0.001);
        elseif scond==3 #strict surface condition
            if bcond==1 #fixed bcond
                sfX = doavgf(X,sneighbors,0); #avg in N+1 neighbors
                lzX = lookznf(X,sneighbors,0); #look for zeros around
                loX = lookonf(X,sneighbors,0); #look for ones around
            elseif bcond==2 #periodic bcond
                sfX = doavgp(X,sneighbors); #average in N+1 neighbors
                lzX = lookznp(X,sneighbors); #look for zeros around
                loX = lookonp(X,sneighbors); #look for ones around
            end
            xid=(sfX.<0.999) .& (sfX.>0.001)
            #xid .= (sfX.<0.999 .& sfX.>0.001);# | (X<0.999 & X>0.001);
            xid[X.<=0.001] .= loX[X.<=0.001];
            xid[X.>=1.000] .= lzX[X.>=1.000];
        end

        
        # UPDATE RULE
        if ctype== "f"
            Xnew[xid] .= efX[xid] + k[1]*sign.(efE[xid]);
        elseif ctype== "p"
            Xnew[xid] .= efX[xid] + k[2]*efE[xid]./tgS[xid];
        elseif ctype==  "i"
            Xnew[xid] .= efX[xid] + k[3]*efEsum[xid]./tgS[xid];
        elseif ctype==  "d"
            Xnew[xid] .= efX[xid] + k[4]*(efE[xid] - efEold[xid])./tgS[xid];
        elseif ctype=="e" #ratio approach

            if ttype==1
                Xnew[xid] = k[5]*efX[xid].*(efS[xid]./tgS[xid]).^(1.0/p);
            elseif ttype==2
                Xnew[xid] = k[5]*efX[xid].*(tgS[xid]./efS[xid]).^(1.0/p-1.0);
            end
        elseif ctype==  "n" #ratio approach variable set point
            if ttype==1
                Xnew[xid] .= k[5].*efX[xid].*(S[xid]./efS[xid]).^(1/p);
            elseif ttype==2
                Xnew[xid] .= k[5].*efX[xid].*(efS[xid]./S[xid]).^(1.0/p-1.0);
            end
        elseif ctype==  "pi"
            Xnew[xid] = efX[xid] + k[2]*efE[xid]./tgS[xid];
            Xnew[xid] = Xnew[xid] + k[3]*efEsum[xid]./tgS[xid];
        elseif ctype==  "pd"
            Xnew[xid] = efX[xid] + k[2]*efE[xid]./tgS[xid];
            Xnew[xid] = Xnew[xid] + k[4]*(efE[xid] - efEold[xid])./tgS[xid];
        elseif ctype==  "id"
            Xnew[xid] = efX[xid] + k[3]*efEsum[xid]./tgS[xid];
            Xnew[xid] = Xnew[xid] + k[4]*(efE[xid] - efEold[xid])./tgS[xid];
        elseif ctype==  "pid"
            Xnew[xid] = efX[xid] + k[2]*efE[xid]./tgS[xid];
            Xnew[xid] = Xnew[xid] + k[3]*efEsum[xid]./tgS[xid];
            Xnew[xid] = Xnew[xid] + k[4]*(efE[xid] - efEold[xid])./tgS[xid];
        end
                
        # FE-ANALYSIS
        SEDnew = fea(copy(Xnew), fxF, fxU, Ymod, nu, p);
        # MECHANICAL STIMULUS
        if ttype==1
            Snew = SEDnew;
        else
            Snew = SEDnew./Xnew;
        end
        # CONVERGENCE
        Mnew = sum(sum(Xnew,dims=1));
        Mtol = (abs(Mnew - M) .+ abs(M - Mold))./2.0;
        # UPDATE
        SED = SEDnew;
        X = Xnew;
        S = Snew;
        Mold = M;
        M = Mnew;
        # FOR DERIVATIVE CONTROL
        efEold = efE;
        # EFFECTIVE RELATIVE MASS
        if xfilter
            if bcond==1 #fixed bcond
                efX = doavgf(X,neighbors,0);
            elseif bcond==2 #periodic bcond
                efX = doavgp(X,neighbors);
            end
        else
            efX = X;
        end
        # EFFECTIVE STIMULUS
        if sfilter
            if bcond==1 #fixed bcond
                efS = doavgf(S,neighbors,0);
            elseif bcond==2 #periodic bcond
                efS = doavgp(S,neighbors);
            end
        else
            efS = S;
        end
        # EFFECTIVE ERROR SIGNAL
        efE = zeros(size(S));
        efE[efS .> tgS] .= efS[efS .> tgS] .- (1+s).*tgS[efS .> tgS];
        efE[efS .< tgS] .= efS[efS .< tgS] .- (1-s).*tgS[efS .< tgS];
        # FOR INTEGRAL CONTROL: SUM OF STORED EFFECTIVE STIMULI EST
        efEsum .= efEsum .+ efE .- efEhis[:,:,1];
        
        efEhis[:,:,T+2] .= efE;
        
        for t in 1:T+1
            efEhis[:,:,t] .= efEhis[:,:,t+1];
        end
        
        # PLOT DENSITIES
        SE = sum(sum(SED));
        
        println("t: $it, SE: $SE, M: $M")

        # STORE FOR CONVERGENCE PLOT

        display(Gray.(Xnew))
        heatmap(clamp.(Xnew,0,1),xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it")
        append!(states,[clamp.(Xnew,0,1)])
        frame(anim)
    end
    
    # SAVE FINAL RESULTS
    cells = Xnew;
    
    return anim,cells,states
    

end

hcastruc1 (generic function with 1 method)

function  doavgf(c,neighbors,val)
    # DOAVGF(X, neighbors,val) sums around each X using fixed boundary
    # conditions. val is the fixed value.
#     if nargin <3
#         val=0;
#     end
    nely,nelx = size(c);
    cext=zeros(nely+2,nelx+2);
    cextm=zeros(nely+4,nelx+4);
    if neighbors > 24 # Global (nelx*ney-1)
        csum = ones(size(c))*sum(sum(c));
    else
        csum = c; # no neighbors (N=0)
        if neighbors > 0 # von Newmann (N=4) > 24 # Global (nelx*ney-1)
            # extends c field for fixed BC
            cext[nely+2, nelx+2] = val;
            cext[2:nely+1, 2:nelx+1] = c;
            x = 2:nelx+1;
            y = 2:nely+1;
            csum = csum + 
            cext[y,x.+1] + cext[y,x.-1] + 
            cext[y.+1,x] +  cext[y.-1,x];
            if neighbors > 4 # Moore (N=8)von Newmann (N=4) > 24 # Global (nelx*ney-1)
                csum = csum + 
                cext[y.+1,x.+1] + cext[y.-1,x.-1] +
                cext[y.+1,x.-1] + cext[y.-1,x.+1];
                
                if neighbors > 8 # MvonN (N=12)Moore (N=8)von Newmann (N=4) > 24 # Global (nelx*ney-1)
                    cextm[nely+4, nelx+4] = val;
                    cextm[3:nely+2, 3:nelx+2] = c;
                    y = 3:nely+2; x = 3:nelx+2;
                    csum = csum + 
                    cextm[y.-2, x] + cextm[y.+2, x] + 
                    cextm[y, x.-2] + cextm[y, x.+2];
                    if neighbors > 12 # Extended Moore (N=24)MvonN (N=12)Moore (N=8)von Newmann (N=4) > 24 # Global (nelx*ney-1)
                        csum = csum + 
                        cextm[y.-2, x.-2] + cextm[y.-2, x.-1] +  
                        cextm[y.-2, x.+1] + cextm[y.-2, x.+2] + 
                        cextm[y.+2, x.-2] + cextm[y.+2, x.-1] +  
                        cextm[y.+2, x.+1] + cextm[y.+2, x.+2] + 
                        cextm[y.-1, x.-2] + cextm[y.+1, x.-2] +  
                        cextm[y.-1, x.+2] + cextm[y.+1, x.+2];
                    end
                end
            end
        end
    end
    cavg = csum./(neighbors+1); # average csum (matrix)
    return cavg
end
#-----------------------------------------------------------------

doavgf (generic function with 1 method)

function doavgp(c,neighbors)
    # DOAVGNP(X, neighbors) sums around each X using fixed periodic
    # conditions
    nely,nelx = size(c);
    x = 1:nelx;
    y = 1:nely;
    if neighbors > 24 # Global (nelx*ney-1)
        csum = ones(size(c))*sum(sum(c))-c;
    else
        csum = zeros(size(c)); # no neighbors (N=0)
        if neighbors > 0 # von Newmann (N=4)
            csum[mod(y,nely)+1, mod(x,nelx)+1] = csum[mod(y,nely)+1, mod(x,nelx)+1] + 
                c[mod(y,nely)+1, mod(x.+1,nelx)+1] + 
                c[mod(y,nely)+1, mod(x.-1,nelx)+1] + 
                c[mod(y.+1,nely)+1, mod(x,nelx)+1] + 
                c[mod(y.-1,nely)+1, mod(x,nelx)+1];
            if neighbors > 4 # Moore (N=8)
                csum[mod(y,nely)+1, mod(x,nelx)+1] = csum[mod(y,nely)+1, mod(x,nelx)+1] + 
                    c[mod(y+1,nely)+1, mod(x+1,nelx)+1] + 
                    c[mod(y-1,nely)+1, mod(x-1,nelx)+1] + 
                    c[mod(y+1,nely)+1, mod(x-1,nelx)+1] + 
                    c[mod(y-1,nely)+1, mod(x+1,nelx)+1];
                if neighbors > 8 # MvonN (N=12)
                    csum[mod(y,nely)+1, mod(x,nelx)+1] = csum[mod(y,nely)+1, mod(x,nelx)+1] + 
                        c[mod(y.-2,nely)+1, mod(x,nelx)+1] + 
                        c[mod(y.+2,nely)+1, mod(x,nelx)+1] + 
                        c[mod(y,nely)+1, mod(x.-2,nelx)+1] + 
                        c[mod(y,nely)+1, mod(x.+2,nelx)+1];
                    
                    if neighbors > 12 # Extended Moore (N=24)
                        csum[mod(y,nely)+1, mod(x,nelx)+1] = csum[mod(y,nely)+1, mod(x,nelx)+1] + 
                            c[mod(y.-2,nely)+1, mod(x.-2,nelx)+1] + 
                            c[mod(y.-2,nely)+1, mod(x.-1,nelx)+1] + 
                            c[mod(y.-2,nely)+1, mod(x.+1,nelx)+1] + 
                            c[mod(y.-2,nely)+1, mod(x.+2,nelx)+1] + 
                            c[mod(y.+2,nely)+1, mod(x.-2,nelx)+1] + 
                            c[mod(y.+2,nely)+1, mod(x.-1,nelx)+1] + 
                            c[mod(y.+2,nely)+1, mod(x.+1,nelx)+1] + 
                            c[mod(y.+2,nely)+1, mod(x.+2,nelx)+1] + 
                            c[mod(y.-1,nely)+1, mod(x.-2,nelx)+1] + 
                            c[mod(y.+1,nely)+1, mod(x.-2,nelx)+1] + 
                            c[mod(y.-1,nely)+1, mod(x.+2,nelx)+1] + 
                            c[mod(y.+1,nely)+1, mod(x.+2,nelx)+1];
                
                    end
                end
            end
        end
    end
    cavg = csum./(neighbors); # average csum (matrix)
    return cavg
end

doavgp (generic function with 1 method)

function  lookznf(c,neighbors,val)
    # LOOKZNF(X, neighbors) products around each X using fixed boundary
    # conditions
#     if nargin <3
#         val=0;
#     end
    c[c.<=0.001].=0;
    nely,nelx = size(c);
    cext=zeros(nely+2,nelx+2);
    cextm=zeros(nely+4,nelx+4);
    cpro = copy(c);
    if neighbors <= 24
        cpro = c; # no neighbors (N=0)
        if neighbors > 0 # von Newmann (N=4)
            # extends c field for fixed BC
            cext[nely+2, nelx+2] = val;
            cext[2:nely+1, 2:nelx+1] = c;
            x = 2:nelx+1;
            y = 2:nely+1;
            cpro = 1*
            cext[y,x.+1] .* cext[y,x.-1] .* 
            cext[y.+1,x] .* cext[y.-1,x];
            if neighbors > 4 # Moore (N=8)
                cpro = cpro .* 
                cext[y.+1,x.+1] .* cext[y.-1,x.-1] .* 
                cext[y.+1,x.-1] .* cext[y.-1,x.+1];
                if neighbors > 8 # MvonN (N=12)
                    cextm[nely+4, nelx+4] = val;
                    cextm[3:nely+2, 3:nelx+2] = c;
                    y = 3:nely+2; x = 3:nelx+2;
                    cpro = cpro .* 
                    cextm[y.-2, x] .* cextm[y.+2, x] .* 
                    cextm[y, x.-2] .* cextm[y, x.+2];
                    if neighbors > 12 # Extended Moore (N=24)
                        csum = csum .* 
                            cextm[y.-2, x.-2] .* cextm[y.-2, x.-1] .* 
                            cextm[y.-2, x.+1] .* cextm[y.-2, x.+2] .* 
                            cextm[y.+2, x.-2] .* cextm[y.+2, x.-1] .* 
                            cextm[y.+2, x.+1] .* cextm[y.+2, x.+2] .* 
                            cextm[y.-1, x.-2] .* cextm[y.+1, x.-2] .* 
                            cextm[y.-1, x.+2] .* cextm[y.+1, x.+2];
                    
                    end
                end
            end
        end
    end
    cproy=copy(cpro)
    cpro[cproy.>0].=0
    cpro[cproy.==0].=1
    return cpro
end
#-----------------------------------------------------------------
function  lookznp(c,neighbors)
    # LOOKZNP(X, neighbors) products around each X using fixed periodic
    # conditions
    c[c.<=0.001].=0;
    nely,nelx = size(c);
    x = 1:nelx;
    y = 1:nely;
    cpro = copy(c);
    if neighbors <= 24
        if neighbors > 0 # von Newmann (N=4)
            cpro[mod.(y,nely).+1, mod.(x,nelx).+1] = 
                c[mod.(y,nely).+1, mod.(x+1,nelx).+1] .* 
                c[mod.(y,nely).+1, mod.(x-1,nelx).+1] .* 
                c[mod.(y.+1,nely).+1, mod.(x,nelx).+1] .* 
                c[mod.(y.-1,nely).+1, mod.(x,nelx).+1];
            if neighbors > 4 # Moore (N=8)
                cpro[mod(y,nely)+1, mod(x,nelx)+1] = 
                    cpro[mod.(y,nely).+1, mod.(x,nelx).+1] .* 
                    c[mod.(y.+1,nely).+1, mod.(x.+1,nelx).+1] .* 
                    c[mod.(y.-1,nely).+1, mod.(x.-1,nelx).+1] .* 
                    c[mod.(y.+1,nely).+1, mod.(x.-1,nelx).+1] .* 
                    c[mod.(y.-1,nely).+1, mod.(x.+1,nelx).+1];
                if neighbors > 8 # MvonN (N=12)
                    cpro[mod(y,nely)+1, mod(x,nelx)+1] = 
                        cpro[mod(y,nely)+1, mod.(x,nelx)+1] .* 
                        c[mod.(y.-2,nely)+1, mod.(x,nelx)+1] .* 
                        c[mod.(y.+2,nely)+1, mod.(x,nelx)+1] .* 
                        c[mod.(y,nely)+1, mod.(x.-2,nelx)+1] .* 
                        c[mod.(y,nely)+1, mod.(x.+2,nelx)+1];
                    if neighbors > 12 # Extended Moore (N=24)
                        cpro[mod(y,nely)+1, mod(x,nelx)+1] = 
                            cpro[mod(y,nely)+1, mod(x,nelx)+1] .* 
                            c[mod.(y.-2,nely).+1, mod.(x.-2,nelx).+1] .* 
                            c[mod.(y.-2,nely).+1, mod.(x.-1,nelx).+1] .* 
                            c[mod.(y.-2,nely).+1, mod.(x.+1,nelx).+1] .* 
                            c[mod.(y.-2,nely).+1, mod.(x.+2,nelx).+1] .* 
                            c[mod.(y.+2,nely).+1, mod.(x.-2,nelx).+1] .* 
                            c[mod.(y.+2,nely).+1, mod.(x.-1,nelx).+1] .* 
                            c[mod.(y.+2,nely).+1, mod.(x.+1,nelx).+1] .* 
                            c[mod.(y.+2,nely).+1, mod.(x.+2,nelx).+1] .* 
                            c[mod.(y.-1,nely).+1, mod.(x.-2,nelx).+1] .* 
                            c[mod.(y.+1,nely).+1, mod.(x.-2,nelx).+1] .* 
                            c[mod.(y.-1,nely).+1, mod.(x.+2,nelx).+1] .* 
                            c[mod.(y.+1,nely).+1, mod.(x.+2,nelx).+1];
                    end
                end
            end
        end
    end
    cproy=copy(cpro)
    cpro[cproy.>0].=0
    cpro[cproy.==0].=1

                                                                return cpro
end

lookznp (generic function with 1 method)

function  lookonf(c,neighbors,val)
    # LOOKONF(X, neighbors) sums around each X using fixed boundary
    # conditions
#     if nargin <3
#         val=0;
#     end
    c[c.>=0.999].=1;
    c[c.<0.999].=0;
    nely,nelx = size(c);
    cext=zeros(nely+2,nelx+2);
    cextm=zeros(nely+4,nelx+4);
    if neighbors <= 24 # Global (nelx*ney-1)
        csum = zeros(size(c)); # no neighbors (N=0)
        if neighbors > 0 # von Newmann (N=4)
            # extends c field for fixed BC
            cext[nely+2, nelx+2] = val;
            cext[2:nely+1, 2:nelx+1] = c;
            x = 2:nelx+1;
            y = 2:nely+1;
            csum = csum + cext[y,x+1] + cext[y,x-1] + 
                cext[y.+1,x] + cext[y.-1,x];
            if neighbors > 4 # Moore (N=8)
                csum = csum + cext[y.+1,x.+1] + cext[y.-1,x.-1] + 
                    cext[y.+1,x.-1] + cext[y.-1,x.+1];
                if neighbors > 8 # MvonN (N=12)
                    cextm[nely+4, nelx+4] = val;
                    cextm[3:nely+2, 3:nelx+2] = c;
                    y = 3:nely+2; x = 3:nelx+2;
                    csum = csum + cextm[y-2, x] + cextm[y+2, x] + 
                        cextm[y, x.-2] + cextm[y, x.+2];
                    if neighbors > 12 # Extended Moore (N=24)
                        csum = csum + 
                            cextm[y.-2, x.-2] + cextm[y.-2, x.-1] + 
                            cextm[y.-2, x.+1] + cextm[y.-2, x.+2] + 
                            cextm[y.+2, x.-2] + cextm[y.+2, x.-1] + 
                            cextm[y.+2, x.+1] + cextm[y.+2, x.+2] + 
                            cextm[y.-1, x.-2] + cextm[y.+1, x.-2] + 
                            cextm[y.-1, x.+2] + cextm[y.+1, x.+2];
                    end
                end
            end
        end
    end
    cavg = csum./(neighbors); # average csum (matrix)
    cout = cavg > 0;        
    return cout
end
#-----------------------------------------------------------------
function  looknp(c,neighbors)
    # LOOKONP(X, neighbors) sums around each X using fixed periodic
    # conditions
    c[c.>=0.999].=1;
    c[c.<0.999].=0;
    nely,nelx = size(c);
    x = 1:nelx;
    y = 1:nely;
    if neighbors <= 24
        csum = zeros(size(c)); # no neighbors (N=0)
        if neighbors > 0 # von Newmann (N=4)
            csum[mod.(y,nely).+1, mod.(x,nelx)+1] = 
                csum[mod(y,nely).+1, mod.(x,nelx).+1] + 
                c[mod.(y,nely).+1, mod.(x.+1,nelx).+1] + 
                c[mod.(y,nely).+1, mod.(x.-1,nelx).+1] + 
                c[mod.(y.+1,nely).+1, mod.(x,nelx).+1] + 
                c[mod.(y.-1,nely).+1, mod.(x,nelx).+1];
            if neighbors > 4 # Moore (N=8)
                csum[mod.(y,nely)+1, mod.(x,nelx).+1] = 
                    csum[mod.(y,nely).+1, mod.(x,nelx).+1] + 
                    c[mod.(y.+1,nely).+1, mod.(x.+1,nelx).+1] + 
                    c[mod.(y.-1,nely).+1, mod.(x.-1,nelx).+1] + 
                    c[mod.(y.+1,nely).+1, mod.(x.-1,nelx).+1] + 
                    c[mod.(y.-1,nely).+1, mod.(x.+1,nelx).+1];
                if neighbors > 8 # MvonN (N=12)
                    csum[mod(y,nely)+1, mod(x,nelx)+1] = 
                        csum[mod.(y,nely).+1, mod.(x,nelx).+1] + 
                        c[mod.(y.-2,nely).+1, mod.(x,nelx).+1] + 
                        c[mod.(y.+2,nely).+1, mod.(x,nelx).+1] + 
                        c[mod.(y,nely).+1, mod.(x.-2,nelx).+1] + 
                        c[mod.(y,nely).+1, mod.(x.+2,nelx).+1];
                    if neighbors > 12 # Extended Moore (N=24)
                        csum[mod(y,nely)+1, mod(x,nelx)+1] = 
                            csum[mod(y,nely)+1, mod(x,nelx)+1] + 
                            c[mod.(y.-2,nely).+1, mod.(x.-2,nelx).+1] + 
                            c[mod.(y.-2,nely).+1, mod.(x.-1,nelx).+1] + 
                            c[mod.(y.-2,nely).+1, mod.(x.+1,nelx).+1] + 
                            c[mod.(y.-2,nely).+1, mod.(x.+2,nelx).+1] + 
                            c[mod.(y.+2,nely).+1, mod.(x.-2,nelx).+1] + 
                            c[mod.(y.+2,nely).+1, mod.(x.-1,nelx).+1] + 
                            c[mod.(y.+2,nely).+1, mod.(x.+1,nelx).+1] + 
                            c[mod.(y.+2,nely).+1, mod.(x.+2,nelx).+1] + 
                            c[mod.(y.-1,nely).+1, mod.(x.-2,nelx).+1] + 
                            c[mod.(y.+1,nely).+1, mod.(x.-2,nelx).+1] + 
                            c[mod.(y.-1,nely).+1, mod.(x.+2,nelx).+1] + 
                            c[mod.(y.+1,nely).+1, mod.(x.+2,nelx).+1];
                    end
                end
            end
        end
    end
    cavg = csum./(neighbors); # average csum (matrix)
    cout = cavg>0;                  
    return cout
end

looknp (generic function with 1 method)

function  fea(x, fxF, fxU, E, nu, p)
    nelx = size(x,2);
    nely = size(x,1);
    lcase = convert(Int,maximum(fxF[:,3]));
    # Element stiffness matrix for a bilinear element
    k=[1/2-nu/6,1/8+nu/8,-1/4-nu/12,-1/8+3*nu/8,-1/4+nu/12,-1/8-nu/8,nu/6,1/8-3*nu/8]
    KE = E/(1-nu^2)*[  k[0+1] k[1+1] k[2+1] k[3+1] k[4+1] k[5+1] k[6+1] k[7+1];
                       k[1+1] k[0+1] k[7+1] k[6+1] k[5+1] k[4+1] k[3+1] k[2+1];
                       k[2+1] k[7+1] k[0+1] k[5+1] k[6+1] k[3+1] k[4+1] k[1+1];
                       k[3+1] k[6+1] k[5+1] k[0+1] k[7+1] k[2+1] k[1+1] k[4+1];
                       k[4+1] k[5+1] k[6+1] k[7+1] k[0+1] k[1+1] k[2+1] k[3+1];
                       k[5+1] k[4+1] k[3+1] k[2+1] k[1+1] k[0+1] k[7+1] k[6+1];
                       k[6+1] k[3+1] k[4+1] k[1+1] k[2+1] k[7+1] k[0+1] k[5+1];
                       k[7+1] k[2+1] k[1+1] k[4+1] k[3+1] k[6+1] k[5+1] k[0+1] ];
 
    # Global stiffness matrix, and force and displacement vectors
    K = spzeros(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
    F = spzeros(2*(nely+1)*(nelx+1), lcase);
    U = spzeros(2*(nely+1)*(nelx+1), lcase);
    # Obtain global stiffness matrix
    for ely = 1:nely
        for elx = 1:nelx
            n1 = (nely+1)*(elx-1)+ely;
            n2 = (nely+1)* elx +ely;
            edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 
                2*n1+1; 2*n1+2];
            K[edof,edof] = K[edof,edof] + x[ely,elx]^p*KE;
        end
    end
    # DEFINE LOADS AND SUPPORTS
    # loads
    for i = 1:size(fxF,1)
        F[fxF[i,1],fxF[i,3]] = F[fxF[i,1],fxF[i,3]] + fxF[i,2];
    end
    # supports
    for lc = 1:lcase
        fxdof = fxU[fxU[:,3].==lc, 1];
        fxdis = fxU[fxU[:,3].==lc, 2];
        if isempty(fxdof)
            fxdof = fxU[fxU[:,3].==1, 1];
            fxdis = fxU[fxU[:,3].==1, 2];
        end
        alldof = 1:2*(nely+1)*(nelx+1);
        freedof = setdiff(alldof,fxdof);
        # solution
        U[freedof,lc] = K[freedof,freedof] \ Array(F[freedof,lc]);
        U[fxdof,lc] = fxdis;
    end
    # STRAIN ENERGY DENSITIES AND SENSITIVITIES
    SEDt = zeros(nely,nelx,lcase);
    SED = zeros(nely,nelx);
    for ely = 1:nely
        for elx = 1:nelx
            n1 = (nely+1)*(elx-1)+ely;
            n2 = (nely+1)* elx +ely;
#             dc=zeros(ely,elx)
#             dc(ely,elx) = 0.;
            for lc = 1:lcase
                Ue = U[[2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2], lc];
                SEDi = 1/2*x[ely,elx]^p*Ue'*KE*Ue;
                xi = x[ely,elx];
                if lcase >1
                    SEDt[ely,elx,lc] = SEDi;
                end
                SED[ely,elx] = SED[ely,elx]  + SEDi;
            end
        end
    end
    SED=SED/lcase;
#     SEDt
    return SED
end

fea (generic function with 1 method)

# SED = fea(X, fxF, fxU, Ymod, nu, p)
# size(SED)

(anim,cells,states)=hcastruc1(X, fxF, fxU, Ymod, nu, p, tgS, s, scond, k, neighbors, bcond, itmax, vtol, ctype, ttype, T, xfilter, sfilter)

running hca...
t: 0, SE: 12454.446433867357, M: 1062.5000000000005

t: 1, SE: 7042.168539704487, M: 1006.6855425152736

t: 2, SE: 3655.305323109011, M: 1090.7786184383322

t: 3, SE: 4657.3674455733, M: 932.1769217592991

t: 4, SE: 3800.076463290977, M: 929.7904231789895

t: 5, SE: 4187.18368700627, M: 844.0868821264332

t: 6, SE: 4362.189878452322, M: 798.7137452456893

t: 7, SE: 4571.71721884907, M: 771.8577299522622

t: 8, SE: 4601.427462909138, M: 762.0504374238097

t: 9, SE: 4696.834127896693, M: 755.1582673656225

t: 10, SE: 4674.898741142548, M: 754.5081323917747

t: 11, SE: 4717.249506603515, M: 751.9140192419256

t: 12, SE: 4711.569617660478, M: 751.7891421122727

t: 13, SE: 4731.186893184453, M: 750.925921508149

t: 14, SE: 4726.244292247353, M: 751.0789675753692

t: 15, SE: 4729.690171794538, M: 750.9476190717361

t: 16, SE: 4727.166810745661, M: 751.0343765179421

t: 17, SE: 4730.997992438792, M: 750.922830999645

t: 18, SE: 4725.995612914041, M: 751.0866367016869

t: 19, SE: 4730.297171705467, M: 750.9254495451745

t: 20, SE: 4725.940990459331, M: 751.0609772498678

t: 21, SE: 4732.727935014125, M: 750.8717271351957

t: 22, SE: 4723.663868957003, M: 751.1823747982032

t: 23, SE: 4731.995769182213, M: 750.8758530926027

t: 24, SE: 4724.2773927003145, M: 751.1440745543048

t: 25, SE: 4733.807509431284, M: 750.850637350602

t: 26, SE: 4722.007461004263, M: 751.2535066750843

t: 27, SE: 4729.761311005515, M: 750.9138150930023

t: 28, SE: 4728.774835472978, M: 751.0007183521186

t: 29, SE: 4728.312445331505, M: 751.046310095217

t: 30, SE: 4728.490738763299, M: 751.0521451588985

t: 31, SE: 4728.247736780054, M: 751.0696601291545

t: 32, SE: 4728.257699072732, M: 751.0700514223231

t: 33, SE: 4728.362908934972, M: 751.0712196029715

(Animation("C:\\Users\\amira\\AppData\\Local\\Temp\\jl_yiLIa5", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000024.png", "000025.png", "000026.png", "000027.png", "000028.png", "000029.png", "000030.png", "000031.png", "000032.png", "000033.png"]), [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 5.874015046735567 9.047153224432313 … 8.577344625427084 5.5688701804425715; 3.509510163302371 5.402276460590305 … 5.1217508861802985 3.3272343841242846], Any[[0.0006233950299422973 0.003089256160619748 … 0.0028240505358137014 0.0005758241678701575; 0.0031078362600568655 0.013053948856162055 … 0.01199245639696917 0.0028424321355615366; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]  …  [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [2.5252343829730185e-6 2.6728064095499655e-5 … 2.382152735854522e-5 2.222818470351869e-6; 2.7076235419536834e-5 0.00023524045226610182 … 0.0002119905618436961 2.410556744740308e-5; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]])

gif(anim,"hca.gif", fps=5)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\MetaVoxels\01_Code\hca.gif
└ @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98

X=clamp.(cells,0,1)

25×50 Array{Float64,2}:
 2.52523e-6   2.67281e-5   0.000150738  …  2.38215e-5   2.22282e-6 
 2.70762e-5   0.00023524   1.1361e-6       0.000211991  2.41056e-5 
 0.00015313   1.32402e-6   1.33371e-5      8.76866e-7   0.00013775 
 7.25021e-7   1.1822e-5    9.69186e-5      7.98278e-6   4.81867e-7 
 5.49212e-6   7.31585e-5   0.000500065     5.03585e-5   3.72855e-6 
 3.05761e-5   0.000341964  1.07893e-5   …  0.000240049  2.12105e-5 
 0.000132913  6.5704e-6    9.48001e-5      2.37919e-6   9.43069e-5 
 1.78465e-6   5.54479e-5   7.44395e-6      2.22294e-5   0.000342706
 1.2539e-5    3.55732e-6   0.00012855      0.000149728  5.74879e-6 
 6.40667e-5   4.98984e-5   2.23903e-5      1.16068e-5   3.27322e-5 
 4.83025e-6   5.2406e-5    0.000145029  …  7.84576e-6   0.000147723
 4.36169e-5   1.64088e-5   0.00369063      3.44817e-5   1.59093e-5 
 1.18509e-5   0.000622605  0.0374176       4.97685e-5   1.06183e-5 
 4.09533e-5   0.00965537   0.170965        0.00195789   1.06828e-5 
 0.000804651  0.0827356    0.498404        0.0316707    0.000139621
 0.0110355    0.264432     0.963807     …  0.154984     0.00385374 
 0.0449524    1.0          1.0             1.0          0.0241523  
 0.13057      1.0          1.0             1.0          0.0834639  
 1.0          1.0          1.0             1.0          1.0        
 1.0          1.0          1.0             1.0          1.0        
 1.0          1.0          1.0          …  1.0          1.0        
 1.0          1.0          1.0             1.0          1.0        
 1.0          1.0          1.0             1.0          1.0        
 1.0          1.0          1.0             1.0          1.0        
 1.0          1.0          1.0             1.0          1.0        

nel=(nelx+1)*(nely+1)
ind=1:nel
function getIndex(i,j,nelx,nely)
    return (i-1)*(nely+1)+(j-1)+1
end
function getIndex(ind,nelx,nely)
    j=mod((ind-1),nely+1)
    i=(ind-1)÷(nely+1)
    return (i+1),(j+1)
end
nel=(nelx)*(nely)
ind=1:nel
function getIndex(i,j,nelx,nely)
    return (i-1)*(nely)+(j-1)+1
end
function getIndex(ind,nelx,nely)
    j=mod((ind-1),nely)
    i=(ind-1)÷(nely)
    return (i+1),(j+1)
end

getIndex (generic function with 2 methods)

#8 neighbours
neigh=[]
indIJ=zeros(nel,2)
XIJ=zeros(nel)
for ii in ind
    n=[]
    (i,j)=getIndex(ii,nelx,nely)
    indIJ[ii,1]=i
    indIJ[ii,2]=j
    XIJ[ii]=X[i,j]
    if i>1
        append!(n,[[i-1,j]])
        if j>1
            append!(n,[[i-1,j-1]])
        end
        if j<=nely
            append!(n,[[i-1,j+1]])
        end
    end
    if i<nelx
        append!(n,[[i+1,j]])
        if j>1
            append!(n,[[i+1,j-1]])
        end
        if j<=nely
            append!(n,[[i+1,j+1]])
        end
    end
    if j>1
        append!(n,[[i,j-1]])
    end
    if j<nely
        append!(n,[[i,j+1]])
    end
    append!(neigh,[n])  
end
neigh;

#4 neighbours
neigh=[]
indIJ=zeros(Int,nel,2)
XIJ=zeros(nel)
for ii in ind
    n=[]
    (i,j)=getIndex(ii,nelx,nely)
    indIJ[ii,1]=i
    indIJ[ii,2]=j
    XIJ[ii]=X[i,j]
    if i>1
        append!(n,[[i-1,j]])
    end
    if i<nelx
        append!(n,[[i+1,j]])
    end
    if j>1
        append!(n,[[i,j-1]])
    end
    if j<nely
        append!(n,[[i,j+1]])
    end
    append!(neigh,[n])  
end
neigh;

links=[]
XIJE=[]

for i in ind
    for j in 1:length(neigh[i])
        append!(links,[[[indIJ[i,2],neigh[i][j][2]],[indIJ[i,1],neigh[i][j][1]]]])
        append!(XIJE,min(X[indIJ[i,1],indIJ[i,2]],X[neigh[i][j][1],neigh[i][j][2]]))
    end
end

XIJEs=[]
for ii in 1:length(states)
    xx=[]
    for i in ind
        for j in 1:length(neigh[i])
            append!(xx,min(reverse(states[ii],dims=1)[indIJ[i,1],indIJ[i,2]],reverse(states[ii],dims=1)[neigh[i][j][1],neigh[i][j][2]]))
        end
    end
    append!(XIJEs,[xx])
end
XIJEs;

linkss=zeros(length(links),2,2)
for i in 1:length(links)
    for j in 1:2
        for k in 1:2
            linkss[i,j,k]=links[i][j][k]
        end
    end
end
linkss[:,1,:]

4850×2 Array{Float64,2}:
  1.0   1.0
  1.0   2.0
  2.0   2.0
  2.0   1.0
  2.0   3.0
  3.0   3.0
  3.0   2.0
  3.0   4.0
  4.0   4.0
  4.0   3.0
  4.0   5.0
  5.0   5.0
  5.0   4.0
  ⋮        
 46.0  47.0
 47.0  47.0
 47.0  46.0
 47.0  48.0
 48.0  48.0
 48.0  47.0
 48.0  49.0
 49.0  49.0
 49.0  48.0
 49.0  50.0
 50.0  50.0
 50.0  49.0

scatter(indIJ[:,2],indIJ[:,1],color=:black,label="",markeralpha = XIJ, aspect_ratio=:equal,xlim=(0,nely+2),ylim=(0,nelx+2))
for i in ind
    for j in 1:length(neigh[i])
        plot!([indIJ[i,2],neigh[i][j][2]],[indIJ[i,1],neigh[i][j][1]],label="",color=:blue,linewidth = 1,linealpha = XIJ[i])
    end
end
plot!()

scatter(indIJ[:,2],indIJ[:,1],color=:black,label="",markeralpha = XIJ, aspect_ratio=:equal,xlim=(0,nely+2),ylim=(0,nelx+2))
plot!(linkss[:,1,:]',linkss[:,2,:]',label="",color=:blue,linealpha = XIJE',linewidth = 1)
# for i in 1:length(XIJE)
#      plot!(links[i][1],links[i][2],label="",color=:blue,linewidth = 1,linealpha = XIJE[i])
# end
# plot!()

anim1=Animation()
######
i=0
display(i)
scatter(indIJ[:,2],indIJ[:,1],color=:black,label="",markeralpha = 0.85, aspect_ratio=:equal,xlim=(0,nely+2),ylim=(0,nelx+2),title="it: $i")
plot!(linkss[:,1,:]',linkss[:,2,:]',label="",color=:blue,linealpha = 1,linewidth = 0.85)
frame(anim1)
####
# for i in 1:length(states)
for i in 1:10
    display(i)
    scatter(indIJ[:,2],indIJ[:,1],color=:black,label="",markeralpha = reverse(states[i],dims=1)'[:], aspect_ratio=:equal,xlim=(0,nely+2),ylim=(0,nelx+2),title="it: $i")
    plot!(linkss[:,1,:]',linkss[:,2,:]',label="",color=:blue,linealpha = XIJEs[i]',linewidth = 1)
    frame(anim1)
end

0

1

2

3

4

5

6

7

8

9

10

gif(anim1,"hca_strut.gif", fps=1)

┌ Info: Saved animation to 
│   fn = C:\Users\amira\Dropbox (MIT)\CBA\MetaVoxels\01_Code\hca_strut.gif
└ @ Plots C:\Users\amira\.julia\packages\Plots\cc8wh\src\animation.jl:98

display(1)

1

i=10
scatter(indIJ[:,2],indIJ[:,1],color=:black,label="",markeralpha = reverse(states[i],dims=1)'[:], aspect_ratio=:equal,xlim=(0,nely+2),ylim=(0,nelx+2),title="it: $i")
plot!(linkss[:,1,:]',linkss[:,2,:]',label="",color=:blue,linealpha = XIJEs[i]',linewidth = 1)

reverse