Commit 62f42901 authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

2.5D optimization and visualization

parent 66ed8c84
Pipeline #19753 passed with stage
in 49 seconds
This diff is collapsed.
......@@ -65,6 +65,7 @@ include("./concurrent3D.jl")
# strain library
include("./microstructureDesignUe.jl")
include("./strain.jl");
include("./strain3D.jl");
# clustering
include("./clustering2D.jl");
......
......@@ -36,7 +36,7 @@ function Uclustering3DStrainLibrary(Macro_struct, Micro_struct,prob,library,verb
if Macro_Vol<1.0
Macro_xPhys,anim=topologyOptimization3d(Macro_nelx,Macro_nely,Macro_nelz,prob,Macro_Vol,Macro_rmin,penal,false)
end
scene= GLMakie.volume(Macro_xPhys, algorithm = :iso, isorange = 0.2, isovalue = 0.8,colormap=:grays)
scene= GLMakie.volume(permutedims(Macro_xPhys, [2, 1, 3]), algorithm = :iso, isorange = 0.4, isovalue = 1.0,colormap=:grays)
display(scene)
KE=lk_H8(nu)
......@@ -358,6 +358,343 @@ function Uclustering3DStrainLibrary(Macro_struct, Micro_struct,prob,library,verb
end
#based on direction not U
function Uclustering3D_new(θ,Macro_struct, Micro_struct,prob)
penal=3.0;Emin=1e-9;nu=0.3;E0=1;
Macro_length = Macro_struct[1]; Macro_width = Macro_struct[2]; Macro_Height = Macro_struct[3];
Micro_length = Micro_struct[1]; Micro_width = Micro_struct[2]; Micro_Height = Micro_struct[3];
Macro_nelx = Int(Macro_struct[4]); Macro_nely = Int(Macro_struct[5]); Macro_nelz = Int(Macro_struct[6]);
Micro_nelx = Int(Micro_struct[4]); Micro_nely = Int(Micro_struct[5]); Micro_nelz = Int(Micro_struct[6]);
Macro_Vol = Macro_struct[7][2];
Macro_Vol_FM = Macro_struct[7][1];
Micro_Vol=Micro_Vol_FM = Micro_struct[7][1];
Micro_rmin=Micro_struct[8]
Macro_rmin=Macro_struct[8]
Macro_max_change=Macro_struct[9];Micro_max_change=Micro_struct[9]
Macro_Elex = Macro_length/Macro_nelx; Macro_Eley = Macro_width/Macro_nely; Macro_Elez = Macro_Height/Macro_nelz;
Macro_nele = Macro_nelx*Macro_nely*Macro_nelz; Micro_nele = Micro_nelx*Micro_nely*Micro_nelz;Micro_nele2D = Micro_nelx*Micro_nely;
Macro_ndof = 3*(Macro_nelx+1)*(Macro_nely+1)*(Macro_nelz+1);
Macro_mgcg=mgcg[1];Micro_mgcg=mgcg[2];
# if length(voxels)==1
voxels=ones(Macro_nely,Macro_nelx,Macro_nelz)
# end
Macro_xPhys = ones(Macro_nely,Macro_nelx,Macro_nelz)
if Macro_Vol<1.0
Macro_xPhys,anim=topologyOptimization3d(Macro_nelx,Macro_nely,Macro_nelz,prob,Macro_Vol,Macro_rmin,penal,false)
end
scene= GLMakie.volume(permutedims(Macro_xPhys, [2, 1, 3]), algorithm = :iso, isorange = 0.4, isovalue = 1.0,colormap=:grays)
display(scene)
KE=lk_H8(nu)
# PREPARE FINITE ELEMENT ANALYSIS
U,F,freedofs=prob(Macro_nelx,Macro_nely,Macro_nelz);
fixeddofs = setdiff(1:Macro_ndof,freedofs)
# edofMat = repeat(3 .*nodeids[:] .+1,1,24).+repeat([0 1 2 3*Macro_nely.+[3 4 5 0 1 2] -3 -2 -1 3*(Macro_nelx+1)*(Macro_nely+1).+[0 1 2 3*Macro_nely.+[3 4 5 0 1 2] -3 -2 -1]], Macro_nele, 1);
nodenrs = reshape(1:(1+Macro_nelx)*(1+Macro_nely)*(1+Macro_nelz),1+Macro_nely,1+Macro_nelx,1+Macro_nelz)
edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1].+1,Macro_nelx*Macro_nely*Macro_nelz,1)
edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*Macro_nely.+[3 4 5 0 1 2] -3 -2 -1 3*(Macro_nely+1)*(Macro_nelx+1).+[0 1 2 3*Macro_nely.+[3 4 5 0 1 2] -3 -2 -1]],Macro_nelx*Macro_nely*Macro_nelz,1)
iK = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*Macro_nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*Macro_nele,1))
### IF MGCG
if Macro_mgcg
Macro_nl=4
Macro_Pu=Array{Any}(undef, Macro_nl-1);
for l = 1:Macro_nl-1
Macro_Pu[l,1] = prepcoarse3(Macro_nelz/2^(l-1),Macro_nely/2^(l-1),Macro_nelx/2^(l-1));
end
fixeddofs = setdiff(1:Macro_ndof,freedofs)
Macro_N = ones(Macro_ndof); Macro_N[fixeddofs] .= 0;
Macro_Null = spdiagm(Macro_ndof,Macro_ndof,0 => Macro_N);
end
### IF MGCG
if Macro_mgcg
K = Array{Any}(undef, Macro_nl);
sK = reshape(Kes.*(Emin.+Macro_xPhys[:]'.^penal*(1 .-Emin)),24*24*Macro_nele,1)
K[1,1] = sparse(vec(iK),vec(jK),vec(sK));
K[1,1] = Macro_Null'*K[1,1]*Macro_Null - (Macro_Null .-sparse(1.0I, Macro_ndof, Macro_ndof));
for l = 1:Macro_nl-1
K[l+1,1] = Macro_Pu[l,1]'*(K[l,1]*Macro_Pu[l,1]);
end
Lfac = factorize(Matrix(Hermitian(K[Macro_nl,1]))).L ;Ufac = Lfac';
cgtol=1e-10;
cgmax=100;
cgiters,cgres,U = MGCG(K,F[:,1],U[:,1],Lfac,Ufac,Macro_Pu,Macro_nl,1,cgtol,cgmax);
else
sK = reshape(KE[:]*(Emin .+Macro_xPhys[:]'.^penal*(1 .-Emin)),24*24*Macro_nele,1);
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2;
U[freedofs,:] = K[freedofs,freedofs]\Array(F[freedofs,:]);
end
ce = reshape(sum((U[edofMat]*KE).*U[edofMat],dims=2),Macro_nely,Macro_nelx,Macro_nelz)
c = sum(sum(sum((Emin.+Macro_xPhys.^penal*(E0-Emin)).*ce)))
display("compliance before $c")
##########clustering#########
Umat=U[edofMat];
ord=reshape(1:(Macro_nely*Macro_nelx*Macro_nelz),Macro_nely,Macro_nelx,Macro_nelz);
#todo plot strain directions
unit1=10
displayDeformation3D(Umat,Macro_nely,Macro_nelx,Macro_nelz,unit1,Macro_xPhys)
# diss=reshape(U_edofMat,Macro_nely, Macro_nelx, Macro_nelz,24);
ny=Macro_nely;nx=Macro_nelx;nz=Macro_nelz;
Θxy=zeros(ny ,nx, nz);ΘxyD=zeros(ny ,nx, nz);
Θz=zeros(ny ,nx, nz);ΘzD=zeros(ny ,nx, nz);
von_mises=zeros(ny , nx, nz);
ratio1=zeros(ny , nx, nz);ratio2=zeros(ny , nx, nz);
xs=ones(nz*ny*nx);ys=ones(nz*ny*nx);zs=ones(nz*ny*nx);vonmises=ones(nz*ny*nx);
count=1;
for k in 1:nz
for i in 1:ny
for j in 1:nx
if Macro_xPhys[i,j,k]>0.6
# Ue = diss[i,j,k,:];
Ue=Umat[ord[i,j,k],:];
Θxy[i,j,k],Θz[i,j,k],von_mises[i,j,k],ratio1[i,j,k],ratio2[i,j,k],E,T =getPrincipalStrains3D(Ue);
xs[count]=i;ys[count]=j;zs[count]=k;
vonmises[count]=von_mises[i,j,k];count+=1;
end
end
end
end
ΘxyD=Θxy.*180/π;
ΘzD=Θz.*180/π;
#plot vonmises
Makie.contour(von_mises, alpha=0.4)
display(current_figure())
rect = Rect(Vec(0.0, 0.0,0.0), Vec(1.0, 1.0,1.0))
mesh = GeometryBasics.mesh(rect)
Makie.meshscatter(nx.-ys .+1, ny.-xs .+1, zs, markersize = 1, marker=mesh,color = vonmises,
transparency=false,opacity=0.8,shading = true,shininess=32.0
,axis = (; type = Axis3, protrusions = (0, 0, 0, 0),
viewmode = :fit,aspect=:data, limits = (1, nx+1, 1, ny+1, 1, nz+1)))
Makie.inline!(true)
display(current_figure())
Makie.inline!(false)
display(current_figure())
# omg=ommh
# X=U[edofMat]'
# X=hcat(Macro_xPhys[:],U[edofMat])'
wd=2;
X=hcat(sin.(2*Θxy[:]).*Macro_xPhys[:],cos.(2*Θxy[:]).*Macro_xPhys[:],sin.(2*Θz[:]).*Macro_xPhys[:],cos.(2*Θz[:]).*Macro_xPhys[:],wd.*Macro_xPhys[:])'
# for i=1:3
# # display(i)
# # XXX=reshape(mean(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
# # XXX=reshape(sum(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
# XXX=reshape(U[edofMat]'[i,:],Macro_nely,Macro_nelx,Macro_nelz)
# # display(minimum(XXX))
# # display(maximum(XXX))
# # scene = Scene(resolution = (200, 200))
# scene= GLMakie.volume(permutedims(XXX, [2, 1, 3]), algorithm = :mip,transparency=true,colorrange=(minimum(XXX), maximum(XXX)),colormap=:lighttest)
# display(scene)
# # save("./img/FM3DU_$(Macro_Vol).png",scene)
# end
# # cluster X into θ clusters using K-means
# # R = kmeans(X, θ; maxiter=200, display=:iter)
# R = Clustering.kmeans(X, θ; maxiter=500)
# @assert Clustering.nclusters(R) == θ # verify the number of clusters
# a = Clustering.assignments(R) # get the assignments of points to clusters
# c = Clustering.counts(R) # get the cluster sizes
# M = R.centers; # get the cluster centers
# D=Distances.pairwise(Euclidean(), M, dims=2)
# hc = hclust(D, linkage=:single)
# display(Plots.plot(hc))
data=X
D=Distances.pairwise(Euclidean(), data, dims=2)
hc = hclust(D, linkage=:ward_presquared) # hc = hclust(D, linkage=:single,branchorder=:optimal)
display(Plots.plot(hc))
gr(size=(400,100))
#for elbow curve
hcElbow=reverse(hc.heights)[1:20]
display(Plots.plot(1:(length(hcElbow)),hcElbow,label=""))
hcElbow=reverse(hc.heights)[2:20]
xy=hcat(1:(length(hcElbow)),hcElbow)
mapp=mapproject(xy, dist2line=(line=[xy[1:1, :]; xy[end:end, :]], unit=:cartesian))
display(Plots.plot(mapp[:,3],label=""))
nclust=(findmax(mapp[:,3])[2]+1)
display("Best Number of clusters is= $nclust")
gr(size=(400,300))
nclust=θ
a=cutree(hc; k=nclust)
#visualize tree with only θ clusters
centers=zeros(size(data)[1],nclust)
for i in 1:nclust
R = Clustering.kmeans(data[:,a.==i], 1); # run K-means just to get the center
M = R.centers; # get the cluster centers
centers[:,i].=M[:]
end
D=Distances.pairwise(Euclidean(), centers, dims=2)
hc = hclust(D, linkage=:ward_presquared,branchorder=:optimal) # hc = hclust(D, linkage=:single,branchorder=:optimal)
display(Plots.plot(hc))
Uclustered=reshape(a,Macro_nely,Macro_nelx,Macro_nelz);
#############################
Macro_masksU=[]
Macro_masks_densU=[]
t=:darktest
cscheme=cgrad(t, θ, categorical = true)
cschemeList=[]
for i=1:θ
append!(cschemeList,[cscheme[i]])
end
cscheme=cgrad(cschemeList)
Macro_xPhys_diagU=copy(Macro_xPhys)
d=1/θ
for i in 1:θ
Macro_xPhys_temp=copy(Macro_xPhys)
out= (Uclustered.==i)
Macro_xPhys_temp[out].=true
Macro_xPhys_temp[.!out].=false
dens=Micro_Vol
println("Mask $i with density $(round(dens,digits=2)) occupies $(round(sum(out)/length(Macro_xPhys)*100,digits=2))% of the macrostructure.")
Macro_xPhys_diagU[out].= i
append!(Macro_masks_densU,[dens])
append!(Macro_masksU,[Macro_xPhys_temp])
# display(heatmap(Macro_xPhys_temp,showaxis = false,xaxis=nothing,yaxis=nothing,legend=nothing,fc=:grays,clims=(0.0, 1.0),aspect_ratio=:equal))
end
# scene = Scene(resolution = (400, 400))
scene= GLMakie.volume(permutedims(Macro_xPhys_diagU, [2, 1, 3]), algorithm = :absorption,transparency=true,colorrange=(1, θ),colormap=cscheme)
display(scene)
# save("./img/FM3DU1_$(Macro_Vol).png",scene)
# scene = Scene(resolution = (400, 400))#todo fix
clu=ones(nz*ny*nx);
count=1;
for k in 1:nz
for i in 1:ny
for j in 1:nx
if Macro_xPhys[i,j,k]>0.6
xs[count]=i;ys[count]=j;zs[count]=k;
clu[count]=Macro_xPhys_diagU[i,j,k];count+=1;
end
end
end
end
Makie.contour(Macro_xPhys_diagU, alpha=0.5,colormap=cscheme,colorrange=(1, θ))
display(current_figure())
# Makie.inline!(true)
# Makie.volume(permutedims(Macro_xPhys_diagU, [2, 1, 3]) , algorithm = :iso,colorrange=(1, θ),isorange = 1.0, isovalue = 1,colormap=cscheme)
# display(current_figure())
rect = Rect(Vec(0.0, 0.0,0.0), Vec(1.0, 1.0,1.0))
mesh = GeometryBasics.mesh(rect)
Makie.meshscatter(nx.-ys .+1, ny.-xs .+1, zs, markersize = 1, marker=mesh,color = clu,
transparency=false,opacity=0.8,shading = true,shininess=32.0,colormap=cscheme,colorrange=(1, θ)
,axis = (; type = Axis3, protrusions = (0, 0, 0, 0),
viewmode = :fit,aspect=:data, limits = (1, nx+1, 1, ny+1, 1, nz+1)))
Makie.inline!(true)
display(current_figure())
Makie.inline!(false)
display(current_figure())
Makie.inline!(true)
t=:darktest
cscheme=cgrad(t, θ, categorical = true)
cschemeList=[]
for i=1:θ
append!(cschemeList,[cscheme[i]])
end
cscheme11=cgrad(cschemeList)
Macro_xPhys_diagUUU=copy(Macro_xPhys_diagU)
Macro_xPhys_diagUUU[Macro_xPhys.<0.6].=NaN
fig = Figure(resolution=(400*ny, 400), fontsize=10)
for i in 1:ny
ax1=Makie.Axis(fig[1, i])
ax1.aspect = DataAspect()
Makie.heatmap!(ax1,Macro_xPhys_diagUUU[i,:,:],colormap=cscheme11,colorrange=(1, θ))
end
display(fig)
fig = Figure(resolution=(400*nx, 400), fontsize=10)
for i in 1:nx
ax1=Makie.Axis(fig[1, i])
ax1.aspect = DataAspect()
Makie.heatmap!(ax1,Macro_xPhys_diagUUU[:,i,:]',colormap=cscheme11,colorrange=(1, θ))
end
display(fig)
fig = Figure(resolution=(400*nz, 400), fontsize=10)
for i in 1:nz
ax1=Makie.Axis(fig[1, i])
ax1.aspect = DataAspect()
Makie.heatmap!(ax1,Macro_xPhys_diagUUU[:,:,i]',colormap=cscheme11,colorrange=(1, θ))
end
display(fig)
# omg=ommh
# for i=2:θ
# Makie.volume!(permutedims(Macro_xPhys_diagU, [2, 1, 3]) , algorithm = :iso,colorrange=(1, θ),isorange = 1.0, isovalue = i,colormap=cscheme)
# end
# display(current_figure())
# scene= volume(Macro_xPhys_diagU, algorithm = :mip,transparency=false,colorrange=(0.0, 1.0),colormap=:lighttest)
# display(scene)
# save("./img/FM3DFU_$(Macro_Vol).png",scene)
return Macro_masksU,Macro_masks_densU,Macro_xPhys_diagU,Macro_xPhys
end
##############################################################
# divide macrostructure into θ parts
function freeMaterial3D(θ,Macro_nely,Macro_nelx,Macro_nelz,Macro_nele,Macro_Vol,iK,jK,Emin,U,F,freedofs,edofMat,voxels=true)
......@@ -494,18 +831,18 @@ function Uclustering3D(θ,Macro_nely,Macro_nelx,Macro_nelz,Macro_nele,Macro_Vol,
# X=U[edofMat]'
X=hcat(Macro_xPhys[:],U[edofMat])'
for i=1:3
# display(i)
# XXX=reshape(mean(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
# XXX=reshape(sum(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
XXX=reshape(U[edofMat]'[i,:],Macro_nely,Macro_nelx,Macro_nelz)
# display(minimum(XXX))
# display(maximum(XXX))
# scene = Scene(resolution = (200, 200))
scene= GLMakie.volume(permutedims(XXX, [2, 1, 3]), algorithm = :mip,transparency=true,colorrange=(minimum(XXX), maximum(XXX)),colormap=:lighttest)
display(scene)
# save("./img/FM3DU_$(Macro_Vol).png",scene)
end
# for i=1:3
# # display(i)
# # XXX=reshape(mean(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
# # XXX=reshape(sum(U[edofMat]',dims=1),Macro_nely,Macro_nelx,Macro_nelz)
# XXX=reshape(U[edofMat]'[i,:],Macro_nely,Macro_nelx,Macro_nelz)
# # display(minimum(XXX))
# # display(maximum(XXX))
# # scene = Scene(resolution = (200, 200))
# scene= GLMakie.volume(permutedims(XXX, [2, 1, 3]), algorithm = :mip,transparency=true,colorrange=(minimum(XXX), maximum(XXX)),colormap=:lighttest)
# display(scene)
# # save("./img/FM3DU_$(Macro_Vol).png",scene)
# end
......@@ -1689,7 +2026,7 @@ function assembleFull3DStructure(θ,Macro_xPhys,Micro_xPhys,DHs,a,aYX,aZY,aZX,si
θTemp=Int(aYX[Int((zz-1)*2+2)][Int(yy),Int(xx)]) ; xPhys3D[:,:,end-thick]=Micro_xPhys[θTemp].*θTemp; xPhys3D_layers[:,:,end-thick,θTemp]=Micro_xPhys[θTemp];
θTemp=Int(aZY[Int((xx-1)*2+1)][Int(zz),Int(yy)]) ; xPhys3D[:,1+thick,:]=Micro_xPhys[θTemp]'.*θTemp; xPhys3D_layers[:,1+thick,:,θTemp]=Micro_xPhys[θTemp]';
θTemp=Int(aZY[Int((xx-1)*2+2)][Int(zz),Int(yy)]) ; xPhys3D[:,end-thick,:]=Micro_xPhys[θTemp]'.*θTemp; xPhys3D_layers[:,end-thick,:,θTemp]=Micro_xPhys[θTemp]';
θTemp=Int(aZX[Int((yy-1)*2+1)][Int(zz),Int(xx)]) ; xPhys3D[1+thick,:,:]=Micro_xPhys[θTemp]'.*θTemp; xPhys3D_layers[1+thick,:,:,θTemp]=Micro_xPhys[θTemp];
θTemp=Int(aZX[Int((yy-1)*2+1)][Int(zz),Int(xx)]) ; xPhys3D[1+thick,:,:]=Micro_xPhys[θTemp]'.*θTemp; xPhys3D_layers[1+thick,:,:,θTemp]=Micro_xPhys[θTemp]';
θTemp=Int(aZX[Int((yy-1)*2+2)][Int(zz),Int(xx)]) ; xPhys3D[end-thick,:,:]=Micro_xPhys[θTemp]'.*θTemp; xPhys3D_layers[end-thick,:,:,θTemp]=Micro_xPhys[θTemp]';
end
xPhys3D1=reshape(xPhys3D,sizeMicro,sizeMicro,sizeMicro,1,1)
......
......@@ -193,6 +193,79 @@ function elementMatVec3D(a, b, c, DH)
return Ke
end
function brick_stiffnessMatrix()
# elastic matrix formulation
nu=0.3;
D = 1.0 /((1+nu)*(1-2*nu))*[1-nu nu nu 0 0 0; nu 1-nu nu 0 0 0;
nu nu 1-nu 0 0 0; 0 0 0 (1-2*nu)/2 0 0; 0 0 0 0 (1-2*nu)/2 0;
0 0 0 0 0 (1-2*nu)/2];
#stiffness matrix formulation
A = [32 6 -8 6 -6 4 3 -6 -10 3 -3 -3 -4 -8;
-48 0 0 -24 24 0 0 0 12 -12 0 12 12 12];
k = 1/144*A'*[1; nu];
K1 = [k[1] k[2] k[2] k[3] k[5] k[5];
k[2] k[1] k[2] k[4] k[6] k[7];
k[2] k[2] k[1] k[4] k[7] k[6];
k[3] k[4] k[4] k[1] k[8] k[8];
k[5] k[6] k[7] k[8] k[1] k[2];
k[5] k[7] k[6] k[8] k[2] k[1]];
K2 = [k[9] k[8] k[12] k[6] k[4] k[7];
k[8] k[9] k[12] k[5] k[3] k[5];
k[10] k[10] k[13] k[7] k[4] k[6];
k[6] k[5] k[11] k[9] k[2] k[10];
k[4] k[3] k[5] k[2] k[9] k[12]
k[11] k[4] k[6] k[12] k[10] k[13]];
K3 = [k[6] k[7] k[4] k[9] k[12] k[8];
k[7] k[6] k[4] k[10] k[13] k[10];
k[5] k[5] k[3] k[8] k[12] k[9];
k[9] k[10] k[2] k[6] k[11] k[5];
k[12] k[13] k[10] k[11] k[6] k[4];
k[2] k[12] k[9] k[4] k[5] k[3]];
K4 = [k[14] k[11] k[11] k[13] k[10] k[10];
k[11] k[14] k[11] k[12] k[9] k[8];
k[11] k[11] k[14] k[12] k[8] k[9];
k[13] k[12] k[12] k[14] k[7] k[7];
k[10] k[9] k[8] k[7] k[14] k[11];
k[10] k[8] k[9] k[7] k[11] k[14]];
K5 = [k[1] k[2] k[8] k[3] k[5] k[4];
k[2] k[1] k[8] k[4] k[6] k[11];
k[8] k[8] k[1] k[5] k[11] k[6];
k[3] k[4] k[5] k[1] k[8] k[2];
k[5] k[6] k[11] k[8] k[1] k[8];
k[4] k[11] k[6] k[2] k[8] k[1]];
K6 = [k[14] k[11] k[7] k[13] k[10] k[12];
k[11] k[14] k[7] k[12] k[9] k[2];
k[7] k[7] k[14] k[10] k[2] k[9];
k[13] k[12] k[10] k[14] k[7] k[11];
k[10] k[9] k[2] k[7] k[14] k[7];
k[12] k[2] k[9] k[11] k[7] k[14]];
KE = 1/((nu+1)*(1-2*nu))*[ K1 K2 K3 K4;
K2' K5 K6 K3';
K3' K6 K5' K2';
K4 K3 K2 K1'];
# strain matrix formulation
B_1=[-0.044658 0 0 0.044658 0 0 0.16667 0;
0 -0.044658 0 0 -0.16667 0 0 0.16667;
0 0 -0.044658 0 0 -0.16667 0 0;
-0.044658 -0.044658 0 -0.16667 0.044658 0 0.16667 0.16667;
0 -0.044658 -0.044658 0 -0.16667 -0.16667 0 -0.62201;
-0.044658 0 -0.044658 -0.16667 0 0.044658 -0.62201 0];
B_2=[0 -0.16667 0 0 -0.16667 0 0 0.16667;
0 0 0.044658 0 0 -0.16667 0 0;
-0.62201 0 0 -0.16667 0 0 0.044658 0;
0 0.044658 -0.16667 0 -0.16667 -0.16667 0 -0.62201;
0.16667 0 -0.16667 0.044658 0 0.044658 -0.16667 0;
0.16667 -0.16667 0 -0.16667 0.044658 0 -0.16667 0.16667];
B_3=[0 0 0.62201 0 0 -0.62201 0 0;
-0.62201 0 0 0.62201 0 0 0.16667 0;
0 0.16667 0 0 0.62201 0 0 0.16667;
0.16667 0 0.62201 0.62201 0 0.16667 -0.62201 0;
0.16667 -0.62201 0 0.62201 0.62201 0 0.16667 0.16667;
0 0.16667 0.62201 0 0.62201 0.16667 0 -0.62201];
B=[B_1 B_2 B_3];
return KE,B,D
end
####################################
## SUB FUNCTION: EBHM2D
function EBHM2D(den, lx, ly, E0, Emin, nu, penal)
......
......@@ -46,53 +46,54 @@ function addFabricationConstraints(ss,Emin,Micro_x,Micro_nelx,Micro_nely,Micro_n
return Micro_x
end
function addFabricationConstraints(ss,Emin,Micro_x,Micro_nelx,Micro_nely,Micro_nelz)
function addFabricationConstraints(ss,Emin,Micro_x,Micro_nelx,Micro_nely,Micro_nelz,discrete=false)
ss2=Int(ss/2)
ss4=1
###
Micro_x[1:ss2,:,:].=0.0;
Micro_x[Micro_nely-ss2+1,:,:].=0.0;
Micro_x[1:ss4,:,:].=0.0;
Micro_x[Micro_nely-ss4+1,:,:].=0.0;
Micro_x[:,1:ss2,:].=0.0;
Micro_x[:,Micro_nelx-ss2+1,:].=0.0;
Micro_x[:,1:ss4,:].=0.0;
Micro_x[:,Micro_nelx-ss4+1,:].=0.0;
Micro_x[:,:,1:ss2].=0.0;
Micro_x[:,:,Micro_nelz-ss2+1].=0.0;
Micro_x[:,:,1:ss4].=0.0;
Micro_x[:,:,Micro_nelz-ss4+1].=0.0;
###############
inds=findall(x-> x>=0, Micro_x)
for ind in inds
# display(ind)
X=ind[1];Y=ind[2];Z=ind[3]*2;
Z1=Micro_nely*2-ind[3]*2;
if min((Micro_nely+1)-Z,min(Z-(0),min((Micro_nely+1)-Y-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(Y-(0)-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(X-(0)-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,(Micro_nely+1)-X-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2))))) ==0
Micro_x[ind]=0.0
elseif min((Micro_nely+1)-Z1,min(Z1-(0),min((Micro_nely+1)-Y-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(Y-(0)-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(X-(0)-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,(Micro_nely+1)-X-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2))))) ==0
if(ind[3]<Micro_nely)
Micro_x[X,Y,ind[3]+1]=0.0
if discrete
inds=findall(x-> x>=0, Micro_x)
for ind in inds
# display(ind)
X=ind[1];Y=ind[2];Z=ind[3]*2;
Z1=Micro_nely*2-ind[3]*2;
if min((Micro_nely+1)-Z,min(Z-(0),min((Micro_nely+1)-Y-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(Y-(0)-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(X-(0)-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,(Micro_nely+1)-X-((Z-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2))))) ==0
Micro_x[ind]=Emin
elseif min((Micro_nely+1)-Z1,min(Z1-(0),min((Micro_nely+1)-Y-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(Y-(0)-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,min(X-(0)-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2,(Micro_nely+1)-X-((Z1-(0))/((Micro_nely+1)-(0)))*((Micro_nely+1)-(0))/2))))) ==0
if(ind[3]<Micro_nely)
Micro_x[X,Y,ind[3]+1]=Emin
end