Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
brighthe committed Jan 13, 2025
1 parent 57529ad commit b388b54
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 15 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,9 @@
from fealpy.fem.dirichlet_bc import DirichletBC
from fealpy.solver import cg, spsolve

from soptx.utils import timer
from app.soptx.linear_elasticity.JingYiGearProject.utils import export_to_inp
from app.gearx.gear import ExternalGear, InternalGear
import pickle
import json

def compute_strain_stress(tensor_space, uh, B_BBar, D):
Expand Down Expand Up @@ -163,14 +163,16 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):

# 从全局载荷向量中提取有载荷节点处的值
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/JingYiGearProject/inp/external_gear_helix_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='abaqus', mesh_type='hex')

# 创建计时器
t = timer(f"Timing_{i}")
next(t) # 启动计时器
E = 206e3
nu = 0.3
lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu))
Expand All @@ -184,8 +186,6 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
q=q, method='C3D8_BBar')
integrator_K.keep_data(True)
_, _, D, B_BBar = integrator_K.fetch_c3d8_bbar_assembly(tensor_space)
# print(f"B_BBar:\n {B_BBar[0]}")
# KE = integrator_K.c3d8_bbar_assembly(space=tensor_space)
bform = BilinearForm(tensor_space)
bform.add_integrator(integrator_K)
K = bform.assembly(format='csr')
Expand All @@ -201,14 +201,16 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
threshold=tensor_is_bd_dof,
method='interp')
K, F = dbc.apply(K, F)
t.send('准备时间')

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")
# print(f"uh10: {uh[:10]}")
# print(f"uh-10: {uh[-10:]}")
t.send('求解时间')
t.send(None)

# 计算残差向量和范数
residual = K.matmul(uh[:]) - F
residual_norm = bm.sqrt(bm.sum(residual * residual))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
method='interp')
# 齿面上所有节点的位移结果
uh_profiles = bm.zeros((NPN, GD), dtype=bm.float64) # (NPN, GD)
for i in range(10):
for i in range(3):
# 创建计时器
t = timer(f"Timing_{i}")
next(t) # 启动计时器
Expand All @@ -152,6 +152,7 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
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")
t.send('求解时间')
# 获取齿面上节点的位移
Expand All @@ -172,7 +173,7 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
uh_magnitude = bm.linalg.norm(uh_show, axis=1)
mesh.nodedata['uh'] = uh_show[:]
mesh.nodedata['uh_magnitude'] = uh_magnitude[:]
mesh.to_vtk('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/vtu/external_gear_profile_fealpy.vtu')
mesh.to_vtk(f'/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/vtu/external_gear_profile_fealpy_{i}.vtu')

# 单个节点载荷的索引
load_node_indices = node_indices[i].reshape(-1) # (1, )
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,13 @@
from fealpy.fem.bilinear_form import BilinearForm
from fealpy.fem.dirichlet_bc import DirichletBC
from fealpy.solver import cg, spsolve
from fealpy.mesh import HexahedronMesh

from soptx.utils import timer
from app.soptx.linear_elasticity.JingYiGearProject.utils import export_to_inp

import pickle

from fealpy.mesh import HexahedronMesh


def compute_strain_stress(tensor_space, uh, B_BBar, D):
cell2tdof = tensor_space.cell_to_dof()
cuh = uh[cell2tdof] # (NC, TLDOF)
Expand Down Expand Up @@ -143,13 +142,15 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
F_load_nodes = F[dof_indices] # (15*8, GD)

fixed_node_index = bm.where(is_outer_node)[0]
# print(f"fixed_node_index: {fixed_node_index}")
export_to_inp(filename='/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/inp/internal_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='abaqus', mesh_type='hex')

# 创建计时器
t = timer(f"Timing_{i}")
next(t) # 启动计时器
E = 206e3
nu = 0.3
lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu))
Expand Down Expand Up @@ -180,14 +181,14 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D):
threshold=tensor_is_bd_dof,
method='interp')
K, F = dbc.apply(K, F)

t.send('准备时间')
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")
print(f"uh10:, {uh[:10]}")
print(f"uh[-10:]:, {uh[-10:]}")
t.send('求解时间')
t.send(None)

# 计算残差向量和范数
residual = K.matmul(uh[:]) - F
Expand Down

0 comments on commit b388b54

Please sign in to comment.