hca.jl

# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020

# based on https://curate.nd.edu/show/dv13zs28316

function hcastruc1(X, E, nu, p, tgS, s, scond, k, neighbors, bcond, itmax, vtol, ctype, ttype, T, xfilter, sfilter)
    states=[]
    anim=Animation()
    # INPUT PARAMETERS
    nelx, nely = 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);
    SED,U = fea(copy(X), p,E,nu)
    
    
    
    # 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")
    
    
    
    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 )
    # 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"
            #display("SED")
            #display(heatmap(SED',xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it"))
            #display("Xnew")
            #display(heatmap(Xnew',xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it"))
            #display("efX")
            #display(heatmap(efX',xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it"))
            #display("efE")
            #display(heatmap(efE',xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it"))
            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
        #println(efE);
        # println(Xnew);
        #if(minimum(Xnew)<0.0)
        #    Xnew.=Xnew.-minimum(Xnew).+s#amira added to remove negative values TODO make sure right ?? //or add xmin?
        #end
        Xnew[Xnew.<=0.0].=s #amira added to remove negative values TODO make sure right ?? //or add xmin?
        
        if(minimum(Xnew)>0)
            SEDnew,U = fea(copy(Xnew), p,E,nu); 
            # 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'))
            display(heatmap(Xnew',xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it"))
            heatmap(clamp.(Xnew',0,1),xaxis=nothing,yaxis=nothing,legend=nothing,c=:greys,title="it: $it")
            append!(states,[clamp.(Xnew',0,1)])
            frame(anim)
        else
            display("ERROR Xnew has negative values!!")
            it=itmax+1
        end
    end
    
    # SAVE FINAL RESULTS
    cells = Xnew;
    
    return anim,cells,states
    

end

#############################################
function  doavgf(c,neighbors,val)
    # DOAVGF(X, neighbors,val) sums around each X using fixed boundary
    # conditions. val is the fixed value.
    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
#-----------------------------------------------------------------
    
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

###############################################
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
#####################################
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
######################################
function fea(x, p,E,nu)
    nelx = size(x,1)+1;
    nely = size(x,2)+1;
    
    latticeSize=2
    setup=getSetup(latticeSize)
    maxNumTimeSteps=5000
    maxNumFiles=200


    saveEvery=round(maxNumTimeSteps/maxNumFiles)
    maxNumFiles=round(maxNumTimeSteps/saveEvery)-2
    setup["maxNumFiles"]=maxNumFiles
    save=false
    metavoxel,displacements=runMetavoxel2DGPULive!(setup,x,nelx,nely,p,E,nu,maxNumTimeSteps,saveEvery,maxNumFiles,save)
    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] ]
    
    
    
    SED = zeros(nelx-1,nely-1)
    U = zeros(nelx,nely,2)
    for ii in 0:(nelx-1)
        for jj in 0:(nely-1)
            i=((ii)*nely +(jj))+1
            
            
            
            if(ii<nelx-1&&jj<nely-1)
                Ue=[displacements[i].x;displacements[i].y; 
                displacements[((ii+1)*nely +(jj))+1].x;displacements[((ii+1)*nely +(jj))+1].y;
                displacements[((ii+1)*nely +(jj+1))+1].x;displacements[((ii+1)*nely +(jj+1))+1].y;
                displacements[((ii)*nely +(jj+1))+1].x;displacements[((ii)*nely +(jj+1))+1].y]
                
                SEDi = 1/2*x[ii+1,jj+1]^p*Ue'*KE*Ue;
                SED[ii+1,jj+1]=SED[ii+1,jj+1]+SEDi
            end
            
            
            
            U[ii+1,jj+1,1]=displacements[i].x
            U[ii+1,jj+1,2]=displacements[i].y
        end
    end
    # return SED
    return SED,U
end
#####################