Commit 711f0ad6 authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

min max top opt

parent 3f025104
Pipeline #19573 passed with stage
in 1 minute and 14 seconds
......@@ -15,4 +15,10 @@ live
*.png
*.mp4
img
res3d
\ No newline at end of file
res3d
*.stlbak
*.dxfbak
*.3mf
*bracket.3dm
*new_face_75.3dm
*.gcode
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -96718,15 +96718,15 @@
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.3",
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.5"
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
"version": "1.6.2"
}
},
"nbformat": 4,
This diff is collapsed.
This diff is collapsed.
......@@ -439,6 +439,8 @@ end
###############################################
# input the homogenized elasticity tensor
# fig=visual(CH)
# save("normal.png", fig)
function visual(CH)
# transform it to 3*3*3*3 tensor
tensor = generate(CH);
......@@ -467,8 +469,18 @@ function visual(CH)
end
x,y,z = sph2cart(a,e,E1);
c = sqrt.(x.^2 .+ y.^2 +z.^2)
display(surface(x,y,z ,color=c))
return s
# display(surface(x,y,z ,color=c))
fig = Figure(resolution=(600, 600), fontsize=10)
ax=Axis3(fig[1, 1],aspect = :data,viewmode =:fit,perspectiveness = 0.0,textsize = 10)
hm =surface!(ax,x,y,z ,color=c,opacity=1.0 ,shading=false,colormap = Reverse(:deep))
Colorbar(fig[1, 2], hm,height=Relative(0.5),ticksize=15,)
display(fig)
# save("normal.png", fig)
return fig
end
function modulus(CH)
......
......@@ -1513,6 +1513,27 @@ function OC2D(x, dc, dv, H, Hs, volfrac, nele, move, beta,fab=false,ss=4)
return x, xPhys, change
end
function OC2D_2(x, dc, dv, H, Hs, volfrac, nele, move, beta,fab=false,ss=4)
l1 = 0; l2 = 1e9;
xTilde=copy(x);xnew=copy(x);xPhys=copy(x)
while (l2-l1)/(l1+l2) > 1e-4
lmid = 0.5*(l2+l1);
xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(abs.(-dc./dv/lmid))))));
xTilde[:] = (H*xnew[:] )./Hs;
xPhys = 1 .-exp.(-beta .*xTilde) .+xTilde .*exp.(-beta);
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
if fab
xnew=addFabConstraint2D(xnew,ss);
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x,xTilde, xPhys, change
end
## SUB FUNCTION: MMA
function OC2DMMA(x, xold1,xold2, dc, dv, H, Hs, volfrac, nele, move, beta,loop,low,upp)
......
......@@ -2135,12 +2135,12 @@ function SlicesUclustering3D(θ,Macro_nely,Macro_nelx,Macro_nelz,Macro_nele,Macr
# println("FMD It.:$loop Obj.:$(round.(c,digits=2)) Macro_Vol.:$(round.(mean(Macro_xPhys[:] ),digits=2)) Macro_ch.:$(round.(Macro_change,digits=2))")
end
AbstractPlotting.inline!(true)
# AbstractPlotting.inline!(true)
scene = Scene(resolution = (200, 200))
scene= volume!(Macro_xPhys, algorithm = :mip,transparency=true,colorrange=(0.0, 1.0),colormap=:lighttest)
display(scene)
AbstractPlotting.inline!(true)
# AbstractPlotting.inline!(true)
scene = Scene(resolution = (200, 200))
# scene= volume!(Macro_xPhys, algorithm = :mip,transparency=true,colorrange=(0.0, 1.0),colormap=:lighttest)
scene= volume!(Macro_xPhys, algorithm = :iso,colorrange=(0.0, 1.0),isorange = 0.2, isovalue = 0.8)
......
......@@ -52,12 +52,16 @@ function addFabConstraint2D(x,ss)
ss2=Int(ss/2)
ss4=2
# sss=Int(nelx/4)
#sides empty
# x[1:ss,1:ss].=0.0;
# x[1:ss,nelx-ss:nelx].=0.0;
# x[nely-ss:nely,1:ss].=0.0;
# x[nely-ss:nely,nelx-ss:nelx].=0.0;
# x[1:sss,1:sss].=0.0;
# x[1:sss,nelx-sss:nelx].=0.0;
# x[nely-sss:nely,1:sss].=0.0;
# x[nely-sss:nely,nelx-sss:nelx].=0.0;
###
x[1:ss4,1:Int(nelx/2)-ss2].=0.0;
x[1:ss4,Int(nelx/2)+ss2:nelx].=0.0;
......
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
#### based on: A 110 LINE TOPOLOGY OPTIMIZATION CODE WITH HEAVISIDE FILTERING Nov, 2010####
# based on: Projection-based two-phase minimum and maximum length scale control in topology optimization
include("./mma/mmasub.jl");
using Images
function top_max_size_mma(alpha,alphamax,beta,betamax,d_beta,d_eta,eta,etamax,nelx,nely,nelz,max_it,rmax,rmin_s,rmin_v,volfrac,xPhys_check,x_check,x_i)
padX=true;
## MATERIAL PROPERTIES
E0 = 100;
Emin = 1e-4;
Emax =100;
nu = 0.3;
## MIN MAX PROPERTIES
rmin_m=rmax+rmin_v;
rmin_xv=rmin_m;
nele = nelx*nely;
dc_s = ones(Float64,nely,nelx);
dc_v = ones(Float64,nely,nelx);
dc_xv = ones(Float64,nely,nelx);
dv_s = ones(Float64,nely,nelx);
dv_v = ones(Float64,nely,nelx);
dv_xv = ones(Float64,nely,nelx);
stop = 0;
## PREPARE FINITE ELEMENT ANALYSIS
# FE: Build the index vectors for the for coo matrix format.
KE=lk()
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx)
edofVec = reshape(2*nodenrs[1:end-1,1:end-1].+1,nelx*nely,1)
edofMat = repeat(edofVec,1,8).+repeat([0 1 2*nely.+[2 3 0 1] -2 -1],nelx*nely,1)
iK = convert(Array{Int},reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1))
jK = convert(Array{Int},reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1))
# DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
nele = nelx*nely;
F = sparse([2],[1],[-1.0],2*(nely+1)*(nelx+1),1)
U = zeros(2*(nely+1)*(nelx+1),1)
fixeddofs = union(1:2:2*(nely+1),2*(nelx+1)*(nely+1))
alldofs = 1:(2*(nely+1)*(nelx+1))
freedofs = setdiff(alldofs,fixeddofs);
## PREPARE FILTER
if padX
H_s,Hs_s= make_filter(nelx+2*rmin_s,nely+2*rmin_s,rmin_s);
H_v,Hs_v= make_filter(nelx+2*rmin_v,nely+2*rmin_v,rmin_v);
H_xv,Hs_xv= make_filter(nelx+2*rmin_xv,nely+2*rmin_xv,rmin_xv);
else
H_s,Hs_s= make_filter(nelx,nely,rmin_s);
H_v,Hs_v= make_filter(nelx,nely,rmin_v);
H_xv,Hs_xv= make_filter(nelx,nely,rmin_xv);
end
## DESIGNATE PASSIVE ELEMENTS
passive = zeros(nely,nelx);
if !padX
for i = 1:rmin_s-1
passive[i,:] .= 1;
end
for i = nely-rmin_s+2:nely
passive[i,:] .= 1;
end
for i = nelx-rmin_s+2:nelx
passive[:,i] .= 1;
end
end
## INITIALIZE ITERATION
x=x_i .* ones(Float64,nely,nelx)
weightx_s,weightx_v,weightx_xv = w_function(alpha,x);
if padX
xTilde_s = parent(padarray(weightx_s, Pad(:reflect,rmin_s,rmin_s)));
xTilde_v = parent(padarray(weightx_v, Pad(:reflect,rmin_v,rmin_v)));
xTilde_xv = parent(padarray(weightx_xv, Pad(:reflect,rmin_xv,rmin_xv)));
else
xTilde_s = weightx_s;
xTilde_v = weightx_v;
xTilde_xv = weightx_xv;
end
xTilde_s[:] = (H_s*xTilde_s[:])./Hs_s;
xTilde_v[:] = (H_v*xTilde_v[:])./Hs_v;
xTilde_xv[:] = (H_xv*xTilde_xv[:])./Hs_xv;
if padX
xTilde_s = xTilde_s[rmin_s+1:end-rmin_s,rmin_s+1:end-rmin_s];
xTilde_v = xTilde_v[rmin_v+1:end-rmin_v,rmin_v+1:end-rmin_v];
xTilde_xv=xTilde_xv[rmin_xv+1:end-rmin_xv,rmin_xv+1:end-rmin_xv];
end
xPhys_s = 1 .-exp.(-beta .*xTilde_s)+xTilde_s.*exp.(-beta);
xPhys_v = 1 .-exp.(-beta .*xTilde_v)+xTilde_v.*exp.(-beta);
xPhys_xv =1 .-exp.(-beta .*xTilde_xv)+xTilde_xv.*exp.(-beta);
xPhys = (xPhys_s+xPhys_xv-xPhys_v)/2;
# INITIALIZE MMA OPTIMIZER
m = 1; # The number of general constraints.
n = nele; # The number of design variables x_j.
xmin = zeros(n); # Column vector with the lower bounds for the variables x_j.
xmax = ones(n); # Column vector with the upper bounds for the variables x_j.
xold1 = copy(x); # xval, one iteration ago (provided that iter>1).
xold2 = copy(x); # xval, two iterations ago (provided that iter>2).
low = copy(x); # Column vector with the lower asymptotes from the previous iteration (provided that iter>1).
upp = copy(x); # Column vector with the upper asymptotes from the previous iteration (provided that iter>1).
a0 = 1; # The constants a_0 in the term a_0*z.
a = zeros(m); # Column vector with the constants a_i in the terms a_i*z.
c_MMA = 100*ones(m); # Column vector with the constants c_i in the terms c_i*y_i. Initially 10000
d = zeros(m); # Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
## START ITERATION
loopbeta = 0;
loop = 0;
change = 1;
while change > 0.001
loopbeta = loopbeta+1;
loop = loop+1;
## FE-ANALYSIS
sK = reshape(KE[:]*(Emin.+xPhys[:]'.^eta*(Emax-Emin)),64*nelx*nely,1)
K = sparse(vec(iK),vec(jK),vec(sK)); K = (K+K')/2
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)
c = sum(sum((Emin.+xPhys.^eta*(Emax-Emin)).*ce))
dc = -eta*(Emax-Emin)*xPhys.^(eta-1).*ce
dv = ones(nely,nelx);
dxTilde_s,dxTilde_v,dxTilde_xv = dw_function(alpha,x);
dx_s =(beta.*exp.(-beta .*xTilde_s).+exp.(-beta));
dx_v =(beta.*exp.(-beta .*xTilde_v).+exp.(-beta));
dx_xv=(beta.*exp.(-beta .*xTilde_xv).+exp.(-beta));
if padX
ddc_s = parent(padarray(dc.*dx_s, Pad(:reflect,rmin_s,rmin_s)));
ddc_v = parent(padarray(dc.*dx_v, Pad(:reflect,rmin_v,rmin_v)));
ddc_xv = parent(padarray(dc.*dx_xv, Pad(:reflect,rmin_xv,rmin_xv)));
ddc_s[:] = H_s* (ddc_s[:]./Hs_s);
ddc_v[:] = H_v* (ddc_v[:]./Hs_v);
ddc_xv[:] = H_xv*(ddc_xv[:]./Hs_xv);
dc_s = ddc_s[rmin_s+1:end-rmin_s,rmin_s+1:end-rmin_s];
dc_v = ddc_v[rmin_v+1:end-rmin_v,rmin_v+1:end-rmin_v];
dc_xv= ddc_xv[rmin_xv+1:end-rmin_xv,rmin_xv+1:end-rmin_xv];
else
dc_s[:] = H_s*(dc[:].*dx_s[:]./Hs_s);
dc_v[:] = H_v*(dc[:].*dx_v[:]./Hs_v);
dc_xv[:] = H_xv*(dc[:].*dx_xv[:]./Hs_xv);
end
dc_s=dc_s .*dxTilde_s;
dc_v=dc_v .*dxTilde_v;
dc_xv=dc_xv .*dxTilde_xv;
dc=0.5*(dc_s+dc_xv-dc_v);
if padX
ddv_s = parent(padarray(dv.*dx_s, Pad(:reflect,rmin_s,rmin_s)));
ddv_v = parent(padarray(dv.*dx_v, Pad(:reflect,rmin_v,rmin_v)));
ddv_xv = parent(padarray(dv.*dx_xv, Pad(:reflect,rmin_xv,rmin_xv)));
ddv_s[:] = H_s* (ddv_s[:]./Hs_s);
ddv_v[:] = H_v* (ddv_v[:]./Hs_v);
ddv_xv[:] = H_xv*(ddv_xv[:]./Hs_xv);
dv_s = ddv_s[rmin_s+1:end-rmin_s,rmin_s+1:end-rmin_s];
dv_v = ddv_v[rmin_v+1:end-rmin_v,rmin_v+1:end-rmin_v];
dv_xv= ddv_xv[rmin_xv+1:end-rmin_xv,rmin_xv+1:end-rmin_xv];
else
dv_s[:] = H_s*(dv[:].*dx_s[:]./Hs_s);
dv_v[:] = H_v*(dv[:].*dx_v[:]./Hs_v);
dv_xv[:] = H_xv*(dv[:].*dx_xv[:]./Hs_xv);
end
dv_s=dv_s .*dxTilde_s;
dv_v=dv_v .*dxTilde_v;
dv_xv=dv_xv .*dxTilde_xv;
dv=0.5*(dv_s+dv_xv-dv_v);
# METHOD OF MOVING ASYMPTOTES
xval = x[:];
f0val = c;
df0dx = dc[:];
fval = sum(xPhys[:])/(volfrac*nele) - 1;
dfdx = dv[:]' / (volfrac*nele);
xmma, ymma, zmma, lam, xsi, eta_, mu, zet, s, low, upp = mmasub(m, n, loop, xval, xmin, xmax[:], xold1[:], xold2[:], f0val,df0dx,fval,dfdx,low[:],upp[:],a0,a,c_MMA,d,beta);
# Update MMA Variables
xnew = reshape(xmma,nely,nelx);
if !padX
xnew[findall(passive==1)] .= 0.5;
xnew[findall(passive==2)] .= 1;
end
xnew_s,xnew_v,xnew_xv = w_function(alpha,xnew); #void weighting function
if padX
xnew_s = parent(padarray(xnew_s, Pad(:reflect,rmin_s,rmin_s)));
xnew_v = parent(padarray(xnew_v, Pad(:reflect,rmin_v,rmin_v)));
xnew_xv = parent(padarray(xnew_xv, Pad(:reflect,rmin_xv,rmin_xv)));
xnew_s[:] = (H_s* xnew_s[:])./Hs_s;
xnew_v[:] = (H_v* xnew_v[:])./Hs_v;
xnew_xv[:] = (H_xv*xnew_xv[:])./Hs_xv;
xTilde_s = xnew_s[rmin_s+1:end-rmin_s,rmin_s+1:end-rmin_s];
xTilde_v = xnew_v[rmin_v+1:end-rmin_v,rmin_v+1:end-rmin_v];
xTilde_xv= xnew_xv[rmin_xv+1:end-rmin_xv,rmin_xv+1:end-rmin_xv];
else
xTilde_s[:] = (H_s*xnew_s[:])./Hs_s;
xTilde_v[:] = (H_v*xnew_v[:])./Hs_v;
xTilde_xv[:] = (H_xv*xnew_xv[:])./Hs_xv;
end
xPhys_s = 1 .-exp.(-beta .*xTilde_s)+xTilde_s.*exp.(-beta);
xPhys_v = 1 .-exp.(-beta .*xTilde_v)+xTilde_v.*exp.(-beta);
xPhys_xv = 1 .-exp.(-beta .*xTilde_xv)+xTilde_xv.*exp.(-beta);
xPhys = (xPhys_s+xPhys_xv-xPhys_v)/2;
xold2 = copy(xold1);
xold1 = copy(x);
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
## PRINT RESULTS
println(" It:$loop Obj:$(round.(c,digits=2)) Vol:$(round.(mean(xPhys[:]),digits=2)) ch:$(round.(change,digits=2)) ")
## PLOT DENSITIES
if mod(loop,5)==0
xxx=1 .- clamp.(xPhys,0,1)
display(Plots.heatmap(xxx,xaxis=nothing,showaxis = false,yaxis=nothing,legend=nothing,fc=:grays,clims=(0.0, 1.0),aspect_ratio=:equal))
Plots.heatmap(xxx,xaxis=nothing,showaxis = false,yaxis=nothing,legend=nothing,fc=:grays,clims=(0.0, 1.0),aspect_ratio=:equal)
end
## UPDATE HEAVISIDE REGULARIZATION PARAMETER
if (loopbeta >= max_it || change <= 0.001) #(beta < betamax || eta < etamax) &&
# beta = min(d_beta*beta,betamax); #beta MULTIPLIES by factor of d_beta every iteration
beta = min(d_beta+beta,betamax); #beta ADDS d_beta every iteration
eta = min(etamax,eta + d_eta);
alpha = min(alphamax,alpha + 2);
loopbeta = 0;
# max_it = 40;
change = min(1,3-floor(beta/betamax)-floor(eta/etamax)-stop);
stop = floor(beta/betamax)*floor(eta/etamax);
println("Parameter beta equal to $beta ");
println("Parameter eta equal to $eta ");
println("Parameter alpha equal to $alpha" );
end
end
return xPhys;
end
function w_function(alpha,x)
######################
# Toy Weighting
# xw_s = x;
# xw_v = 1-x;
# xw_m = 1-x;
# ######################
# # Normal Weighting
# alpha_m = alpha+1;
# xw_s = 2.^(-(alpha*x-alpha).^2);
# xw_v = 2.^(-(alpha*x ).^2); #void weighting function
# xw_m = 2.^(-(alpha_m*x ).^2); #void weighting function
######################
# BoxCar Weighting
step_s = 3/12;
step_v = 3/12;
step_m = .3;
alpha_m = alpha;
xw_s = 1 ./(1 .+exp.(-(x .-(1 .-step_s))*alpha)) .- 1 ./(1 .+exp.(-(x .-(1 +step_s))*alpha));
xw_v = 1 ./(1 .+exp.(-(x .-(0 .-step_v))*alpha)) .- 1 ./(1 .+exp.(-(x .-(0 +step_v))*alpha));
xw_m = 1 ./(1 .+exp.(-(x .-(0 .-step_m))*alpha_m)) .- 1 ./(1 .+exp.(-(x .-(0 +step_m))*alpha_m));
# #####################
# # Heaviside Weighting
#
# as = 0.002;
# av = 0.0005;
# am = 0.0001;
#
# ns = -2*log(as);
# nv = -2*log(av);
# nm = -2*log(am);
#
# xw_s = (1+as)./(1+as*exp.(2*ns*(1-x)));
# xw_v = (1+av)./(1+av*exp.(2*nv*(x )));
# xw_m = (1+am)./(1+am*exp.(2*nm*(x )));
return xw_s,xw_v,xw_m
end
function dw_function(alpha,x)
######################
# Toy Weighting
# xdw_s = 1;
# xdw_v = -1;
# #####################
# # Normal Weighting
# alpha_m = alpha+1;
# xdw_s = -(2.^(1-(-alpha+alpha *x).^2).*alpha .*(-alpha+alpha .*x).*log(2)); #derivative of solid weighting function
# xdw_v = -(2.^(1-( alpha *x).^2).*alpha .*( alpha .*x).*log(2)); #derivative of void weighting function
# xdw_m = -(2.^(1-( alpha_m*x).^2).*alpha_m.*( alpha_m.*x).*log(2)); #derivative of void weighting function
######################
# BoxCar Weighting
step_s = 3/12;
step_v = 3/12;
step_m = .3;
alpha_m = alpha;
xdw_s = (alpha*exp.(alpha*(-1*(x) .-step_s .+1 )))./(1 .+exp.(alpha*(-1*(x) .-step_s .+1))).^2 .- (alpha*exp.(alpha*(-1*(x) .+step_s .+1 )))./(1 .+exp.(alpha*(-1*(x) .+step_s .+1))).^2;
xdw_v = (alpha*exp.(alpha*(-1*(x) .-step_v .+0 )))./(1 .+exp.(alpha*(-1*(x) .-step_v .+0))).^2 .- (alpha*exp.(alpha*(-1*(x) .+step_v .+0 )))./(1 .+exp.(alpha*(-1*(x) .+step_v .+0))).^2;
xdw_m = (alpha_m*exp.(alpha_m*(-1*(x) .-step_m .+0 )))./(1 .+exp.(alpha_m*(-1*(x) .-step_m .+0))).^2 .- (alpha_m*exp.(alpha_m*(-1*(x) .+step_m .+ 0 )))./(1 .+exp.(alpha_m*(-1*(x) .+step_m .+0))).^2;
######################
# ######################
# # Heaviside Weighting
#
# as = 0.002;
# av = 0.0005;
# am = 0.0001;
#
# # ns = -2*log(as)/(phi_max-phi_min);
# # nv = -2*log(av)/(phi_max-phi_min);
# ns = -2*log(as);
# nv = -2*log(av);
# nm = -2*log(am);
#
# # tmp = exp.(2*ns*(phi_max-phi));
# tmp = exp.(2*ns*(-x));
# xdw_s = (1+as)*( 2*ns*as*tmp)./(1+as*tmp).^2;
#
# # tmp = exp.(2*nv*(phi-phi_min));
# tmp = exp.(2*nv*( x));
# xdw_v = (1+av)*(-2*nv*av*tmp)./(1+av*tmp).^2;
#
# tmp = exp.(2*nv*( x));
# xdw_m = (1+am)*(-2*nm*am*tmp)./(1+am*tmp).^2;
# ######################
return xdw_s,xdw_v,xdw_m
end
\ No newline at end of file
#-------------------------------------------------------------
#
# Copyright (C) 2007, 2008 Krister Svanberg
#
# This file, mmasub.m, is part of GCMMA-MMA-code.
#
# GCMMA-MMA-code is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3 of
# the License, or (at your option) any later version.
#
# This code is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.