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

🎯要点

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

📜皂膜用例：Python计算物理粒子及拉格朗日和哈密顿动力学 | Python和MATLAB粘性力接触力动态模型半隐式欧拉算法

🍇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;

👉参阅：计算思维 | 亚图跨际