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

clustering and inverse top

parent 711f0ad6
Pipeline #19643 passed with stage
in 36 seconds
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
# (c) Massachusetts Institute of Technology 2021
#======================================================================================================================#
# Asymptotic Homogenization Implementation in Julia
# Based on: https://asmedigitalcollection.asme.org/materialstechnology/article-abstract/141/1/011005/368579
......
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
#======================================================================================================================#
##############################################################
include("./mmaSolver.jl");
## SUB FUNCTION: OC
function OC2D(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, 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)
# display("min dc $(minimum(dc)),max dc $(maximum(dc))")
xTilde=copy(x);
nelx=size(x)[1]
nely=size(x)[2]
# METHOD OF MOVING ASYMPTOTES
xval = x[:];
# f0val = c;
f0val = 0;
df0dx = dc[:];
fval = sum(x[:])/(volfrac*nele) - 1;
fval=[fval,1-sum(x[:])/(volfrac*nele)];
dfdx = dv[:]' / (volfrac*nele);
dfdx = [dfdx;-dv[:]' / (volfrac*nele)];
## Solving with MMA
m=length(fval);
# m = 2; # The number of general constraints.
mdof = 1:m;
a0 = 1;
a = zeros(m);
c_MMA = ones(m)*1000;
c_MMA = ones(m)*10000;
d = zeros(m);
n = nele;
# xmin=max.(0.0, x[:].-move);
# xmax=min.(1, x[:].+move);
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.
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
## SUB FUNCTION: MMA with connectivity
function OC2DMMACON(x, xold1,xold2, dc, dv, H, Hs, volfrac, nele, move, beta,loop,low,upp,L)
# display("min dc $(minimum(dc)),max dc $(maximum(dc))")
xTilde=copy(x);
nelx=size(x)[1]
nely=size(x)[2]
ss=6
x=addFabConstraint2D(x,ss);
# METHOD OF MOVING ASYMPTOTES
xval = x[:];
# f0val = c;
f0val = 0;
df0dx = dc[:];
fval = sum(x[:])/(volfrac*nele) - 1 ;
dfdx = dv[:]' / (volfrac*nele);
fval,dfdx=connectivityConstraint(x,L,H,Hs,fval,dfdx);
## Solving with MMA
# m = 1;
m=length(fval) # The number of general constraints.
mdof = 1:m;
a0 = 1;
a = zeros(m);
c_MMA = ones(m)*1000;
c_MMA = ones(m)*10000;
d = zeros(m);
n = nele;
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.
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
## SUB FUNCTION: OC
function CompliantOC(x, dc, dv, H, Hs, volfrac, nele, move, beta)
l1 = 0; l2 = 1e9;
l1 = 0; l2 = 1e9;
xTilde=copy(x);xnew=copy(x);xPhys=copy(x)
# while (l2-l1)/(l1+l2) > 1e-4
while (l2-l1)/(l2+l1) > 1e-4 && l2>1e-40
lmid = 0.5*(l2+l1);
xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(max.(1e-10,-dc./dv./lmid))))))
# xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(-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
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
## SUB FUNCTION: OC_reg (no filter or HEAVISIDE)
function OC_reg(x, dc, dv, volfrac, nele, move)
l1 = 0; l2 = 1e9;
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.(-dc./dv/lmid)))));
xPhys[:] = xnew[:]
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
function CompliantOC_reg(x, dc, dv, volfrac, nele, move)
l1 = 0; l2 = 1e9;
xnew=copy(x);
xPhys=copy(x)
# while (l2-l1)/(l1+l2) > 1e-4
while (l2-l1)/(l2+l1) > 1e-4 && l2>1e-40
lmid = 0.5*(l2+l1);
xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(max.(1e-10,-dc./dv./lmid))))))
# xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(-dc./dv/lmid)))));
xPhys[:] = xnew[:]
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
##############################################################
## SUB FUNCTION: OC #todo remove as same in 2D
function OC(x, dc, dv, H, Hs, volfrac, nele, move, beta,voxels=true,fabric=false,ss=4)
nely=size(x)[1]
nelx=size(x)[2]
nelz=size(x)[3]
if length(voxels)==1
fabric=false
voxels=ones(size(x)[1],size(x)[2],size(x)[3])
end
l1 = 0; l2 = 1e9;
xTilde=copy(x);xnew=copy(x);xPhys=copy(x)
while (l2-l1)/(l1+l2) > 1e-3
lmid = 0.5*(l2+l1);
if !fabric
xnew[Bool.(1 .-voxels)] .=0
end
if fabric
xnew.=addFabricationConstraints(ss,0,xnew,nelx,nely,nelz)
end
# xnew = max.(0,max.(x .-move,min.(1,min.(x .+move,x.*sqrt.((-dc./dv/lmid))))));
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 !fabric
xPhys[Bool.(1 .-voxels)] .=0
end
if fabric
xPhys.=addFabricationConstraints(ss,0,xPhys,nelx,nely,nelz)
end
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
## SUB FUNCTION: OC_reg3 (no filter or HEAVISIDE) #todo remove as same in 2D
function OC_reg3(x, dc, dv, volfrac, nele, move,voxels=true)
if length(voxels)==1
voxels=ones(size(x)[1],size(x)[2],size(x)[3])
end
l1 = 0; l2 = 1e9;
xnew=copy(x);
xPhys=copy(x)
while (l2-l1)/(l1+l2) > 1e-4
lmid = 0.5*(l2+l1);
xnew[Bool.(1 .-voxels)] .=0
xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(-dc./dv/lmid)))));
xnew[Bool.(1 .-voxels)] .=0
xPhys[:] = xnew[:]
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
## SUB FUNCTION: OC
function CompliantOC3(x, dc, dv, H, Hs, volfrac, nele, move, beta,voxels=true)
if length(voxels)==1
voxels=ones(size(x)[1],size(x)[2],size(x)[3])
end
l1 = 0; l2 = 1e9;
l1 = 0; l2 = 1e9;
xTilde=copy(x);xnew=copy(x);xPhys=copy(x)
# while (l2-l1)/(l1+l2) > 1e-4
while (l2-l1)/(l2+l1) > 1e-4 && l2>1e-40
lmid = 0.5*(l2+l1);
# xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(max.(1e-10,-dc./dv./lmid))))))
# xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(-dc./dv/lmid)))));
xnew[Bool.(1 .-voxels)] .=0
xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*(max.(1e-10,-dc./dv./lmid)).^0.3))));
xTilde[:] = (H*xnew[:] )./Hs;
xPhys = 1 .-exp.(-beta .*xTilde) .+xTilde .*exp.(-beta);
xPhys[Bool.(1 .-voxels)] .=0
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
function CompliantOC_reg3(x, dc, dv, volfrac, nele, move,voxels=true)
if length(voxels)==1
voxels=ones(size(x)[1],size(x)[2],size(x)[3])
end
l1 = 0; l2 = 1e9;
xnew=copy(x);
xPhys=copy(x)
# while (l2-l1)/(l1+l2) > 1e-4
while (l2-l1)/(l2+l1) > 1e-4 && l2>1e-40
lmid = 0.5*(l2+l1);
xnew[Bool.(1 .-voxels)] .=0
# xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(max.(1e-10,-dc./dv./lmid))))))
xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*(max.(1e-10,-dc./dv./lmid)).^0.3))));
# xnew = max.(0,max.(x.-move,min.(1,min.(x.+move,x.*sqrt.(-dc./dv/lmid)))));
xPhys[:] = xnew[:]
xPhys[Bool.(1 .-voxels)] .=0
if sum(xPhys[:] ) > volfrac*nele
l1 = lmid;
else
l2 = lmid;
end
end
change = maximum(abs.(xnew[:]-x[:] )); x = xnew;
return x, xPhys, change
end
\ No newline at end of file
using StaticArrays: minimum
using Base: return_types
using Plots: display
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
#======================================================================================================================#
#########################################################################################################
using StaticArrays, BenchmarkTools, SparseArrays
using LinearAlgebra, Statistics
using StaticArrays: minimum
using Base: return_types
using Plots: display
using StaticArrays, BenchmarkTools, SparseArrays
using LinearAlgebra, Statistics
# import JSON
using NLopt
using Plots
# using Images
using VectorizedRoutines #matlab.meshgrid
#in case multistructre compliant with clustering
# using Clustering
# using StatsPlots
# using ParallelKMeans
# using Distances
using JLD #load and save data
gr(size=(400,300))
......@@ -34,33 +31,51 @@ gr(size=(400,300))
#########################################################################################################
####### utilities
include("./utils.jl")
include("./filter.jl")
include("./element.jl")
include("./OC.jl")
include("./plot.jl")
####### testcases
include("./problems.jl")
####### optimization
include("./top2D.jl")
include("./top3D.jl")
include("./topmgcg.jl")
# Asymp_homogenization
# include("./timeSpace.jl")
# multiscale and homogenization
include("./asymp_homogenization.jl")
include("./multimaterial.jl")
include("./microstructure_topX.jl")
# include("./multiscale.jl")
include("./problems.jl")
# concurrent
include("./concurrent2D.jl")
include("./concurrent3D.jl")
# strain library
include("./microstructureDesignUe.jl")
include("./strain.jl");
include("./fabrication.jl")
include("./topmgcg.jl")
# clustering
include("./clustering2D.jl");
include("./clustering3D.jl");
# fabrication constraints
include("./fabrication.jl")
# include("./maxlength2D.jl")
#########################################################################################################
println("Loaded Topology Optimization Library!")
#########################################################################################################
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
......@@ -135,4 +135,96 @@ function exportFacesToObj(Micro_xPhys,scale,Micro_rmin)
end
#######################################################
function connectivityConstraint(x,L,H,Hs,fval,dfdx)
nelx=size(x)[1]
nely=size(x)[2]
# tPhys=xPhys;
xPhys=copy(x)
xPhys[x.>=0.7].=1.0;
xPhys[x.<0.7].=0.0;
# display(Plots.heatmap(xPhys))
seg = fast_scanning(xPhys, 0.5);
tPhys=Float32.(labels_map(seg));
ii=1
for i=1:length(segment_labels(seg))
if segment_mean(seg, segment_labels(seg)[i]) ==0
tPhys[tPhys.==segment_labels(seg)[i]].=0
else
tPhys[tPhys.==segment_labels(seg)[i]].=ii
ii+=1
end
end
# display(Plots.heatmap(tPhys))
# seg = fast_scanning(x, 0.5);
# tPhys=labels_map(seg);
# tPhys=tPhys./length(segment_labels(seg))
# println(segment_labels(seg))
# print("seg: $(length(segment_labels(seg))),")
# display(Plots.heatmap(tPhys)) # returns the assigned labels map
# segment_labels(seg) # returns a list of all assigned labels
# segment_mean(seg, 1) # returns the mean intensity of label 1
# segment_pixel_count(seg, 1) # returns the pixel count of label 1
## Continuity constraint
Nei = 1:nely;
LL = L;
kk = 2*(nely*nelx); # controlling the smoothness of the time field
# kk = sum(x); # controlling the smoothness of the time field
A = LL*tPhys[:];
# A =A .-minimum(A);
B = A.^2/(nely*nelx);
fval = [fval;-(kk*(sum(B)-1.0e-6))];
# print_out = fval[end]/kk;
dft = kk*2*LL'*A;
dft = H*(dft./Hs)/(nely*nelx);
dfdx=[dfdx;-dft']
# fval = [fval; (tPhys[Nei]'[:] .- 1.0e-9)];
# ss = zeros(length(Nei), nely*nelx);
# for ii = 1 : length(Nei)
# ss[ii, Nei[ii]] = 1;
# end
# dfdx = [dfdx; (H*(ss'./repeat(Hs, 1, length(Nei))))'];
return fval,dfdx
end
function connectivityMatrix(nelx,nely)
lrmin = 2;
iH = ones(convert(Int,nelx*nely*(2*(ceil(lrmin)-1)+1)^2),1)
jH = ones(Int,size(iH))
sH = zeros(size(iH))
k = 0;
for i1 = 1:nelx
for j1 = 1:nely
e1 = (i1-1)*nely+j1
for i2 = max(i1-(ceil(lrmin)-1),1):min(i1+(ceil(lrmin)-1),nelx)
for j2 = max(j1-(ceil(lrmin)-1),1):min(j1+(ceil(lrmin)-1),nely)
e2 = (i2-1)*nely+j2
if e1 == e2
continue;
end
k = k+1
iH[k] = e1
jH[k] = e2
sH[k] = 1
end
end
end
end
L = sparse(vec(iH),vec(jH),vec(sH));
M = sparse(repeat(sum(L, dims=2), 1, size(L)[2]));
E = sparse(1.0I, size(L)[1], size(L)[2]);
L = E - L./M;
return L
end
#######################################################
\ No newline at end of file
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
#======================================================================================================================#
#### based on: A 110 LINE TOPOLOGY OPTIMIZATION CODE WITH HEAVISIDE FILTERING Nov, 2010####
......
#======================================================================================================================#
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2021
#======================================================================================================================#