25

我有一组定义 3D 轮廓的 3D 点。我想要做的是获得与这个轮廓相对应的最小表面表示（参见维基百科中的最小表面）。基本上这需要求解非线性偏微分方程。

在 Matlab 中，使用该pdenonlin函数几乎很简单（参见 Matlab 的文档）。可以在此处找到其用于解决最小曲面问题的示例：单元盘上的最小曲面问题。

我需要在 Python 中进行这样的实现，但据我所知，我还没有找到任何关于如何实现的网络资源。

谁能指出这种实现的任何资源/示例？

谢谢，米格尔。

更新

我想要找到的 3D 表面（理想情况下是三角形网格表示）由这组 3D 点限定（如图所示，这些点位于最佳拟合平面中）：

在此处输入图像描述

好的，所以做了一些研究，我发现这个最小表面问题与双调和方程的解有关，我还发现薄板样条是这个方程的基本解。

所以我认为该方法是尝试使用薄板样条曲线拟合表面的这种稀疏表示（由点的 3D 轮廓给出）。我在 scipy.interpolate 中找到了这个示例，其中使用薄板样条插值分散数据（x，y，z 格式）以获得均匀网格（XI，YI）上的 ZI 坐标。

出现了两个问题：（1）薄板样条插值是否是从 3D 轮廓点集计算表面问题的正确方法？（2）如果是这样，如何使用非均匀网格对scipy进行薄板插值？

再次感谢！米格尔

更新：在 MATLAB 中的实现（但它不适用于 SCIPY PYTHON）

我使用 Matlab 的函数遵循这个例子tpaps，并在均匀的网格上获得了适合我的轮廓的最小表面。这是 Matlab 中的结果（看起来很棒！）： 在此处输入图像描述

但是我需要在 Python 中实现它，所以我使用包 scipy.interpolate.Rbf 和thin-plate函数。这是python中的代码（XYZ包含轮廓中每个点的3D坐标）：

GRID_POINTS = 25
x_min = XYZ[:,0].min()
x_max = XYZ[:,0].max()
y_min = XYZ[:,1].min()
y_max = XYZ[:,1].max()
xi = np.linspace(x_min, x_max, GRID_POINTS)
yi = np.linspace(y_min, y_max, GRID_POINTS)
XI, YI = np.meshgrid(xi, yi)

from scipy.interpolate import Rbf
rbf = Rbf(XYZ[:,0],XYZ[:,1],XYZ[:,2],function='thin-plate',smooth=0.0)
ZI = rbf(XI,YI)

然而，这是结果（与在 Matlab 中获得的结果大不相同）：

在此处输入图像描述

很明显，scipy 的结果并不对应于最小表面。

scipy.interpolate.Rbf + thin-plate 是否按预期进行，为什么它与 Matlab 的结果不同？

4

3 回答 3

1

问题表明我们需要求解非线性偏微分方程。然而，维基百科指出“它们很难研究：几乎没有适用于所有这些方程的通用技术，通常每个单独的方程都必须作为一个单独的问题来研究。” 但是，您没有给出方程式！Matlab 有时会使用遗传算法来得出它的表面吗？也就是说，它是否使用经验法则做出最佳猜测，然后尝试在组件方格中进行小的变化，直到找不到更小的表面。实施这种解决方案会很费力，但在概念上并不困难（假设你喜欢那种东西）。还要记住，连续函数的微积分只是函数的所有线性逼近的微积分的一个特例（增量设置为零而不是某个有限值）。通过阅读 JL Bell 关于平滑无穷小分析的书籍，我清楚地知道了这一点——只需使用具有有限增量的代数，并将结果因子留在推导中，而不是“忽略”它们。

于 2013-05-20T10:56:19.253 回答
0

显然，Matlab 和 SciPy 以不同的方式理解 TPS。Matlab 实现看起来是正确的。SciPy 以与其他 RBF 相同的方式处理 TPS，因此您可以自己在 Python 中正确实现它 - 形成相关线性方程组的矩阵并求解它以接收 TPS 的系数就足够了。

于 2013-09-29T06:56:51.183 回答
0

您可以使用 FEniCS：

from fenics import (
    UnitSquareMesh,
    FunctionSpace,
    Expression,
    interpolate,
    assemble,
    sqrt,
    inner,
    grad,
    dx,
    TrialFunction,
    TestFunction,
    Function,
    solve,
    DirichletBC,
    DomainBoundary,
    MPI,
    XDMFFile,
)

# Create mesh and define function space
mesh = UnitSquareMesh(100, 100)
V = FunctionSpace(mesh, "Lagrange", 2)

# initial guess (its boundary values specify the Dirichlet boundary conditions)
# (larger coefficient in front of the sin term makes the problem "more nonlinear")
u0 = Expression("a*sin(2.5*pi*x[1])*x[0]", a=0.2, degree=5)
u = interpolate(u0, V)
print(
    "initial surface area: {}".format(assemble(sqrt(1 + inner(grad(u), grad(u))) * dx))
)

# Define the linearized weak formulation for the Newton iteration
du = TrialFunction(V)
v = TestFunction(V)
q = (1 + inner(grad(u), grad(u))) ** (-0.5)
a = (
    q * inner(grad(du), grad(v)) * dx
    - q ** 3 * inner(grad(u), grad(du)) * inner(grad(u), grad(v)) * dx
)
L = -q * inner(grad(u), grad(v)) * dx
du = Function(V)

# Newton iteration
tol = 1.0e-5
maxiter = 30
for iter in range(maxiter):
    # compute the Newton increment by solving the linearized problem;
    # note that the increment has *homogeneous* Dirichlet boundary conditions
    solve(a == L, du, DirichletBC(V, 0.0, DomainBoundary()))
    u.vector()[:] += du.vector()  # update the solution
    eps = sqrt(
        abs(assemble(inner(grad(du), grad(du)) * dx))
    )  # check increment size as convergence test
    area = assemble(sqrt(1 + inner(grad(u), grad(u))) * dx)
    print(
        f"iteration{iter + 1:3d}  H1 seminorm of delta: {eps:10.2e}  area: {area:13.5e}"
    )
    if eps < tol:
        break

if eps > tol:
    print("no convergence after {} Newton iterations".format(iter + 1))
else:
    print("convergence after {} Newton iterations".format(iter + 1))

with XDMFFile(MPI.comm_world, "out.xdmf") as xdmf_file:
    xdmf_file.write(u)

（修改自http://www-users.math.umn.edu/~arnold/8445/programs/minimalsurf-newton.py。）

于 2019-12-19T13:36:42.640 回答