204 lines
8.1 KiB
Python
204 lines
8.1 KiB
Python
# -*- coding: utf-8 -*-
|
|
'''
|
|
@ Copyright (c) 2022 by Zeyu Zhang, All Rights Reserved.
|
|
@ Author : Zeyu Zhang
|
|
@ Email : zhangzeyu_work@outlook.com
|
|
@ Date : 2022-10-25 09:30:27
|
|
@ LastEditTime : 2022-10-25 09:37:58
|
|
@ FilePath : /ZZY_CODE/Env_JAX/IDRL/Linear_TO/TO_Define.py
|
|
@
|
|
@ Description :
|
|
@ Reference :
|
|
'''
|
|
|
|
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
|
|
edofMat_3D = TO_edofMat_3DTensor(nelx, nely)
|
|
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):
|
|
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):
|
|
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):
|
|
return matrixA.T.flatten()
|
|
|
|
|
|
def matrix_array(matrixA):
|
|
return matrixA.T.reshape(-1, 1)
|
|
|
|
|
|
def matrix_order_arrange(numbegin, numend, ydim, xdim):
|
|
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):
|
|
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
|
|
<<<<<<< HEAD
|
|
=======
|
|
|
|
>>>>>>> 96d2b2f07c0a482615e28909e1bda8fb9d15edb2
|