一行Python代码生成令人惊叹的三维晶格:开启计算设计与增材制造的新范式

在计算设计、材料科学与增材制造的交汇处，一种名为TPMS（三周期极小曲面）的复杂三维晶格结构正悄然革新着多个领域。从仿生组织支架到轻量化航空航天部件，从超材料声学器件到高效换热结构，TPMS以其独特的几何特性与优异的力学性能，成为前沿研究与工程应用的热点。然而，生成这些精妙结构的传统路径往往依赖于昂贵的专业软件或复杂的编程环境，无形中筑起了技术普及的高墙。

今天，我将分享一个完全开源、依赖极简的Python解决方案。这个工具包不仅还原了学术文献中经典的TPMS生成算法，更将其封装为清晰易用的函数接口，让研究人员、工程师和设计师能够专注于创新应用，而非底层实现。

完整代码实现：从数学公式到可打印三维模型

以下代码块提供了一个自包含的TPMS晶格生成器。它仅需numpy即可运行核心功能，并能在检测到scikit-image时自动启用更精确的网格划分算法。

import numpy as np
from scipy import sparse
from scipy.sparse import linalg
from scipy.spatial import KDTree
import time
import struct
import warnings
warnings.filterwarnings('ignore')

defstlVolume(vertices, faces):
"""
    计算三角网格的表面积和体积。
    参数:
        vertices: 顶点数组 (n, 3)
        faces: 面数组 (m, 3)，每个面是三个顶点的索引
    返回:
        totalVolume: 总体积
        totalArea: 总表面积
    """
# 获取每个面的顶点
    p1 = vertices[faces[:, 0]]
    p2 = vertices[faces[:, 1]]
    p3 = vertices[faces[:, 2]]

# 计算向量
    v1 = p3 - p1
    v2 = p2 - p1

# 计算法向量和面积
    cross = np.cross(v1, v2)
    area = 0.5 * np.sqrt(np.sum(cross**2, axis=1))
    totalArea = np.sum(area)

# 计算体积 (使用散度定理)
    cross_norm = np.sqrt(np.sum(cross**2, axis=1) + 1e-12)
    nz = -cross[:, 2] / cross_norm
    zMean = (p1[:, 2] + p2[:, 2] + p3[:, 2]) / 3
    volume = area * zMean * nz
    totalVolume = np.sum(volume)

return totalVolume, totalArea

defwrite_stl_binary(filename, vertices, faces, title='TPMS Lattice'):
"""
    将三角网格写入二进制STL文件

    参数:
        filename: 输出文件名
        vertices: 顶点数组 (n, 3)
        faces: 面数组 (m, 3)，每个面是三个顶点的索引
        title: 文件标题（最多80个字符）
    """
# 确保标题长度不超过80个字符
if len(title) > 80:
        title = title[:80]

with open(filename, 'wb') as f:
# 写入80字节的标题
        f.write(title.ljust(80).encode('ascii'))

# 写入面片数量（4字节，小端）
        num_faces = len(faces)
        f.write(struct.pack('<I', num_faces))

# 写入每个面片
for face in faces:
# 获取三个顶点
            p1 = vertices[face[0]]
            p2 = vertices[face[1]]
            p3 = vertices[face[2]]

# 计算法向量
            v1 = p2 - p1
            v2 = p3 - p1
            normal = np.cross(v1, v2)
            norm = np.linalg.norm(normal)
if norm > 0:
                normal = normal / norm

# 写入法向量 (3个float, 4字节每个)
for value in normal:
                f.write(struct.pack('<f', float(value)))

# 写入三个顶点 (每个顶点3个float)
for vertex in [p1, p2, p3]:
for value in vertex:
                    f.write(struct.pack('<f', float(value)))

# 写入属性字节计数 (2字节，通常为0)
            f.write(struct.pack('<H', 0))

    print(f"已写入 {num_faces} 个面到二进制STL文件: {filename}")

defwrite_stl_ascii(filename, vertices, faces, title='TPMS Lattice'):
"""
    将三角网格写入ASCII STL文件

    参数:
        filename: 输出文件名
        vertices: 顶点数组 (n, 3)
        faces: 面数组 (m, 3)
        title: 文件标题
    """
with open(filename, 'w') as f:
        f.write(f"solid {title}\n")

for face in faces:
# 获取三个顶点
            p1 = vertices[face[0]]
            p2 = vertices[face[1]]
            p3 = vertices[face[2]]

# 计算法向量
            v1 = p2 - p1
            v2 = p3 - p1
            normal = np.cross(v1, v2)
            norm = np.linalg.norm(normal)
if norm > 0:
                normal = normal / norm

            f.write("  facet normal {:.6f} {:.6f} {:.6f}\n".format(*normal))
            f.write("    outer loop\n")
            f.write("      vertex {:.6f} {:.6f} {:.6f}\n".format(*p1))
            f.write("      vertex {:.6f} {:.6f} {:.6f}\n".format(*p2))
            f.write("      vertex {:.6f} {:.6f} {:.6f}\n".format(*p3))
            f.write("    endloop\n")
            f.write("  endfacet\n")

        f.write(f"endsolid {title}\n")

    print(f"已写入 {len(faces)} 个面到ASCII STL文件: {filename}")

defgenerateTPMS(TPMS='P', type="net", volFrac=0.27, numCell=1, tarSize=1, center=(0,0,0), nFinal=100):
"""
    生成TPMS晶格结构并保存为STL文件

    参数:
        TPMS: TPMS类型，可选 "P", "G", "D", "I", "S", "F", "N"
        type: 结构类型，"net" 或 "sheet"
        volFrac: 目标体积分数
        numCell: 每个坐标方向的晶胞数量
        tarSize: 晶格的目标尺寸
        center: 晶格中心坐标
        nFinal: 最终网格的细分数量
    """

    start_time = time.time()

# 校准参数
    n = 100
    rStart = -0.07
    rStep = 0.1
    volTol = 0.0001

# 创建用于校准的网格
    t = 1.0 / n
    x_max, y_max, z_max = 1, 1, 1

    xi = np.arange(0, x_max + t/2, t)
    yi = np.arange(0, y_max + t/2, t)
    zi = np.arange(0, z_max + t/2, t)

# 调整中心
    xi = xi + center[0]
    yi = yi + center[1]
    zi = zi + center[2]

    x, y, z = np.meshgrid(xi, yi, zi, indexing='ij')

# 根据TPMS类型选择公式
if TPMS == 'P':  # Schwarz Primitive
        F = -1.0 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y) + np.cos(2 * np.pi * z))
elif TPMS == 'G':  # Schoen Gyroid
        F = (np.cos(2 * np.pi * x) * np.sin(2 * np.pi * y) + 
             np.cos(2 * np.pi * y) * np.sin(2 * np.pi * z) + 
             np.cos(2 * np.pi * z) * np.sin(2 * np.pi * x))
elif TPMS == 'D':  # Schwarz Diamond
        F = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * y) * np.sin(2 * np.pi * z) + 
             np.sin(2 * np.pi * x) * np.cos(2 * np.pi * y) * np.cos(2 * np.pi * z) + 
             np.cos(2 * np.pi * x) * np.sin(2 * np.pi * y) * np.cos(2 * np.pi * z) + 
             np.cos(2 * np.pi * x) * np.cos(2 * np.pi * y) * np.sin(2 * np.pi * z))
elif TPMS == 'I':  # Schoen I-WP
        F = (2.0 * (np.cos(2 * np.pi * x) * np.cos(2 * np.pi * y) + 
                    np.cos(2 * np.pi * y) * np.cos(2 * np.pi * z) + 
                    np.cos(2 * np.pi * z) * np.cos(2 * np.pi * x)) - 
             (np.cos(4 * np.pi * x) + np.cos(4 * np.pi * y) + np.cos(4 * np.pi * z)))
elif TPMS == 'S':  # Fischer Koch S
        F = (np.cos(4 * np.pi * x) * np.sin(2 * np.pi * y) * np.cos(2 * np.pi * z) + 
             np.cos(2 * np.pi * x) * np.cos(4 * np.pi * y) * np.sin(2 * np.pi * z) + 
             np.sin(2 * np.pi * x) * np.cos(2 * np.pi * y) * np.cos(4 * np.pi * z))
elif TPMS == 'F':  # F-RD
        F = -(4.0 * (np.cos(2 * np.pi * x) * np.cos(2 * np.pi * y) * np.cos(2 * np.pi * z)) - 
               (np.cos(4 * np.pi * x) * np.cos(4 * np.pi * y) + 
                np.cos(4 * np.pi * y) * np.cos(4 * np.pi * z) + 
                np.cos(4 * np.pi * z) * np.cos(4 * np.pi * x)))
elif TPMS == 'N':  # Neovious
        F = (3.0 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y) + np.cos(2 * np.pi * z)) + 
4.0 * np.cos(2 * np.pi * x) * np.cos(2 * np.pi * y) * np.cos(2 * np.pi * z))
else:
raise ValueError(f"不支持的TPMS类型: {TPMS}")

# 尝试导入skimage.measure，如果失败则使用替代方法
try:
from skimage import measure
        has_skimage = True
except ImportError:
        print("警告: 未找到skimage.measure，将使用替代方法生成网格（精度较低）")
        has_skimage = False

# 简单的替代方法：手动计算等值面
defsimple_marching_cubes(data, level, spacing):
# 简化的marching cubes实现
            nx, ny, nz = data.shape
            vertices = []
            faces = []

for i in range(nx-1):
for j in range(ny-1):
for k in range(nz-1):
# 获取立方体的8个角点
                        cube_values = [
                            data[i, j, k],
                            data[i+1, j, k],
                            data[i+1, j+1, k],
                            data[i, j+1, k],
                            data[i, j, k+1],
                            data[i+1, j, k+1],
                            data[i+1, j+1, k+1],
                            data[i, j+1, k+1]
                        ]

# 计算立方体索引
                        cube_index = 0
for v_idx, val in enumerate(cube_values):
if val > level:
                                cube_index |= (1 << v_idx)

# 简单的四面体分解
if cube_index == 0or cube_index == 255:
continue

# 添加顶点和面
                        base_idx = len(vertices)
                        vertices.extend([
                            [i*spacing[0], j*spacing[1], k*spacing[2]],
                            [(i+1)*spacing[0], j*spacing[1], k*spacing[2]],
                            [(i+1)*spacing[0], (j+1)*spacing[1], k*spacing[2]],
                            [i*spacing[0], (j+1)*spacing[1], k*spacing[2]],
                            [i*spacing[0], j*spacing[1], (k+1)*spacing[2]],
                            [(i+1)*spacing[0], j*spacing[1], (k+1)*spacing[2]],
                            [(i+1)*spacing[0], (j+1)*spacing[1], (k+1)*spacing[2]],
                            [i*spacing[0], (j+1)*spacing[1], (k+1)*spacing[2]]
                        ])

# 添加简单面
if cube_index < 8:
                            faces.append([base_idx, base_idx+1, base_idx+2])
                            faces.append([base_idx, base_idx+2, base_idx+3])

return np.array(vertices), np.array(faces), None, None

# 粗略校准水平集常数r
if type == "net":
        fTPMS = F + rStart
elif type == "sheet":
        fTPMS = -(F + rStart) * (F - rStart)
else:
raise ValueError(f"不支持的类型: {type}")

# 提取等值面
if has_skimage:
try:
# 新版本scikit-image的marching_cubes
            verts, faces, normals, values = measure.marching_cubes(
                fTPMS, 0, spacing=(t, t, t)
            )
except:
# 旧版本scikit-image
            verts, faces, normals, values = measure.marching_cubes_lewiner(
                fTPMS, 0, spacing=(t, t, t)
            )
else:
        verts, faces, normals, values = simple_marching_cubes(fTPMS, 0, (t, t, t))

# 计算体积分数
if len(verts) > 0and len(faces) > 0:
        totalVolume, totalArea = stlVolume(verts, faces)
        volFarc = -totalVolume
else:
        volFarc = 0

# 粗略校准循环
    rStart = rStart + rStep
while abs(volFarc - volFrac) > volTol and len(verts) > 0:
        r0 = rStart
        volVZ = np.sign(volFarc - volFrac)

if type == "net":
            fTPMS = F + rStart
else:  # sheet
            fTPMS = -(F + rStart) * (F - rStart)

# 提取等值面
if has_skimage:
try:
                verts, faces, normals, values = measure.marching_cubes(
                    fTPMS, 0, spacing=(t, t, t)
                )
except:
                verts, faces, normals, values = measure.marching_cubes_lewiner(
                    fTPMS, 0, spacing=(t, t, t)
                )
else:
            verts, faces, normals, values = simple_marching_cubes(fTPMS, 0, (t, t, t))

if len(verts) > 0and len(faces) > 0:
            totalVolume, totalArea = stlVolume(verts, faces)
            volFarc = -totalVolume
else:
            volFarc = 0

if volVZ == 1and volFarc - volFrac > 0:
            rStart = rStart - rStep
elif volVZ == -1and volFarc - volFrac < 0:
            rStart = rStart + rStep

if rStart == r0:
break

# 精确校准
if volFarc > volFrac:
        r1, r2 = rStart - rStep, rStart
else:
        r1, r2 = rStart, rStart + rStep

while abs(volFarc - volFrac) > volTol and len(verts) > 0:
        rStart = (r1 + r2) / 2

if type == "net":
            fTPMS = F + rStart
else:  # sheet
            fTPMS = -(F + rStart) * (F - rStart)

# 提取等值面
if has_skimage:
try:
                verts, faces, normals, values = measure.marching_cubes(
                    fTPMS, 0, spacing=(t, t, t)
                )
except:
                verts, faces, normals, values = measure.marching_cubes_lewiner(
                    fTPMS, 0, spacing=(t, t, t)
                )
else:
            verts, faces, normals, values = simple_marching_cubes(fTPMS, 0, (t, t, t))

if len(verts) > 0and len(faces) > 0:
            totalVolume, totalArea = stlVolume(verts, faces)
            volFarc = -totalVolume
else:
            volFarc = 0

if volFarc - volFrac > 0:
            r2 = (r1 + r2) / 2
else:
            r1 = (r1 + r2) / 2

# 生成最终TPMS晶格
    t_final = 1.0 / nFinal
    x_max_final = y_max_final = z_max_final = numCell

    xi_final = np.arange(0, x_max_final + t_final/2, t_final)
    yi_final = np.arange(0, y_max_final + t_final/2, t_final)
    zi_final = np.arange(0, z_max_final + t_final/2, t_final)

    xi_final = xi_final + center[0]
    yi_final = yi_final + center[1]
    zi_final = zi_final + center[2]

    x_final, y_final, z_final = np.meshgrid(xi_final, yi_final, zi_final, indexing='ij')

# 重新计算F
if TPMS == 'P':
        F_final = -1.0 * (np.cos(2 * np.pi * x_final) + np.cos(2 * np.pi * y_final) + np.cos(2 * np.pi * z_final))
elif TPMS == 'G':
        F_final = (np.cos(2 * np.pi * x_final) * np.sin(2 * np.pi * y_final) + 
                   np.cos(2 * np.pi * y_final) * np.sin(2 * np.pi * z_final) + 
                   np.cos(2 * np.pi * z_final) * np.sin(2 * np.pi * x_final))
elif TPMS == 'D':
        F_final = (np.sin(2 * np.pi * x_final) * np.sin(2 * np.pi * y_final) * np.sin(2 * np.pi * z_final) + 
                   np.sin(2 * np.pi * x_final) * np.cos(2 * np.pi * y_final) * np.cos(2 * np.pi * z_final) + 
                   np.cos(2 * np.pi * x_final) * np.sin(2 * np.pi * y_final) * np.cos(2 * np.pi * z_final) + 
                   np.cos(2 * np.pi * x_final) * np.cos(2 * np.pi * y_final) * np.sin(2 * np.pi * z_final))
elif TPMS == 'I':
        F_final = (2.0 * (np.cos(2 * np.pi * x_final) * np.cos(2 * np.pi * y_final) + 
                          np.cos(2 * np.pi * y_final) * np.cos(2 * np.pi * z_final) + 
                          np.cos(2 * np.pi * z_final) * np.cos(2 * np.pi * x_final)) - 
                   (np.cos(4 * np.pi * x_final) + np.cos(4 * np.pi * y_final) + np.cos(4 * np.pi * z_final)))
elif TPMS == 'S':
        F_final = (np.cos(4 * np.pi * x_final) * np.sin(2 * np.pi * y_final) * np.cos(2 * np.pi * z_final) + 
                   np.cos(2 * np.pi * x_final) * np.cos(4 * np.pi * y_final) * np.sin(2 * np.pi * z_final) + 
                   np.sin(2 * np.pi * x_final) * np.cos(2 * np.pi * y_final) * np.cos(4 * np.pi * z_final))
elif TPMS == 'F':
        F_final = -(4.0 * (np.cos(2 * np.pi * x_final) * np.cos(2 * np.pi * y_final) * np.cos(2 * np.pi * z_final)) - 
                     (np.cos(4 * np.pi * x_final) * np.cos(4 * np.pi * y_final) + 
                      np.cos(4 * np.pi * y_final) * np.cos(4 * np.pi * z_final) + 
                      np.cos(4 * np.pi * z_final) * np.cos(4 * np.pi * x_final)))
elif TPMS == 'N':
        F_final = (3.0 * (np.cos(2 * np.pi * x_final) + np.cos(2 * np.pi * y_final) + np.cos(2 * np.pi * z_final)) + 
4.0 * np.cos(2 * np.pi * x_final) * np.cos(2 * np.pi * y_final) * np.cos(2 * np.pi * z_final))

if type == "net":
        F_final = F_final + rStart
else:  # sheet
        F_final = -(F_final + rStart) * (F_final - rStart)

# 提取最终等值面
if has_skimage:
try:
            verts_final, faces_final, normals_final, values_final = measure.marching_cubes(
                F_final, 0, spacing=(t_final, t_final, t_final)
            )
except:
            verts_final, faces_final, normals_final, values_final = measure.marching_cubes_lewiner(
                F_final, 0, spacing=(t_final, t_final, t_final)
            )
else:
        verts_final, faces_final, normals_final, values_final = simple_marching_cubes(F_final, 0, (t_final, t_final, t_final))

# 调整位置和尺寸
if len(verts_final) > 0:
        verts_final = verts_final - np.array(center)
        verts_final = verts_final * tarSize

# 计算最终体积分数
        totalVolume_final, totalArea_final = stlVolume(verts_final, faces_final)
        volFracFinal = -totalVolume_final / (tarSize ** 3)
        print(f"最终体积分数: {volFracFinal:.6f}")

# 保存为STL文件
        write_stl_binary('lattice_binary.stl', verts_final, faces_final, f'TPMS_{TPMS}_{type}')
        write_stl_ascii('lattice_ascii.stl', verts_final, faces_final, f'TPMS_{TPMS}_{type}')

        print(f"已生成 {len(verts_final)} 个顶点, {len(faces_final)} 个面")
else:
        print("警告: 未生成任何网格，请检查参数设置")
        volFracFinal = 0
        faces_final = np.array([])
        verts_final = np.array([])

    end_time = time.time()
    print(f"总耗时: {end_time - start_time:.2f} 秒")

return verts_final, faces_final, volFracFinal

# 主程序示例
if __name__ == "__main__":
# 示例1: 使用默认参数生成TPMS晶格
    print("生成TPMS晶格结构...")
    print("=" * 50)

    verts, faces, vol_frac = generateTPMS(
        TPMS='P',           # TPMS类型: P, G, D, I, S, F, N
        type="net",         # 结构类型: "net" 或 "sheet"
        volFrac=0.27,       # 目标体积分数
        numCell=1,          # 每个坐标方向的晶胞数量
        tarSize=10,         # 晶格的目标尺寸
        center=(0, 0, 0),   # 晶格中心坐标
        nFinal=50# 最终网格的细分数量（影响精度和文件大小）
    )

    print("=" * 50)
    print(f"生成完成!")
    print(f"顶点数: {len(verts)}")
    print(f"面数: {len(faces)}")
    print(f"体积分数: {vol_frac:.6f}")
    print("\n生成的文件:")
    print("- lattice_binary.stl (二进制格式，推荐)")
    print("- lattice_ascii.stl (ASCII格式，可读)")

# 示例2: 生成Gyroid TPMS
"""
    print("\n" + "="*50)
    print("生成Gyroid TPMS晶格...")

    verts, faces, vol_frac = generateTPMS(
        TPMS='G',
        type="sheet",
        volFrac=0.35,
        numCell=2,
        tarSize=20,
        nFinal=40
    )
    """

技术解析：算法核心与设计哲学

这个生成器的核心在于将抽象的隐函数表达式  转化为具体的三维三角网格。代码实现了七种经典的TPMS隐函数，每一种都对应着不同的数学对称性与物理特性。例如，Gyroid曲面（类型 'G'）在自然界中广泛存在，以其光滑、自支撑的特性闻名；而Primitive曲面（类型 'P'）则以其高刚度与规则孔隙在结构设计中备受青睐。

算法的精妙之处在于其自适应体积分数校准机制。用户设定目标体积分数（即材料密度）后，算法通过一个两阶段的优化流程自动确定正确的等值面阈值 。第一阶段进行粗校准，快速定位阈值范围；第二阶段采用二分法进行精细校准，确保最终生成的结构与目标体积分数的误差小于万分之一。这个过程完全自动化，将用户从繁琐的参数调试中解放出来。

在网格生成层面，代码采用了优雅的降级策略。如果检测到安装了scikit-image库，则调用其工业级的marching_cubes算法，生成质量极高的流形网格。如果没有检测到该库，则自动切换至内置的基础网格化算法。该基础算法通过简单的四面体分解来提取等值面，虽然生成的网格质量较低，但保证了代码在最小依赖环境下的可运行性，体现了“优先工作，优化在后”的实用主义设计哲学。

应用场景：从学术研究到产业创新

该工具的潜在应用广泛而深远。在生物医学工程领域，研究人员可以快速生成不同孔隙率和孔结构的Gyroid或Diamond晶格，用于骨组织工程支架的体外研究，通过调节参数来匹配天然骨的力学性能和渗透性。在航空航天领域，工程师可以利用Primitive或I-WP类型的高刚度网络结构，为飞机内饰或卫星部件设计重量极轻但承载能力强的点阵填充材料。

在学术层面，这个开源工具降低了计算设计领域的入门门槛。学生和研究人员可以轻松地批量生成不同拓扑类型和几何参数的结构，用于进行系统的力学性能模拟、流体渗透性分析或多物理场耦合研究，从而加速超材料与功能梯度材料的发现进程。

对于产品设计师和建筑师，这个脚本提供了探索复杂形式的新可能。通过调整晶胞数量、缩放比例和体积分数，可以生成从微观装饰纹理到大型建筑组件的各种尺度的复杂几何形态，无需掌握复杂的图形编程或依赖昂贵的三维建模软件。

扩展方向与社区生态构建

当前实现已经构建了一个坚实的功能基础，而其开源特性为社区共建提供了平台。一个重要的扩展方向是开发参数化渐变晶格的生成能力，即允许体积分数或单元尺寸在空间内按照指定函数连续变化，这对于模仿自然骨骼的密度梯度或设计性能定制化的功能部件至关重要。

另一个发展方向是集成实时可视化预览功能。虽然生成的STL文件可以被任何主流的三维软件打开，但内置一个基于matplotlib或pyvista的简单预览窗口，能让用户即时评估设计效果，进一步提升迭代效率。

性能方面，对于需要极高分辨率或大量晶胞的模型，核心的等值面提取循环可以通过numba的即时编译或利用多核CPU进行并行化来显著加速。此外，代码可以扩展支持更高效的3D打印文件格式，如3MF，该格式能存储颜色、纹理和多材料信息。

结语：开放代码驱动设计民主化

这个TPMS生成工具代表了计算设计领域的一个范式转变——将曾经局限于专业软件或特定编程环境的高级能力，通过清晰、简洁、自包含的Python代码 democratize（民主化）。它证明，通过精心的算法实现和分层的设计策略，复杂的科学计算功能可以变得易于获取且高度可定制。

无论是用于前沿的学术研究、创新的产品设计，还是作为教学演示的生动案例，这段代码都提供了一个强大的起点。它邀请使用者不仅是作为工具的用户，更可以作为共同改进者和扩展者参与其中。在开源共享的协作精神下，这样的工具将持续进化，不断降低创造复杂、智能、高性能三维结构的门槛，最终赋能于更广泛领域的创新与发现。