Merge pull request #1409 from brighthe/develop

完善 SOptx 中的程序

app/soptx/linear_elasticity/jy_gear15_api.py → 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]

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("----------------------")
 

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("---------------------")

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

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