🍍C++(Python)肥皂泡沫普拉托边界膜曲面模型算法

Python | C++ | MATLAB | 数学 | 算法 | 物理 | 流体力学 | 泡沫 | 皂膜 | 曲面模型 | 厚度

C++(Python)肥皂泡沫普拉托边界膜曲面模型算法 | 亚图跨际亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

🎯肥皂泡二维流体模拟 | 🎯泡沫普拉托边界膜曲面模型算法演化厚度变化 | 🎯螺旋曲面三周期最小结构生成

📜皂膜用例：Python计算物理粒子及拉格朗日和哈密顿动力学

🍇Python(MATLAB)三维曲面图

普拉托定律描述了肥皂膜的结构。这些定律是由比利时物理学家普拉托在 19 世纪根据他的实验观察制定的。自然界中的许多图案都是基于遵守这些定律的泡沫。此定律描述了皂膜的形状和构造如下：

• 肥皂膜由整个（完整的）光滑表面制成

• 肥皂膜的一部分的平均曲率在同一块肥皂膜上的任何点处处处恒定

• 肥皂膜总是沿着称为普拉托边界的边缘成三块相交，并且它们以 arccos(− 1 / 2 ) = 120°‘

• 这些普拉托边界在一个顶点处以四面体角 arccos(− 1 / 3 ) ≈ 109.47°。

除普拉托定律之外的结构是不稳定的，并且薄膜将很快倾向于重新排列自身以符合这些定律。美国数学家吉恩·艾伦·泰勒使用几何测量理论在数学上证明了这些定律适用于最小曲面。

肥皂膜是被空气包围的薄层液体（通常为水基）。例如，如果两个肥皂泡接触，它们会合并并在其间形成一层薄膜。因此，泡沫由通过普拉托边界连接的薄膜网络组成。肥皂膜可用作极小曲面的模型系统，极小曲面在数学中被广泛使用。

从数学角度来看，肥皂膜是最小表面。表面张力是单位面积产生表面所需的能量。薄膜——与任何物体或结构一样——倾向于以最小势能状态存在。为了最小化其能量，自由空间中的液滴自然呈现球形，对于给定的体积，其表面积最小。水坑和薄膜可以在其他力的存在下存在，例如重力和对基质原子的分子间吸引力。后一种现象称为润湿：基质原子和薄膜原子之间的结合力会导致总能量降低。在这种情况下，物体的最低能量配置是尽可能多的薄膜原子尽可能靠近基质。这将导致无限薄的薄膜，无限广泛地分布在基质上。实际上，粘附润湿效应（导致表面最大化）和表面张力效应（导致表面最小化）会相互平衡：稳定的结构可以是液滴、水坑或薄膜，具体取决于作用于身体的力。

💦Python绘制三维曲面：

Matplotlib 的 mpl_toolkits.mplot3d 工具包中的 axis3d 提供了用于创建三维曲面图的必要函数。曲面图是通过使用 ax.plot_surface() 函数创建的。

语法：

 ax.plot_surface(X, Y, Z)

其中 X 和 Y 是 x 和 y 点的二维数组，而 Z 是高度的二维数组。

示例：

 # Import libraries
 from mpl_toolkits import mplot3d
 import numpy as np
 import matplotlib.pyplot as plt
 ​
 x = np.outer(np.linspace(-3, 3, 32), np.ones(32))
 y = x.copy().T # transpose
 z = (np.sin(x **2) + np.cos(y **2) )
 ​
 fig = plt.figure(figsize =(14, 9))
 ax = plt.axes(projection ='3d')
 ​
 ax.plot_surface(x, y, z)
 plt.show()

梯度曲面图是三维曲面图与二维轮廓图的组合。在此图中，三维曲面的颜色与二维轮廓图相同。曲面较高的部分与曲面较低的部分颜色不同。

语法：

 surf = ax.plot_surface(X, Y, Z, cmap=, linewidth=0, antialiased=False)

属性 cmap= 设置表面的颜色。还可以通过调用Fig.colorbar来添加颜色条。下面的代码创建一个梯度曲面图：

 from mpl_toolkits import mplot3d
 import numpy as np
 import matplotlib.pyplot as plt
 ​
 x = np.outer(np.linspace(-3, 3, 32), np.ones(32))
 y = x.copy().T # transpose
 z = (np.sin(x **2) + np.cos(y **2) )
 ​
 fig = plt.figure(figsize =(14, 9))
 ax = plt.axes(projection ='3d')
 my_cmap = plt.get_cmap('hot')
 ​
 surf = ax.plot_surface(x, y, z,
                     cmap = my_cmap,
                     edgecolor ='none')
 ​
 fig.colorbar(surf, ax = ax,
             shrink = 0.5, aspect = 5)
 ​
 ax.set_title('Surface plot')
 plt.show()
 ​

使用 Matplotlib 绘制的三维曲面图可以投影到二维曲面上。下面的代码创建三维图并可视化其在二维等高线图上的投影：

 from mpl_toolkits import mplot3d
 import numpy as np
 import matplotlib.pyplot as plt
 ​
 x = np.outer(np.linspace(-3, 3, 32), np.ones(32))
 y = x.copy().T # transpose
 z = (np.sin(x **2) + np.cos(y **2) )
 ​
 fig = plt.figure(figsize =(14, 9))
 ax = plt.axes(projection ='3d')
 ​
 my_cmap = plt.get_cmap('hot')
 ​
 surf = ax.plot_surface(x, y, z, 
                     rstride = 8,
                     cstride = 8,
                     alpha = 0.8,
                     cmap = my_cmap)
 cset = ax.contourf(x, y, z,
                 zdir ='z',
                 offset = np.min(z),
                 cmap = my_cmap)
 cset = ax.contourf(x, y, z,
                 zdir ='x',
                 offset =-5,
                 cmap = my_cmap)
 cset = ax.contourf(x, y, z, 
                 zdir ='y',
                 offset = 5,
                 cmap = my_cmap)
 fig.colorbar(surf, ax = ax, 
             shrink = 0.5,
             aspect = 5)
 ​
 ​
 ax.set_xlabel('X-axis')
 ax.set_xlim(-5, 5)
 ax.set_ylabel('Y-axis')
 ax.set_ylim(-5, 5)
 ax.set_zlabel('Z-axis')
 ax.set_zlim(np.min(z), np.max(z))
 ax.set_title('3D surface having 2D contour plot projections')
 ​
 plt.show()

💦MATLAB最小曲面

构造最小曲面

 n=35;
 [X,Y,Z]=meshgrid(linspace(-pi,pi,n));
 [F,V] = isosurface(X,Y,Z,S,0);

曲面可视化

 cFigure;
 hold on;
 title('Schwarz P-surface','FontSize',fontSize);
 gpatch(F,V,'kw','k',1);
 axisGeom;
 camlight headlight;
 drawnow;

可视化曲面所有变体：

 cFigure;
 ​
 subplot(2,3,1);
 title('Schwarz P-surface','FontSize',fontSize);
 hold on;
 gpatch(F,V,pColors(1,:),'none',1);
 axisGeom;
 camlight headlight; lighting gouraud;
 view(-50,30);
 ​
 [X,Y,Z]=meshgrid(linspace(-pi,pi,n));
 [F,V] = isosurface(X,Y,Z,S,0.1);
 ​
 subplot(2,3,2);
 title('d','FontSize',fontSize);
 hold on;
 gpatch(F,V,pColors(2,:),'none',1);
 axisGeom;
 camlight headlight; lighting gouraud;
 view(-50,30);
 ​
 [X,Y,Z]=meshgrid(linspace(-2*pi,2*pi,n));
 [F,V] = isosurface(X,Y,Z,S,0.6);
 ​
 subplot(2,3,3);
 title('Gyroid','FontSize',fontSize);
 hold on;
 gpatch(F,V,pColors(3,:),'none',1);
 axisGeom;
 camlight headlight; lighting gouraud;
 view(-50,30);
 ​
 [X,Y,Z]=meshgrid(linspace(-pi,pi,n));
 [F,V] = isosurface(X,Y,Z,S,0);
 ​
 subplot(2,3,4);
 title('Neovius','FontSize',fontSize);
 hold on;
 gpatch(F,V,pColors(4,:),'none',1);
 axisGeom;
 camlight headlight; lighting gouraud;
 view(-50,30);
 ​
 [X,Y,Z]=meshgrid(linspace(-pi,pi,n));
 [F,V] = isosurface(X,Y,Z,S,0);
 ​
 subplot(2,3,5);
 title('w','FontSize',fontSize);
 hold on;
 gpatch(F,V,pColors(5,:),'none',1);
 axisGeom;
 camlight headlight; lighting gouraud;
 view(-50,30);
 ​
 [X,Y,Z]=meshgrid(linspace(-pi,pi,n));
 [F,V] = isosurface(X,Y,Z,S,0.5);
 ​
 subplot(2,3,6);
 title('pw','FontSize',fontSize);
 hold on;
 gpatch(F,V,pColors(6,:),'none',1);
 axisGeom;
 camlight headlight; lighting gouraud;
 view(-50,30);
 drawnow;