From 10fcf8023ab6fdfdbda22e6d25aff97413ba8d23 Mon Sep 17 00:00:00 2001 From: "brighthe98@gmail.com" Date: Wed, 18 Dec 2024 09:30:43 +0800 Subject: [PATCH] update --- ...ear15_api.py => jy_external_gear15_api.py} | 8 +- .../linear_elasticity/jy_linear_exact_api.py | 190 ++++++++--------- app/soptx/linear_elasticity/test_sri.py | 40 ++++ app/soptx/soptx/filter/matrix.py | 19 +- app/soptx/soptx/opt/compliance.py | 14 +- app/soptx/soptx/opt/mma.py | 193 +++++++++--------- app/soptx/soptx/opt/utils.py | 110 +++++----- app/soptx/soptx/opt/volume.py | 4 +- app/soptx/soptx/pde/cantilever_3d.py | 11 +- app/soptx/soptx/solver/elastic_fem_solver.py | 43 ++-- app/soptx/soptx/tests/test_assemble.py | 4 +- app/soptx/soptx/tests/test_mma.py | 66 +++--- fealpy/fem/linear_elastic_integrator.py | 24 ++- fealpy/material/elastic_material.py | 1 - 14 files changed, 400 insertions(+), 327 deletions(-) rename app/soptx/linear_elasticity/{jy_gear15_api.py => jy_external_gear15_api.py} (99%) create mode 100644 app/soptx/linear_elasticity/test_sri.py diff --git a/app/soptx/linear_elasticity/jy_gear15_api.py b/app/soptx/linear_elasticity/jy_external_gear15_api.py similarity index 99% rename from app/soptx/linear_elasticity/jy_gear15_api.py rename to app/soptx/linear_elasticity/jy_external_gear15_api.py index 643408417..c3eb2fd88 100644 --- a/app/soptx/linear_elasticity/jy_gear15_api.py +++ b/app/soptx/linear_elasticity/jy_external_gear15_api.py @@ -1,3 +1,6 @@ +""" +外齿轮 +""" from fealpy.backend import backend_manager as bm from fealpy.functionspace import LagrangeFESpace, TensorFunctionSpace from fealpy.sparse import COOTensor @@ -490,11 +493,11 @@ def compute_equivalent_stress(stress, nu): F_load_nodes = F[dof_indices] # (15*8, GD) fixed_node_index = bm.where(is_inner_node)[0] -export_to_inp(filename='/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/inp/external_gear_analys.inp', +export_to_inp(filename='/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/inp/external_gear_abaqus.inp', nodes=node, elements=cell, fixed_nodes=fixed_node_index, load_nodes=load_node_indices, loads=F_load_nodes, young_modulus=206e3, poisson_ratio=0.3, density=7.85e-9, - used_app='ansys', mesh_type='hex') + used_app='abaqus', mesh_type='hex') E = 206e3 nu = 0.3 @@ -564,6 +567,7 @@ def compute_equivalent_stress(stress, nu): uh_show = uh.reshape(GD, NN).T else: uh_show = uh.reshape(NN, GD) + uh_x = uh_show[:, 0] uh_y = uh_show[:, 1] uh_z = uh_show[:, 2] diff --git a/app/soptx/linear_elasticity/jy_linear_exact_api.py b/app/soptx/linear_elasticity/jy_linear_exact_api.py index 7fe7a4975..5c212f15e 100644 --- a/app/soptx/linear_elasticity/jy_linear_exact_api.py +++ b/app/soptx/linear_elasticity/jy_linear_exact_api.py @@ -617,99 +617,99 @@ def dirichlet(self, points: TensorLike) -> TensorLike: strain5, stress5, nstrain5, nstress5 = compute_strain_stress_at_quadpoints3(tensor_space, uh, B_BBar, D) -# 单元积分点处的位移梯度 -q = 2 -qf = mesh.quadrature_formula(q) -bcs_quadrature, ws = qf.get_quadrature_points_and_weights() -tgphi = tensor_space.grad_basis(bcs_quadrature) # (NC, NQ, tldof, GD, GD) -tgrad = bm.einsum('cqimn, ci -> cqmn', tgphi, uh_cell) # (NC, NQ, GD, GD) - -# 应变张量 -strain = (tgrad + bm.transpose(tgrad, (0, 1, 3, 2))) / 2 # (NC, NQ, GD, GD) -strain_xx = strain[..., 0, 0] # (NC, NQ) -strain_yy = strain[..., 1, 1] # (NC, NQ) -strain_zz = strain[..., 2, 2] # (NC, NQ) -strain_xy = strain[..., 0, 1] # (NC, NQ) -strain_yz = strain[..., 1, 2] # (NC, NQ) -strain_xz = strain[..., 0, 2] # (NC, NQ) - -# 应力张量 -trace_e = bm.einsum("...ii", strain) # (NC, NQ) -I = bm.eye(GD, dtype=bm.float64) -stress = 2 * mu * strain + lam * trace_e[..., None, None] * I # (NC, NQ, GD, GD) -stress_xx = stress[..., 0, 0] # (NC, 8) -stress_yy = stress[..., 1, 1] # (NC, 8) -stress_zz = stress[..., 2, 2] # (NC, 8) -stress_xy = stress[..., 0, 1] # (NC, 8) -stress_yz = stress[..., 1, 2] # (NC, 8) -stress_xz = stress[..., 0, 2] # (NC, 8) - - -mesh.nodedata['uh'] = uh_show[:] -mesh.nodedata['uh_magnitude'] = uh_magnitude[:] -mesh.nodedata['strain_vertices'] = nstrain1 -mesh.nodedata['stress_vertices'] = nstress1 -mesh.nodedata['strian_centers'] = nstrain2 -mesh.nodedata['stress_centers'] = nstress2 -mesh.nodedata['strain_quadpoints'] = nstrain3 -mesh.nodedata['stress_quadpoints'] = nstress3 -mesh.to_vtk('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/box_hex_linear_exact_fealpy.vtu') - - -# 验证收敛阶 -maxit = 4 -errorType = ['$|| u_{bd} - uh_{bd}$', '$|| u - u_h ||_{L2}$', '$|| u - u_h||_{l2}$'] -errorMatrix = bm.zeros((len(errorType), maxit), dtype=bm.float64) -NDof = bm.zeros(maxit, dtype=bm.int32) -for i in range(maxit): - space = LagrangeFESpace(mesh, p=p, ctype='C') - # tensor_space = TensorFunctionSpace(space, shape=(3, -1)) # dof_priority - tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority - - NDof[i] = tensor_space.number_of_global_dofs() - - linear_elastic_material = LinearElasticMaterial(name='E_nu', - elastic_modulus=E, poisson_ratio=nu, - hypo='3D', device=bm.get_device(mesh)) - - integrator_K = LinearElasticIntegrator(material=linear_elastic_material, q=q, method='voigt') - bform = BilinearForm(tensor_space) - bform.add_integrator(integrator_K) - K = bform.assembly(format='csr') - - integrator_F = VectorSourceIntegrator(source=pde.source, q=q) - lform = LinearForm(tensor_space) - lform.add_integrator(integrator_F) - F = lform.assembly() - - dbc = DirichletBC(space=tensor_space, - gd=pde.dirichlet, - threshold=None, - method='interp') - uh_bd = bm.zeros(tensor_space.number_of_global_dofs(), - dtype=bm.float64, device=bm.get_device(mesh)) - uh_bd, isDDof = tensor_space.boundary_interpolate(gd=pde.dirichlet, - uh=uh_bd, threshold=None, method='interp') - F = F - K.matmul(uh_bd) - F = bm.set_at(F, isDDof, uh_bd[isDDof]) - - K = dbc.apply_matrix(matrix=K, check=True) - - uh = tensor_space.function() - # uh[:] = cg(K, F, maxiter=1000, atol=1e-8, rtol=1e-8) - uh[:] = spsolve(K, F, solver="mumps") - - u_exact = tensor_space.interpolate(pde.solution) - errorMatrix[0, i] = bm.sqrt(bm.sum(bm.abs(uh[:] - u_exact)**2 * (1 / NDof[i]))) - errorMatrix[1, i] = mesh.error(u=uh, v=pde.solution, q=tensor_space.p+3, power=2) - errorMatrix[2, i] = bm.sqrt(bm.sum(bm.abs(uh[isDDof] - u_exact[isDDof])**2 * (1 / NDof[i]))) - - if i < maxit-1: - mesh.uniform_refine() - -print("errorMatrix:\n", errorType, "\n", errorMatrix) -print("NDof:", NDof) -print("order_l2:\n", bm.log2(errorMatrix[0, :-1] / errorMatrix[0, 1:])) -print("order_L2:\n ", bm.log2(errorMatrix[1, :-1] / errorMatrix[1, 1:])) -print("----------------------") +# # 单元积分点处的位移梯度 +# q = 2 +# qf = mesh.quadrature_formula(q) +# bcs_quadrature, ws = qf.get_quadrature_points_and_weights() +# tgphi = tensor_space.grad_basis(bcs_quadrature) # (NC, NQ, tldof, GD, GD) +# tgrad = bm.einsum('cqimn, ci -> cqmn', tgphi, uh_cell) # (NC, NQ, GD, GD) + +# # 应变张量 +# strain = (tgrad + bm.transpose(tgrad, (0, 1, 3, 2))) / 2 # (NC, NQ, GD, GD) +# strain_xx = strain[..., 0, 0] # (NC, NQ) +# strain_yy = strain[..., 1, 1] # (NC, NQ) +# strain_zz = strain[..., 2, 2] # (NC, NQ) +# strain_xy = strain[..., 0, 1] # (NC, NQ) +# strain_yz = strain[..., 1, 2] # (NC, NQ) +# strain_xz = strain[..., 0, 2] # (NC, NQ) + +# # 应力张量 +# trace_e = bm.einsum("...ii", strain) # (NC, NQ) +# I = bm.eye(GD, dtype=bm.float64) +# stress = 2 * mu * strain + lam * trace_e[..., None, None] * I # (NC, NQ, GD, GD) +# stress_xx = stress[..., 0, 0] # (NC, 8) +# stress_yy = stress[..., 1, 1] # (NC, 8) +# stress_zz = stress[..., 2, 2] # (NC, 8) +# stress_xy = stress[..., 0, 1] # (NC, 8) +# stress_yz = stress[..., 1, 2] # (NC, 8) +# stress_xz = stress[..., 0, 2] # (NC, 8) + + +# mesh.nodedata['uh'] = uh_show[:] +# mesh.nodedata['uh_magnitude'] = uh_magnitude[:] +# mesh.nodedata['strain_vertices'] = nstrain1 +# mesh.nodedata['stress_vertices'] = nstress1 +# mesh.nodedata['strian_centers'] = nstrain2 +# mesh.nodedata['stress_centers'] = nstress2 +# mesh.nodedata['strain_quadpoints'] = nstrain3 +# mesh.nodedata['stress_quadpoints'] = nstress3 +# mesh.to_vtk('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/box_hex_linear_exact_fealpy.vtu') + + +# # 验证收敛阶 +# maxit = 4 +# errorType = ['$|| u_{bd} - uh_{bd}$', '$|| u - u_h ||_{L2}$', '$|| u - u_h||_{l2}$'] +# errorMatrix = bm.zeros((len(errorType), maxit), dtype=bm.float64) +# NDof = bm.zeros(maxit, dtype=bm.int32) +# for i in range(maxit): +# space = LagrangeFESpace(mesh, p=p, ctype='C') +# # tensor_space = TensorFunctionSpace(space, shape=(3, -1)) # dof_priority +# tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority + +# NDof[i] = tensor_space.number_of_global_dofs() + +# linear_elastic_material = LinearElasticMaterial(name='E_nu', +# elastic_modulus=E, poisson_ratio=nu, +# hypo='3D', device=bm.get_device(mesh)) + +# integrator_K = LinearElasticIntegrator(material=linear_elastic_material, q=q, method='voigt') +# bform = BilinearForm(tensor_space) +# bform.add_integrator(integrator_K) +# K = bform.assembly(format='csr') + +# integrator_F = VectorSourceIntegrator(source=pde.source, q=q) +# lform = LinearForm(tensor_space) +# lform.add_integrator(integrator_F) +# F = lform.assembly() + +# dbc = DirichletBC(space=tensor_space, +# gd=pde.dirichlet, +# threshold=None, +# method='interp') +# uh_bd = bm.zeros(tensor_space.number_of_global_dofs(), +# dtype=bm.float64, device=bm.get_device(mesh)) +# uh_bd, isDDof = tensor_space.boundary_interpolate(gd=pde.dirichlet, +# uh=uh_bd, threshold=None, method='interp') +# F = F - K.matmul(uh_bd) +# F = bm.set_at(F, isDDof, uh_bd[isDDof]) + +# K = dbc.apply_matrix(matrix=K, check=True) + +# uh = tensor_space.function() +# # uh[:] = cg(K, F, maxiter=1000, atol=1e-8, rtol=1e-8) +# uh[:] = spsolve(K, F, solver="mumps") + +# u_exact = tensor_space.interpolate(pde.solution) +# errorMatrix[0, i] = bm.sqrt(bm.sum(bm.abs(uh[:] - u_exact)**2 * (1 / NDof[i]))) +# errorMatrix[1, i] = mesh.error(u=uh, v=pde.solution, q=tensor_space.p+3, power=2) +# errorMatrix[2, i] = bm.sqrt(bm.sum(bm.abs(uh[isDDof] - u_exact[isDDof])**2 * (1 / NDof[i]))) + +# if i < maxit-1: +# mesh.uniform_refine() + +# print("errorMatrix:\n", errorType, "\n", errorMatrix) +# print("NDof:", NDof) +# print("order_l2:\n", bm.log2(errorMatrix[0, :-1] / errorMatrix[0, 1:])) +# print("order_L2:\n ", bm.log2(errorMatrix[1, :-1] / errorMatrix[1, 1:])) +# print("----------------------") diff --git a/app/soptx/linear_elasticity/test_sri.py b/app/soptx/linear_elasticity/test_sri.py new file mode 100644 index 000000000..102f8b704 --- /dev/null +++ b/app/soptx/linear_elasticity/test_sri.py @@ -0,0 +1,40 @@ +from fealpy.backend import backend_manager as bm + +from fealpy.mesh import HexahedronMesh, TetrahedronMesh +from fealpy.material.elastic_material import LinearElasticMaterial +from fealpy.fem.linear_elastic_integrator import LinearElasticIntegrator +from fealpy.fem.vector_source_integrator import VectorSourceIntegrator +from fealpy.fem.bilinear_form import BilinearForm +from fealpy.fem.linear_form import LinearForm +from fealpy.fem.dirichlet_bc import DirichletBC +from fealpy.functionspace import LagrangeFESpace, TensorFunctionSpace +from fealpy.typing import TensorLike +from fealpy.decorator import cartesian +from fealpy.sparse import COOTensor, CSRTensor +from fealpy.solver import cg, spsolve +# 刚度矩阵 +E = 2.1e5 +nu = 0.3 +# E = 1e6 +# nu = 0.25 +lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)) +mu = E / (2.0 * (1.0 + nu)) +mesh = HexahedronMesh.from_box(box=[0, 1, 0, 1, 0, 1], nx=1, ny=1, nz=1) +p = 1 +q = p+1 +space = LagrangeFESpace(mesh, p=p, ctype='C') +sgdof = space.number_of_global_dofs() +print(f"sgdof: {sgdof}") +tensor_space = TensorFunctionSpace(space, shape=(3, -1)) # dof_priority +linear_elastic_material = LinearElasticMaterial(name='E_nu', + elastic_modulus=E, poisson_ratio=nu, + hypo='3D', device=bm.get_device(mesh)) +integrator_K_sri = LinearElasticIntegrator(material=linear_elastic_material, + q=q, method='C3D8_SRI') +qf2 = mesh.quadrature_formula(q) +bcs2, ws2 = qf2.get_quadrature_points_and_weights() +gphi2 = space.grad_basis(bcs2, variable='x') + +# B0_q1 = linear_elastic_material._normal_strain_sri(gphi=gphi1) +KE_sri_yz_xz_xy = integrator_K_sri.c3d8_sri_assembly(space=tensor_space) +print("---------------------") \ No newline at end of file diff --git a/app/soptx/soptx/filter/matrix.py b/app/soptx/soptx/filter/matrix.py index a5182c707..5989f08e0 100644 --- a/app/soptx/soptx/filter/matrix.py +++ b/app/soptx/soptx/filter/matrix.py @@ -44,6 +44,7 @@ def _compute_filter_2d(nx: int, ny: int, rmin: float) -> Tuple[COOTensor, Tensor for i in range(nx): for j in range(ny): + # 单元的编号顺序: y->x row = i * ny + j kk1 = int(max(i - (ceil(rmin) - 1), 0)) kk2 = int(min(i + ceil(rmin), nx)) @@ -82,17 +83,19 @@ def _compute_filter_3d(nx: int, ny: int, nz: int, rmin: float) -> Tuple[TensorLi for i in range(nx): for j in range(ny): for k in range(nz): - row = i * ny * nz + j * nz + k + # 单元的编号顺序: z -> y -> x + row = k + j * nz + i * ny * nz ii1 = max(i - (ceil_rmin - 1), 0) - ii2 = min(i + ceil_rmin, nx) + ii2 = min(i + (ceil_rmin - 1), nx - 1) jj1 = max(j - (ceil_rmin - 1), 0) - jj2 = min(j + ceil_rmin, ny) + jj2 = min(j + (ceil_rmin - 1), ny - 1) kk1 = max(k - (ceil_rmin - 1), 0) - kk2 = min(k + ceil_rmin, nz) - for ii in range(ii1, ii2): - for jj in range(jj1, jj2): - for kk in range(kk1, kk2): - col = ii * ny * nz + jj * nz + kk + kk2 = min(k + (ceil_rmin - 1), nz - 1) + for ii in range(ii1, ii2 + 1): + for jj in range(jj1, jj2 + 1): + for kk in range(kk1, kk2 + 1): + # 单元的编号顺序: z -> y -> x + col = kk + jj * nz + ii * ny * nz fac = rmin - sqrt((i - ii)**2 + (j - jj)**2 + (k - kk)**2) iH[cc] = row jH[cc] = col diff --git a/app/soptx/soptx/opt/compliance.py b/app/soptx/soptx/opt/compliance.py index eba7fb8f2..1c68eccce 100644 --- a/app/soptx/soptx/opt/compliance.py +++ b/app/soptx/soptx/opt/compliance.py @@ -228,25 +228,25 @@ def jac(self, dc : 目标函数对密度的梯度 """ - # 创建计时器 - t = timer(f"Gradient Computation ({diff_mode} mode)") - next(t) + # # 创建计时器 + # t = timer(f"Gradient Computation ({diff_mode} mode)") + # next(t) # 选择计算方法 if diff_mode == "manual": dc = self._compute_gradient_manual(rho, u) - t.send('Manual gradient computed') + # t.send('Manual gradient computed') elif diff_mode == "auto": dc = self._compute_gradient_auto(rho) - t.send('Automatic gradient computed') + # t.send('Automatic gradient computed') # 应用滤波（如果需要） if self.filter is not None: dc = self.filter.filter_sensitivity(dc, rho, 'objective', filter_params) - t.send("Sensitivity filter applied") + # t.send("Sensitivity filter applied") # 结束计时 - t.send(None) + # t.send(None) return dc diff --git a/app/soptx/soptx/opt/mma.py b/app/soptx/soptx/opt/mma.py index 531d68106..642d6db33 100644 --- a/app/soptx/soptx/opt/mma.py +++ b/app/soptx/soptx/opt/mma.py @@ -19,7 +19,7 @@ class MMAOptions: n: int # 设计变量的数量 # 算法控制参数 max_iterations: int = 100 # 最大迭代次数 - tolerance: float = 0.01 # 收敛容差 + tolerance: float = 0.001 # 收敛容差 # MMA 子问题参数 a0: float = 1.0 # a_0*z 项的常数系数 a_0 a: Optional[TensorLike] = None # a_i*z 项的线性系数 a_i @@ -41,9 +41,9 @@ def log_iteration(self, iter_idx: int, obj_val: float, volume: float, self.densities.append(density.copy()) print(f"Iteration: {iter_idx + 1}, " - f"Objective: {obj_val:.3f}, " - f"Volume: {volume:.12f}, " - f"Change: {change:.3f}, " + f"Objective: {obj_val:.6f}, " + f"Volume: {volume:.6f}, " + f"Change: {change:.6f}, " f"Time: {time:.3f} sec") class MMAOptimizer(OptimizerBase): @@ -54,21 +54,21 @@ class MMAOptimizer(OptimizerBase): """ # MMA 算法的固定参数 - _MOVE_LIMIT = 0.01 # 移动限制 - _ASYMP_INIT = 0.01 # 渐近初始系数 - _ASYMP_INCR = 1.2 # 渐近递增系数 - _ASYMP_DECR = 0.4 # 渐近递减系数 - _ALBEFA = 0.1 # 渐近线移动系数 - _RAA0 = 1e-5 # 正则化参数 + _ASYMP_INIT = 0.5 # 渐近线初始距离的因子 + _ASYMP_INCR = 1.2 # 渐近线矩阵减小的因子 + _ASYMP_DECR = 0.7 # 渐近线矩阵增加的因子 + _MOVE_LIMIT = 0.2 # 移动限制 + _ALBEFA = 0.1 # 计算边界 alpha 和 beta 的因子 + _RAA0 = 1e-5 # 函数近似精度的参数 _EPSILON_MIN = 1e-7 # 最小容差 def __init__(self, - objective: ObjectiveBase, - constraint: ConstraintBase, - m: int, - n: int, - filter: Optional[Filter] = None, - options: Optional[Dict[str, Any]] = None): + objective: ObjectiveBase, + constraint: ConstraintBase, + m: int, + n: int, + filter: Optional[Filter] = None, + options: Optional[Dict[str, Any]] = None): """初始化 MMA 优化器 """ self.objective = objective self.constraint = constraint @@ -90,18 +90,22 @@ def __init__(self, self.options.c = 1e4 * bm.ones((m, 1)) if self.options.d is None: self.options.d = bm.zeros((m, 1)) + + # 初始化 xmin 和 xmax + self.xmin = bm.zeros((n, 1)) # 下界 + self.xmax = bm.ones((n, 1)) # 上界 # MMA 内部状态 self._epoch = 0 - self._xold1 = None - self._xold2 = None self._low = None self._upp = None def _update_asymptotes(self, xval: TensorLike, xmin: TensorLike, - xmax: TensorLike) -> Tuple[TensorLike, TensorLike]: + xmax: TensorLike, + xold1: TensorLike, + xold2: TensorLike) -> Tuple[TensorLike, TensorLike]: """更新渐近线位置""" asyinit = self._ASYMP_INIT asyincr = self._ASYMP_INCR @@ -110,22 +114,18 @@ def _update_asymptotes(self, xmami = xmax - xmin if self._epoch <= 2: - # 初始化渐近线 self._low = xval - asyinit * xmami self._upp = xval + asyinit * xmami else: - # 基于历史信息调整渐近线 factor = bm.ones((xval.shape[0], 1)) - xxx = (xval - self._xold1) * (self._xold1 - self._xold2) - # 根据变化趋势调整系数 + xxx = (xval - xold1) * (xold1 - xold2) factor[xxx > 0] = asyincr factor[xxx < 0] = asydecr + factor[xxx == 0] = 1.0 - # 更新渐近线位置 - self._low = xval - factor * (self._xold1 - self._low) - self._upp = xval + factor * (self._upp - self._xold1) + self._low = xval - factor * (xold1 - self._low) + self._upp = xval + factor * (self._upp - xold1) - # 限制渐近线范围 lowmin = xval - 10 * xmami lowmax = xval - 0.01 * xmami uppmin = xval + 0.01 * xmami @@ -143,8 +143,10 @@ def _solve_subproblem(self, fval: TensorLike, df0dx: TensorLike, dfdx: TensorLike, - xmin: TensorLike, - xmax: TensorLike) -> TensorLike: + low: TensorLike, + upp: TensorLike, + xold1: TensorLike, + xold2: TensorLike) -> TensorLike: """求解 MMA 子问题""" m = self.options.m # 使用配置的约束数量 n = self.options.n # 使用配置的设计变量数量 @@ -153,27 +155,35 @@ def _solve_subproblem(self, c = self.options.c d = self.options.d + move = self._MOVE_LIMIT + albefa = self._ALBEFA raa0 = self._RAA0 epsimin = self._EPSILON_MIN - move = self._MOVE_LIMIT + + xmin = self.xmin + xmax = self.xmax + + eeen = bm.ones((n, 1), dtype=bm.float64) + eeem = bm.ones((m, 1), dtype=bm.float64) # 更新渐近线 - low, upp = self._update_asymptotes(xval, xmin, xmax) + low, upp = self._update_asymptotes(xval, xmin, xmax, xold1, xold2) - # 计算移动限制 - alpha = bm.maximum(low + 0.1 * (xval - low), xval - move * (xmax - xmin)) - beta = bm.minimum(upp - 0.1 * (upp - xval), xval + move * (xmax - xmin)) - - # 一些辅助量 - eeen = bm.ones(n) - eeem = bm.ones((m, 1)) + # 计算边界 alpha, beta + xxx1 = low + albefa * (xval - low) + xxx2 = xval - move * (xmax - xmin) + xxx = bm.maximum(xxx1, xxx2) + alpha = bm.maximum(xmin, xxx) + xxx1 = upp - albefa * (upp - xval) + xxx2 = xval + move * (xmax - xmin) + xxx = bm.minimum(xxx1, xxx2) + beta = bm.minimum(xmax, xxx) - # 计算 xmami, xmamiinv 等参数 + # 计算 p0, q0, P, Q 和 b xmami = xmax - xmin - xmamieps = raa0 * eeen - xmami = bm.maximum(xmami, xmamieps) - xmamiinv = eeen / xmami - + xmami_eps = raa0 * eeen + xmami = bm.maximum(xmami, xmami_eps) + xmami_inv = eeen / xmami # 定义当前设计点 ux1 = upp - xval xl1 = xval - low @@ -181,63 +191,44 @@ def _solve_subproblem(self, xl2 = xl1 * xl1 uxinv = eeen / ux1 xlinv = eeen / xl1 - # 构建 p0, q0 - p0 = bm.maximum(df0dx, 0) # (NC, ) - q0 = bm.maximum(-df0dx, 0) # (NC, ) - pq0 = 0.001 * (p0 + q0) + raa0 * xmamiinv + p0 = bm.maximum(df0dx, 0) + q0 = bm.maximum(-df0dx, 0) + pq0 = 0.001 * (p0 + q0) + raa0 * xmami_inv p0 = p0 + pq0 q0 = q0 + pq0 p0 = p0 * ux2 q0 = q0 * xl2 - # 构建 P, Q - P = csr_matrix((m, n)) - Q = csr_matrix((m, n)) - P_data = bm.maximum(dfdx.reshape(m, n), 0) - Q_data = bm.maximum(-dfdx.reshape(m, n), 0) - P = csr_matrix(P_data) - Q = csr_matrix(Q_data) - - PQ = 0.001 * (P + Q) + raa0*(eeem @ xmamiinv[None,:]) + P = bm.zeros((m, n), dtype=bm.float64) + Q = bm.zeros((m, n), dtype=bm.float64) + P = bm.maximum(dfdx, 0) + Q = bm.maximum(-dfdx, 0) + PQ = 0.001 * (P + Q) + raa0 * bm.dot(eeem, xmami_inv.T) P = P + PQ Q = Q + PQ - - P = P @ spdiags(ux2, 0, n, n) - Q = Q @ spdiags(xl2, 0, n, n) - - # 计算 b - b = (P @ uxinv + Q @ xlinv - fval) + from numpy import diag as diags + P = (diags(ux2.flatten(), 0) @ P.T).T + Q = (diags(xl2.flatten(), 0) @ Q.T).T + b = bm.dot(P, uxinv) + bm.dot(Q, xlinv) - fval # 求解子问题 xmma, ymma, zmma, lam, xsi, eta, mu, zet, s = solve_mma_subproblem( - m, n, epsimin, low, upp, alpha, beta, - p0, q0, P, Q, - a0, a, b, c, d) + m, n, epsimin, low, upp, alpha, beta, + p0, q0, P, Q, + a0, a, b, c, d + ) - return xmma + return xmma.reshape(-1), low, upp def optimize(self, rho: TensorLike, **kwargs) -> Tuple[TensorLike, OptimizationHistory]: - """运行 MMA 优化算法 - - Parameters - ---------- - rho-(NC, ): 初始密度场 - **kwargs : 其他参数，例如： - - beta: Heaviside 投影参数 - - Returns - ------- - rho : 最优密度场 - history : 优化历史记录 - """ + """运行 MMA 优化算法""" # 获取参数 max_iters = self.options.max_iterations tol = self.options.tolerance - - # 设置变量界限 - xmin = bm.zeros_like(rho) - xmax = bm.ones_like(rho) + + low = bm.ones_like(rho) + upp = bm.ones_like(rho) # 准备 Heaviside 投影的参数 filter_params = {'beta': kwargs.get('beta')} if 'beta' in kwargs else None @@ -248,6 +239,9 @@ def optimize(self, rho: TensorLike, **kwargs) -> Tuple[TensorLike, OptimizationH # 初始化历史记录 history = OptimizationHistory() + + xold1 = bm.copy(rho) # 当前的设计变量 + xold2 = bm.copy(rho) # 初始化为当前的设计变量 # 优化主循环 for iter_idx in range(max_iters): @@ -256,36 +250,37 @@ def optimize(self, rho: TensorLike, **kwargs) -> Tuple[TensorLike, OptimizationH # 更新迭代计数 self._epoch = iter_idx + 1 - # 保存历史信息 - if self._xold1 is None: - self._xold1 = rho.copy() - self._xold2 = rho.copy() - else: - self._xold2 = self._xold1.copy() - self._xold1 = rho.copy() - # 计算目标函数值和梯度 obj_val = self.objective.fun(rho_phys) obj_grad = self.objective.jac(rho_phys) - # 计算约束值和梯度 + # 计算约束值和约束值梯度 con_val = self.constraint.fun(rho_phys) con_grad = self.constraint.jac(rho_phys) # 求解 MMA 子问题 - rho_new = self._solve_subproblem(rho, con_val, obj_grad, con_grad, xmin, xmax) + volfrac = self.constraint.volume_fraction + dfdx = con_grad[:, None].T / (volfrac * con_grad.shape[0]) + rho_new, low, upp = self._solve_subproblem(rho[:, None], + con_val, obj_grad[:, None], + dfdx, + low, upp, + xold1[:, None], xold2[..., None]) + + # 更新物理密度 + if self.filter is not None: + rho_phys = self.filter.filter_density(rho_new, filter_params) + else: + rho_phys = rho_new + + xold2 = xold1 + xold1 = rho # 计算收敛性 change = bm.max(bm.abs(rho_new - rho)) # 更新密度场 rho = rho_new - - # 更新物理密度 - if self.filter is not None: - rho_phys = self.filter.filter_density(rho, filter_params) - else: - rho_phys = rho # 记录当前迭代信息 iteration_time = time() - start_time diff --git a/app/soptx/soptx/opt/utils.py b/app/soptx/soptx/opt/utils.py index 6bbe9911b..0d25c34ae 100644 --- a/app/soptx/soptx/opt/utils.py +++ b/app/soptx/soptx/opt/utils.py @@ -6,8 +6,8 @@ """ from typing import Tuple -from scipy.sparse import diags -from scipy.linalg import solve +from numpy import diag as diags +from numpy.linalg import solve from fealpy.backend import backend_manager as bm from fealpy.typing import TensorLike @@ -70,27 +70,27 @@ def solve_mma_subproblem(m: int, epsvecn = epsi * een epsvecm = epsi * eem x = 0.5 * (alfa + beta) - y = eem.copy() + y = bm.copy(eem) z = bm.array([[1.0]]) - lam = eem.copy() + lam = bm.copy(eem) xsi = een / (x - alfa) xsi = bm.maximum(xsi, een) eta = een / (beta - x) eta = bm.maximum(eta, een) mu = bm.maximum(eem, 0.5*c) zet = bm.array([[1.0]]) - s = eem.copy() + s = bm.copy(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 + epsvecn = epsi * een + epsvecm = epsi * eem + ux1 = upp - x + xl1 = x - low + ux2 = ux1 * ux1 + xl2 = xl1 * xl1 + uxinv1 = een / ux1 + xlinv1 = een / xl1 # 计算梯度和残差 plam = p0 + bm.dot(P.T, lam) @@ -123,13 +123,13 @@ def solve_mma_subproblem(m: int, ux1 = upp-x xl1 = x-low ux2 = ux1*ux1 - xl2 = xl1*xl2 + xl2 = xl1*xl1 ux3 = ux1*ux2 xl3 = xl1*xl2 - uxinv1 = een/ux1 - xlinv1 = een/xl1 - uxinv2 = een/ux2 - xlinv2 = een/xl2 + uxinv1 = een / ux1 + xlinv1 = een / xl1 + uxinv2 = een / ux2 + xlinv2 = een / xl2 # 计算梯度矩阵和向量 plam = p0 + bm.dot(P.T, lam) @@ -142,7 +142,7 @@ def solve_mma_subproblem(m: int, dpsidx = plam/ux2 - qlam/xl2 delx = dpsidx - epsvecn/(x-alfa) + epsvecn/(beta-x) dely = c + d*y - lam - epsvecm/y - delz = a0 - np.dot(a.T, lam) - epsi/z + delz = a0 - bm.dot(a.T, lam) - epsi/z dellam = gvec - a*z - y - b + epsvecm/lam # 计算 Hessian 对角线 @@ -156,17 +156,17 @@ def solve_mma_subproblem(m: int, # 求解线性系统 if m < n: - blam = dellam + dely/diagy - bm.dot(GG, (delx/diagx)) + blam = dellam + dely / diagy - bm.dot(GG, (delx / diagx)) bb = bm.concatenate((blam, delz), axis=0) Alam = bm.asarray(diags(diaglamyi.flatten(), 0) + \ - (diags(diagxinv.flatten(), 0).dot(GG.T).T).dot(GG.T)) + (diags(diagxinv.flatten(), 0).dot(GG.T).T).dot(GG.T)) AAr1 = bm.concatenate((Alam, a), axis=1) AAr2 = bm.concatenate((a, -zet/z), axis=0).T AA = bm.concatenate((AAr1, AAr2), axis=0) solut = solve(AA, bb) dlam = solut[0:m] dz = solut[m:m+1] - dx = -delx/diagx - bm.dot(GG.T, dlam)/diagx + dx = -delx / diagx - bm.dot(GG.T, dlam) / diagx else: diaglamyiinv = eem/diaglamyi dellamyi = dellam + dely/diagy @@ -193,32 +193,31 @@ def solve_mma_subproblem(m: int, dmu = -mu + epsvecm/y - (mu*dy)/y dzet = -zet + epsi/z - zet*dz/z ds = -s + epsvecm/lam - (s*dlam)/lam - - # 计算步长 xx = bm.concatenate((y, z, lam, xsi, eta, mu, zet, s), axis=0) dxx = bm.concatenate((dy, dz, dlam, dxsi, deta, dmu, dzet, ds), axis=0) - stepxx = -1.01*dxx/xx + # 步长确定 + stepxx = -1.01 * dxx / xx stmxx = bm.max(stepxx) - stepalfa = -1.01*dx/(x-alfa) + stepalfa = -1.01 * dx / (x - alfa) stmalfa = bm.max(stepalfa) - stepbeta = 1.01*dx/(beta-x) + stepbeta = 1.01 * dx / (beta - x) stmbeta = bm.max(stepbeta) stmalbe = max(stmalfa, stmbeta) stmalbexx = max(stmalbe, stmxx) stminv = max(stmalbexx, 1.0) - steg = 1.0/stminv + steg = 1.0 / stminv # 保存旧值 - xold = x.copy() - yold = y.copy() - zold = z.copy() - lamold = lam.copy() - xsiold = xsi.copy() - etaold = eta.copy() - muold = mu.copy() - zetold = zet.copy() - sold = s.copy() + xold = bm.copy(x) + yold = bm.copy(y) + zold = bm.copy(z) + lamold = bm.copy(lam) + xsiold = bm.copy(xsi) + etaold = bm.copy(eta) + muold = bm.copy(mu) + zetold = bm.copy(zet) + sold = bm.copy(s) # 线搜索 itto = 0 @@ -234,14 +233,12 @@ def solve_mma_subproblem(m: int, 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 + ux1 = upp - x + xl1 = x - low + ux2 = ux1 * ux1 + xl2 = xl1 * xl1 + uxinv1 = een / ux1 + xlinv1 = een / xl1 plam = p0 + bm.dot(P.T, lam) qlam = q0 + bm.dot(Q.T, lam) gvec = bm.dot(P, uxinv1) + bm.dot(Q, xlinv1) @@ -252,10 +249,9 @@ def solve_mma_subproblem(m: int, 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 - + remu = mu * y - epsvecm + rezet = zet * z - epsi + res = lam * s - epsvecm residu1 = bm.concatenate((rex, rey, rez), axis=0) residu2 = bm.concatenate((relam, rexsi, reeta, remu, rezet, res), axis=0) residu = bm.concatenate((residu1, residu2), axis=0) @@ -264,15 +260,19 @@ def solve_mma_subproblem(m: int, residunorm = resinew.copy() residumax = bm.max(bm.abs(residu)) - steg = 2*steg + steg = 2 * steg - epsi = 0.1*epsi + epsi = 0.1 * epsi # 返回最优解 - xmma = x.copy() - ymma = y.copy() - zmma = z.copy() + xmma = bm.copy(x) + ymma = bm.copy(y) + zmma = bm.copy(z) lamma = lam xsimma = xsi etamma = eta - mumma = mu \ No newline at end of file + mumma = mu + zetmma = zet + smma = s + + return xmma, ymma, zmma, lamma, xsimma, etamma, mumma, zetmma, smma \ No newline at end of file diff --git a/app/soptx/soptx/opt/volume.py b/app/soptx/soptx/opt/volume.py index a88c7a299..7c8935d7b 100644 --- a/app/soptx/soptx/opt/volume.py +++ b/app/soptx/soptx/opt/volume.py @@ -78,7 +78,9 @@ def jac(self, return gradient # 明确指定这是约束函数的梯度 - return self.filter.filter_sensitivity(gradient, rho, 'constraint', filter_params) + filtered_gradient = self.filter.filter_sensitivity(gradient, rho, 'constraint', filter_params) + + return filtered_gradient def hess(self, rho: TensorLike, lambda_: Dict[str, Any]) -> TensorLike: """计算体积约束 Hessian 矩阵（未实现） diff --git a/app/soptx/soptx/pde/cantilever_3d.py b/app/soptx/soptx/pde/cantilever_3d.py index 53eb1c85a..918b98ef2 100644 --- a/app/soptx/soptx/pde/cantilever_3d.py +++ b/app/soptx/soptx/pde/cantilever_3d.py @@ -12,14 +12,15 @@ def __init__(self, ymin: float=0, ymax: float=20, zmin: float=0, zmax: float=4): """ - flip_direction = 'y' - 1------- 5 + 3------- 7 / | /| - 3 ------- 7 | + 1 ------- 5 | | | | | - | 0------|-4 + | 2------|-6 | / |/ - 2 ------- 6 + 0 ------- 4 + 位移边界条件: x 坐标为 0 的节点全部固定 + 载荷: x 坐标为 xmax, y 坐标为 ymin 的节点施加载荷 """ self.xmin, self.xmax = xmin, xmax self.ymin, self.ymax = ymin, ymax diff --git a/app/soptx/soptx/solver/elastic_fem_solver.py b/app/soptx/soptx/solver/elastic_fem_solver.py index 9712f8b75..3b1933ce7 100644 --- a/app/soptx/soptx/solver/elastic_fem_solver.py +++ b/app/soptx/soptx/solver/elastic_fem_solver.py @@ -26,6 +26,7 @@ class AssemblyMethod(Enum): """矩阵组装方法的枚举类""" STANDARD = auto() # 标准组装方法 VOIGT = auto() # Voigt 格式组装 + VOIGT_UNIFORM = auto # Voigt 格式一致网格组装 FAST_STRAIN = auto() # 应变快速组装 FAST_STRESS = auto() # 应力快速组装 FAST_3D = auto() # 3D 快速组装 @@ -33,7 +34,7 @@ class AssemblyMethod(Enum): @dataclass class AssemblyConfig: """矩阵组装的配置类""" - method: AssemblyMethod = AssemblyMethod.VOIGT + method: AssemblyMethod = AssemblyMethod.VOIGT_UNIFORM quadrature_degree_increase: int = 3 class ElasticFEMSolver: @@ -93,12 +94,7 @@ def current_density(self) -> Optional[TensorLike]: # 状态管理相关方法 #--------------------------------------------------------------------------- def update_density(self, density: TensorLike) -> None: - """更新密度并更新相关状态 - - Parameters - ---------- - density : 新的密度场 - """ + """更新密度并更新相关状态""" if density is None: raise ValueError("'density' cannot be None") @@ -134,9 +130,10 @@ def get_base_local_stiffness_matrix(self) -> TensorLike: base_material = self.material_properties.get_base_material() integrator = LinearElasticIntegrator( material=base_material, - q=self.tensor_space.p + 3 + q=self.tensor_space.p + 3, + method='voigt_uniform' ) - self._base_local_stiffness_matrix = integrator.assembly(space=self.tensor_space) + self._base_local_stiffness_matrix = integrator.voigt_uniform_assembly(space=self.tensor_space) return self._base_local_stiffness_matrix def compute_local_stiffness_matrix(self) -> TensorLike: @@ -151,6 +148,7 @@ def compute_local_stiffness_matrix(self) -> TensorLike: method_map = { AssemblyMethod.STANDARD: integrator.assembly, AssemblyMethod.VOIGT: integrator.voigt_assembly, + AssemblyMethod.VOIGT_UNIFORM: integrator.voigt_uniform_assembly, AssemblyMethod.FAST_STRAIN: integrator.fast_assembly_strain, AssemblyMethod.FAST_STRESS: integrator.fast_assembly_stress, AssemblyMethod.FAST_3D: integrator.fast_assembly @@ -178,6 +176,7 @@ def _create_integrator(self) -> LinearElasticIntegrator: method_map = { AssemblyMethod.STANDARD: 'assembly', AssemblyMethod.VOIGT: 'voigt', + AssemblyMethod.VOIGT_UNIFORM: 'voigt_uniform', AssemblyMethod.FAST_STRAIN: 'fast_strain', AssemblyMethod.FAST_STRESS: 'fast_stress', AssemblyMethod.FAST_3D: 'fast_3d' @@ -188,11 +187,10 @@ def _create_integrator(self) -> LinearElasticIntegrator: # 创建积分器 q = self.tensor_space.p + self.assembly_config.quadrature_degree_increase - integrator = LinearElasticIntegrator( - material=self._current_material, - q=q, - method=method - ) + integrator = LinearElasticIntegrator(material=self._current_material, + q=q, + method=method + ) return integrator @@ -208,7 +206,6 @@ def _assemble_global_stiffness_matrix(self) -> CSRTensor: # 创建适当的积分器 integrator = self._create_integrator() - bform = BilinearForm(self.tensor_space) bform.add_integrator(integrator) # t.send('Local Assembly') @@ -295,18 +292,18 @@ def solve_cg(self, if self._current_density is None: raise ValueError("Density not set. Call update_density first.") - # 创建计时器 - t = timer("CG Solver Timing") - next(t) # 启动计时器 + # # 创建计时器 + # t = timer("CG Solver Timing") + # next(t) # 启动计时器 K0 = self._assemble_global_stiffness_matrix() - t.send('Stiffness Matrix Assembly') + # t.send('Stiffness Matrix Assembly') F0 = self._assemble_global_force_vector() - t.send('Force Vector Assembly') + # t.send('Force Vector Assembly') K, F = self._apply_boundary_conditions(K0, F0) - t.send('Boundary Conditions') + # t.send('Boundary Conditions') uh = self.tensor_space.function() @@ -315,10 +312,10 @@ def solve_cg(self, # uh[:], info = cg(K, F[:], x0=x0, maxiter=maxiter, atol=atol, rtol=rtol) # TODO 目前 FEALPy 中的 cg 只能通过 logger 获取迭代步数，无法直接返回 uh[:] = cg(K, F[:], x0=x0, maxiter=maxiter, atol=atol, rtol=rtol) - t.send('Solving Phase') # 记录求解阶段时间 + # t.send('Solving Phase') # 记录求解阶段时间 # 结束计时 - t.send(None) + # t.send(None) except Exception as e: raise RuntimeError(f"CG solver failed: {str(e)}") diff --git a/app/soptx/soptx/tests/test_assemble.py b/app/soptx/soptx/tests/test_assemble.py index 4cb92dc64..339826b33 100644 --- a/app/soptx/soptx/tests/test_assemble.py +++ b/app/soptx/soptx/tests/test_assemble.py @@ -126,8 +126,8 @@ def create_base_components(config: TestConfig): tensor_space=tensor_space_C, pde=pde, assembly_config=assembly_config, - solver_type=config.solver_type, # 添加默认求解器类型 - solver_params=solver_params # 添加求解器参数 + solver_type=config.solver_type, + solver_params=solver_params ) # Initialize density field diff --git a/app/soptx/soptx/tests/test_mma.py b/app/soptx/soptx/tests/test_mma.py index 5a4ee3c70..b8ef44074 100644 --- a/app/soptx/soptx/tests/test_mma.py +++ b/app/soptx/soptx/tests/test_mma.py @@ -14,7 +14,7 @@ SIMPInterpolation ) from soptx.pde import MBBBeam2dData1, Cantilever2dData1, Cantilever3dData1 -from soptx.solver import ElasticFEMSolver +from soptx.solver import ElasticFEMSolver, AssemblyMethod, AssemblyConfig from soptx.opt import ComplianceObjective, VolumeConstraint from soptx.filter import Filter, FilterConfig from soptx.opt import MMAOptimizer @@ -43,6 +43,11 @@ class TestConfig: tolerance: float = 0.01 initial_lambda: float = 1e9 bisection_tol: float = 1e-3 + + assembly_method: AssemblyMethod = AssemblyMethod.VOIGT_UNIFORM # 矩阵组装方法 + quadrature_degree_increase: int = 3 # 积分阶数增量 + solver_type: Literal['cg', 'direct'] = 'cg' # 求解器类型 + solver_params: Optional[Dict[str, Any]] = None # 求解器参数 def create_base_components(config: TestConfig): @@ -56,7 +61,7 @@ def create_base_components(config: TestConfig): origin = [0.0, 0.0, 0.0] mesh = UniformMesh3d( extent=extent, h=h, origin=origin, - ipoints_ordering='zyx', flip_direction='y', + ipoints_ordering='zyx', flip_direction=None, device='cpu' ) dimension = 3 @@ -66,7 +71,7 @@ def create_base_components(config: TestConfig): origin = [0.0, 0.0] mesh = UniformMesh2d( extent=extent, h=h, origin=origin, - ipoints_ordering='yx', flip_direction='y', + ipoints_ordering='yx', flip_direction=None, device='cpu' ) dimension = 2 @@ -107,15 +112,26 @@ def create_base_components(config: TestConfig): ymin=0, ymax=config.ny*h[1], zmin=0, zmax=config.nz*h[2] ) + assembly_config = AssemblyConfig( + method=config.assembly_method, + quadrature_degree_increase=config.quadrature_degree_increase + ) - # Create solver + # 设置默认的求解器参数 + default_solver_params = { + 'cg': {'maxiter': 5000, 'atol': 1e-12, 'rtol': 1e-12}, + 'direct': {'solver_type': 'mumps'} + } + solver_params = config.solver_params or default_solver_params[config.solver_type] + solver = ElasticFEMSolver( - material_properties=material_properties, - tensor_space=tensor_space_C, - pde=pde, - solver_type='cg', # 添加默认求解器类型 - solver_params={'maxiter': 5000, 'atol': 1e-12, 'rtol': 1e-12} # 添加求解器参数 - ) + material_properties=material_properties, + tensor_space=tensor_space_C, + pde=pde, + assembly_config=assembly_config, + solver_type=config.solver_type, + solver_params=solver_params + ) # Initialize density field array = config.volume_fraction * bm.ones(mesh.number_of_cells(), dtype=bm.float64) @@ -130,9 +146,9 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]: # Create filter based on configuration filter_config = FilterConfig( - filter_type={'sensitivity': 0, 'density': 1, 'heaviside': 2}[config.filter_type], - filter_radius=config.filter_radius - ) + filter_type={'sensitivity': 0, 'density': 1, 'heaviside': 2}[config.filter_type], + filter_radius=config.filter_radius + ) filter_obj = Filter(filter_config) filter_obj.initialize(mesh) @@ -150,16 +166,16 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]: # Create optimizer optimizer = MMAOptimizer( - objective=objective, - constraint=constraint, - m=1, - n=mesh.number_of_cells(), - filter=filter_obj, - options={ - 'max_iterations': config.max_iterations, - 'tolerance': config.tolerance - } - ) + objective=objective, + constraint=constraint, + m=1, + n=mesh.number_of_cells(), + filter=filter_obj, + options={ + 'max_iterations': config.max_iterations, + 'tolerance': config.tolerance + } + ) # Prepare optimization parameters opt_params = {} @@ -167,7 +183,7 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]: opt_params['beta'] = config.projection_beta # Run optimization - rho_opt, history = optimizer.optimize(rho=rho[:], **opt_params) + rho_opt, history = optimizer.optimize(rho[:], **opt_params) return { 'optimal_density': rho_opt, @@ -198,7 +214,7 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]: config3 = TestConfig( problem_type='cantilever_3d', - nx=6, ny=2, nz=4, + nx=60, ny=20, nz=4, volume_fraction=0.3, filter_radius=1.5, filter_type='density', diff --git a/fealpy/fem/linear_elastic_integrator.py b/fealpy/fem/linear_elastic_integrator.py index 5bb34d6c2..9513bcc78 100644 --- a/fealpy/fem/linear_elastic_integrator.py +++ b/fealpy/fem/linear_elastic_integrator.py @@ -307,8 +307,7 @@ def voigt_uniform_assembly(self, space: _TS) -> TensorLike: B = self.material.strain_matrix(dof_priority=space.dof_priority, gphi=gphi) - KK = bm.einsum('q, c, cqki, cqkl, cqlj, cq -> cij', - ws, cm, B, D, B) + KK = bm.einsum('q, c, cqki, cqkl, cqlj -> cij', ws, cm, B, D, B) return KK @@ -646,8 +645,25 @@ def c3d8_bbar_assembly(self, space: _TS) -> TensorLike: @assemblymethod('C3D8_SRI') def c3d8_sri_assembly(self, space: _TS) -> TensorLike: scalar_space = space.scalar_space - bcs1, ws1, gphi1, detJ1, bcs2, ws2, gphi2, detJ2 = \ - self.fetch_c3d8_sri_assembly(scalar_space) + # bcs1, ws1, gphi1, detJ1, bcs2, ws2, gphi2, detJ2 = \ + # self.fetch_c3d8_sri_assembly(scalar_space) + + index = self.index + mesh = getattr(space, 'mesh', None) + + q1 = 1 + qf1 = mesh.quadrature_formula(q1) + bcs1, ws1 = qf1.get_quadrature_points_and_weights() + gphi1 = scalar_space.grad_basis(bcs1, index=index, variable='x') + J1 = mesh.jacobi_matrix(bcs1) + detJ1 = bm.linalg.det(J1) + + q2 = 2 + qf2 = mesh.quadrature_formula(q2) + bcs2, ws2 = qf2.get_quadrature_points_and_weights() + gphi2 = scalar_space.grad_basis(bcs2, index=index, variable='x') + J2 = mesh.jacobi_matrix(bcs2) + detJ2 = bm.linalg.det(J2) D_q1 = self.material.elastic_matrix(bcs1) D0_q1 = D_q1[..., :3, :3] # (1, 1, 3, 3) diff --git a/fealpy/material/elastic_material.py b/fealpy/material/elastic_material.py index 43e58d718..312d23876 100644 --- a/fealpy/material/elastic_material.py +++ b/fealpy/material/elastic_material.py @@ -286,7 +286,6 @@ def _normal_strain_sri(self, gphi: TensorLike, for i in range(GD): for j in range(GD): out = bm.set_at(out, (..., i, indices[:, j]), gphi[..., :, j]) - return out def _normal_strain_bbar(self, gphi: TensorLike,