Skip to content

Commit

Permalink
Merge pull request #1407 from brighthe/develop
Browse files Browse the repository at this point in the history 
精益项目
  • Loading branch information
@AlbertZyy
AlbertZyy authored Dec 16, 2024
2 parents 4e96fc1 + 2177908 commit d1a073b
Show file tree
Hide file tree
Showing 3 changed files with 171 additions and 83 deletions.
56 changes: 41 additions & 15 deletions app/soptx/linear_elasticity/jy_linear_exact_api.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -443,22 +443,28 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
cell = mesh.entity('cell')

p = 1
q = p+2
q = p+3
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
# tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority
tldof = tensor_space.number_of_local_dofs()
tgdof = tensor_space.number_of_global_dofs()
print(f"tgdof: {tgdof}")
pde = BoxDomainLinear3d()

filename = "/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/box_linear_exact_fealpy_STIF2.mtx"
matrices = read_mtx_file(filename)
# matrices1 = matrices + bm.transpose(matrices, axes=(0, 2, 1))
KE0_true = matrices[0]
KE0_true_sum = bm.sum(bm.abs(KE0_true[0, :]))
KE0_true = matrices[0].round(4)
print(f"KE0_true_max: {bm.max(KE0_true)}")
abs_value0 = bm.abs(KE0_true)
print(f"KE0_abs_min: {bm.min(abs_value0[abs_value0 > 0])}")
print(f"KE0_true_min: {bm.min(KE0_true)}")
KE_true = matrices.round(4)
print(f"KE_true_max: {bm.max(KE_true)}")
abs_value = bm.abs(KE_true)
print(f"KE_abs_min: {bm.min(abs_value[abs_value > 0])}")
print(f"KE_true_min: {bm.min(KE_true)}")

# 刚度矩阵
E = 2.1e5
Expand All @@ -470,18 +476,36 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
linear_elastic_material = LinearElasticMaterial(name='E_nu',
elastic_modulus=E, poisson_ratio=nu,
hypo='3D', device=bm.get_device(mesh))
integrator_K_bbar = LinearElasticIntegrator(material=linear_elastic_material,
q=q, method='C3D8_BBar')
KE_bbar_yz_xz_xy = integrator_K_bbar.c3d8_bbar_assembly(space=tensor_space)
KE0_bbar_yz_xz_xy = KE_bbar_yz_xz_xy[0].round(4)
print(f"KE0_bbar_yz_xz_xy_max: {bm.max(KE0_bbar_yz_xz_xy)}")
print(f"KE0_bbar_yz_xz_xy_abs_min: {bm.min(bm.abs(KE0_bbar_yz_xz_xy))}")
print(f"KE0_bbar_yz_xz_xy_min: {bm.min(KE0_bbar_yz_xz_xy)}")
KE_bbar_yz_xz_xy = KE_bbar_yz_xz_xy.round(4)
print(f"KE_bbar_yz_xz_xy_max: {bm.max(KE_bbar_yz_xz_xy)}")
print(f"KE_bbar_yz_xz_xy_abs_min: {bm.min(bm.abs(KE_bbar_yz_xz_xy))}")
print(f"KE_bbar_yz_xz_xy_min: {bm.min(KE_bbar_yz_xz_xy)}")


integrator_K_sri = LinearElasticIntegrator(material=linear_elastic_material,
q=q, method='C3D8_SRI')
KE_sri_yz_xz_xy = integrator_K_sri.c3d8_sri_assembly(space=tensor_space)
KE0_sri_yz_xz_xy = KE_sri_yz_xz_xy[0].round(4)
print(f"KE0_sri_yz_xz_xy_max: {bm.max(KE0_sri_yz_xz_xy)}")
print(f"KE0_sri_yz_xz_xy_min: {bm.min(KE0_sri_yz_xz_xy)}")

integrator_K = LinearElasticIntegrator(material=linear_elastic_material,
q=q, method='voigt_bbar')
KE = integrator_K.voigt_bbar_assembly(space=tensor_space)
KE0 = KE[0]
KE0_sum = bm.sum(bm.abs(KE0[0, :]))
integrator_K1 = LinearElasticIntegrator(material=linear_elastic_material,
q=q, method='voigt')
KE1 = integrator_K.voigt_assembly(space=tensor_space)
KE10 = KE1[0]
KE10_sum = bm.sum(bm.abs(KE10[0, :]))
q=q, method='voigt_tensor')
KE_yz_xz_xy = integrator_K.voigt_tensor_assembly(space=tensor_space)
KE0_yz_xz_xy = KE_yz_xz_xy[0].round(4)
print(f"KE0_yz_xz_xy_max: {bm.max(KE0_yz_xz_xy)}")
print(f"KE0_yz_xz_xy_min: {bm.min(KE0_yz_xz_xy)}")


bform = BilinearForm(tensor_space)
bform.add_integrator(integrator_K)
bform.add_integrator(integrator_K_sri)
K = bform.assembly(format='csr')
Kdense = K.to_dense()
values = K.values()
Expand Down Expand Up @@ -536,6 +560,8 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
print(f"uh: {uh.shape}:\n {uh[:]}")
u_exact = tensor_space.interpolate(pde.solution)
print(f"u_exact: {u_exact[:].shape}:\n {u_exact[:]}")
error = bm.sum(bm.abs(uh - u_exact))
print(f"error: {error:.6f}")

if tensor_space.dof_priority:
uh_show = uh.reshape(GD, NN).T
Expand Down
133 changes: 92 additions & 41 deletions fealpy/fem/linear_elastic_integrator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,18 +55,29 @@ def fetch_assembly(self, space: _FS):
is_constant_J = True
else:
J = mesh.jacobi_matrix(bcs)
is_constant_J = bm.allclose(
J[:, 0:1, :, :], J[:, 1:, :, :],
rtol=1e-10, atol=1e-12
)

if is_constant_J:
detJ = None
else:
detJ = bm.linalg.det(J)
detJ = bm.linalg.det(J)

return cm, bcs, ws, gphi, detJ, is_constant_J

@enable_cache
def fetch_tensor_assembly(self, space: _FS):
index = self.index
mesh = getattr(space, 'mesh', None)
if not isinstance(mesh, HomogeneousMesh):
raise RuntimeError("The LinearElasticIntegrator only support spaces on"
f"homogeneous meshes, but {type(mesh).__name__} is"
"not a subclass of HomoMesh.")

cm = mesh.entity_measure('cell', index=index)
q = space.p+3 if self.q is None else self.q
qf = mesh.quadrature_formula(q)
bcs, ws = qf.get_quadrature_points_and_weights()
gphi = space.grad_basis(bcs, index=index, variable='x')
J = mesh.jacobi_matrix(bcs)
detJ = bm.linalg.det(J)

return cm, bcs, ws, gphi, detJ

@enable_cache
def fetch_fast_assembly(self, space: _FS):
index = self.index
Expand Down Expand Up @@ -112,6 +123,31 @@ def fetch_symbolic_assembly(self, space: _TS) -> TensorLike:
M = bm.tensor(symbolic_int.gphi_gphi_matrix())

return cm, mesh, glambda_x, bm.asarray(M, dtype=bm.float64)

@enable_cache
def fetch_c3d8_sri_assembly(self, space: _FS):
index = self.index
mesh = getattr(space, 'mesh', None)
if not isinstance(mesh, HomogeneousMesh):
raise RuntimeError("The LinearElasticIntegrator only support spaces on"
f"homogeneous meshes, but {type(mesh).__name__} is"
"not a subclass of HomoMesh.")

q1 = 1
qf1 = mesh.quadrature_formula(q1)
bcs1, ws1 = qf1.get_quadrature_points_and_weights()
gphi1 = 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 = space.grad_basis(bcs2, index=index, variable='x')
J2 = mesh.jacobi_matrix(bcs2)
detJ2 = bm.linalg.det(J2)

return bcs1, ws1, gphi1, detJ1, bcs2, ws2, gphi2, detJ2

def assembly(self, space: _TS) -> TensorLike:
scalar_space = space.scalar_space
Expand Down Expand Up @@ -236,41 +272,15 @@ def assembly(self, space: _TS) -> TensorLike:

return KK

@assemblymethod('voigt')
def voigt_assembly(self, space: _TS) -> TensorLike:
@assemblymethod('voigt_tensor')
def voigt_tensor_assembly(self, space: _TS) -> TensorLike:
scalar_space = space.scalar_space
cm, bcs, ws, gphi, detJ, is_constant_J = self.fetch_assembly(scalar_space)
_, bcs, ws, gphi, detJ = self.fetch_tensor_assembly(scalar_space)

D = self.material.elastic_matrix(bcs)
B = self.material.strain_matrix(dof_priority=space.dof_priority, gphi=gphi)

if is_constant_J:
# 常数 Jacobian：使用单元测度
KK = bm.einsum('q, c, cqki, cqkl, cqlj -> cij',
ws, cm, B, D, B)
else:
# 非常数 Jacobian：使用行列式
KK = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
ws, detJ, B, D, B)

return KK

@assemblymethod('voigt_bbar')
def voigt_bbar_assembly(self, space: _TS) -> TensorLike:
scalar_space = space.scalar_space
cm, bcs, ws, gphi, detJ, is_constant_J = self.fetch_assembly(scalar_space)

D = self.material.elastic_matrix(bcs)
B = self.material.strain_matrix(dof_priority=space.dof_priority, gphi=gphi,
bbar=True, cm=cm, ws=ws, detJ=detJ)

if is_constant_J:
# 常数 Jacobian：使用单元测度
KK = bm.einsum('q, c, cqki, cqkl, cqlj -> cij',
ws, cm, B, D, B)
else:
# 非常数 Jacobian：使用行列式
KK = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
B = self.material.strain_matrix(dof_priority=space.dof_priority,
gphi=gphi)
KK = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
ws, detJ, B, D, B)

return KK
Expand Down Expand Up @@ -596,3 +606,44 @@ def fast_assembly(self, space: _TS) -> TensorLike:
# KK[:, 2:KK.shape[1]:GD, 1:KK.shape[2]:GD] = D01 * A_zy + D55 * A_yz

return KK

@assemblymethod('C3D8_BBar')
def c3d8_bbar_assembly(self, space: _TS) -> TensorLike:
scalar_space = space.scalar_space
cm, bcs, ws, gphi, detJ = self.fetch_tensor_assembly(scalar_space)

D = self.material.elastic_matrix(bcs)
B = self.material.strain_matrix(dof_priority=space.dof_priority,
gphi=gphi,
bbar=True, cm=cm, ws=ws, detJ=detJ)

KK = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
ws, detJ, B, D, B)

return KK

@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)

D_q1 = self.material.elastic_matrix(bcs1)
D0_q1 = D_q1[..., :3, :3] # (1, 1, 3, 3)
B_q1 = self.material.strain_matrix(dof_priority=space.dof_priority,
gphi=gphi1,
bbar=None)
B0_q1 = B_q1[..., :3, :] # (NC, NQ, 3, TLDOF)

D_q2 = self.material.elastic_matrix(bcs2)
B_q2 = self.material.strain_matrix(dof_priority=space.dof_priority,
gphi=gphi2,
bbar=None)

KK = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
ws2, detJ2, B_q2, D_q2, B_q2)
KK0 = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
ws1, detJ1, B0_q1, D0_q1, B0_q1)
KK += KK0

return KK
65 changes: 38 additions & 27 deletions fealpy/material/elastic_material.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,34 +171,45 @@ def elastic_matrix(self, bcs: Optional[TensorLike] = None) -> TensorLike:

def strain_matrix(self, dof_priority: bool,
gphi: TensorLike,
shear_order: List[str]=['xy', 'zx', 'yz'],
shear_order: List[str]=['yz', 'xz', 'xy'],
# shear_order: List[str]=['xy', 'yz', 'xz'],
bbar: bool = False,
cm: TensorLike = None, ws: TensorLike = None, detJ: TensorLike = None) -> TensorLike:
'''
Constructs the strain-displacement matrix B for the material based on the gradient of the shape functions.
Constructs the strain-displacement matrix B for the material \n
based on the gradient of the shape functions.
B = [∂Ni/∂x 0 0 ]
[0 ∂Ni/∂y 0 ]
[0 0 ∂Ni/∂z]
[0 ∂Ni/∂z ∂Ni/∂y]
[∂Ni/∂z 0 ∂Ni/∂x]
[∂Ni/∂y ∂Ni/∂x 0 ]
B = [∂Ni/∂x 0 0 ]
[0 ∂Ni/∂y 0 ]
[0 0 ∂Ni/∂z]
[∂Ni/∂y ∂Ni/∂x 0 ]
[∂Ni/∂z 0 ∂Ni/∂x]
[0 ∂Ni/∂z ∂Ni/∂y]
B = [∂Ni/∂x 0 0 ]
[0 ∂Ni/∂y 0 ]
[0 0 ∂Ni/∂z]
[∂Ni/∂y ∂Ni/∂x 0 ]
[0 ∂Ni/∂z ∂Ni/∂y]
[∂Ni/∂z 0 ∂Ni/∂x]
[0 ∂Ni/∂z ∂Ni/∂y]
Parameters:
dof_priority (bool): A flag that determines the ordering of DOFs.
If True, the priority is given to the first dimension of degrees of freedom.
gphi (TensorLike): A tensor representing the gradient of the shape functions. Its shape
typically includes the number of local degrees of freedom and the geometric
dimension (GD).
gphi - (NC, NQ, LDOF, GD).
shear_order (List[str], optional): Specifies the order of shear strain components for GD=3.
Valid options are permutations of {'xy', 'yz', 'zx'}.
Defaults to ['xy', 'zx', 'yz']
Valid options are permutations of {'xy', 'yz', 'xz'}.
Returns:
TensorLike: The strain-displacement matrix `B`, which is a tensor with shape:
- For 2D problems (GD=2): (NC, NQ, 3, tldof)
- For 3D problems (GD=3): (NC, NQ, 6, tldof)
Here, NC is the number of cells, NQ is the number of quadrature points,
and tldof is the number of local degrees of freedom.
- For 2D problems (GD=2): (NC, NQ, 3, TLDOF)
- For 3D problems (GD=3): (NC, NQ, 6, TLDOF)
'''
ldof, GD = gphi.shape[-2:]
if dof_priority:
Expand All @@ -222,16 +233,16 @@ def _normal_strain(self, gphi: TensorLike,
"""Assembly normal strain tensor.
Parameters:
gphi - (NC, NQ, ldof, GD).\n
indices - (ldof, GD): Indices of DoF components in the flattened DoF, shaped .\n
gphi - (NC, NQ, LDOF, GD).
indices - (LDOF, GD): Indices of DoF components in the flattened DoF, shaped .
out - (TensorLike | None, optional): Output tensor. Defaults to None.
Returns:
out - Normal strain shaped (NC, NQ, GD, GD*ldof):
out - Normal strain shaped (NC, NQ, GD, GD*LDOF):
"""
kwargs = bm.context(gphi)
ldof, GD = gphi.shape[-2:]
new_shape = gphi.shape[:-2] + (GD, GD*ldof) # (NC, NQ, GD, GD*ldof)
new_shape = gphi.shape[:-2] + (GD, GD*ldof) # (NC, NQ, GD, GD*LDOF)

if out is None:
out = bm.zeros(new_shape, **kwargs)
Expand Down Expand Up @@ -269,16 +280,16 @@ def _normal_strain_bbar(self, gphi: TensorLike,
if out.shape != new_shape:
raise ValueError(f'out.shape={out.shape} != {new_shape}')

average_gphi = bm.einsum('cqid, c, q -> cid', gphi, cm, ws) # (NC, LDOF, GD)
# average_gphi = bm.einsum('cqid, cq, q -> cid', gphi, detJ, ws) # (NC, LDOF, GD)

average_gphi = bm.einsum('cqid, cq, q -> cid', gphi, detJ, ws) # (NC, LDOF, GD)
for i in range(GD):
for j in range(GD):
if i == j:
corrected_phi = (2.0 / 3.0) * gphi[..., :, i] \
+ (1.0 / (3.0 * cm[:, None, None]) ) * average_gphi[..., None, :, i] # (NC, NQ, LDOF)
+ (1.0 / (3.0 * cm[:, None, None]) ) * average_gphi[..., None, :, i] # (NC, NQ, LDOF)
else:
corrected_phi = (-1.0 / 3.0) * gphi[..., :, j] \
+ (1.0 / (3.0 * cm[:, None, None]) ) * average_gphi[..., None, :, j] # (NC, NQ, LDOF)
+ (1.0 / (3.0 * cm[:, None, None]) ) * average_gphi[..., None, :, j] # (NC, NQ, LDOF)

out = bm.set_at(out, (..., i, indices[:, j]), corrected_phi)

Expand All @@ -292,31 +303,31 @@ def _shear_strain(self, gphi: TensorLike,
"""Assembly shear strain tensor.
Parameters:
gphi - (NC, NQ, ldof, GD).\n
indices (bool, optional): Indices of DoF components in the flattened DoF, shaped (ldof, GD).\n
gphi - (NC, NQ, 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:
out - Shear strain shaped (NC, NQ, NNZ, GD*ldof) where NNZ = (GD + (GD+1))//2: .
out - Shear strain shaped (NC, NQ, 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)
NNZ = (GD * (GD-1))//2 # 剪切应变分量的数量
new_shape = gphi.shape[:-2] + (NNZ, GD*ldof) # (NC, NQ, NNZ, GD*LDOF)

if GD == 2:
shear_indices = [(0, 1)] # Corresponds to 'xy'
elif GD == 3:
valid_pairs = {'xy', 'zx', 'yz'}
valid_pairs = {'xy', 'yz', 'xz'}
if not set(shear_order).issubset(valid_pairs):
raise ValueError(f"Invalid shear_order: {shear_order}. Valid options are {valid_pairs}")

index_map = {
'xy': (0, 1),
'yz': (1, 2),
'zx': (2, 0),
'xz': (2, 0),
}
shear_indices = [index_map[pair] for pair in shear_order]
else:
Expand Down

0 comments on commit d1a073b

Please sign in to comment.