You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
816 lines
28 KiB
816 lines
28 KiB
from __future__ import division
|
|
import numpy as np
|
|
from scipy.sparse import coo_matrix
|
|
from scipy.sparse.linalg import spsolve
|
|
from matplotlib import colors
|
|
import matplotlib.pyplot as plt
|
|
import time
|
|
from localSMRT import LocalSMRT
|
|
from p2voxel import drawVoxel_3
|
|
from p2voxel import setPara
|
|
from getDH import getDH
|
|
from p2voxel import p2voxel
|
|
|
|
|
|
def lk(nu):
|
|
A = np.array([[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]])
|
|
b = np.array([[1],[nu]])
|
|
k = 1/float(144)*np.dot(A.T,b).flatten()
|
|
|
|
K1 = np.array([[k[0],k[1],k[1],k[2],k[4],k[4]],
|
|
[k[1],k[0],k[1],k[3],k[5],k[6]],
|
|
[k[1],k[1],k[0],k[3],k[6],k[5]],
|
|
[k[2],k[3],k[3],k[0],k[7],k[7]],
|
|
[k[4],k[5],k[6],k[7],k[0],k[1]],
|
|
[k[4],k[6],k[5],k[7],k[1],k[0]]])
|
|
|
|
K2 = np.array([[k[8],k[7],k[11],k[5],k[3],k[6]],
|
|
[k[7],k[8],k[11],k[4],k[2],k[4]],
|
|
[k[9],k[9],k[12],k[6],k[3],k[5]],
|
|
[k[5],k[4],k[10],k[8],k[1],k[9]],
|
|
[k[3],k[2],k[4],k[1],k[8],k[11]],
|
|
[k[10],k[3],k[5],k[11],k[9],k[12]]])
|
|
|
|
K3 = np.array([[k[5],k[6],k[3],k[8],k[11],k[7]],
|
|
[k[6],k[5],k[3],k[9],k[12],k[9]],
|
|
[k[4],k[4],k[2],k[7],k[11],k[8]],
|
|
[k[8],k[9],k[1],k[5],k[10],k[4]],
|
|
[k[11],k[12],k[9],k[10],k[5],k[3]],
|
|
[k[1],k[11],k[8],k[3],k[4],k[2]]])
|
|
|
|
K4 = np.array([[k[13],k[10],k[10],k[12],k[9],k[9]],
|
|
[k[10],k[13],k[10],k[11],k[8],k[7]],
|
|
[k[10],k[10],k[13],k[11],k[7],k[8]],
|
|
[k[12],k[11],k[11],k[13],k[6],k[6]],
|
|
[k[9],k[8],k[7],k[6],k[13],k[10]],
|
|
[k[9],k[7],k[8],k[6],k[10],k[13]]])
|
|
|
|
K5 = np.array([[k[0],k[1],k[7],k[2],k[4],k[3]],
|
|
[k[1],k[0],k[7],k[3],k[5],k[10]],
|
|
[k[7],k[7],k[0],k[4],k[10],k[5]],
|
|
[k[2],k[3],k[4],k[0],k[7],k[1]],
|
|
[k[4],k[5],k[10],k[7],k[0],k[7]],
|
|
[k[3],k[10],k[5],k[1],k[7],k[0]]])
|
|
|
|
K6 = np.array([[k[13],k[10],k[6],k[12],k[9],k[11]],
|
|
[k[10],k[13],k[6],k[11],k[8],k[1]],
|
|
[k[6],k[6],k[13],k[9],k[1],k[8]],
|
|
[k[12],k[11],k[9],k[13],k[6],k[10]],
|
|
[k[9],k[8],k[1],k[6],k[13],k[6]],
|
|
[k[11],k[1],k[8],k[10],k[6],k[13]]])
|
|
|
|
KE1=np.hstack((K1,K2,K3,K4))
|
|
KE2=np.hstack((K2.T,K5,K6,K3.T))
|
|
KE3=np.hstack((K3.T,K6,K5.T,K2.T))
|
|
KE4=np.hstack((K4,K3,K2,K1.T))
|
|
KE = 1/float(((nu+1)*(1-2*nu)))*np.vstack((KE1,KE2,KE3,KE4))
|
|
|
|
return(KE)
|
|
|
|
def save_x(data, pre):
|
|
dirr = 'data'
|
|
# pre='test'
|
|
np.savetxt(f'./{dirr}/{pre}_x.txt', data.reshape(1, -1), fmt="%.4f", delimiter=",")
|
|
|
|
|
|
# MAIN DRIVER
|
|
def main(nelx, nely, nelz, x, penal, rmin, ft):
|
|
print("Minimum compliance problem with OC")
|
|
print("ndes: " + str(nelx) + " x " + str(nely) + " x " + str(nelz))
|
|
print( "rmin: " + str(rmin) + ", penal: " + str(penal))
|
|
print("Filter method: Sensitivity based")
|
|
|
|
|
|
deltax=0.005
|
|
|
|
# Max and min stiffness
|
|
Emin = 1e-9
|
|
Emax = 1.0
|
|
nu = 0.3
|
|
|
|
# dofs
|
|
ndof = 3 * (nelx + 1) * (nely + 1) * (nelz + 1)
|
|
|
|
# USER - DEFINED LOAD DOFs
|
|
kl = np.arange(nelz + 1)
|
|
loadnid = kl * (nelx + 1) * (nely + 1) + (nely + 1) * (nelx + 1) - 1 # Node IDs
|
|
loaddof = 3 * loadnid + 1 # DOFs
|
|
|
|
# USER - DEFINED SUPPORT FIXED DOFs
|
|
[jf, kf] = np.meshgrid(np.arange(nely + 1), np.arange(nelz + 1)) # Coordinates
|
|
fixednid = (kf) * (nely + 1) * (nelx + 1) + jf # Node IDs
|
|
fixeddof = np.array([3 * fixednid, 3 * fixednid + 1, 3 * fixednid + 2]).flatten() # DOFs
|
|
|
|
# Allocate design variables (as array), initialize and allocate sens.
|
|
|
|
xold = x.copy()
|
|
xPhys = x.copy()
|
|
|
|
# FE: Build the index vectors for the for coo matrix format.
|
|
KE = np.zeros((24,24),dtype=float)
|
|
DH=getDH(x)
|
|
|
|
edofMat = np.zeros((nelx * nely * nelz, 24), dtype=int)
|
|
for elz in range(nelz):
|
|
for elx in range(nelx):
|
|
for ely in range(nely):
|
|
el = ely + (elx * nely) + elz * (nelx * nely)
|
|
n1 = elz * (nelx + 1) * (nely + 1) + (nely + 1) * elx + ely
|
|
n2 = elz * (nelx + 1) * (nely + 1) + (nely + 1) * (elx + 1) + ely
|
|
n3 = (elz + 1) * (nelx + 1) * (nely + 1) + (nely + 1) * elx + ely
|
|
n4 = (elz + 1) * (nelx + 1) * (nely + 1) + (nely + 1) * (elx + 1) + ely
|
|
edofMat[el, :] = np.array(
|
|
[3 * n1 + 3, 3 * n1 + 4, 3 * n1 + 5, 3 * n2 + 3, 3 * n2 + 4, 3 * n2 + 5, \
|
|
3 * n2, 3 * n2 + 1, 3 * n2 + 2, 3 * n1, 3 * n1 + 1, 3 * n1 + 2, \
|
|
3 * n3 + 3, 3 * n3 + 4, 3 * n3 + 5, 3 * n4 + 3, 3 * n4 + 4, 3 * n4 + 5, \
|
|
3 * n4, 3 * n4 + 1, 3 * n4 + 2, 3 * n3, 3 * n3 + 1, 3 * n3 + 2])
|
|
|
|
|
|
|
|
Z=np.asarray([0,0,0,0,1,1,1,1])
|
|
Y=np.asarray([0,0,1,1,0,0,1,1])
|
|
X=np.asarray([0,1,0,1,0,1,0,1])
|
|
|
|
|
|
# Construct the index pointers for the coo format
|
|
iK = np.kron(edofMat, np.ones((24, 1))).flatten()
|
|
jK = np.kron(edofMat, np.ones((1, 24))).flatten()
|
|
|
|
# Filter: Build (and assemble) the index+data vectors for the coo matrix format
|
|
nfilter = nelx * nely * nelz * ((2 * (np.ceil(rmin) - 1) + 1) ** 3)
|
|
iH = np.zeros(int(nfilter))
|
|
jH = np.zeros(int(nfilter))
|
|
sH = np.zeros(int(nfilter))
|
|
cc = 0
|
|
for z in range(nelz):
|
|
for i in range(nelx):
|
|
for j in range(nely):
|
|
row = i * nely + j + z * (nelx * nely)
|
|
kk1 = int(np.maximum(i - (np.ceil(rmin) - 1), 0))
|
|
kk2 = int(np.minimum(i + np.ceil(rmin), nelx))
|
|
ll1 = int(np.maximum(j - (np.ceil(rmin) - 1), 0))
|
|
ll2 = int(np.minimum(j + np.ceil(rmin), nely))
|
|
mm1 = int(np.maximum(z - (np.ceil(rmin) - 1), 0))
|
|
mm2 = int(np.minimum(z + np.ceil(rmin), nelz))
|
|
for m in range(mm1, mm2):
|
|
for k in range(kk1, kk2):
|
|
for l in range(ll1, ll2):
|
|
col = k * nely + l + m * (nelx * nely)
|
|
fac = rmin - np.sqrt((i - k) * (i - k) + (j - l) * (j - l) + (z - m) * (z - m))
|
|
iH[cc] = row
|
|
jH[cc] = col
|
|
sH[cc] = np.maximum(0.0, fac)
|
|
cc = cc + 1
|
|
# Finalize assembly and convert to csc format
|
|
H = coo_matrix((sH, (iH, jH)), shape=(nelx * nely * nelz, nelx * nely * nelz)).tocsc()
|
|
Hs = H.sum(1)
|
|
|
|
# BC's and support
|
|
dofs = np.arange(3 * (nelx + 1) * (nely + 1) * (nelz + 1))
|
|
free = np.setdiff1d(dofs, fixeddof)
|
|
# Solution and RHS vectors
|
|
f = np.zeros((ndof, 1), dtype=float)
|
|
u = np.zeros((ndof, 1), dtype=float)
|
|
uold=np.zeros((ndof, 1), dtype=float)
|
|
# Set load
|
|
f[loaddof, 0] = -1
|
|
# Initialize plot and plot the initial design
|
|
|
|
# Set loop counter and gradient vectors
|
|
loop = 0
|
|
change = 1
|
|
|
|
|
|
ce = np.ones(nely * nelx * nelz, dtype=float)
|
|
cold=0
|
|
c=0
|
|
while change > 0.001 and loop < 2000:
|
|
dc = np.zeros(len(x), dtype=float)
|
|
t0 = time.perf_counter()
|
|
loop = loop + 1
|
|
print (f'loop:{loop}')
|
|
# Setup and solve FE problem
|
|
DH = getDH(x)
|
|
KE=LocalSMRT(X,Y,Z,DH)
|
|
KE=KE.reshape(24,24)
|
|
|
|
sK = ((KE.flatten()[np.newaxis]).T * (Emin + np.ones(nelx*nely*nelz) ** penal * (Emax - Emin))).flatten(order='F')
|
|
#print (len(sK),len(iK),len(jK))
|
|
K = coo_matrix((sK, (iK, jK)), shape=(ndof, ndof)).tocsc()
|
|
#print ("K = coo_matrix((sK, (iK, jK)), shape=(ndof, ndof)).tocsc()")
|
|
# Remove constrained dofs from matrix
|
|
K = K[free, :][:, free]
|
|
#print ("K = K[free, :][:, free]")
|
|
# Solve system
|
|
u[free, 0] = spsolve(K, f[free, 0])
|
|
#print ("u[free, 0] = spsolve(K, f[free, 0])")
|
|
# Objective and sensitivity
|
|
|
|
#part1=np.dot(u[edofMat].reshape(nelx * nely * nelz, 24), KE)
|
|
#part2=u[edofMat].reshape(nelx * nely * nelz,24)
|
|
#part3=(part1*part2)
|
|
#part3=part3.sum(1)
|
|
#ce[:] = part3
|
|
|
|
cold = c
|
|
c = (np.dot(uold[edofMat].reshape(nelx * nely * nelz, 24), KE) * uold[edofMat].reshape(nelx * nely * nelz,
|
|
24)).sum()
|
|
|
|
uold[:]=u
|
|
|
|
#print ("ce[:] = (np.dot(u[edofMat].reshape(nelx * nely * nelz, 24), KE) * u[edofMat].reshape(nelx * nely * nelz,24)).sum(1)")
|
|
#obj = ((Emin + xPhys ** penal * (Emax - Emin)) * ce).sum()
|
|
|
|
|
|
#print ("ue=np.asarray(u[edofMat[0]])")
|
|
|
|
#print (f'ue:{ue}')
|
|
voxel=get_voxel(x)
|
|
#print (f'DH:{DH}')
|
|
print (f'c:%e deltac:{c - cold}'%c)
|
|
|
|
for i in range(0,len(x)):
|
|
x[i]+=deltax
|
|
DH1=getDH(x)
|
|
x[i] -= 2 * deltax
|
|
DH2=getDH(x)
|
|
x[i]+=deltax
|
|
|
|
#voxel1 = get_voxel(x)
|
|
K1=LocalSMRT(X,Y,Z,DH1)
|
|
K2=LocalSMRT(X,Y,Z,DH2)
|
|
#print (f'x[{i}]: sum:{(voxel-voxel1).sum()}')
|
|
#print (x)
|
|
#print(DH1)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#print (ue.shape)
|
|
|
|
|
|
#for id in range(0,nelx*nely*nelz):
|
|
#ue=np.asarray(u[edofMat[id]])
|
|
#temp=np.matmul(np.matmul(ue.T,K2),ue)
|
|
#dc[i]+=temp.reshape(1)
|
|
#temp = np.matmul(np.matmul(ue.T, K1), ue)
|
|
#T1+=temp.reshape(1)
|
|
|
|
T2=(np.dot(u[edofMat].reshape(nelx * nely * nelz, 24), K2) * u[edofMat].reshape(nelx * nely * nelz,
|
|
24)).sum()
|
|
T1 = (np.dot(u[edofMat].reshape(nelx * nely * nelz, 24), K1) * u[edofMat].reshape(nelx * nely * nelz,
|
|
24)).sum()
|
|
dc[i]=T2-T1
|
|
|
|
|
|
if T1 > c and T2 > c:
|
|
dc[i]=0
|
|
else:
|
|
if dc[i]>0:
|
|
1
|
|
#print (f'++newc=%e'%T1)
|
|
if dc[i]<0:
|
|
1
|
|
#print (f'--newc=%e'%T2)
|
|
break
|
|
#print (KE-K1)
|
|
#print (temp)
|
|
#print(temp.shape)
|
|
|
|
|
|
#print (f'dc:{dc}')
|
|
|
|
|
|
|
|
# Optimality criteria
|
|
xold[:] = x
|
|
x[:] = MMA(x, dc)
|
|
# Filter design variables
|
|
if ft == 0:
|
|
xPhys[:] = x
|
|
elif ft == 1:
|
|
xPhys[:] = np.asarray(H * x[np.newaxis].T / Hs)[:, 0]
|
|
|
|
# Compute the change by the inf. norm
|
|
change = np.linalg.norm(x - xold, np.inf)
|
|
t1 = time.perf_counter() - t0
|
|
# Plot to screen
|
|
|
|
|
|
print(f'change:{change}')
|
|
#print(f'cost: {t1}')
|
|
save_x(xPhys, f"{loop}")
|
|
#np.save(f"data/ori_{volfrac}_{loop}epoch_{nelx}_{nely}_{nelz}",xPhys)
|
|
print (f'x:{xPhys}')
|
|
|
|
|
|
|
|
'''plt.show()'''
|
|
|
|
draw_x(xPhys)
|
|
|
|
def draw_x(x):
|
|
|
|
mtype = 387 - 1
|
|
|
|
in_parameter_values = x[:7]
|
|
|
|
out_parameter_values = x[7:]
|
|
|
|
parameters = setPara(mtype, in_parameter_values, out_parameter_values)
|
|
|
|
voxel = p2voxel(mtype, parameters, resolution=39)
|
|
#print (voxel.shape)
|
|
drawVoxel_3(voxel[0:20,0:20,0:20])
|
|
drawVoxel_3(voxel)
|
|
|
|
def get_voxel(x):
|
|
mtype = 387 - 1
|
|
|
|
in_parameter_values = x[:7]
|
|
|
|
out_parameter_values = x[7:]
|
|
|
|
parameters = setPara(mtype, in_parameter_values, out_parameter_values)
|
|
|
|
voxel = p2voxel(mtype, parameters, resolution=39)
|
|
|
|
return voxel
|
|
|
|
|
|
def MMA( x, dc):
|
|
|
|
alpha=1
|
|
#print (dv)
|
|
move = 0.005
|
|
# reshape to perform vector operations
|
|
xnew = np.zeros(len(x), dtype=float)
|
|
|
|
xnew[:] = np.maximum(0.0,np.maximum(x - move, np.minimum(0.5, np.minimum(x + move, np.maximum(x-move,x+alpha*dc)))))
|
|
|
|
return xnew
|
|
|
|
|
|
import numpy as np
|
|
import scipy.sparse
|
|
import scipy.sparse.linalg
|
|
|
|
|
|
def mmasub(m, n, iter, xval, xmin, xmax, xold1, xold2,
|
|
f0val, df0dx, fval, dfdx, low, upp, a0, a, c, d):
|
|
# Version September 2007 (and a small change August 2008)
|
|
#
|
|
# Krister Svanberg <krille@math.kth.se>
|
|
# Department of Mathematics, SE-10044 Stockholm, Sweden.
|
|
#
|
|
# This function mmasub performs one MMA-iteration, aimed at
|
|
# solving the nonlinear programming problem:
|
|
#
|
|
# Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 )
|
|
# subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m
|
|
# xmin_j <= x_j <= xmax_j, j = 1,...,n
|
|
# z >= 0, y_i >= 0, i = 1,...,m
|
|
# *** INPUT:
|
|
#
|
|
# m = The number of general constraints.
|
|
# n = The number of variables x_j.
|
|
# iter = Current iteration number ( =1 the first time mmasub is called).
|
|
# xval = Column vector with the current values of the variables x_j.
|
|
# xmin = Column vector with the lower bounds for the variables x_j.
|
|
# xmax = Column vector with the upper bounds for the variables x_j.
|
|
# xold1 = xval, one iteration ago (provided that iter>1).
|
|
# xold2 = xval, two iterations ago (provided that iter>2).
|
|
# f0val = The value of the objective function f_0 at xval.
|
|
# df0dx = Column vector with the derivatives of the objective function
|
|
# f_0 with respect to the variables x_j, calculated at xval.
|
|
# fval = Column vector with the values of the constraint functions f_i,
|
|
# calculated at xval.
|
|
# dfdx = (m x n)-matrix with the derivatives of the constraint functions
|
|
# f_i with respect to the variables x_j, calculated at xval.
|
|
# dfdx(i,j) = the derivative of f_i with respect to x_j.
|
|
# low = Column vector with the lower asymptotes from the previous
|
|
# iteration (provided that iter>1).
|
|
# upp = Column vector with the upper asymptotes from the previous
|
|
# iteration (provided that iter>1).
|
|
# a0 = The constants a_0 in the term a_0*z.
|
|
# a = Column vector with the constants a_i in the terms a_i*z.
|
|
# c = Column vector with the constants c_i in the terms c_i*y_i.
|
|
# d = Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
#
|
|
# *** OUTPUT:
|
|
#
|
|
# xmma = Column vector with the optimal values of the variables x_j
|
|
# in the current MMA subproblem.
|
|
# ymma = Column vector with the optimal values of the variables y_i
|
|
# in the current MMA subproblem.
|
|
# zmma = Scalar with the optimal value of the variable z
|
|
# in the current MMA subproblem.
|
|
# lam = Lagrange multipliers for the m general MMA constraints.
|
|
# xsi = Lagrange multipliers for the n constraints alfa_j - x_j <= 0.
|
|
# eta = Lagrange multipliers for the n constraints x_j - beta_j <= 0.
|
|
# mu = Lagrange multipliers for the m constraints -y_i <= 0.
|
|
# zet = Lagrange multiplier for the single constraint -z <= 0.
|
|
# s = Slack variables for the m general MMA constraints.
|
|
# low = Column vector with the lower asymptotes, calculated and used
|
|
# in the current MMA subproblem.
|
|
# upp = Column vector with the upper asymptotes, calculated and used
|
|
# in the current MMA subproblem.
|
|
|
|
epsimin = 1e-10
|
|
raa0 = 0.01
|
|
albefa = 0.4
|
|
asyinit = 0.1
|
|
asyincr = 0.8
|
|
asydecr = 0.6
|
|
|
|
# epsimin = sqrt(m+n)*10^(-9);
|
|
# epsimin = 1e-7
|
|
# raa0 = 1e-5
|
|
move = 1.0
|
|
# move = 0.05 # MODIFICATION BY DK
|
|
# albefa = 0.1
|
|
# asyinit = 0.5
|
|
# asyincr = 1.2
|
|
# asydecr = 0.7
|
|
# asyincr = 1.05 # MODIFICATION BY DK
|
|
# asydecr = 0.65 # MODIFICATION BY DK
|
|
eeen = np.ones((n, 1))
|
|
eeem = np.ones((m, 1))
|
|
zeron = np.zeros((n, 1))
|
|
|
|
# Calculation of the asymptotes low and upp :
|
|
if iter < 2.5:
|
|
low = xval - asyinit * (xmax - xmin)
|
|
upp = xval + asyinit * (xmax - xmin)
|
|
else:
|
|
zzz = (xval - xold1) * (xold1 - xold2)
|
|
factor = eeen
|
|
factor[zzz > 0] = asyincr
|
|
factor[zzz < 0] = asydecr
|
|
low = xval - factor * (xold1 - low)
|
|
upp = xval + factor * (upp - xold1)
|
|
lowmin = xval - 10 * (xmax - xmin)
|
|
lowmax = xval - 0.01 * (xmax - xmin)
|
|
uppmin = xval + 0.01 * (xmax - xmin)
|
|
uppmax = xval + 10 * (xmax - xmin)
|
|
low = np.maximum(low, lowmin)
|
|
low = np.minimum(low, lowmax)
|
|
upp = np.minimum(upp, uppmax)
|
|
upp = np.maximum(upp, uppmin)
|
|
|
|
# Calculation of the bounds alfa and beta :
|
|
zzz1 = low + albefa * (xval - low)
|
|
zzz2 = xval - move * (xmax - xmin)
|
|
zzz = np.maximum(zzz1, zzz2)
|
|
alfa = np.maximum(zzz, xmin)
|
|
zzz1 = upp - albefa * (upp - xval)
|
|
zzz2 = xval + move * (xmax - xmin)
|
|
zzz = np.minimum(zzz1, zzz2)
|
|
beta = np.minimum(zzz, xmax)
|
|
|
|
# Calculations of p0, q0, P, Q and b.
|
|
xmami = xmax - xmin
|
|
xmamieps = 0.00001 * eeen
|
|
xmami = np.maximum(xmami, xmamieps)
|
|
xmamiinv = eeen / xmami
|
|
# print(xmamiinv.shape) #112 x 1
|
|
ux1 = upp - xval
|
|
ux2 = ux1 * ux1
|
|
xl1 = xval - low
|
|
xl2 = xl1 * xl1
|
|
uxinv = eeen / ux1
|
|
xlinv = eeen / xl1
|
|
|
|
# p0 = zeron
|
|
# q0 = zeron
|
|
|
|
p0 = np.maximum(df0dx, 0)
|
|
q0 = np.maximum(-df0dx, 0)
|
|
|
|
# p0(find(df0dx > 0)) = df0dx(find(df0dx > 0))
|
|
# q0(find(df0dx < 0)) = -df0dx(find(df0dx < 0))
|
|
|
|
pq0 = 0.001 * (p0 + q0) + raa0 * xmamiinv
|
|
p0 = p0 + pq0
|
|
q0 = q0 + pq0
|
|
p0 = p0 * ux2
|
|
q0 = q0 * xl2
|
|
|
|
# P = coo_matrix((m,n))
|
|
# Q = coo_matrix((m,n))
|
|
|
|
P = np.maximum(dfdx, 0)
|
|
Q = np.maximum(-dfdx, 0)
|
|
|
|
# P(find(dfdx > 0)) = dfdx(find(dfdx > 0))
|
|
# Q(find(dfdx < 0)) = -dfdx(find(dfdx < 0))
|
|
# print("--",(eeem @ xmamiinv.T).shape)
|
|
PQ = 0.001 * (P + Q) + raa0 * eeem @ xmamiinv.T
|
|
# print(PQ.shape)
|
|
P = P + PQ
|
|
Q = Q + PQ
|
|
P = P @ scipy.sparse.spdiags(ux2.T, 0, n, n)
|
|
Q = Q @ scipy.sparse.spdiags(xl2.T, 0, n, n)
|
|
b = P @ uxinv + Q @ xlinv - fval
|
|
# print("--",P.shape)
|
|
|
|
# %%% Solving the subproblem by a primal-dual Newton method
|
|
xmma, ymma, zmma, lam, xsi, eta, mu, zet, s = subsolve(m, n, epsimin, low, upp, alfa,
|
|
beta, p0, q0, P, Q, a0, a, b, c, d)
|
|
|
|
return xmma, ymma, zmma, lam, xsi, eta, mu, zet, s, low, upp
|
|
|
|
|
|
def subsolve(m, n, epsimin, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d):
|
|
# -------------------------------------------------------------
|
|
# This is the file subsolv.m
|
|
#
|
|
# Version Dec 2006.
|
|
# Krister Svanberg <krille@math.kth.se>
|
|
# Department of Mathematics, KTH,
|
|
# SE-10044 Stockholm, Sweden.
|
|
#
|
|
#
|
|
# This function subsolv solves the MMA subproblem:
|
|
#
|
|
# minimize SUM[ p0j/(uppj-xj) + q0j/(xj-lowj) ] + a0*z +
|
|
# + SUM[ ci*yi + 0.5*di*(yi)^2 ],
|
|
#
|
|
# subject to SUM[ pij/(uppj-xj) + qij/(xj-lowj) ] - ai*z - yi <= bi,
|
|
# alfaj <= xj <= betaj, yi >= 0, z >= 0.
|
|
#
|
|
# Input: m, n, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d.
|
|
# Output: xmma,ymma,zmma, slack variables and Lagrange multiplers.
|
|
#
|
|
een = np.ones((n, 1))
|
|
eem = np.ones((m, 1))
|
|
epsi = 1
|
|
# epsvecn = epsi*een
|
|
# epsvecm = epsi*eem
|
|
x = 0.5 * (alfa + beta)
|
|
y = eem
|
|
z = 1
|
|
lam = eem
|
|
xsi = een / (x - alfa)
|
|
xsi = np.maximum(xsi, een)
|
|
eta = een / (beta - x)
|
|
eta = np.maximum(eta, een)
|
|
mu = np.maximum(eem, 0.5 * c)
|
|
zet = 1
|
|
s = eem
|
|
itera = 0
|
|
|
|
while epsi > epsimin:
|
|
|
|
epsvecn = epsi * een
|
|
epsvecm = epsi * eem
|
|
|
|
ux1 = upp - x
|
|
xl1 = x - low
|
|
ux2 = ux1 * ux1
|
|
xl2 = xl1 * xl1
|
|
|
|
uxinv1 = een / ux1
|
|
xlinv1 = een / xl1
|
|
|
|
# print(p0.shape)
|
|
# print(P.T.shape)
|
|
# print(lam.shape)
|
|
plam = p0 + P.T @ lam
|
|
qlam = q0 + Q.T @ lam
|
|
|
|
gvec = P @ uxinv1 + Q @ xlinv1
|
|
dpsidx = plam / ux2 - qlam / xl2
|
|
|
|
rex = dpsidx - xsi + eta
|
|
rey = c + d * y - mu - lam
|
|
rez = a0 - zet - a.T @ lam
|
|
|
|
relam = gvec - a * z - y + s - b
|
|
rexsi = xsi * (x - alfa) - epsvecn
|
|
reeta = eta * (beta - x) - epsvecn
|
|
remu = mu * y - epsvecm
|
|
rezet = zet * z - epsi
|
|
res = lam * s - epsvecm
|
|
|
|
residu1 = np.concatenate((rex, rey, rez), axis=0)
|
|
residu2 = np.concatenate((relam, rexsi, reeta, remu, np.array([[rezet]]), res), axis=0)
|
|
residu = np.concatenate((residu1, residu2), axis=0)
|
|
residunorm = np.linalg.norm(residu)
|
|
residumax = np.amax(np.absolute(residu))
|
|
|
|
ittt = 0
|
|
|
|
while (residumax > 0.9 * epsi) and (ittt < 200):
|
|
|
|
ittt = ittt + 1
|
|
itera = itera + 1
|
|
|
|
ux1 = upp - x
|
|
xl1 = x - low
|
|
ux2 = ux1 * ux1
|
|
xl2 = xl1 * xl1
|
|
ux3 = ux1 * ux2
|
|
xl3 = xl1 * xl2
|
|
|
|
uxinv1 = een / ux1
|
|
xlinv1 = een / xl1
|
|
uxinv2 = een / ux2
|
|
xlinv2 = een / xl2
|
|
|
|
plam = p0 + P.T @ lam
|
|
qlam = q0 + Q.T @ lam
|
|
|
|
gvec = P @ uxinv1 + Q @ xlinv1
|
|
GG = P @ scipy.sparse.spdiags(uxinv2.T, 0, n, n) - Q @ scipy.sparse.spdiags(xlinv2.T, 0, n, n)
|
|
dpsidx = plam / ux2 - qlam / xl2
|
|
# print('line 261, ittt=',ittt)
|
|
# print(np.amin(np.absolute(x-alfa)))
|
|
# print(np.amin(np.absolute(beta-x)))
|
|
# print()
|
|
|
|
delx = dpsidx - epsvecn / (x - alfa) + epsvecn / (beta - x)
|
|
dely = c + d * y - lam - epsvecm / y
|
|
delz = a0 - a.T @ lam - epsi / z
|
|
dellam = gvec - a * z - y - b + epsvecm / lam
|
|
diagx = plam / ux3 + qlam / xl3
|
|
diagx = 2 * diagx + xsi / (x - alfa) + eta / (beta - x)
|
|
diagxinv = een / diagx
|
|
diagy = d + mu / y
|
|
diagyinv = eem / diagy
|
|
diaglam = s / lam
|
|
diaglamyi = diaglam + diagyinv
|
|
|
|
if m < n:
|
|
blam = dellam + dely / diagy - GG @ (delx / diagx)
|
|
bb = np.concatenate((blam, delz), axis=0)
|
|
Alam = scipy.sparse.spdiags(diaglamyi.T, 0, m, m) + GG @ scipy.sparse.spdiags(diagxinv.T, 0, n,
|
|
n) @ GG.T
|
|
AA = np.concatenate((np.concatenate((Alam, a), axis=1),
|
|
np.concatenate((a.T, np.array([[-zet / z]])), axis=1)), axis=0)
|
|
solut = scipy.sparse.linalg.spsolve(AA, bb)[:, np.newaxis]
|
|
dlam = solut[0:m]
|
|
dz = solut[m, 0]
|
|
dx = -delx / diagx - (GG.T @ dlam) / diagx
|
|
|
|
|
|
else:
|
|
diaglamyiinv = eem / diaglamyi
|
|
dellamyi = dellam + dely / diagy
|
|
Axx = scipy.sparse.spdiags(diagx.T, 0, n, n) + GG.T @ scipy.sparse.spdiags(diaglamyiinv.T, 0, m, m) @ GG
|
|
azz = zet / z + a.T @ (a / diaglamyi)
|
|
axz = -GG.T @ (a / diaglamyi)
|
|
bx = delx + GG.T @ (dellamyi / diaglamyi)
|
|
bz = delz - a.T @ (dellamyi / diaglamyi)
|
|
AA = np.concatenate((np.concatenate((Axx, axz), axis=1),
|
|
np.concatenate((axz.T, azz), axis=1)), axis=0)
|
|
bb = np.concatenate((-bx, -bz), axis=0)
|
|
solut = scipy.sparse.linalg.spsolve(AA, bb)
|
|
dx = solut[0:n]
|
|
dz = solut[n, 0]
|
|
dlam = (GG @ dx) / diaglamyi - dz * (a / diaglamyi) + dellamyi / diaglamyi
|
|
|
|
dy = -dely / diagy + dlam / diagy
|
|
dxsi = -xsi + epsvecn / (x - alfa) - (xsi * dx) / (x - alfa)
|
|
deta = -eta + epsvecn / (beta - x) + (eta * dx) / (beta - x)
|
|
dmu = -mu + epsvecm / y - (mu * dy) / y
|
|
dzet = -zet + epsi / z - zet * dz / z
|
|
ds = -s + epsvecm / lam - (s * dlam) / lam
|
|
xx = np.concatenate((y, np.array([[z]]), lam, xsi, eta, mu, np.array([[zet]]), s), axis=0)
|
|
dxx = np.concatenate((dy, np.array([[dz]]), dlam, dxsi, deta, dmu, np.array([[dzet]]), ds), axis=0)
|
|
|
|
stepxx = -1.01 * dxx / xx
|
|
stmxx = np.amax(stepxx)
|
|
stepalfa = -1.01 * dx / (x - alfa)
|
|
stmalfa = np.amax(stepalfa)
|
|
stepbeta = 1.01 * dx / (beta - x)
|
|
stmbeta = np.amax(stepbeta)
|
|
stmalbe = np.maximum(stmalfa, stmbeta)
|
|
stmalbexx = np.maximum(stmalbe, stmxx)
|
|
stminv = np.maximum(stmalbexx, 1)
|
|
steg = 1 / stminv
|
|
|
|
xold = x
|
|
yold = y
|
|
zold = z
|
|
lamold = lam
|
|
xsiold = xsi
|
|
etaold = eta
|
|
muold = mu
|
|
zetold = zet
|
|
sold = s
|
|
|
|
itto = 0
|
|
resinew = 2 * residunorm
|
|
|
|
while (resinew > residunorm) and (itto < 50):
|
|
itto = itto + 1
|
|
x = xold + steg * dx
|
|
y = yold + steg * dy
|
|
z = zold + steg * dz
|
|
lam = lamold + steg * dlam
|
|
xsi = xsiold + steg * dxsi
|
|
eta = etaold + steg * deta
|
|
mu = muold + steg * dmu
|
|
zet = zetold + steg * dzet
|
|
s = sold + steg * ds
|
|
|
|
ux1 = upp - x
|
|
xl1 = x - low
|
|
ux2 = ux1 * ux1
|
|
xl2 = xl1 * xl1
|
|
|
|
uxinv1 = een / ux1
|
|
xlinv1 = een / xl1
|
|
plam = p0 + P.T @ lam
|
|
qlam = q0 + Q.T @ lam
|
|
gvec = P @ uxinv1 + Q @ xlinv1
|
|
|
|
dpsidx = plam / ux2 - qlam / xl2
|
|
rex = dpsidx - xsi + eta
|
|
rey = c + d * y - mu - lam
|
|
rez = a0 - zet - a.T @ lam
|
|
relam = gvec - a * z - y + s - b
|
|
rexsi = xsi * (x - alfa) - epsvecn
|
|
reeta = eta * (beta - x) - epsvecn
|
|
remu = mu * y - epsvecm
|
|
rezet = zet * z - epsi
|
|
res = lam * s - epsvecm
|
|
|
|
residu1 = np.concatenate((rex, rey, rez), axis=0)
|
|
residu2 = np.concatenate((relam, rexsi, reeta, remu, np.array([[rezet]]), res), axis=0)
|
|
residu = np.concatenate((residu1, residu2), axis=0)
|
|
resinew = np.linalg.norm(residu)
|
|
|
|
steg = steg / 2
|
|
|
|
residunorm = resinew
|
|
residumax = np.amax(np.absolute(residu))
|
|
steg = 2 * steg
|
|
|
|
if ittt > 198:
|
|
print(epsi)
|
|
print(ittt)
|
|
|
|
epsi = 0.1 * epsi
|
|
|
|
# store final output
|
|
xmma = x
|
|
ymma = y
|
|
zmma = z
|
|
lamma = lam
|
|
xsimma = xsi
|
|
etamma = eta
|
|
mumma = mu
|
|
zetmma = zet
|
|
smma = s
|
|
|
|
return xmma, ymma, zmma, lamma, xsimma, etamma, mumma, zetmma, smma
|
|
|
|
|
|
# The real main driver
|
|
if __name__ == "__main__":
|
|
# Default input parameters
|
|
|
|
#draw_x(x)
|
|
|
|
#初值0.1
|
|
x=[0.07,0.12,0.15,0.03,0.11,0.175,0.1,0.03,0.03,0.18]
|
|
#draw_x(x)
|
|
|
|
#初值0.36
|
|
x=[0.166,0.183,0.166,0.158,0.213,0.138,0.153,0.223,0.078,0.168]
|
|
|
|
nelx = 10
|
|
nely = 10
|
|
nelz = 10
|
|
|
|
rmin = 1.2
|
|
penal = 3.0
|
|
ft = 0 # ft==0 -> sens, ft==1 -> dens
|
|
default = 0.123
|
|
x = [default] * 10
|
|
x[0]=0.036
|
|
x[2]=0.036
|
|
#draw_x(x)
|
|
|
|
default = 0.1
|
|
x = [default] * 10
|
|
|
|
import sys
|
|
|
|
if len(sys.argv) > 1: nelx = int(sys.argv[1])
|
|
if len(sys.argv) > 2: nely = int(sys.argv[2])
|
|
if len(sys.argv) > 3: nelz = int(sys.argv[3])
|
|
if len(sys.argv) > 4: volfrac = float(sys.argv[4])
|
|
if len(sys.argv) > 5: rmin = float(sys.argv[5])
|
|
if len(sys.argv) > 6: penal = float(sys.argv[6])
|
|
if len(sys.argv) > 7: ft = int(sys.argv[7])
|
|
|
|
x=np.asarray(x)
|
|
main(nelx, nely, nelz, x,penal, rmin, ft)
|
|
print ("FINISH")
|
|
|
|
|
|
|