Commit 4f23fda7 authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

fixed 3d

parent 33030734
Pipeline #19659 passed with stage
in 31 seconds
This diff is collapsed.
This diff is collapsed.
......@@ -412,8 +412,7 @@ end
# multiple microstructures (θ)
function MultiConTop3DU(θ,Macro_struct, Micro_struct,prob, penal,saveItr=5,maxloop = 90,fabric=false,mgcg=[false,false],voxels=true)
# USER-DEFINED LOOP PARAMETERS
E0 = 1; Emin = 1e-9; nu = 0.3;
......
......@@ -448,4 +448,6 @@ function evaluateCH(CH,dens)
return SH
end
\ No newline at end of file
end
#####################################
\ No newline at end of file
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
# (c) Massachusetts Institute of Technology 2022
#======================================================================================================================#
###################################################################
......@@ -944,7 +944,122 @@ function MicroTop2DU_max(Ue,Micro_struct, maxParameters,saveItr=5,maxloop = 200,
end
###################################################################
# cuboct vs strain library vs full 3D vs vs 2D3D
function MicroTop3DU(Ue, Micro_struct,library, penal=3,mgcg=[false,true],size3D=20,saveItr=5,maxloop = 200,fab=false)
unit=10;
displayElementDeformation3D(Ue,unit)
Macro_length = 1; Macro_width = 1; Macro_Height = 1;
Micro_length = Micro_struct[1]; Micro_width = Micro_struct[2]; Micro_Height = Micro_struct[3];
Macro_nelx = 1; Macro_nely = 1; Macro_nelz = 1;
Micro_nelx = Int(Micro_struct[4]); Micro_nely = Int(Micro_struct[5]); Micro_nelz = Int(Micro_struct[6]);
Micro_Vol=Micro_Vol_FM = Micro_struct[7][1];
Micro_rmin=Micro_struct[8]
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];
Macro_xPhys = ones(Macro_nely,Macro_nelx,Macro_nelz)
E0 = 1; Emin = 1e-9; nu = 0.3;penal=3;
###################################################################
#cuboct
cuboct=load("./img/library/1/cuboct_vol_$((Micro_Vol)).jld")["data"];
cuboctDH, dDH = EBHM2D(cuboct, Micro_length, Micro_width, E0, Emin, nu, penal);
# cuboct3D ,cuboctDH3=visualize3D2DMicrostructure(cuboct,false)
###################################################################
#library
Umat=zeros(1,8*3)
Umat[1,:].=Ue
dissYX,dissZY,dissZX,X,Macro_xPhys_2D=reshapeU2DSlices_new(Umat,Macro_nely, Macro_nelx, Macro_nelz,Macro_xPhys,true);
order=[1,2,3,4,5,6,7,8];
n=size(X)[1];
Θ=zeros(n);ΘD=zeros(n);
von_mises=zeros(n);
ratio=zeros(n);
for i=1:n
Ue = X[i];
Θ[i] , von_mises[i], ratio[i],α1,α2 =getPrincipalStrains(Ue);
end
ΘD=Θ.*180/π;
E0 = 1; Emin = 1e-9; nu = 0.3;penal=3;
θ=Int((num-1)*4+1)+1
Micro_xPhys,DHs=loadVisualizeLibrary(library)
Micro_x=Micro_xPhys
t=:darktest
cscheme=cgrad(t, θ, categorical = true)
cschemeList=[]
append!(cschemeList,[:white])
for i=1:θ
append!(cschemeList,[cscheme[i]])
end
# cscheme=ColorGradient(cschemeList)
cscheme1=cgrad(cschemeList)
function findNearestList(theta)
findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2]
values=[]
for i in range(0.0, 180, length=θ)
append!(values,[i])
end
return values[findnearest(values,theta)]
end
function findNearestListIndex(theta)
findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2]
values=[]
for i in range(0.0, 180, length=θ)
append!(values,[i])
end
return findnearest(values,theta)
end
ΘDNearest=findNearestList.(ΘD)
Uclustered=findNearestListIndex.(ΘD)
a=Uclustered
display(a[1])
display(θ)
for i=1:n
Ue = X[i];
if sum(Ue)==0
a[i]=θ #if no load cuboct
end
end
######
aYX,aZY,aZX=reshape_back_reshapeU2DSlices_new(θ,a,Macro_nely, Macro_nelx, Macro_nelz,true)
###################################################################
θTempList="";DH="";
yy=1;xx=1;zz=1
xPhys3D=zeros(Micro_nely,Micro_nelx,Micro_nelz);
θTemp=Int(aYX[Int((zz-1)*2+1)][Int(yy),Int(xx)]) ;xPhys3D[:,:,1:thickness]=Micro_x[θTemp]';
θTemp=Int(aYX[Int((zz-1)*2+2)][Int(yy),Int(xx)]) ;xPhys3D[:,:,end-thickness+1:end]=Micro_x[θTemp]';
θTemp=Int(aZY[Int((xx-1)*2+1)][Int(zz),Int(yy)]) ;xPhys3D[:,1:thickness,:]=Micro_x[θTemp]';
θTemp=Int(aZY[Int((xx-1)*2+2)][Int(zz),Int(yy)]) ;xPhys3D[:,end-thickness+1:end,:]=Micro_x[θTemp]';
θTemp=Int(aZX[Int((yy-1)*2+1)][Int(zz),Int(xx)]) ;xPhys3D[1:thickness,:,:]=Micro_x[θTemp]';
θTemp=Int(aZX[Int((yy-1)*2+2)][Int(zz),Int(xx)]) ;xPhys3D[end-thickness+1:end,:,:]=Micro_x[θTemp]';
xPhys3D1=reshape(xPhys3D,Micro_nely,Micro_nelx,Micro_nelz,1,1)
xPhys3D1=reshape(Flux.AdaptiveMaxPool((size3D, size3D,size3D))(xPhys3D1),size3D,size3D,size3D)
scene= GLMakie.volume( permutedims(xPhys3D1, [2, 1, 3]), algorithm = :iso, isorange = 0.3, isovalue = 1.0,colormap=:grays)
display(scene)
DH, dDH = EBHM3D(xPhys3D1, Micro_length, Micro_width, Micro_Height, E0, Emin, nu, penal,Micro_mgcg);
dispFlag = 1;plotFlag = 1;outOption = "struct";
dens=sum(xPhys3D1)/(size(xPhys3D1)[1]*size(xPhys3D1)[2]*size(xPhys3D1)[3])
props, SH = evaluateCH(DH, dens, outOption, dispFlag);
fig=visual(DH);
# Kes[:,yy,xx,zz].=elementMatVec3D(Macro_Elex/2, Macro_Eley/2, Macro_Elez/2, DH)[:]
#
return xPhys3D1
end
###################################################################
\ No newline at end of file
......@@ -52,6 +52,94 @@ function plotDisplacement(nelx,nely,nelz,xPhysC,getProblem=inverter,factor=0.04,
end
function plotDisplacement3D(nelx,nely,nelz,prob,volfrac,rmin,penal,factor=0.04,cameraX=30,cameraY=60)
# Max and min stiffness
Emin=1e-9
Emax=1.0
nu=0.3
border=4
ndof = 3*(nelx+1)*(nely+1)*(nelz+1)
# Allocate design variables (as array), initialize and allocate sens.
x= ones(Float64,nely,nelx,nelz)
xold=copy(x)
xPhys=copy(x)
g=0 # must be initialized to use the NGuyen/Paulino OC approach
dc=zeros(Float64,(nely,nelx,nelz))
# FE: Build the index vectors for the for coo matrix format.
KE=lk_H8(nu)
nele = nelx*nely*nelz
nodenrs = reshape(1:(1+nelx)*(1+nely)*(1+nelz),1+nely,1+nelx,1+nelz)
edofVec = reshape(3*nodenrs[1:end-1,1:end-1,1:end-1].+1,nelx*nely*nelz,1)
edofMat = repeat(edofVec,1,24).+repeat([0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1 3*(nely+1)*(nelx+1).+[0 1 2 3*nely.+[3 4 5 0 1 2] -3 -2 -1]],nelx*nely*nelz,1)
iK = convert(Array{Int},reshape(kron(edofMat,ones(24,1))',24*24*nele,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,24))',24*24*nele,1))
U,F,freedofs=prob(nelx,nely,nelz)
# FE-ANALYSIS
sK = reshape(KE[:]*(Emin.+xPhys[:]'.^penal*(Emax-Emin)),24*24*nelx*nely*nelz,1)
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
@timed U[freedofs] = K[freedofs,freedofs] \ Array(F[freedofs])
# Objective function and sensitivity analysis
ce = reshape(sum((U[edofMat]*KE).*U[edofMat],dims=2),nely,nelx,nelz)
c = sum(sum(sum((Emin.+xPhys.^penal*(Emax-Emin)).*ce)))
Umat=U[edofMat];
ord=reshape(1:(nely*nelx*nelz),nely,nelx,nelz);
gr(size=(800,600))
function getYX3(el3D)
el3D=reshape(el3D,3,8)'
Ue1=[el3D[1,:] el3D[2,:] el3D[3,:] el3D[4,:]]
Ue2=[el3D[5,:] el3D[6,:] el3D[7,:] el3D[8,:]]
return reshape(Ue1,3,4)',reshape(Ue2,3,4)'
end
# cameraX=60 # (azimuthal, elevation)
# cameraY=30
# border=4
# cameraX=90 # (azimuthal,
# cameraY=0 # elevation)
anim1=Animation()
ix=[];iy=[];iz=[];
for step in 0:0.5:10
ix=[];iy=[];iz=[];
exg=factor*step
for k in 1:nelz
for i in 1:nely
for j in 1:nelx
Ue1,Ue2=getYX3(Umat[ord[i,j,k],:])
append!(iy,-i+exg*Ue1[4,2]); append!(ix,j+exg*Ue1[4,1]);append!(iz,k+exg*Ue1[4,3]);
append!(iy,-(i+1)+exg*Ue1[2,2]);append!(ix,(j+1)+exg*Ue1[2,1]);append!(iz,k+exg*Ue1[2,3]);
append!(iy,-(i)+exg*Ue1[3,2]); append!(ix,(j+1)+exg*Ue1[3,1]);append!(iz,k+exg*Ue1[3,3]);
append!(iy,-(i+1)+exg*Ue1[1,2]);append!(ix,(j)+exg*Ue1[1,1]);append!(iz,k+exg*Ue1[1,3]);
append!(iy,-i+exg*Ue2[4,2]); append!(ix,j+exg*Ue2[4,1]);append!(iz,(k+1)+exg*Ue2[4,3]);
append!(iy,-(i+1)+exg*Ue2[2,2]);append!(ix,(j+1)+exg*Ue2[2,1]);append!(iz,(k+1)+exg*Ue2[2,3]);
append!(iy,-(i)+exg*Ue2[3,2]); append!(ix,(j+1)+exg*Ue2[3,1]);append!(iz,(k+1)+exg*Ue2[3,3]);
append!(iy,-(i+1)+exg*Ue2[1,2]);append!(ix,(j)+exg*Ue2[1,1]);append!(iz,(k+1)+exg*Ue2[1,3]);
end
end
end
Plots.scatter(iy,ix,iz,color="black",label="",markersize =4, yaxis ="xxx", xaxis ="yyy", zaxis ="zzz",aspect_ratio=:equal,markerstrokealpha = 0.2,markeralpha = 0.6,markershape = :square,camera = (cameraX,cameraY),xlim=(-border-nely,border),ylim=(-border,nelx+border),zlim=(-border,nelz+border))#,markershape = :square
frame(anim1)
end
return anim1
gif(anim1, "anim.gif", fps = 10)
end
function make_bitmap(p,nx,ny,alpha)
color = [1 0 0; 0 0 .45; 0 1 0; 0 0 0; 1 1 1];
I = zeros(nx*ny,3);
......
......@@ -172,8 +172,6 @@ function microStructureConnection_normal2(nelx,nely)
end
# cuboct
function microStructureConnection_shear1(nelx,nely)
ndof= Int(2*(nelx+1)*(nely+1));
......@@ -566,6 +564,39 @@ function stool(nelx,nely,nelz)
return U,F,freedofs
end
function stressTest3D(nelx,nely,nelz)
ndof=3*(nelx+1)*(nely+1)*(nelz+1);
il = [Int(nelx/2.)];
jl = [ 0];
kl = [Int(nelz/2.)];
loadnid = kl.*(nelx+1).*(nely+1).+il.*(nely+1).+(nely+1 .-jl)
loaddof = 3 .*loadnid[:] .- 1;
val=-1;
iif = [0 0 nelx nelx];
jf = [0 0 0 0];
kf = [0 nelz 0 nelz];
fixednid = kf.*(nelx+1).*(nely+1) .+iif .*(nely .+1) .+(nely .+1 .-jf) # Node IDs
fixeddof = [3 .*fixednid[:]; 3 .*fixednid[:].-1; 3 .*fixednid[:].-2]
F= sparse(loaddof,fill(1,length(loaddof)),fill(val,length(loaddof)),ndof,1);
U = zeros(ndof,1)
alldofs = 1:ndof
freedofs = sort(setdiff(1:ndof,fixeddof))
return U,F,freedofs
end
function stressTest3D(nelx,nely,nelz)
ndof=3*(nelx+1)*(nely+1)*(nelz+1);
......@@ -1345,10 +1376,10 @@ function seat3DM(nelx,nely,nelz)
fixeddof_s = 3*fixednid[:];
# # four pins
# iif = [0;0;0;0];
# jf = [0;1;nely-1 ;nely];
# kf = [nelz;nelz;nelz;nelz];
# # four pins
# iif = [0;0;0;0];
# jf = [0;1;nely-1 ;nely];
# kf = [nelz;nelz;nelz;nelz];
m = Matlab.meshgrid([0],0:nely)
iif=m[1]
......
......@@ -342,8 +342,6 @@ function getShapeStrainDisplacement(Ue,x,y)
# end
end
###################################################################
## SUB FUNCTION: display_3D
function display_3D(rho)
......@@ -369,7 +367,116 @@ function display_3D(rho)
# axis equal; axis tight; axis off; box on; view( [30,30] ); pause(1e-6);
end
function displayElementDeformation3D(Ue,unit)
Ue1 =[
unit,0,0, #1
unit,unit,0, #2
0,unit,0, #3
0,0,0, #4
unit,0,unit, #5
unit,unit,unit, #6
0,unit,unit, #7
0,0,unit, #8
]
Uee1=reshape(Ue1,3,8)
fig = Figure(resolution=(300, 300), fontsize=10)
pointsX=[];pointsY=[];pointsZ=[];
for i =1:4
append!(pointsY,[Uee1[1,i]]);append!(pointsX,[Uee1[2,i]]);append!(pointsZ,[Uee1[3,i]]);
end
append!(pointsY,[Uee1[1,1]]);append!(pointsX,[Uee1[2,1]]);append!(pointsZ,[Uee1[3,1]]);
# Makie.lines(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ))
for i =5:8
append!(pointsY,[Uee1[1,i]]);append!(pointsX,[Uee1[2,i]]);append!(pointsZ,[Uee1[3,i]])
end
append!(pointsY,[Uee1[1,5]]);append!(pointsX,[Uee1[2,5]]);append!(pointsZ,[Uee1[3,5]]);
Makie.lines(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:black, linestyle = :dash,figure=(resolution=(400, 400), fontsize=10),axis = (; type = Axis3, protrusions = (0, 0, 0, 0),
viewmode = :fit, limits = (-unit*0.5, unit+unit*0.5, -unit*0.5, unit+unit*0.5, -unit*0.5, unit+unit*0.5)))
pointsX=[];pointsY=[];pointsZ=[];
append!(pointsY,[Uee1[1,4]]);append!(pointsX,[Uee1[2,4]]);append!(pointsZ,[Uee1[3,4]]);
append!(pointsY,[Uee1[1,8]]);append!(pointsX,[Uee1[2,8]]);append!(pointsZ,[Uee1[3,8]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:black, linestyle = :dash)
pointsX=[];pointsY=[];pointsZ=[];
append!(pointsY,[Uee1[1,2]]);append!(pointsX,[Uee1[2,2]]);append!(pointsZ,[Uee1[3,2]]);
append!(pointsY,[Uee1[1,6]]);append!(pointsX,[Uee1[2,6]]);append!(pointsZ,[Uee1[3,6]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:black, linestyle = :dash)
pointsX=[];pointsY=[];pointsZ=[];
append!(pointsY,[Uee1[1,3]]);append!(pointsX,[Uee1[2,3]]);append!(pointsZ,[Uee1[3,3]]);
append!(pointsY,[Uee1[1,7]]);append!(pointsX,[Uee1[2,7]]);append!(pointsZ,[Uee1[3,7]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:black, linestyle = :dash)
Ueee=reshape(Ue,3,8)
Uee1=reshape(Ue1,3,8)
Uee1.=Uee1.+Ueee[[2,1,3],:]
pointsX=[];pointsY=[];pointsZ=[];
for i =1:4
append!(pointsY,[Uee1[1,i]]);append!(pointsX,[Uee1[2,i]]);append!(pointsZ,[Uee1[3,i]]);
end
append!(pointsY,[Uee1[1,1]]);append!(pointsX,[Uee1[2,1]]);append!(pointsZ,[Uee1[3,1]]);
# Makie.lines(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ))
for i =5:8
append!(pointsY,[Uee1[1,i]]);append!(pointsX,[Uee1[2,i]]);append!(pointsZ,[Uee1[3,i]])
end
append!(pointsY,[Uee1[1,5]]);append!(pointsX,[Uee1[2,5]]);append!(pointsZ,[Uee1[3,5]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:blue)
pointsX=[];pointsY=[];pointsZ=[];
append!(pointsY,[Uee1[1,4]]);append!(pointsX,[Uee1[2,4]]);append!(pointsZ,[Uee1[3,4]]);
append!(pointsY,[Uee1[1,8]]);append!(pointsX,[Uee1[2,8]]);append!(pointsZ,[Uee1[3,8]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:blue)
pointsX=[];pointsY=[];pointsZ=[];
append!(pointsY,[Uee1[1,2]]);append!(pointsX,[Uee1[2,2]]);append!(pointsZ,[Uee1[3,2]]);
append!(pointsY,[Uee1[1,6]]);append!(pointsX,[Uee1[2,6]]);append!(pointsZ,[Uee1[3,6]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:blue)
pointsX=[];pointsY=[];pointsZ=[];
append!(pointsY,[Uee1[1,3]]);append!(pointsX,[Uee1[2,3]]);append!(pointsZ,[Uee1[3,3]]);
append!(pointsY,[Uee1[1,7]]);append!(pointsX,[Uee1[2,7]]);append!(pointsZ,[Uee1[3,7]]);
Makie.lines!(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ),color =:blue)
display(current_figure())
end
function displayElementDeformation2D(Ue,unit)
Ue1 =[
unit,1, #1
unit,unit, #2
0,unit, #3
0,0 #4
]
Uee1=reshape(Ue1,2,4)
fig = Figure(resolution=(300, 300), fontsize=10)
pointsX=[];pointsY=[];
for i =1:4
append!(pointsY,[Uee1[1,i]]);append!(pointsX,[Uee1[2,i]]);
end
append!(pointsY,[Uee1[1,1]]);append!(pointsX,[Uee1[2,1]]);
Makie.lines(Float16.(pointsX), Float16.(pointsY),color =:black, linestyle = :dash,figure=(resolution=(300, 300), fontsize=10),axis = (; limits = (-unit, 2*unit, -unit, 2*unit)))
Ueee=reshape(Ue,2,4)
Uee1=reshape(Ue1,2,4)
Uee1.=Uee1.+Ueee[[2,1],:]
pointsX=[];pointsY=[];
for i =1:4
append!(pointsY,[Uee1[1,i]]);append!(pointsX,[Uee1[2,i]]);
end
append!(pointsY,[Uee1[1,1]]);append!(pointsX,[Uee1[2,1]]);
# Makie.lines(Float16.(pointsX), Float16.(pointsY),Float16.(pointsZ))
Makie.lines!(Float16.(pointsX), Float16.(pointsY),color =:blue)
display(current_figure())
end
###################################################################
......@@ -433,11 +540,19 @@ function loadVisualizeLibrary(library)
append!(DHs,[[]])
append!(Micro_xPhys,[fetchMicrostructureConnect(i,library)])
end
cuboct=load("./img/library/1/cuboct_vol_$((Micro_Vol)).jld")["data"];
append!(DHs,[[]])
append!(Micro_xPhys,[cuboct])
θ+=1
E0 = 1; Emin = 1e-9; nu = 0.3;penal=3;Micro_length=1.0;Micro_width=1.0;
for i=1:θ
DH, dDH = EBHM2D(Micro_xPhys[i], Micro_length, Micro_width, E0, Emin, nu, penal);
DHs[i]=DH;
end
# t=:dark
t=:darktest
# # theme(t);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment