246 lines
9.8 KiB
Python
246 lines
9.8 KiB
Python
|
# -*- coding: utf-8 -*-
|
||
|
|
"""
|
||
|
|
Created on Wed Apr 22 17:40:00 2020
|
||
|
|
|
||
|
|
@author: zzy
|
||
|
|
|
||
|
|
Object: 对TO问题进行定义
|
||
|
|
1. 对于部分确定的static variable 应采用numpy定义
|
||
|
|
"""
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
import jax.numpy as jnp
|
||
|
|
import jax.ops as jops
|
||
|
|
import TO_Cal
|
||
|
|
import TO_FEA
|
||
|
|
|
||
|
|
|
||
|
|
def Toparams(problem):
|
||
|
|
young = 1.0
|
||
|
|
rou = 1.0
|
||
|
|
nelx, nely, ele_length, ele_width = problem.nelx, problem.nely, problem.ele_length, problem.ele_width
|
||
|
|
Total_ele, Total_nod, Total_dof = nelx * \
|
||
|
|
nely, (nelx+1)*(nely+1), 2*(nelx+1)*(nely+1)
|
||
|
|
gobalcoords0xy, M_elenod, M_edofMat, M_iK, M_jK = coordinate_and_preFEA(
|
||
|
|
nelx, nely, ele_length, ele_width, Total_ele, Total_nod)
|
||
|
|
elenod, edofMat, iK, jK = M_elenod-1, M_edofMat-1, M_iK-1, M_jK - \
|
||
|
|
1 # gobalcoords0xy是绝对坐标无需更改 其他根据索引方式均-1 纯粹的数值处理 后面索引用的IE也要从0起
|
||
|
|
edofMat_3D = TO_edofMat_3DTensor(nelx, nely) # 这是一个Tensor版本edofMat
|
||
|
|
design_map_bool, num_design_ele, num_nondesign_ele, design_area_indices, nondesign_area_indices = handle_design_region(
|
||
|
|
nelx, nely, Total_ele, problem.design_map)
|
||
|
|
params = {
|
||
|
|
'young': young,
|
||
|
|
'young_min': 1e-9*young,
|
||
|
|
'rou': rou,
|
||
|
|
'poisson': 0.3,
|
||
|
|
'g': 0,
|
||
|
|
'ele_length': ele_length,
|
||
|
|
'ele_width': ele_width,
|
||
|
|
'ele_thickness': problem.ele_thickness,
|
||
|
|
'nelx': nelx,
|
||
|
|
'nely': nely,
|
||
|
|
'Total_ele': Total_ele,
|
||
|
|
'Total_nod': Total_nod,
|
||
|
|
'Total_dof': Total_dof,
|
||
|
|
'volfrac': problem.volfrac,
|
||
|
|
'xmin': 0.001,
|
||
|
|
'xmax': 1.0,
|
||
|
|
'design_map': problem.design_map,
|
||
|
|
# 'design_map_bool': design_map_bool,
|
||
|
|
# 'num_design_ele': num_design_ele,
|
||
|
|
# 'num_nondesign_ele': num_nondesign_ele,
|
||
|
|
'design_area_indices': design_area_indices,
|
||
|
|
'nondesign_area_indices': nondesign_area_indices,
|
||
|
|
'freedofs': problem.freedofs,
|
||
|
|
'fixeddofs': problem.fixeddofs,
|
||
|
|
'fext': problem.fext,
|
||
|
|
'penal': 3.0,
|
||
|
|
'filter_width': problem.filter_width,
|
||
|
|
'gobalcoords0xy': gobalcoords0xy,
|
||
|
|
'M_elenod': M_elenod,
|
||
|
|
'M_edofMat': M_edofMat,
|
||
|
|
'M_iK': M_iK,
|
||
|
|
'M_jK': M_jK,
|
||
|
|
'elenod': elenod,
|
||
|
|
'edofMat': edofMat,
|
||
|
|
'edofMat_3D': edofMat_3D,
|
||
|
|
'iK': iK,
|
||
|
|
'jK': jK,
|
||
|
|
'design_condition': problem.design_condition,
|
||
|
|
}
|
||
|
|
return params
|
||
|
|
|
||
|
|
|
||
|
|
class Topology_Optimization:
|
||
|
|
|
||
|
|
def __init__(self, args):
|
||
|
|
self.args = args
|
||
|
|
self.young, self.young_min, self.poisson = args['young'], args['young_min'], args['poisson']
|
||
|
|
self.Total_dof, self.Total_ele, self.nelx, self.nely, self.ele_thickness = args[
|
||
|
|
'Total_dof'], args['Total_ele'], args['nelx'], args['nely'], args['ele_thickness']
|
||
|
|
self.freedofs, self.fext, self.iK, self.jK, self.edofMat = args[
|
||
|
|
'freedofs'], args['fext'], args['iK'], args['jK'], args['edofMat']
|
||
|
|
self.filter_width, self.design_condition, self.volfrac = args[
|
||
|
|
'filter_width'], args['design_condition'], args['volfrac']
|
||
|
|
self.design_map, self.design_area_indices, self.nondesign_area_indices = args[
|
||
|
|
'design_map'], args['design_area_indices'], args['nondesign_area_indices']
|
||
|
|
self.ke = TO_FEA.linear_stiffness_matrix(
|
||
|
|
self.young, self.poisson, self.ele_thickness)
|
||
|
|
self.filter_kernal, self.filter_weight = TO_Cal.filter_define(
|
||
|
|
self.design_map, self.filter_width)
|
||
|
|
self.s_K_predefined, self.dispTD_predefined, self.idx = jnp.zeros(
|
||
|
|
[self.Total_dof, self.Total_dof]), jnp.zeros([self.Total_dof]), jops.index[self.iK, self.jK]
|
||
|
|
self.loop_1, self.iCont = 0, 0
|
||
|
|
self.penal_use, self.beta_x_use = 3.0, 1.0
|
||
|
|
self.beta_x_use_max, self.beta_x_deltath1, self.beta_x_deltath2 = 16, 50, 30,
|
||
|
|
self.xPhys = self.volfrac * jnp.ones([self.nely, self.nelx])
|
||
|
|
|
||
|
|
def xPhys_transform(self, params, projection=False):
|
||
|
|
x = params.reshape(self.nely, self.nelx)
|
||
|
|
beta_x = self.penal_and_betax()
|
||
|
|
xPhys_1D = TO_Cal.design_variable(x, self.design_area_indices, self.nondesign_area_indices,
|
||
|
|
self.filter_kernal, self.filter_weight, beta_x, self.volfrac, projection=projection)
|
||
|
|
return xPhys_1D
|
||
|
|
|
||
|
|
def objective(self, params, projection=False):
|
||
|
|
self.loop_1 = self.loop_1 + 1
|
||
|
|
xPhys_1D = self.xPhys_transform(params, projection=projection)
|
||
|
|
Obj = TO_Cal.objective(xPhys_1D, self.young, self.young_min, self.freedofs, self.fext,
|
||
|
|
self.s_K_predefined, self.dispTD_predefined, self.idx, self.edofMat, self.ke, self.penal_use)
|
||
|
|
self.xPhys = xPhys_1D.reshape(self.nely, self.nelx, order='F')
|
||
|
|
return Obj
|
||
|
|
|
||
|
|
def penal_and_betax(self, *args):
|
||
|
|
self.iCont = self.iCont + 1
|
||
|
|
if (self.iCont > self.beta_x_deltath1 and self.beta_x_use < 2):
|
||
|
|
self.beta_x_use = self.beta_x_use * 2
|
||
|
|
self.iCont = 1
|
||
|
|
print(' beta_x increased to :%7.3f\n' % (self.beta_x_use))
|
||
|
|
elif (self.iCont > self.beta_x_deltath2 and self.beta_x_use < self.beta_x_use_max and self.beta_x_use >= 2):
|
||
|
|
self.beta_x_use = self.beta_x_use * 2
|
||
|
|
self.iCont = 1
|
||
|
|
print(' beta_x increased to :%7.3f\n' % (self.beta_x_use))
|
||
|
|
return self.beta_x_use
|
||
|
|
|
||
|
|
|
||
|
|
def coordinate_and_preFEA(nelx, nely, ele_length, ele_width, Total_ele, Total_nod):
|
||
|
|
"""
|
||
|
|
有限元计算相关量 单元及节点编号顺序为从左至右 从上至下
|
||
|
|
采用Matlab计数法 对单一单元内物理量 按左下逆时针顺序
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
nelx : 横向网格划分
|
||
|
|
nely : 纵向网格划分
|
||
|
|
ele_length : 单元长度
|
||
|
|
ele_width : 单元宽度
|
||
|
|
Total_ele : 单元总数
|
||
|
|
Total_nod : 节点总数
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
gobalcoords0xy : shape:[Total_nod, 2](2D) 按照节点顺序存储每一个节点的物理坐标(x, y)
|
||
|
|
elenod : shape:[Total_ele, 4](2D) 按照单元顺序存储每一个单元的节点编号
|
||
|
|
edofMat : shape:[Total_ele, 8](2D) 按照单元顺序存储每一个单元的自由度编号
|
||
|
|
iK_int : shape:[Total_ele * 64](2D) 用于刚度阵组装 按照特定顺序储存自由度编号
|
||
|
|
jK_int : shape:[Total_ele * 64](2D) 用于刚度阵组装 按照特定顺序储存自由度编号
|
||
|
|
"""
|
||
|
|
i0, j0 = jnp.meshgrid(jnp.arange(nelx+1), jnp.arange(nely+1))
|
||
|
|
gobalcoords0x = i0*ele_length
|
||
|
|
gobalcoords0y = ele_width*nely-j0*ele_width
|
||
|
|
gobalcoords0xy = jnp.hstack((matrix_array(gobalcoords0x), matrix_array(
|
||
|
|
gobalcoords0y)))
|
||
|
|
gobalcoords0xyT = gobalcoords0xy.T
|
||
|
|
gobalcoords0 = matrix_array(gobalcoords0xyT)
|
||
|
|
nodegrd = matrix_order_arrange(1, Total_nod, nely+1, nelx+1)
|
||
|
|
nodelast = Matlab_reshape(Matlab_matrix_extract(
|
||
|
|
nodegrd, 1, -1, 1, -1), Total_ele, 1)
|
||
|
|
edofVec = 2*matrix_array(nodelast)+1
|
||
|
|
elenod = jnp.tile(nodelast, (1, 4)) + \
|
||
|
|
jnp.tile(jnp.array([1, nely+2, nely+1, 0]), (Total_ele, 1))
|
||
|
|
edofMat = jnp.tile(edofVec, (1, 8))+jnp.tile(jnp.array(
|
||
|
|
[0, 1, 2*nely+2, 2*nely+3, 2*nely+0, 2*nely+1, -2, -1]), (Total_ele, 1))
|
||
|
|
iK = jnp.kron(edofMat, jnp.ones((8, 1))).flatten()
|
||
|
|
jK = jnp.kron(edofMat, jnp.ones((1, 8))).flatten()
|
||
|
|
iK_int = iK.astype(jnp.int32)
|
||
|
|
jK_int = jK.astype(jnp.int32)
|
||
|
|
return gobalcoords0xy, elenod, edofMat, iK_int, jK_int
|
||
|
|
|
||
|
|
|
||
|
|
def handle_design_region(nelx, nely, Total_ele, design_map):
|
||
|
|
"""
|
||
|
|
对设计域/非设计域相关参数进行预定义
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
nelx : 横向网格划分
|
||
|
|
nely : 纵向网格划分
|
||
|
|
Total_ele : 单元总数
|
||
|
|
design_map : 设计域和非设计域定义矩阵
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
design_map_bool : bool形式的设计域和非设计域定义矩阵
|
||
|
|
num_design_ele : 可设计单元总数
|
||
|
|
num_nondesign_ele : 非设计单元总数
|
||
|
|
design_area_indices : 设计域的单元编号索引
|
||
|
|
nondesign_area_indices : 非设计域的单元编号索引
|
||
|
|
"""
|
||
|
|
shape = (nely, nelx)
|
||
|
|
design_map_bool = np.broadcast_to(design_map, shape) > 0
|
||
|
|
num_design_ele = np.sum(design_map)
|
||
|
|
num_nondesign_ele = Total_ele - num_design_ele
|
||
|
|
|
||
|
|
design_map_bool_1D = design_map_bool.flatten(order='F')
|
||
|
|
all_indices = np.arange(Total_ele, dtype=jnp.int64)
|
||
|
|
design_area_indices = jnp.flatnonzero(
|
||
|
|
design_map_bool_1D, size=num_design_ele)
|
||
|
|
nondesign_area_indices = jnp.setdiff1d(
|
||
|
|
all_indices, design_area_indices, size=num_nondesign_ele, assume_unique=True) # setdiff1d通常不支持jit 需要指定size
|
||
|
|
return design_map_bool, num_design_ele, num_nondesign_ele, design_area_indices, nondesign_area_indices
|
||
|
|
|
||
|
|
|
||
|
|
def Matlab_matrix_array1D(matrixA):
|
||
|
|
"""
|
||
|
|
shape:[*] 1D array 按照matlab的排列
|
||
|
|
"""
|
||
|
|
return matrixA.T.flatten()
|
||
|
|
|
||
|
|
|
||
|
|
def matrix_array(matrixA):
|
||
|
|
"""
|
||
|
|
shape:[*, 1] 2D array 按照matlab的排列 matrixA(:)
|
||
|
|
"""
|
||
|
|
return matrixA.T.reshape(-1, 1)
|
||
|
|
|
||
|
|
|
||
|
|
def matrix_order_arrange(numbegin, numend, ydim, xdim):
|
||
|
|
'''
|
||
|
|
shape:[ydim, xdim] 2D array 实现一个元素按列顺序排列的矩阵 注意索引编号和matlab的差1
|
||
|
|
'''
|
||
|
|
matrixA = jnp.array([[i for i in range(numbegin, numend+1)]])
|
||
|
|
matrixA = matrixA.reshape(xdim, ydim)
|
||
|
|
matrixA = matrixA.T
|
||
|
|
return matrixA
|
||
|
|
|
||
|
|
|
||
|
|
def Matlab_reshape(matrixA, ydim, xdim):
|
||
|
|
'''
|
||
|
|
shape:[ydim, xdim] 2D array matlab的reshape效果
|
||
|
|
'''
|
||
|
|
return matrixA.reshape((ydim, xdim), order='F')
|
||
|
|
|
||
|
|
|
||
|
|
def Matlab_matrix_extract(matrixA, xbegin, xend, ybegin, yend):
|
||
|
|
matrixA = matrixA[ybegin-1:yend, xbegin-1:xend]
|
||
|
|
return matrixA
|
||
|
|
|
||
|
|
|
||
|
|
def TO_edofMat_3DTensor(nelx, nely):
|
||
|
|
ely, elx = jnp.meshgrid(jnp.arange(nely), jnp.arange(nelx)) # x, y coords
|
||
|
|
n1 = (nely+1)*(elx+0) + (ely+0)
|
||
|
|
n2 = (nely+1)*(elx+1) + (ely+0)
|
||
|
|
n3 = (nely+1)*(elx+1) + (ely+1)
|
||
|
|
n4 = (nely+1)*(elx+0) + (ely+1)
|
||
|
|
edofMat_3D = jnp.array(
|
||
|
|
[2*n4, 2*n4+1, 2*n3, 2*n3+1, 2*n2, 2*n2+1, 2*n1, 2*n1+1])
|
||
|
|
return edofMat_3D
|