Commit 69a31706 authored by Amira Abdel-Rahman's avatar Amira Abdel-Rahman
Browse files

max min thickness

parent 18a19ccd
Pipeline #15881 passed with stage
in 24 seconds
This diff is collapsed.
......@@ -84,11 +84,15 @@ function mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2,f0val,df0dx,fval,dfdx,low,up
#epsimin = sqrt(m+n)*10^(-9);
epsimin = 10^(-7);
raa0 = 0.00001;
move = 0.5;
move = 0.5;# FOR COMPLIANT
# move = 1.0;
albefa = 0.1;
asyinit = 0.5;
asyinit = 0.5;# FOR COMPLIANT
# asyinit = 0.01;
asyincr = 1.2;
asydecr = 0.7;
eeen = ones(n);
......
......@@ -213,6 +213,132 @@ function topologyOptimizationMMA(nelx,nely,volfrac,rmin,penal,maxEval)
return xPhys
end
include("./mmaSolver.jl");
function topologyOptimizationMMA1(nelx,nely,prob,volfrac,rmin,penal,maxIter)
# Max and min stiffness
Emin=1e-2;Emax=1.0;
# dofs:
ndof = 2*(nelx+1)*(nely+1)
# Allocate design variables (as array), initialize and allocate sens.
x=volfrac .* ones(Float64,nely,nelx)
xold=copy(x)
xPhys=copy(x)
beta = 1;
xPhys = 1 .- exp.(-beta .*x) .+ x * exp.(-beta);
g=0 # must be initialized to use the NGuyen/Paulino OC approach
dc=zeros(Float64,(nely,nelx))
# 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)
U,F,freedofs=prob(nelx,nely)
H,Hs=make_filter(nelx,nely,rmin)
nel=nely*nelx
xold1 = copy(x); # xval, one iteration ago (provided that iter>1).
xold2 = copy(x); # xval, two iterations ago (provided that iter>2).
low = ones(nel); # Column vector with the lower asymptotes from the previous iteration (provided that iter>1).
upp = ones(nel); # Column vector with the upper asymptotes from the previous iteration (provided that iter>1).
change=1;loop=0;
for i =1:maxIter
if (change > 0.01)
# Start iteration
loop += 1
# FE-ANALYSIS
sK = reshape(KE[:]*(Emin.+xPhys[:]'.^penal*(Emax-Emin)),64*nelx*nely,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)
c = sum(sum((Emin.+xPhys.^penal*(Emax-Emin)).*ce))
dc = -penal*(Emax-Emin)*xPhys.^(penal-1).*ce
dv = ones(nely,nelx)
dc[:] = H*(dc[:]./Hs)
dv[:] = H*(dv[:]./Hs)
move = 0.2;
x, xPhys, change, low, upp, xold1, xold2 = OCMMA(x, xold1,xold2,c,dc, dv, H, Hs, volfrac, nel, move, beta,loop,low,upp);
m=mean(xPhys[:])
println(" It:$loop Obj:$(round.(c,digits=2)) Vol:$(round.(m,digits=2)) ch:$(round.(change,digits=2)) ")
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
xPhys = copy(x)
end
end
display(Plots.heatmap(xPhys, aspect_ratio=:equal, legend=false, axis=nothing,fc=:grays))
return xPhys
end
## SUB FUNCTION: MMA
function OCMMA(x, xold1,xold2,c, dc, dv, H, Hs, volfrac, nele, move, beta,loop,low,upp)
# display("min dc $(minimum(dc)),max dc $(maximum(dc))")
xTilde=copy(x);nelx=size(x)[2];nely=size(x)[1];
m = 1; # The number of general constraints.
n = nele; # The number of design variables x_j.
xmin = ones(n).*Emin; # Column vector with the lower bounds for the variables x_j.
xmax = ones(n).*Emax; # Column vector with the upper bounds for the variables x_j.
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 = 10000*ones(m); # Column vector with the constants c_i in the terms c_i*y_i.
d = zeros(m); # Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
# METHOD OF MOVING ASYMPTOTES
xval = x[:];
f0val = c;
df0dx = dc[:];
fval = sum(x[:])/(volfrac*nele) .- 1.0;
dfdx = dv[:]' / (volfrac*nele);
# xmin=max.(0.0, x[:].-move);
# xmax=min.(1, x[:].+move);
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);
# Update MMA Variables
xnew = reshape(xmma,nely,nelx);
xTilde[:] = (H*xnew[:] )./Hs;
xPhys = 1 .-exp.(-beta .*xTilde) .+xTilde .*exp.(-beta);
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
xold2 = copy(xold1);
xold1 = copy(x);
return x, xPhys, change,low,upp,xold1,xold2
end
#########################################################################################################
function CompliantTopologyOptimization(nelx,nely,volfrac,rmin,penal,maxIter,Load,Support,Spring,DOut)
......
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