# 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
#####################