主题014：水平集方法（Level Set Method）

1. 引言

1.1 什么是水平集方法

水平集方法（Level Set Method）是由Osher和Sethian于1988年提出的一种用于追踪界面演化的数值技术。在结构优化领域，Allaire、Wang等人于2000年代初将其引入拓扑优化，开创了基于几何描述的优化新范式。

核心思想：
使用隐式函数 $\varphi (\mathbf{x})$ 描述结构边界：

• $\varphi (\mathbf{x})>0$：实体材料区域
• $\varphi (\mathbf{x})<0$：空腔区域
• $\varphi (\mathbf{x})=0$：结构边界（零水平集）
  

1.2 水平集 vs 传统方法

特性	水平集方法	SIMP	ESO/BESO
边界描述	隐式（光滑）	像素化	像素化
灰度单元	无	有	无
边界光滑性	高	低	中
拓扑变化	自动处理	需要过滤	离散变化
实现复杂度	高	低	中
计算成本	高	低	中

1.3 本章学习目标

通过本章学习，你将掌握：

1. 水平集方法的数学基础
2. Hamilton-Jacobi方程与边界演化
3. 形状灵敏度分析
4. 水平集的重新初始化技术
5. Python实现与工程应用

---

2. 理论基础

2.1 水平集函数

2.1.1 定义与性质

水平集函数 $\varphi (\mathbf{x},t)$ 是一个定义在设计域 $\Omega$ 上的标量场，满足：

Ωsolid={x:ϕ(x)>0}\Omega_{solid} = \{\mathbf{x} : \phi(\mathbf{x}) > 0\}Ωsolid​={x:ϕ(x)>0}
Ωvoid={x:ϕ(x)<0}\Omega_{void} = \{\mathbf{x} : \phi(\mathbf{x}) < 0\}Ωvoid​={x:ϕ(x)<0}
Γ={x:ϕ(x)=0}\Gamma = \{\mathbf{x} : \phi(\mathbf{x}) = 0\}Γ={x:ϕ(x)=0}

理想性质：

1. 符号距离函数：$|\nabla \varphi |=1$
2. 光滑性：$\varphi \in C^1$
3. 正则性：避免尖锐特征

2.1.2 初始化方法

带孔洞的初始设计：

# 创建多个圆形孔洞
phi = np.ones((nely + 1, nelx + 1))
for hole in holes:
    dist = np.sqrt((X - cx)**2 + (Y - cy)**2)
    phi = np.minimum(phi, dist - r)  # 孔洞内部phi < 0

符号距离函数：
$$\varphi (\mathbf{x})=\pm \underset{\mathbf{y}\in \Gamma }{\min }\Vert \mathbf{x}-\mathbf{y}\Vert$$

2.2 Hamilton-Jacobi方程

2.2.1 边界演化方程

结构边界随时间演化遵循Hamilton-Jacobi方程：
$$\frac{\partial \varphi }{\partial t}+V_n|\nabla \varphi |=0$$

其中：

• $V_n$：边界的法向速度
• $|\nabla \varphi |$：水平集梯度模长

2.2.2 速度场的确定

速度场由形状灵敏度决定：
$$V_n=-\frac{\partial J}{\partial \Omega }$$

对于柔度最小化问题：
$$V_n=-\mathbf{u}:\mathbf{\sigma }$$

其中：

• $\mathbf{u}$：位移场
• $\mathbf{\sigma }$：应力张量

2.2.3 数值求解

迎风格式（Upwind Scheme）：
根据速度方向选择差分方向：

$$\frac{\partial \varphi }{\partial x}\approx \{\begin{array}{ll}\frac{{\varphi }_{i,j}-{\varphi }_{i-1,j}}{\Delta x}& V_x>0\\ \frac{{\varphi }_{i+1,j}-{\varphi }_{i,j}}{\Delta x}& V_x<0\end{array}$$

时间步进：
$${\varphi }^{n+1}={\varphi }^n-\Delta t\cdot V_n|\nabla \varphi |$$

稳定性条件（CFL条件）：
$$\Delta t\le \frac{\Delta x}{\max |V_n|}$$

2.3 形状灵敏度分析

2.3.1 形状导数

目标函数 $J$ 对边界形状的导数：
$$\frac{dJ}{d\Omega }={\int }_{\Gamma }G(\mathbf{u},\mathbf{p})V_{n}d\Gamma$$

其中：

• $G$：形状梯度
• $\mathbf{p}$：伴随变量

2.3.2 柔度最小化

对于柔度 $C={\mathbf{F}}^T\mathbf{U}$：
$$G=-\mathbf{u}:\mathbf{\sigma }$$

物理意义：

• 高应变能区域：应保留材料（$V_n<0$，边界向内移动）
• 低应变能区域：可删除材料（$V_n>0$，边界向外移动）

2.3.3 体积约束处理

使用Lagrange乘子法：
$$\mathcal{L}=C+\lambda (V-V^{\ast })$$

速度场修正：
$$V_n^{eff}=V_n+\lambda$$

其中 $\lambda$ 通过二分法确定，满足体积约束。

2.4 重新初始化

2.4.1 为什么需要重新初始化

演化过程中水平集函数会偏离符号距离函数：

• 梯度模长 $|\nabla \varphi |$ 不再等于1
• 导致数值不稳定
• 影响边界精度

2.4.2 重新初始化方程

求解到稳态：
$$\frac{\partial \varphi }{\partial \tau }+S({\varphi }_0)(|\nabla \varphi |-1)=0$$

其中：

• $\tau$：伪时间
• $S({\varphi }_0)$：符号函数

2.4.3 简化方法

高斯平滑法：

phi_smooth = gaussian_filter(phi, sigma=1.0)
phi_new = np.sign(phi) * np.abs(phi_smooth)

快速行进法（Fast Marching Method）：
精确求解距离函数，计算复杂度 $O(N\log N)$。

---

3. Python实现

3.1 水平集优化器类

class LevelSetTopologyOptimizer:
    """
    水平集拓扑优化器
    
    基于Hamilton-Jacobi方程的边界演化
    """
    
    def __init__(self, nelx, nely, volfrac, tau=0.1):
        """
        初始化水平集优化器
        
        参数:
            nelx: x方向单元数
            nely: y方向单元数
            volfrac: 目标体积分数
            tau: 时间步长
        """
        self.nelx = nelx
        self.nely = nely
        self.volfrac = volfrac
        self.tau = tau
        self.nu = 0.3
        
        # 创建网格
        self.x = np.linspace(0, nelx, nelx + 1)
        self.y = np.linspace(0, nely, nely + 1)
        self.X, self.Y = np.meshgrid(self.x, self.y)
        
        # 初始化水平集函数
        self.phi = self._initialize_level_set()

3.2 水平集初始化

def _initialize_level_set(self):
    """初始化水平集函数"""
    # 创建带孔洞的初始设计
    phi = np.ones((self.nely + 1, self.nelx + 1))
    
    # 添加一些初始孔洞
    n_holes = 5
    np.random.seed(42)
    for _ in range(n_holes):
        cx = np.random.uniform(10, self.nelx - 10)
        cy = np.random.uniform(5, self.nely - 5)
        r = np.random.uniform(3, 6)
        
        dist = np.sqrt((self.X - cx)**2 + (self.Y - cy)**2)
        phi = np.minimum(phi, dist - r)
    
    # 确保边界有材料
    phi[0, :] = 1.0
    phi[-1, :] = 1.0
    phi[:, 0] = 1.0
    phi[:, -1] = 1.0
    
    return phi

3.3 密度转换

def _phi_to_density(self, phi):
    """将水平集函数转换为单元密度"""
    x_elem = np.zeros((self.nely, self.nelx))
    
    for elx in range(self.nelx):
        for ely in range(self.nely):
            # 单元四个节点的水平集值
            phi_nodes = [
                phi[ely, elx],
                phi[ely+1, elx],
                phi[ely+1, elx+1],
                phi[ely, elx+1]
            ]
            # 平均值作为单元密度
            phi_avg = np.mean(phi_nodes)
            # 平滑Heaviside函数
            epsilon = 0.5
            x_elem[ely, elx] = 1.0 / (1.0 + np.exp(-phi_avg / epsilon))
    
    return x_elem

3.4 水平集演化

def _evolve_level_set(self, phi, velocity, dt):
    """
    演化水平集函数
    
    使用简单的迎风差分求解 Hamilton-Jacobi 方程
    """
    dx = 1.0
    dy = 1.0
    
    phi_new = phi.copy()
    
    # 内部节点更新
    for i in range(1, self.nelx):
        for j in range(1, self.nely):
            # 计算梯度（迎风格式）
            if velocity[j, i] > 0:
                phi_x = (phi[j, i] - phi[j, i-1]) / dx
            else:
                phi_x = (phi[j, i+1] - phi[j, i]) / dx
            
            if velocity[j, i] > 0:
                phi_y = (phi[j, i] - phi[j-1, i]) / dy
            else:
                phi_y = (phi[j+1, i] - phi[j, i]) / dy
            
            # 更新水平集
            phi_new[j, i] = phi[j, i] - dt * velocity[j, i] * np.sqrt(phi_x**2 + phi_y**2)
    
    return phi_new

3.5 重新初始化

def _reinitialize(self, phi, n_iter=10):
    """
    重新初始化水平集函数为符号距离函数
    """
    # 简化的重新初始化：使用高斯平滑
    phi_smooth = gaussian_filter(phi, sigma=1.0)
    
    # 保持符号
    phi_new = np.sign(phi) * np.abs(phi_smooth)
    
    return phi_new

---

4. 案例研究

4.1 悬臂梁优化

问题设置：

• 设计域：60×20悬臂梁
• 目标体积分数：40%
• 初始设计：带5个随机孔洞

优化过程：

• 每代演化水平集函数
• 每5代重新初始化
• 共30代收敛

结果分析：

• 最终柔度：约940
• 体积分数：81%（未达到目标，需调整参数）
• 边界光滑清晰

4.2 参数敏感性分析

时间步长 $\tau$：

• 过大：数值不稳定
• 过小：收敛缓慢
• 推荐：0.05-0.2

重新初始化频率：

• 每代：计算成本高
• 每5-10代：平衡选择
• 从不：数值不稳定

Heaviside参数 $ϵ$：

• 控制边界过渡区宽度
• 推荐：0.5-1.0

4.3 水平集 vs SIMP对比

指标	水平集	SIMP
边界质量	光滑	锯齿状
灰度单元	无	有
计算时间	长	短
实现难度	高	低
拓扑变化	自动	依赖过滤

选择建议：

• 工程应用：SIMP（成熟稳定）
• 学术研究：水平集（边界精确）
• 混合方法：SIMP+水平集（两阶段）

---

5. 进阶主题

5.1 多相水平集

使用多个水平集函数描述多材料：
$${\varphi }_1,{\varphi }_2,...,{\varphi }_n$$

材料区域：
Ωi={x:ϕi(x)>0 且 ϕj(x)<0,∀j≠i}\Omega_i = \{\mathbf{x} : \phi_i(\mathbf{x}) > 0 \text{ 且 } \phi_j(\mathbf{x}) < 0, \forall j \neq i\}Ωi​={x:ϕi​(x)>0 且 ϕj​(x)<0,∀j=i}

5.2 参数化水平集

使用基函数展开：
$$\varphi (\mathbf{x})=\sum _{i=1}^N{\alpha }_i{\psi }_i(\mathbf{x})$$

其中 ${\psi }_i$ 可以是：

• 径向基函数（RBF）
• 紧支径向基函数（CS-RBF）
• 小波基函数

优势：

• 降低设计变量维度
• 提高光滑性
• 便于制造约束处理

5.3 与等几何分析结合

等几何分析（IGA）：

• 使用NURBS基函数
• 精确几何描述
• 与水平集天然契合

优势：

• 消除几何误差
• 高阶连续性
• 适合薄壳结构

5.4 制造约束集成

最小尺寸约束：
通过速度场修正：
$$V_n^{eff}=V_n+V_n^{size}$$

悬空约束：
检测局部梯度方向，限制边界演化方向。

---

"""
水平集方法（Level Set Method）拓扑优化仿真
=============================================
基于隐式边界描述的拓扑优化方法

核心思想：
- 使用隐式函数 φ(x) 描述结构边界
- φ > 0: 实体材料
- φ < 0: 空腔
- φ = 0: 边界（零水平集）
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
from scipy.ndimage import gaussian_filter
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.size'] = 10


def lk():
    """计算4节点四边形单元的单元刚度矩阵"""
    E = 1.0
    nu = 0.3
    k = np.array([1/2-nu/6, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8,
                  -1/4+nu/12, -1/8-nu/8, nu/6, 1/8-3*nu/8])
    KE = E / (1-nu**2) * np.array([
        [k[0], k[1], k[2], k[3], k[4], k[5], k[6], k[7]],
        [k[1], k[0], k[7], k[6], k[5], k[4], k[3], k[2]],
        [k[2], k[7], k[0], k[5], k[6], k[3], k[4], k[1]],
        [k[3], k[6], k[5], k[0], k[7], k[2], k[1], k[4]],
        [k[4], k[5], k[6], k[7], k[0], k[1], k[2], k[3]],
        [k[5], k[4], k[3], k[2], k[1], k[0], k[7], k[6]],
        [k[6], k[3], k[4], k[1], k[2], k[7], k[0], k[5]],
        [k[7], k[2], k[1], k[4], k[3], k[6], k[5], k[0]]
    ])
    return KE


class LevelSetTopologyOptimizer:
    """
    水平集拓扑优化器
    
    基于Hamilton-Jacobi方程的边界演化
    """
    
    def __init__(self, nelx, nely, volfrac, tau=0.1):
        """
        初始化水平集优化器
        
        参数:
            nelx: x方向单元数
            nely: y方向单元数
            volfrac: 目标体积分数
            tau: 时间步长
        """
        self.nelx = nelx
        self.nely = nely
        self.volfrac = volfrac
        self.tau = tau
        self.nu = 0.3
        
        # 创建网格
        self.x = np.linspace(0, nelx, nelx + 1)
        self.y = np.linspace(0, nely, nely + 1)
        self.X, self.Y = np.meshgrid(self.x, self.y)
        
        # 初始化水平集函数（带孔洞的初始设计）
        self.phi = self._initialize_level_set()
        
        # 历史记录
        self.history = {
            'phi': [],
            'compliance': [],
            'volume': [],
            'boundary_length': []
        }
        
        self._prepare_fem()
    
    def _initialize_level_set(self):
        """初始化水平集函数"""
        # 创建带孔洞的初始设计
        phi = np.ones((self.nely + 1, self.nelx + 1))
        
        # 添加一些初始孔洞
        n_holes = 5
        np.random.seed(42)
        for _ in range(n_holes):
            cx = np.random.uniform(10, self.nelx - 10)
            cy = np.random.uniform(5, self.nely - 5)
            r = np.random.uniform(3, 6)
            
            dist = np.sqrt((self.X - cx)**2 + (self.Y - cy)**2)
            phi = np.minimum(phi, dist - r)
        
        # 确保边界有材料
        phi[0, :] = 1.0
        phi[-1, :] = 1.0
        phi[:, 0] = 1.0
        phi[:, -1] = 1.0
        
        return phi
    
    def _prepare_fem(self):
        """准备有限元分析"""
        self.ndof = 2 * (self.nelx + 1) * (self.nely + 1)
        self.KE = lk()
        
        # 构建自由度映射
        self.edofMat = np.zeros((self.nelx * self.nely, 8), dtype=int)
        for elx in range(self.nelx):
            for ely in range(self.nely):
                el = ely + elx * self.nely
                n1 = (self.nely + 1) * elx + ely
                n2 = (self.nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = [2*n1+2, 2*n1+3, 2*n2+2, 2*n2+3,
                                        2*n2, 2*n2+1, 2*n1, 2*n1+1]
        
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten()
        
        # 边界条件：左端固定
        fixeddofs = np.union1d(np.arange(0, 2*(self.nely+1), 2),
                                np.arange(1, 2*(self.nely+1), 2))
        alldofs = np.arange(self.ndof)
        self.freedofs = np.setdiff1d(alldofs, fixeddofs)
    
    def _phi_to_density(self, phi):
        """将水平集函数转换为单元密度"""
        # 使用节点水平集的平均值作为单元密度
        x_elem = np.zeros((self.nely, self.nelx))
        
        for elx in range(self.nelx):
            for ely in range(self.nely):
                # 单元四个节点的水平集值
                phi_nodes = [
                    phi[ely, elx],
                    phi[ely+1, elx],
                    phi[ely+1, elx+1],
                    phi[ely, elx+1]
                ]
                # 平均值作为单元密度（使用平滑过渡）
                phi_avg = np.mean(phi_nodes)
                # 平滑Heaviside函数
                epsilon = 0.5
                x_elem[ely, elx] = 1.0 / (1.0 + np.exp(-phi_avg / epsilon))
        
        return x_elem
    
    def _fem_analysis(self, x):
        """执行有限元分析"""
        E = 1e-9 + x * (1.0 - 1e-9)
        
        sK = ((self.KE.flatten()[np.newaxis]).T * E.flatten('F')).flatten(order='F')
        K = coo_matrix((sK, (self.iK, self.jK)),
                      shape=(self.ndof, self.ndof)).tocsc()
        K = (K + K.T) / 2
        
        F = np.zeros(self.ndof)
        F[-1] = -1.0
        
        U = np.zeros(self.ndof)
        try:
            U[self.freedofs] = spsolve(
                K[self.freedofs, :][:, self.freedofs],
                F[self.freedofs])
        except:
            U[self.freedofs] = 0
        
        return U, E
    
    def _compute_sensitivity(self, U, x):
        """计算形状灵敏度"""
        ce = np.zeros((self.nely, self.nelx))
        
        for el in range(self.nelx * self.nely):
            edofs = self.edofMat[el, :]
            Ue = U[edofs]
            ce.flat[el] = Ue @ self.KE @ Ue
        
        # 将单元灵敏度映射到节点
        sensitivity = np.zeros((self.nely + 1, self.nelx + 1))
        
        for elx in range(self.nelx):
            for ely in range(self.nely):
                val = ce[ely, elx]
                # 分配到四个节点
                sensitivity[ely, elx] += val / 4
                sensitivity[ely+1, elx] += val / 4
                sensitivity[ely+1, elx+1] += val / 4
                sensitivity[ely, elx+1] += val / 4
        
        return sensitivity
    
    def _compute_curvature(self, phi):
        """计算水平集的曲率（用于正则化）"""
        # 使用中心差分计算梯度
        dx = 1.0
        dy = 1.0
        
        phi_x = (phi[:, 2:] - phi[:, :-2]) / (2 * dx)
        phi_y = (phi[2:, :] - phi[:-2, :]) / (2 * dy)
        
        phi_xx = (phi[:, 2:] - 2*phi[:, 1:-1] + phi[:, :-2]) / (dx**2)
        phi_yy = (phi[2:, :] - 2*phi[1:-1, :] + phi[:-2, :]) / (dy**2)
        
        phi_xy = (phi[2:, 2:] - phi[2:, :-2] - phi[:-2, 2:] + phi[:-2, :-2]) / (4 * dx * dy)
        
        # 计算曲率（在内部点）
        curvature = np.zeros_like(phi)
        phi_x_pad = np.zeros_like(phi)
        phi_y_pad = np.zeros_like(phi)
        phi_x_pad[:, 1:-1] = phi_x
        phi_y_pad[1:-1, :] = phi_y
        
        grad_mag = np.sqrt(phi_x_pad**2 + phi_y_pad**2 + 1e-10)
        
        # 简化的曲率计算
        curvature[1:-1, 1:-1] = (phi_xx + phi_yy) / 2.0
        
        return curvature
    
    def _evolve_level_set(self, phi, velocity, dt):
        """
        演化水平集函数
        
        使用简单的迎风差分求解 Hamilton-Jacobi 方程
        """
        dx = 1.0
        dy = 1.0
        
        phi_new = phi.copy()
        
        # 内部节点更新
        for i in range(1, self.nelx):
            for j in range(1, self.nely):
                # 计算梯度（迎风格式）
                if velocity[j, i] > 0:
                    phi_x = (phi[j, i] - phi[j, i-1]) / dx
                else:
                    phi_x = (phi[j, i+1] - phi[j, i]) / dx
                
                if velocity[j, i] > 0:
                    phi_y = (phi[j, i] - phi[j-1, i]) / dy
                else:
                    phi_y = (phi[j+1, i] - phi[j, i]) / dy
                
                # 更新水平集
                phi_new[j, i] = phi[j, i] - dt * velocity[j, i] * np.sqrt(phi_x**2 + phi_y**2)
        
        return phi_new
    
    def _reinitialize(self, phi, n_iter=10):
        """
        重新初始化水平集函数为符号距离函数
        
        保持零水平集不变，将φ恢复为距离函数
        """
        # 简化的重新初始化：使用高斯平滑
        phi_smooth = gaussian_filter(phi, sigma=1.0)
        
        # 保持符号
        phi_new = np.sign(phi) * np.abs(phi_smooth)
        
        return phi_new
    
    def optimize(self, max_iter=50):
        """执行水平集优化"""
        print("="*70)
        print("水平集方法拓扑优化")
        print("="*70)
        print(f"\n优化参数:")
        print(f"  网格: {self.nelx} x {self.nely}")
        print(f"  目标体积分数: {self.volfrac}")
        print(f"  时间步长: {self.tau}")
        
        phi = self.phi.copy()
        loop = 0
        
        print("\n开始优化...")
        print(f"{'迭代':>6} {'柔度':>12} {'体积分数':>12} {'边界长度':>12}")
        print("-"*50)
        
        while loop < max_iter:
            loop += 1
            
            # 转换密度场
            x = self._phi_to_density(phi)
            
            # 有限元分析
            U, E = self._fem_analysis(x)
            
            # 计算灵敏度
            sensitivity = self._compute_sensitivity(U, x)
            
            # 总柔度
            c = np.sum(E * x * sensitivity[:-1, :-1])
            
            # 当前体积
            vol = np.mean(x)
            
            # 计算速度场（基于灵敏度和体积约束）
            # Lagrange乘子处理体积约束
            lambda_vol = 0.0
            if vol > self.volfrac:
                lambda_vol = 0.1  # 惩罚超体积
            
            # 速度 = -灵敏度 + 体积惩罚
            velocity = -sensitivity + lambda_vol
            
            # 演化水平集
            phi = self._evolve_level_set(phi, velocity, self.tau)
            
            # 定期重新初始化
            if loop % 5 == 0:
                phi = self._reinitialize(phi)
            
            # 计算边界长度（近似）
            boundary = np.sum(np.abs(phi[1:-1, 1:-1]) < 0.5)
            
            # 记录历史
            self.history['phi'].append(phi.copy())
            self.history['compliance'].append(c)
            self.history['volume'].append(vol)
            self.history['boundary_length'].append(boundary)
            
            if loop % 5 == 0 or loop == 1:
                print(f"{loop:>6} {c:>12.4f} {vol:>12.4f} {boundary:>12}")
        
        print("-"*50)
        print(f"\n优化完成!")
        print(f"  总迭代: {loop}")
        print(f"  最终柔度: {c:.4f}")
        print(f"  最终体积分数: {vol:.4f}")
        
        self.phi_opt = phi
        return phi
    
    def visualize_result(self, filename='level_set_result.png'):
        """可视化优化结果"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        
        # 最终水平集函数
        ax = axes[0, 0]
        im = ax.contourf(self.X, self.Y, self.phi_opt, levels=20, cmap='RdBu')
        ax.contour(self.X, self.Y, self.phi_opt, levels=[0], colors='k', linewidths=2)
        ax.set_title('Level Set Function φ')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
        
        # 最终拓扑
        ax = axes[0, 1]
        x_final = self._phi_to_density(self.phi_opt)
        im = ax.imshow(x_final, cmap='gray_r', origin='lower',
                      extent=[0, self.nelx, 0, self.nely], vmin=0, vmax=1)
        ax.set_title('Final Topology')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax, label='Density')
        
        # 收敛曲线
        ax = axes[1, 0]
        ax.plot(self.history['compliance'], 'b-', linewidth=2)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Compliance')
        ax.set_title('Compliance Convergence')
        ax.grid(True)
        
        # 体积变化
        ax = axes[1, 1]
        ax.plot(self.history['volume'], 'r-', linewidth=2, label='Current')
        ax.axhline(y=self.volfrac, color='k', linestyle='--', label='Target')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Volume Fraction')
        ax.set_title('Volume Evolution')
        ax.legend()
        ax.grid(True)
        
        plt.tight_layout()
        plt.savefig(filename, dpi=150, bbox_inches='tight')
        print(f"\n结果已保存: {filename}")


class SimplifiedLevelSet:
    """
    简化版水平集方法（基于密度梯度）
    
    适用于教学演示的简化实现
    """
    
    def __init__(self, nelx, nely, volfrac):
        self.nelx = nelx
        self.nely = nely
        self.volfrac = volfrac
        
        # 初始化密度场（带孔洞）
        self.x = np.ones((nely, nelx))
        
        # 创建初始孔洞
        n_holes = 4
        np.random.seed(42)
        for _ in range(n_holes):
            cx = int(np.random.uniform(nelx * 0.2, nelx * 0.8))
            cy = int(np.random.uniform(nely * 0.2, nely * 0.8))
            r = int(np.random.uniform(2, 5))
            
            for i in range(max(0, cx-r), min(nelx, cx+r)):
                for j in range(max(0, cy-r), min(nely, cy+r)):
                    if (i-cx)**2 + (j-cy)**2 < r**2:
                        self.x[j, i] = 0.001
        
        self.history = {'compliance': [], 'volume': []}
        self._prepare_fem()
    
    def _prepare_fem(self):
        """准备有限元分析"""
        self.ndof = 2 * (self.nelx + 1) * (self.nely + 1)
        self.KE = lk()
        
        self.edofMat = np.zeros((self.nelx * self.nely, 8), dtype=int)
        for elx in range(self.nelx):
            for ely in range(self.nely):
                el = ely + elx * self.nely
                n1 = (self.nely + 1) * elx + ely
                n2 = (self.nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = [2*n1+2, 2*n1+3, 2*n2+2, 2*n2+3,
                                        2*n2, 2*n2+1, 2*n1, 2*n1+1]
        
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten()
        
        fixeddofs = np.union1d(np.arange(0, 2*(self.nely+1), 2),
                                np.arange(1, 2*(self.nely+1), 2))
        alldofs = np.arange(self.ndof)
        self.freedofs = np.setdiff1d(alldofs, fixeddofs)
    
    def _fem_analysis(self, x):
        """执行有限元分析"""
        E = 1e-9 + x * (1.0 - 1e-9)
        
        sK = ((self.KE.flatten()[np.newaxis]).T * E.flatten('F')).flatten(order='F')
        K = coo_matrix((sK, (self.iK, self.jK)),
                      shape=(self.ndof, self.ndof)).tocsc()
        K = (K + K.T) / 2
        
        F = np.zeros(self.ndof)
        F[-1] = -1.0
        
        U = np.zeros(self.ndof)
        try:
            U[self.freedofs] = spsolve(
                K[self.freedofs, :][:, self.freedofs],
                F[self.freedofs])
        except:
            U[self.freedofs] = 0
        
        return U, E
    
    def optimize(self, max_iter=50):
        """执行简化水平集优化"""
        print("="*70)
        print("简化水平集方法")
        print("="*70)
        
        x = self.x.copy()
        
        print(f"\n开始优化...")
        print(f"{'迭代':>6} {'柔度':>12} {'体积分数':>12}")
        print("-"*35)
        
        for loop in range(1, max_iter + 1):
            # 有限元分析
            U, E = self._fem_analysis(x)
            
            # 计算灵敏度
            ce = np.zeros((self.nely, self.nelx))
            for el in range(self.nelx * self.nely):
                edofs = self.edofMat[el, :]
                Ue = U[edofs]
                ce.flat[el] = Ue @ self.KE @ Ue
            
            c = np.sum(E * ce)
            vol = np.mean(x > 0.5)
            
            # 基于灵敏度的边界演化
            # 找到边界单元（密度在0.1-0.9之间）
            boundary = (x > 0.1) & (x < 0.9)
            solid = x > 0.9
            void = x < 0.1
            
            # 演化规则
            # 实体单元中灵敏度低的变为边界
            low_sens = solid & (ce < np.percentile(ce[solid], 20))
            x[low_sens] = 0.5
            
            # 边界单元根据灵敏度决定方向
            high_sens_boundary = boundary & (ce > np.median(ce))
            x[high_sens_boundary] = 1.0
            
            low_sens_boundary = boundary & (ce <= np.median(ce))
            x[low_sens_boundary] = 0.001
            
            # 体积控制
            if vol > self.volfrac * 1.05:
                # 删除更多
                solid_indices = np.where(solid)
                if len(solid_indices[0]) > 0:
                    n_remove = int(0.02 * len(solid_indices[0]))
                    remove_idx = np.argsort(ce[solid])[:n_remove]
                    x[solid_indices[0][remove_idx], solid_indices[1][remove_idx]] = 0.001
            
            self.history['compliance'].append(c)
            self.history['volume'].append(vol)
            
            if loop % 10 == 0:
                print(f"{loop:>6} {c:>12.4f} {vol:>12.4f}")
        
        self.x_opt = x
        return x


def compare_methods():
    """比较水平集方法与SIMP"""
    print("\n" + "="*70)
    print("水平集方法与SIMP对比")
    print("="*70)
    
    nelx, nely = 60, 20
    volfrac = 0.4
    
    # 水平集方法
    print("\n【方法1】水平集方法")
    lsm = LevelSetTopologyOptimizer(nelx, nely, volfrac, tau=0.1)
    phi = lsm.optimize(max_iter=30)
    lsm.visualize_result('level_set_result.png')
    
    print("\n" + "="*70)
    print("对比总结")
    print("="*70)
    print("\n水平集方法特点:")
    print("  - 边界光滑清晰，无灰度单元")
    print("  - 基于几何描述，物理意义明确")
    print("  - 需要定期重新初始化")
    print("  - 实现复杂度较高")
    print("\n与SIMP对比:")
    print("  - SIMP: 连续密度，易实现，有灰度")
    print("  - LSM: 隐式边界，边界清晰，计算复杂")


def main():
    """主函数"""
    print("="*70)
    print("水平集方法（Level Set Method）拓扑优化仿真")
    print("="*70)
    
    # 方法对比
    compare_methods()
    
    print("\n所有仿真完成!")
    print("="*70)


if __name__ == "__main__":
    main()