update

app/soptx/linear_elasticity/jy_linear_exact_api.py

@@ -509,12 +509,12 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
 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_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')
@@ -525,7 +525,7 @@ def dirichlet(self, points: TensorLike) -> TensorLike:
 
 
 bform = BilinearForm(tensor_space)
-bform.add_integrator(integrator_K_bbar)
+bform.add_integrator(integrator_K_sri)
 K = bform.assembly(format='csr')
 Kdense = K.to_dense()
 values = K.values()

fealpy/fem/linear_elastic_integrator.py

@@ -127,7 +127,8 @@ def fetch_c3d8_bbar_assembly(self, space: _FS):
         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)
+                                        correction='BBar', 
+                                        cm=cm, ws=ws, detJ=detJ)
 
         return ws, detJ, D, B
 
@@ -652,13 +653,13 @@ def c3d8_sri_assembly(self, space: _TS) -> TensorLike:
         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)
+                                        correction='SRI')
         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)
+                                        correction=None)
 
         KK = bm.einsum('q, cq, cqki, cqkl, cqlj -> cij',
                         ws2, detJ2, B_q2, D_q2, B_q2)

fealpy/material/elastic_material.py

@@ -173,7 +173,7 @@ def strain_matrix(self, dof_priority: bool,
                     gphi: TensorLike, 
                     shear_order: List[str]=['yz', 'xz', 'xy'],
                     # shear_order: List[str]=['xy', 'yz', 'xz'],
-                    bbar: bool = False, 
+                    correction: Optional[str] = None,  # 'None', 'BBar' 或 'SRI'
                     cm: TensorLike = None, ws: TensorLike = None, detJ: TensorLike = None) -> TensorLike:
         '''
         Constructs the strain-displacement matrix B for the material \n
@@ -216,8 +216,12 @@ def strain_matrix(self, dof_priority: bool,
             indices = flatten_indices((ldof, GD), (1, 0))
         else:
             indices = flatten_indices((ldof, GD), (0, 1))
-        if bbar:
+        if correction == 'BBar':
+            if any(param is None for param in (cm, ws, detJ)):  
+                raise ValueError("BBar correction requires cm, ws, and detJ parameters")
             normal_B = self._normal_strain_bbar(gphi, cm, ws, detJ, indices)
+        elif correction == 'SRI':
+            normal_B = self._normal_strain_sri(gphi, indices)
         else:
             normal_B = self._normal_strain(gphi, indices)
 
@@ -255,20 +259,49 @@ def _normal_strain(self, gphi: TensorLike,
 
         return out
 
+    def _normal_strain_sri(self, gphi: TensorLike,
+                    indices: TensorLike, *,
+                    out: Optional[TensorLike]=None) -> TensorLike:
+        """Assembly normal strain tensor with SRI correction.
+
+        Parameters:
+            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: SRI corrected normal strain shaped (NC, NQ, GD, GD*LDOF).
+        """
+        kwargs = bm.context(gphi)
+        ldof, GD = gphi.shape[-2:]
+        new_shape = gphi.shape[:-2] + (3, 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}')
+
+        # 构建修正的正应变矩阵，所有分量都使用完整积分
+        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,
                         cm, ws, detJ,
                         indices: TensorLike, *,
                         out: Optional[TensorLike] = None) -> TensorLike:
-        """
-        Assembly normal strain tensor with B-Bar correction.
+        """Assembly normal strain tensor with B-Bar correction.
 
         Parameters:
             gphi - (NC, NQ, LDOF, GD).
             indices (TensorLike): Indices of DoF components in the flattened DoF shaped (LDOF, GD).
             out (TensorLike | None, optional): Output tensor. Defaults to None.
 
         Returns:
-            TensorLike: Corrected normal strain shaped (NC, NQ, GD, GD*LDOF).
+            out: B-Bar corrected normal strain shaped (NC, NQ, GD, GD*LDOF).
         """
         kwargs = bm.context(gphi)
         ldof, GD = gphi.shape[-2:]