From 94586c19eafbd86ce40d0eab5e8711e461e82981 Mon Sep 17 00:00:00 2001 From: "brighthe98@gmail.com" Date: Tue, 3 Dec 2024 10:20:04 +0800 Subject: [PATCH 1/2] update --- app/soptx/linear_elasticity/test_abaqus.py | 99 ++++++++++++++++--- .../linear_elasticity/test_local_assembly.py | 78 ++++++++++----- 2 files changed, 140 insertions(+), 37 deletions(-) diff --git a/app/soptx/linear_elasticity/test_abaqus.py b/app/soptx/linear_elasticity/test_abaqus.py index ed0fd667d..7ce518b7a 100644 --- a/app/soptx/linear_elasticity/test_abaqus.py +++ b/app/soptx/linear_elasticity/test_abaqus.py @@ -11,6 +11,17 @@ from fealpy.solver import cg class BoxDomainPolyLoaded3d(): + def __init__(self): + """ + flip_direction = True + 0 ------- 3 ------- 6 + | 0 | 2 | + 1 ------- 4 ------- 7 + | 1 | 3 | + 2 ------- 5 ------- 8 + """ + self.eps = 1e-12 + def domain(self): return [0, 1, 0, 1, 0, 1] @@ -58,14 +69,77 @@ def dirichlet(self, points: TensorLike) -> TensorLike: return bm.zeros(points.shape, dtype=points.dtype, device=bm.get_device(points)) + + @cartesian + def is_dirichlet_boundary_dof_x(self, points: TensorLike) -> TensorLike: + domain = self.domain() + + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + + tag1 = bm.abs(x - domain[0]) < self.eps + tag2 = bm.abs(y - domain[0]) < self.eps + tag3 = bm.abs(z - domain[0]) < self.eps + + return tag1 & tag2 &tag3 + + @cartesian + def is_dirichlet_boundary_dof_y(self, points: TensorLike) -> TensorLike: + domain = self.domain() + + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + + tag1 = bm.abs(x - domain[0]) < self.eps + tag2 = bm.abs(y - domain[0]) < self.eps + tag3 = bm.abs(z - domain[0]) < self.eps + + return tag1 & tag2 &tag3 + + @cartesian + def is_dirichlet_boundary_dof_z(self, points: TensorLike) -> TensorLike: + domain = self.domain() + + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + + tag1 = bm.abs(x - domain[0]) < self.eps + tag2 = bm.abs(y - domain[0]) < self.eps + tag3 = bm.abs(z - domain[0]) < self.eps + + return tag1 & tag2 &tag3 + + @cartesian + def threshold(self): + + temp = (self.is_dirichlet_boundary_dof_x, + self.is_dirichlet_boundary_dof_y, + self.is_dirichlet_boundary_dof_z) + + return temp + + @cartesian + def is_bc(self, p): + x = p[..., 0] + y = p[..., 1] + z = p[..., 2] + + tag1 = bm.abs(x - 0) TensorLike: phi = space.basis(bcs) # (1, NQ, ldof) gphi = space.grad_basis(bc=bcs) # (NC, NQ, ldof, GD) -J = mesh.jacobi_matrix(bcs) -detJ = bm.linalg.det(J) - -tensor_space = TensorFunctionSpace(space, shape=(3, -1)) +tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority +# tensor_space = TensorFunctionSpace(space, shape=(3, -1)) # dof_priority tgdof = tensor_space.number_of_global_dofs() phi_tensor = tensor_space.basis(bcs) # (1, NQ, tldof, GD) cell2tdof = tensor_space.cell_to_dof() # (NC, tldof) @@ -90,7 +162,7 @@ def dirichlet(self, points: TensorLike) -> TensorLike: nu = 0.3 lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)) mu = E / (2.0 * (1.0 + nu)) -linear_elastic_material = LinearElasticMaterial(name='lam1_mu1', +linear_elastic_material = LinearElasticMaterial(name='E_nu', lame_lambda=lam, shear_modulus=mu, hypo='3D', device=bm.get_device(mesh)) @@ -98,10 +170,6 @@ def dirichlet(self, points: TensorLike) -> TensorLike: gphi=gphi, shear_order=['xy', 'yz', 'zx']) D = linear_elastic_material.elastic_matrix(bcs) KE = bm.einsum('q, c, cqki, cqkl, cqlj -> cij', ws, cm, B, D, B) -KE = bm.einsum('c, cqki, cqkl, cqlj -> cij', cm, B, D, B) -KE2 = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij', ws, detJ, B, D, B) - -error = bm.max(bm.abs(KE[0] - KE2[0])) integrator_K = LinearElasticIntegrator(material=linear_elastic_material, q=q) KE_maual = integrator_K.assembly(space=tensor_space) @@ -138,9 +206,11 @@ def dirichlet(self, points: TensorLike) -> TensorLike: 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) +isDDof = tensor_space.is_boundary_dof(threshold=pde.threshold(), method='interp') +# isDDof2 = tensor_space.is_boundary_dof(threshold=pde.is_bc, method='interp') +# isDDof = bm.zeros((24, ), dtype=bm.bool) +# isDDof[0::NN] = True - -isDDof = tensor_space.is_boundary_dof(threshold=None, method='interp') kwargs = K.values_context() # 1. 移除边界自由度相关的非零元素 indices = K.indices() @@ -161,8 +231,7 @@ def dirichlet(self, points: TensorLike) -> TensorLike: # 1. 边界插值 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') +uh_bd, _ = tensor_space.boundary_interpolate(gd=pde.dirichlet, uh=uh_bd, threshold=pde.threshold(), method='interp') # 2. 修改右端向量 F = F - K.matmul(uh_bd) F = bm.set_at(F, isDDof, uh_bd[isDDof]) diff --git a/app/soptx/linear_elasticity/test_local_assembly.py b/app/soptx/linear_elasticity/test_local_assembly.py index 0228718e4..d67182c8d 100644 --- a/app/soptx/linear_elasticity/test_local_assembly.py +++ b/app/soptx/linear_elasticity/test_local_assembly.py @@ -4,7 +4,9 @@ from fealpy.decorator import cartesian -from fealpy.mesh import TriangleMesh, TetrahedronMesh, QuadrangleMesh +from fealpy.mesh import (TriangleMesh, TetrahedronMesh, + QuadrangleMesh, HexahedronMesh, + UniformMesh2d, UniformMesh3d) from fealpy.functionspace import LagrangeFESpace, TensorFunctionSpace @@ -16,43 +18,74 @@ bm.set_backend('numpy') -nx, ny = 10, 10 -mesh_simplex = TriangleMesh.from_box(box=[0, 1, 0, 1], nx=nx, ny=ny) -mesh_tensor = QuadrangleMesh.from_box(box=[0, 1, 0, 1], nx=nx, ny=ny) +box_2d = [0, 1, 0, 1] +box_3d = [0, 1, 0, 1, 0, 1] +nx, ny, nz = 4, 4, 4 +extent_2d = [0, nx, 0, ny] +extent_3d = [0, nx, 0, ny, 0, nz] +h_2d = [(box_2d[1]-box_2d[0])/nx, (box_2d[3]-box_2d[2])/ny] +h_3d = [(box_3d[1]-box_3d[0])/nx, (box_3d[3]-box_3d[2])/ny, (box_3d[5]-box_3d[4])/nz] +origin_2d = [0, 0] +origin_3d = [0, 0, 0] +mesh_2d_simplex = TriangleMesh.from_box(box=box_2d, nx=nx, ny=ny) +mesh_3d_simplex = TetrahedronMesh.from_box(box=box_3d, nx=nx, ny=ny, nz=nx) +mesh_2d_tensor = QuadrangleMesh.from_box(box=box_2d, nx=nx, ny=ny) +mesh_3d_tensor = HexahedronMesh.from_box(box=box_3d, nx=nx, ny=ny, nz=nx) +mesh_2d_struct = UniformMesh2d(extent=extent_2d, h=h_2d, origin=origin_2d) +mesh_3d_struct = UniformMesh3d(extent=extent_3d, h=h_3d, origin=origin_3d) p = 4 -q = p+1 -qf_simplex = mesh_simplex.quadrature_formula(q) -bcs_simplex, _ = qf_simplex.get_quadrature_points_and_weights() -phi_simplex = mesh_simplex.shape_function(bcs=bcs_simplex, p=p) # (NQ, ldof) -qf_tensor = mesh_tensor.quadrature_formula(q) -bcs_tensor, _ = qf_tensor.get_quadrature_points_and_weights() -phi_tensor = mesh_tensor.shape_function(bcs=bcs_tensor, p=p) # (NQ, ldof) +space_2d = LagrangeFESpace(mesh_2d_simplex, p=p, ctype='C') +tensor_space_2d = TensorFunctionSpace(space_2d, shape=(-1, 2)) -space = LagrangeFESpace(mesh_simplex, p=p, ctype='C') -tensor_space = TensorFunctionSpace(space, shape=(-1, 2)) -tldof = tensor_space.number_of_global_dofs() -print(f"tldof: {tldof}") -linear_elastic_material = LinearElasticMaterial(name='lam1_mu1', +tldof_2d_simplex = tensor_space_2d_simplex.number_of_global_dofs() +print(f"tldof_2d_simplex: {tldof_2d_simplex}") + +linear_elastic_material_2d = LinearElasticMaterial(name='lam1_mu1', lame_lambda=1, shear_modulus=1, hypo='plane_stress') -integrator0 = LinearElasticIntegrator(material=linear_elastic_material, - q=tensor_space.p+3) +linear_elastic_material_3d = LinearElasticMaterial(name='lam1_mu1', + lame_lambda=1, shear_modulus=1, + hypo='3D') + +integrator_standard = LinearElasticIntegrator(material=linear_elastic_material_2d, q=p+3) + integrator1 = LinearElasticIntegrator(material=linear_elastic_material, - q=tensor_space.p+3, method='fast_stress') + q=p+3, method='fast_stress') integrator2 = LinearElasticIntegrator(material=linear_elastic_material, - q=tensor_space.p+3, method='symbolic_stress') + q=p+3, method='symbolic_stress') integrator1.keep_data() integrator2.keep_data() # 保留中间数据 # integrator2.keep_result() # 保留积分结果 # 创建计时器 -t = timer("Local Assembly Timing") +t = timer("2d Local Assembly Timing") +next(t) # 启动计时器 + +KE_2d_standard = integrator_standard.assembly(space=tensor_space_2d) +t.send('Assembly') + +KE_2d_fast1 = integrator1.fast_assembly_stress(space=tensor_space) +t.send('Fast Assembly1') +KE12 = integrator1.fast_assembly_stress(space=tensor_space) +t.send('Fast Assembly2') + +KE21 = integrator2.symbolic_assembly_stress(space=tensor_space) +t.send('Symbolic Assembly1') +KE22 = integrator2.symbolic_assembly_stress(space=tensor_space) +t.send('Symbolic Assembly2') +KE23 = integrator2.symbolic_assembly_stress(space=tensor_space) +t.send('Symbolic Assembly3') +# 结束计时 +t.send(None) + +t = timer("3d Local Assembly Timing") next(t) # 启动计时器 -KE0 = integrator0.assembly(space=tensor_space) +KE_2d_standard = integrator_standard.assembly(space=tensor_space) t.send('Assembly') + KE11 = integrator1.fast_assembly_stress(space=tensor_space) t.send('Fast Assembly1') KE12 = integrator1.fast_assembly_stress(space=tensor_space) @@ -66,4 +99,5 @@ t.send('Symbolic Assembly3') # 结束计时 t.send(None) + print("-------------------------------") \ No newline at end of file From af4c5a0d4e64a077d6f21b63511a05c90c91f5f4 Mon Sep 17 00:00:00 2001 From: "brighthe98@gmail.com" Date: Tue, 3 Dec 2024 21:15:50 +0800 Subject: [PATCH 2/2] update --- ...p_detail_linear_element_hexahedron_mesh.py | 3 + app/soptx/linear_elasticity/jy_box_api.py | 228 ++++++++++++++++++ .../{test_abaqus.py => jy_box_details.py} | 26 +- app/soptx/linear_elasticity/jy_gear_api.py | 152 ++++++++++++ fealpy/functionspace/tensor_space.py | 30 ++- fealpy/material/elastic_material.py | 41 +--- 6 files changed, 422 insertions(+), 58 deletions(-) create mode 100644 app/soptx/linear_elasticity/jy_box_api.py rename app/soptx/linear_elasticity/{test_abaqus.py => jy_box_details.py} (93%) create mode 100644 app/soptx/linear_elasticity/jy_gear_api.py diff --git a/app/soptx/linear_elasticity/exp_detail_linear_element_hexahedron_mesh.py b/app/soptx/linear_elasticity/exp_detail_linear_element_hexahedron_mesh.py index fa774a905..db1817be6 100644 --- a/app/soptx/linear_elasticity/exp_detail_linear_element_hexahedron_mesh.py +++ b/app/soptx/linear_elasticity/exp_detail_linear_element_hexahedron_mesh.py @@ -141,6 +141,9 @@ def dirichlet(self, points: TensorLike) -> TensorLike: uh = tensor_space.function() uh[:] = cg(K, F, maxiter=1000, atol=1e-14, rtol=1e-14) + + + u_exact = tensor_space.interpolate(pde.solution) error = mesh.error(u=uh, v=pde.solution, q=tensor_space.p+3, power=2) print("----------------------") diff --git a/app/soptx/linear_elasticity/jy_box_api.py b/app/soptx/linear_elasticity/jy_box_api.py new file mode 100644 index 000000000..920dbc328 --- /dev/null +++ b/app/soptx/linear_elasticity/jy_box_api.py @@ -0,0 +1,228 @@ +from fealpy.backend import backend_manager as bm + +from fealpy.mesh import HexahedronMesh +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 +from fealpy.solver import cg + +class BoxDomainPolyLoaded3d(): + def __init__(self): + """ + flip_direction = True + 0 ------- 3 ------- 6 + | 0 | 2 | + 1 ------- 4 ------- 7 + | 1 | 3 | + 2 ------- 5 ------- 8 + """ + self.eps = 1e-12 + + def domain(self): + return [0, 1, 0, 1, 0, 1] + + @cartesian + def source(self, points: TensorLike): + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + val = bm.zeros(points.shape, + dtype=points.dtype, device=bm.get_device(points)) + mu = 1 + factor1 = -400 * mu * (2 * y - 1) * (2 * z - 1) + term1 = 3 * (x ** 2 - x) ** 2 * (y ** 2 - y + z ** 2 - z) + term2 = (1 - 6 * x + 6 * x ** 2) * (y ** 2 - y) * (z ** 2 - z) + val[..., 0] = factor1 * (term1 + term2) + + factor2 = 200 * mu * (2 * x - 1) * (2 * z - 1) + term1 = 3 * (y ** 2 - y) ** 2 * (x ** 2 - x + z ** 2 - z) + term2 = (1 - 6 * y + 6 * y ** 2) * (x ** 2 - x) * (z ** 2 - z) + val[..., 1] = factor2 * (term1 + term2) + + factor3 = 200 * mu * (2 * x - 1) * (2 * y - 1) + term1 = 3 * (z ** 2 - z) ** 2 * (x ** 2 - x + y ** 2 - y) + term2 = (1 - 6 * z + 6 * z ** 2) * (x ** 2 - x) * (y ** 2 - y) + val[..., 2] = factor3 * (term1 + term2) + + return val + + @cartesian + def solution(self, points: TensorLike): + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + val = bm.zeros(points.shape, + dtype=points.dtype, device=bm.get_device(points)) + + mu = 1 + val[..., 0] = 200*mu*(x-x**2)**2 * (2*y**3-3*y**2+y) * (2*z**3-3*z**2+z) + val[..., 1] = -100*mu*(y-y**2)**2 * (2*x**3-3*x**2+x) * (2*z**3-3*z**2+z) + val[..., 2] = -100*mu*(z-z**2)**2 * (2*y**3-3*y**2+y) * (2*x**3-3*x**2+x) + + return val + + def dirichlet(self, points: TensorLike) -> TensorLike: + + result = bm.zeros(points.shape, + dtype=points.dtype, device=bm.get_device(points)) + + return result + + @cartesian + def is_dirichlet_boundary_dof_x(self, points: TensorLike) -> TensorLike: + domain = self.domain() + + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + + tag1 = bm.abs(x - domain[0]) < self.eps + tag2 = bm.abs(y - domain[0]) < self.eps + tag3 = bm.abs(z - domain[0]) < self.eps + + return tag1 & tag2 & tag3 + + @cartesian + def is_dirichlet_boundary_dof_y(self, points: TensorLike) -> TensorLike: + domain = self.domain() + + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + + tag1 = bm.abs(x - domain[0]) < self.eps + tag2 = bm.abs(y - domain[0]) < self.eps + tag3 = bm.abs(z - domain[0]) < self.eps + + return tag1 & tag2 & tag3 + @cartesian + def is_dirichlet_boundary_dof_z(self, points: TensorLike) -> TensorLike: + domain = self.domain() + + x = points[..., 0] + y = points[..., 1] + z = points[..., 2] + + tag1 = bm.abs(x - domain[0]) < self.eps + tag2 = bm.abs(y - domain[0]) < self.eps + tag3 = bm.abs(z - domain[0]) < self.eps + + return tag1 & tag2 & tag3 + + @cartesian + def threshold(self): + + temp = (self.is_dirichlet_boundary_dof_x, + self.is_dirichlet_boundary_dof_y, + self.is_dirichlet_boundary_dof_z) + + return temp + + +bm.set_backend('numpy') +nx, ny, nz = 1, 1, 1 +mesh = HexahedronMesh.from_box(box=[0, 1, 0, 1, 0, 1], + nx=nx, ny=ny, nz=nz, device=bm.get_device('cpu')) +GD = mesh.geo_dimension() +NN = mesh.number_of_nodes() +NC = mesh.number_of_cells() +cm = mesh.cell_volume() +node = mesh.entity('node') +cell = mesh.entity('cell') + +space = LagrangeFESpace(mesh, p=1, ctype='C') +cell2dof = space.cell_to_dof() + +q = 2 + +# tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority +tensor_space = TensorFunctionSpace(space, shape=(3, -1)) # dof_priority +print(f"dofs = {tensor_space.dof_priority}") +E = 206e3 +nu = 0.3 +lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)) +mu = E / (2.0 * (1.0 + nu)) +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) +KE = integrator_K.assembly(space=tensor_space) +bform = BilinearForm(tensor_space) +bform.add_integrator(integrator_K) +K = bform.assembly(format='csr') +# print(f"K.shape = {K.shape}:\n {K.to_dense()}, ") + +pde = BoxDomainPolyLoaded3d() + +ip1 = mesh.interpolation_points(p=1) +integrator_F = VectorSourceIntegrator(source=pde.source, q=q) +FE = integrator_F.assembly(space=tensor_space) +lform = LinearForm(tensor_space) +lform.add_integrator(integrator_F) +F = lform.assembly() +print(f"F.shape = {F.shape}:\n {F}, ") + +from app.gearx.utils import * + +if tensor_space.dof_priority == True: + F_load_nodes = bm.transpose(F.reshape(3, -1)) +else: + F_load_nodes = F.reshape(NN, GD) +print(f"F_load_nodes.shape = {F_load_nodes.shape}:\n {F_load_nodes}, ") + +load_node_indices = cell[0] +fixed_node_index = bm.tensor([0]) +export_to_inp(filename='/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/box.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) + + +dbc = DirichletBC(space=tensor_space, + gd=pde.dirichlet, + threshold=pde.threshold(), + method='interp') +K = dbc.apply_matrix(matrix=K, check=True) +# print(f"K.shape = {K.shape}:\n {K.to_dense()}, ") + +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=pde.threshold(), method='interp') + +# 2. 修改右端向量 +F = F - K.matmul(uh_bd) +F = bm.set_at(F, isDDof, uh_bd[isDDof]) +print(f"F.shape = {F.shape}:\n {F}, ") + +uh = tensor_space.function() +uh[:] = cg(K, F, maxiter=1000, atol=1e-14, rtol=1e-14) +uh_dof_show = uh.reshape(GD, NN).T +print(f"uh_dof_show.shape = {uh_dof_show.shape}:\n {uh_dof_show}, ") +uh_magnitude = bm.linalg.norm(uh_dof_show, axis=1) +print(f"uh_magnitude = {uh_magnitude}") + + +qf = mesh.quadrature_formula(1) +bcs, ws = qf.get_quadrature_points_and_weights() +phi = space.basis(bcs) # (1, 1, ldof) +gphi = space.grad_basis(bc=bcs) # (NC, 1, ldof, GD) +B = linear_elastic_material.strain_matrix(dof_priority=True, + gphi=gphi, shear_order=['xy', 'yz', 'zx']) # (NC, 1, 6, tldof) +print(f"B.shape = {B.shape}:\n {B}, ") +cell2tdof = tensor_space.cell_to_dof() # (NC, tldof) +tldof = tensor_space.number_of_local_dofs() +uh_cell = bm.zeros((NC, tldof)) +for c in range(NC): + uh_cell[c] = uh[cell2tdof[c]] +print(f"uh_cell.shape = {uh_cell.shape}:\n {uh_cell}, ") +strain = bm.einsum('cqij, cj -> ci', B, uh_cell) # (NC, 6) +print(f"strain.shape = {strain.shape}:\n {strain}, ") + + +print("----------------------") \ No newline at end of file diff --git a/app/soptx/linear_elasticity/test_abaqus.py b/app/soptx/linear_elasticity/jy_box_details.py similarity index 93% rename from app/soptx/linear_elasticity/test_abaqus.py rename to app/soptx/linear_elasticity/jy_box_details.py index 7ce518b7a..5c97ba290 100644 --- a/app/soptx/linear_elasticity/test_abaqus.py +++ b/app/soptx/linear_elasticity/jy_box_details.py @@ -120,25 +120,12 @@ def threshold(self): self.is_dirichlet_boundary_dof_z) return temp - - @cartesian - def is_bc(self, p): - x = p[..., 0] - y = p[..., 1] - z = p[..., 2] - - tag1 = bm.abs(x - 0) pd', load_values, target_cellnorm) # (15, GD) + +u = parameters[..., 0] +v = parameters[..., 1] +w = parameters[..., 2] + +bcs_list = [ + ( + bm.tensor([[u, 1 - u]]), + bm.tensor([[v, 1 - v]]), + bm.tensor([[w, 1 - w]]) + ) + for u, v, w in zip(u, v, w) +] + +space = LagrangeFESpace(mesh, p=1, ctype='C') + +q = 2 + +# tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority +tensor_space = TensorFunctionSpace(space, shape=(3, -1)) # dof_priority + +tgdof = tensor_space.number_of_global_dofs() +tldof = tensor_space.number_of_local_dofs() +cell2tdof = tensor_space.cell_to_dof() + +phi_loads = [] +for bcs in bcs_list: + phi = tensor_space.basis(bcs) + phi_loads.append(phi) + +phi_loads_array = bm.concatenate(phi_loads, axis=1) # (1, NP, tldof, GD) + +FE_load = bm.einsum('pd, cpld -> pl', P, phi_loads_array) # (NP, 24) + +FE = bm.zeros((NC, tldof), dtype=bm.float64) +FE[target_cell_idx, :] = FE_load[:, :] # (NC, tldof) + +F = COOTensor(indices = bm.empty((1, 0), dtype=bm.int32, device=bm.get_device(space)), + values = bm.empty((0, ), dtype=bm.float64, device=bm.get_device(space)), + spshape = (tgdof, )) +indices = cell2tdof.reshape(1, -1) +F = F.add(COOTensor(indices, FE.reshape(-1), (tgdof, ))).to_dense() # (tgdof, ) + +E = 206e3 +nu = 0.3 +lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)) +mu = E / (2.0 * (1.0 + nu)) +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) + +KE = integrator_K.assembly(space=tensor_space) +bform = BilinearForm(tensor_space) +bform.add_integrator(integrator_K) +K = bform.assembly(format='csr') + +from app.gearx.utils import * +F_load_nodes = bm.transpose(F.reshape(3, -1)) + +@cartesian +def dirichlet(points: TensorLike) -> TensorLike: + result = bm.zeros(points.shape, + dtype=points.dtype, device=bm.get_device(points)) + + return result + +scalar_is_bd_dof = is_inner_node +scalar_num = bm.sum(scalar_is_bd_dof) +tensor_is_bd_dof = tensor_space.is_boundary_dof( + threshold=(scalar_is_bd_dof, scalar_is_bd_dof, scalar_is_bd_dof), + method='interp') + +dbc = DirichletBC(space=tensor_space, + gd=dirichlet, + threshold=tensor_is_bd_dof, + method='interp') + +K = dbc.apply_matrix(matrix=K, check=True) + +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=dirichlet, + uh=uh_bd, + threshold=tensor_is_bd_dof, + method='interp') + +F = F - K.matmul(uh_bd) +F = bm.set_at(F, isDDof, uh_bd[isDDof]) + +from fealpy import logger +logger.setLevel('INFO') + +uh = tensor_space.function() +uh[:] = cg(K, F, maxiter=10000, atol=1e-8, rtol=1e-8) +# uh[:] = spsolve(K, F, solver='mumps') + +# 计算残差向量和范数 +residual = K.matmul(uh[:]) - F # 使用 CSRTensor 的 matmul 方法 +residual_norm = bm.sqrt(bm.sum(residual * residual)) +print(f"Final residual norm: {residual_norm:.6e}") + +uh = uh.reshape(GD, NN).T + +mesh.nodedata['deform'] = uh[:] +mesh.to_vtk('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/gear.vtu') +print("-----------") \ No newline at end of file diff --git a/fealpy/functionspace/tensor_space.py b/fealpy/functionspace/tensor_space.py index 5d4646010..888618c7e 100644 --- a/fealpy/functionspace/tensor_space.py +++ b/fealpy/functionspace/tensor_space.py @@ -332,8 +332,9 @@ def boundary_interpolate(self, dof_threshold)) elif method == 'interp': assert len(threshold) == self.dof_numel - isScalarBDof = [scalar_space.is_boundary_dof(i, method=method) for i in threshold] - gd_tensor = [gd(ipoints[isScalarBDof[i]])[...,i] for i in range(self.dof_numel)] + isScalarBDof = [scalar_space.is_boundary_dof(i, method=method) for i in threshold] + + gd_tensor = [gd(ipoints[isScalarBDof[i]])[..., i] for i in range(self.dof_numel)] isTensorBDof = self.is_boundary_dof(threshold=threshold, method=method) if uh is None: uh = self.function() @@ -341,9 +342,32 @@ def boundary_interpolate(self, gd_tensor = bm.concatenate(gd_tensor) else: scalar_gdof = scalar_space.number_of_global_dofs() - gd_tensor = bm.concatenate([bm.array([j[i] for j in isScalarBDof]) for i in range(scalar_gdof)]) + # 使用列表推导式重新排列数据 + gd_values = [] + for i in range(scalar_gdof): + for j in range(self.dof_numel): + # 从 gd_tensor[j] 中获取第 i 个标量基函数的值 + gd_values.append(gd_tensor[j][i] if i < len(gd_tensor[j]) else 0.0) + gd_tensor = bm.array(gd_values) + gd_tensor = gd_tensor[isTensorBDof] + uh[:] = bm.set_at(uh[:], isTensorBDof, gd_tensor) + return uh, isTensorBDof + + # gd_tensor = [gd(ipoints[isScalarBDof[i]]) for i in range(self.dof_numel)] + # gd_tensor = bm.concatenate(gd_tensor, axis=0) + + # isTensorBDof = self.is_boundary_dof(threshold=threshold, method=method) + # index = bm.where(isTensorBDof) + # if uh is None: + # uh = self.function() + # if self.dof_priority: + # gd_tensor = gd_tensor.T.reshape(-1) + # else: + # gd_tensor = gd_tensor.reshape(-1) + # uh[index] = gd_tensor[index] + # return uh, isTensorBDof else: raise ValueError(f"Unknown method: {method}") else: diff --git a/fealpy/material/elastic_material.py b/fealpy/material/elastic_material.py index 9be1cdfe1..b8dc7a4dd 100644 --- a/fealpy/material/elastic_material.py +++ b/fealpy/material/elastic_material.py @@ -142,6 +142,11 @@ def shear_modulus(self) -> float: """获取剪切模量 μ""" return self.mu + @property + def bulk_modulus(self) -> float: + """获取体积模量 K""" + return self.calculate_bulk_modulus() + @property def hypothesis(self) -> str: """获取平面假设""" @@ -279,42 +284,6 @@ def _shear_strain(self, gphi: TensorLike, return out - # def shear_strain(self, gphi: TensorLike, - # indices: TensorLike, *, - # out: Optional[TensorLike]=None) -> TensorLike: - # """Assembly shear strain tensor. - - # Parameters: - # gphi (TensorLike): Gradient of the scalar basis functions shaped (..., ldof, GD).\n - # indices (bool, optional): Indices of DoF components in the flattened DoF, shaped (ldof, GD).\n - # out (TensorLike | None, optional): Output tensor. Defaults to None. - - # Returns: - # TensorLike: Sheared strain shaped (..., NNZ, GD*ldof) where NNZ = (GD + (GD+1))//2. - # """ - # kwargs = bm.context(gphi) - # ldof, GD = gphi.shape[-2:] - # if GD < 2: - # raise ValueError(f"The shear strain requires GD >= 2, but GD = {GD}") - # NNZ = (GD * (GD-1))//2 # 剪切应变分量的数量 - # new_shape = gphi.shape[:-2] + (NNZ, GD*ldof) # (NC, NQ, NNZ, GD*ldof) - - # if out is None: - # out = bm.zeros(new_shape, **kwargs) - # else: - # if out.shape != new_shape: - # raise ValueError(f'out.shape={out.shape} != {new_shape}') - - # cursor = 0 - - # for i in range(0, GD-1): - # for j in range(i+1, GD): - # out = bm.set_at(out, (..., cursor, indices[:, i]), gphi[..., :, j]) - # out = bm.set_at(out, (..., cursor, indices[:, j]), gphi[..., :, i]) - # cursor += 1 - - # return out -