Merge pull request #1367 from brighthe/develop

modify app

app/soptx/soptx/solver/test.py

@@ -0,0 +1,130 @@
+import sympy as sp
+
+from fealpy.backend import backend_manager as bm
+
+from fealpy.functionspace import LagrangeFESpace
+from fealpy.fem import ScalarMassIntegrator
+from fealpy.mesh import TriangleMesh
+
+class SymbolicIntegration:
+    def __init__(self, space1, space2=None):
+        self.space1 = space1
+        self.mesh = space1.mesh
+        self.p1 = space1.p  # 第一个空间的多项式次数
+        self.GD = self.mesh.geo_dimension()  # 几何维度
+        self.ldof1 = space1.number_of_local_dofs()  # 第一个空间的局部自由度数量
+        self.mi1 = self.mesh.multi_index_matrix(p=self.p1, etype=self.GD)  # 第一个空间的多重指标矩阵
+
+        # 如果没有提供第二个空间，则假设两个空间相同
+        if space2 is None:
+            self.space2 = space1
+            self.p2 = self.p1
+            self.ldof2 = self.ldof1
+            self.mi2 = self.mi1
+        else:
+            self.space2 = space2
+            self.p2 = space2.p  # 第二个空间的多项式次数
+            self.ldof2 = space2.number_of_local_dofs()  # 第二个空间的局部自由度数量
+            self.mi2 = self.mesh.multi_index_matrix(p=self.p2, etype=self.GD)  # 第二个空间的多重指标矩阵
+
+        # 定义符号重心坐标 λ_i
+        self.l = sp.symbols('l0:%d' % (self.GD+1), real=True)
+
+    def basis(self, p, mi):
+        GD = self.GD
+        ldof = mi.shape[0]
+        l = self.l
+
+        # 初始化矩阵 A，大小为 (p+1) x (GD+1)
+        A = sp.ones(p+1, GD+1)
+
+        # 计算 A 矩阵的元素
+        for i in range(1, p+1):
+            for j in range(GD+1):
+                A[i, j] = (p*l[j] - (i-1))*A[i-1, j]
+            # 除以阶乘 i!
+            A[i, :] /= sp.factorial(i)
+
+        # 初始化基函数数组
+        phi = sp.ones(ldof)
+
+        # 计算基函数
+        for i in range(ldof):
+            phi_i = 1
+            for j in range(GD+1):
+                mi_ij = mi[i, j]
+                phi_i *= A[mi_ij, j]
+            phi[i] = phi_i
+
+        return phi
+
+    def multi_index(self, monomial):
+        l = self.l
+        GD = self.GD
+
+        m = monomial.as_powers_dict()
+        a = bm.zeros(GD+1, dtype=bm.int32) # 返回幂指标
+        for i in range(GD+1):
+            a[i] = int(m.get(l[i]) or 0)
+
+        return a
+
+    def integrate(self, f):
+        GD = self.GD
+
+        f = f.expand()
+        r = 0    # 积分值
+        for m in f.as_coeff_add()[1]:
+            c = m.as_coeff_mul()[0] # 返回系数
+            a = self.multi_index(m) # 返回单项式的幂指标
+            temp = 1
+            for i in range(GD+1):
+                temp *= sp.factorial(a[i])
+            r += sp.factorial(GD) * c * temp / sp.factorial(sum(a) + GD)
+
+        return r + f.as_coeff_add()[0]
+
+    def phi_phi_matrix(self):
+        phi1 = self.basis(self.p1, self.mi1)
+        phi2 = self.basis(self.p2, self.mi2)
+        ldof1 = self.ldof1
+        ldof2 = self.ldof2
+
+        # 初始化矩阵 M
+        M = sp.zeros(ldof1, ldof2)
+
+        # 计算单元面积
+        area = self.mesh.entity_measure('cell')[0]  # 获取第一个单元的面积
+
+        for i in range(ldof1):
+            for j in range(ldof2):
+                integrand = phi1[i] * phi2[j]
+                M[i, j] = self.integrate(integrand)
+                integral_value = self.integrate(integrand)
+                M[i, j] = integral_value * area  # 将面积乘到积分结果上
+
+        return M
+
+# 初始化网格和空间
+mesh = TriangleMesh.from_box(box=[0, 1, 0, 1], nx=2, ny=2, device='cpu')
+# 定义两个不同的多项式次数
+p1 = 1
+p2 = 1
+
+# 创建两个有限元空间
+space1 = LagrangeFESpace(mesh=mesh, p=p1, ctype='C')
+space2 = LagrangeFESpace(mesh=mesh, p=p2, ctype='C')
+
+# 创建符号积分类的实例，传入两个空间
+symbolic_int = SymbolicIntegration(space1, space2)
+
+# 计算质量矩阵
+M = symbolic_int.phi_phi_matrix()
+
+integrator = ScalarMassIntegrator(q=5)
+M1 = integrator.assembly(space=space1)
+print("M1:", M1[0])
+
+# 输出质量矩阵
+print("Mass Matrix M:")
+sp.pprint(M)

app/soptx/soptx/tests/test_assemble.py

@@ -0,0 +1,210 @@
+"""Test cases for compliance objective and volume constraint."""
+
+from dataclasses import dataclass
+import dataclasses
+from typing import Literal, Optional, Union, Dict, Any
+
+from fealpy.backend import backend_manager as bm
+from fealpy.mesh import UniformMesh2d, UniformMesh3d, TriangleMesh
+from fealpy.functionspace import LagrangeFESpace, TensorFunctionSpace
+
+from soptx.material import (
+    ElasticMaterialConfig,
+    ElasticMaterialProperties,
+    SIMPInterpolation
+)
+from soptx.pde import MBBBeam2dData1, Cantilever2dData1, Cantilever3dData1
+from soptx.solver import ElasticFEMSolver, AssemblyMethod, AssemblyConfig
+
+bm.set_backend('numpy')
+
+@dataclass
+class TestConfig:
+    """Configuration for topology optimization test cases."""
+    # Required parameters (no default values)
+    problem_type: Literal['mbb_2d', 'cantilever_2d', 'cantilever_3d']
+    nx: int
+    ny: int
+    volume_fraction: float
+    filter_radius: float
+    filter_type: Literal['sensitivity', 'density', 'heaviside']
+
+    # Optional parameters (with default values)
+    nz: Optional[int] = None  # Only for 3D problems
+    elastic_modulus: float = 1.0
+    poisson_ratio: float = 0.3
+    minimal_modulus: float = 1e-9
+    penalty_factor: float = 3.0
+    projection_beta: Optional[float] = None  # Only for Heaviside projection
+    move_limit: float = 0.2
+    tolerance: float = 0.01
+    initial_lambda: float = 1e9
+    bisection_tol: float = 1e-3
+
+    assembly_method: AssemblyMethod = AssemblyMethod.STANDARD  # 矩阵组装方法
+    quadrature_degree_increase: int = 3  # 积分阶数增量
+    solver_type: Literal['cg', 'direct'] = 'cg'  # 求解器类型
+    solver_params: Optional[Dict[str, Any]] = None  # 求解器参数
+
+def create_base_components(config: TestConfig):
+    """Create basic components needed for topology optimization based on configuration."""
+    # Create mesh and determine dimensionality
+    if config.problem_type == 'cantilever_3d':
+        extent = [0, config.nx, 0, config.ny, 0, config.nz]
+        h = [1.0, 1.0, 1.0]
+        origin = [0.0, 0.0, 0.0]
+        mesh = UniformMesh3d(
+            extent=extent, h=h, origin=origin,
+            ipoints_ordering='zyx', flip_direction='y',
+            device='cpu'
+        )
+        dimension = 3
+    else:
+        extent = [0, config.nx, 0, config.ny]
+        h = [1.0, 1.0]
+        origin = [0.0, 0.0]
+        # mesh = UniformMesh2d(
+        #     extent=extent, h=h, origin=origin,
+        #     ipoints_ordering='yx', flip_direction='y',
+        #     device='cpu'
+        # )
+        mesh = TriangleMesh.from_box(box=[0, config.nx*h[0], 0, config.ny*h[1]], 
+                                    nx=config.nx, ny=config.ny, device='cpu')
+        dimension = 2
+
+    # Create function spaces
+    p = 1
+    space_C = LagrangeFESpace(mesh=mesh, p=p, ctype='C')
+    tensor_space_C = TensorFunctionSpace(space_C, (-1, dimension))
+    space_D = LagrangeFESpace(mesh=mesh, p=p-1, ctype='D')
+
+    # Create material properties
+    material_config = ElasticMaterialConfig(
+        elastic_modulus=config.elastic_modulus,
+        minimal_modulus=config.minimal_modulus,
+        poisson_ratio=config.poisson_ratio,
+        plane_assumption="3D" if dimension == 3 else "plane_stress"
+    )
+    interpolation_model = SIMPInterpolation(penalty_factor=config.penalty_factor)
+    material_properties = ElasticMaterialProperties(
+        config=material_config,
+        interpolation_model=interpolation_model
+    )
+
+    # Create PDE problem based on problem type
+    if config.problem_type == 'mbb_2d':
+        pde = MBBBeam2dData1(
+            xmin=0, xmax=config.nx*h[0],
+            ymin=0, ymax=config.ny*h[1]
+        )
+    elif config.problem_type == 'cantilever_2d':
+        pde = Cantilever2dData1(
+            xmin=0, xmax=config.nx*h[0],
+            ymin=0, ymax=config.ny*h[1]
+        )
+    elif config.problem_type == 'cantilever_3d':  
+        pde = Cantilever3dData1(
+            xmin=0, xmax=config.nx*h[0],
+            ymin=0, ymax=config.ny*h[1],
+            zmin=0, zmax=config.nz*h[2]
+        )
+
+    assembly_config = AssemblyConfig(
+        method=config.assembly_method,
+        quadrature_degree_increase=config.quadrature_degree_increase
+    )
+
+    # 设置默认的求解器参数
+    default_solver_params = {
+        'cg': {'maxiter': 5000, 'atol': 1e-12, 'rtol': 1e-12},
+        'direct': {'solver_type': 'mumps'}
+    }
+    solver_params = config.solver_params or default_solver_params[config.solver_type]
+
+    solver = ElasticFEMSolver(
+       material_properties=material_properties,
+       tensor_space=tensor_space_C,
+       pde=pde,
+       assembly_config=assembly_config,
+       solver_type=config.solver_type,  # 添加默认求解器类型
+       solver_params=solver_params  # 添加求解器参数
+   )
+
+    # Initialize density field
+    array = config.volume_fraction * bm.ones(mesh.number_of_cells(), dtype=bm.float64)
+    rho = space_D.function(array)
+
+    return mesh, space_D, material_properties, solver, rho
+
+def test_solver(config: TestConfig):
+    """测试不同矩阵组装方法的 solver."""
+    print(f"\n=== Testing Solver with {config.assembly_method} ===")
+
+    # Create base components
+    mesh, space_D, material_properties, solver, rho = create_base_components(config)
+
+    # Test solver
+    solver.update_density(rho[:])
+    solver_result = solver.solve_cg()
+    displacement = solver_result.displacement
+    # print(f"\nSolver information:")
+    # print(f"- Displacement shape: {displacement.shape}:\n {displacement[:]}")
+
+    # base_local_K = solver.get_base_local_stiffness_matrix()
+    # print("\n=== 基础局部刚度矩阵信息 ===")
+    # print(f"基础局部刚度矩阵 - {base_local_K.shape}:\n {base_local_K[0]}")
+    # local_K = solver.compute_local_stiffness_matrix()
+    # print("\n=== 当前材料局部刚度矩阵 ===")
+    # print(f"局部刚度矩阵 - {local_K.shape}:\n {local_K.round(4)}")
+
+    # K = solver.get_global_stiffness_matrix()
+    # F = solver.get_global_force_vector()
+    # print("\n=== 全局矩阵和载荷向量信息 ===")
+    # print(f"全局刚度矩阵 - {K.shape}:\n {K.to_dense().round(4)}")
+    # print(f"全局刚度矩阵最大值: {bm.max(bm.abs(K.to_dense()))}")
+    # print(f"全局载荷向量 -  {F.shape}:\n {F[:]}")
+
+
+if __name__ == "__main__":
+    # Test 3D case with density filter
+    config_3d = TestConfig(
+        problem_type='cantilever_3d',
+        nx=60, ny=20, nz=4,
+        volume_fraction=0.3,
+        filter_radius=1.5,
+        filter_type='density'
+    )
+    # results_3d = test_compliance_objective(config_3d)
+
+    # Test 2D case with sensitivity filter
+    config_2d = TestConfig(
+        problem_type='mbb_2d',
+        nx=60, ny=20,
+        volume_fraction=0.5,
+        filter_radius=2.4,
+        filter_type='sensitivity'
+    )
+    # results_2d = test_compliance_objective(config_2d)
+
+    # Test 2D case with sensitivity filter
+    config_cantilever_2d = TestConfig(
+        problem_type='cantilever_2d',
+        nx=160, ny=100,
+        # nx=8, ny=5,
+        volume_fraction=0.4,
+        filter_radius=6,
+        filter_type='sensitivity'
+    )
+    # 测试标准组装方法
+    config_standard = dataclasses.replace(
+        config_cantilever_2d, 
+        assembly_method=AssemblyMethod.STANDARD
+    )
+    results_standard = test_solver(config_standard)
+
+    # 测试快速应力组装方法
+    config_fast_stress = dataclasses.replace(
+        config_cantilever_2d, 
+        assembly_method=AssemblyMethod.FAST_STRESS
+    )
+    results_fast_stress = test_solver(config_fast_stress)