计算科学 - Python中的最小曲面解决方案 - 吾爱随笔录

Python中的最小曲面解决方案

计算科学计算几何插值网格生成几何学

2021-12-20 03:28:45

注意：这个问题也发布在 StackOverflow 和 math.stackexchange 中。

我有一组定义 3D 轮廓的 3D 点，如下所示。此轮廓中的点位于其最佳拟合平面中，我想获得此轮廓内表面的 3D 三角形网格表示。

在此处输入图像描述

做了一些研究，我发现这基本上是一个最小表面问题，它的解决方案与Biharmonic Equation相关。我还发现薄板样条是这个方程的基本解。

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

出现两个问题：

薄板样条插值是否是从 3D 轮廓点集计算表面问题的正确方法？
如果是这样，如何在具有非均匀网格的 scipy 上执行薄板插值？

2个回答

我不认为你可以插值。你可能知道控制方程的基本解，但你没有关于问题内部点的数据，我怀疑远离边界的解会很差。

此外，已经有一段时间了，但我通常将以下方程与最小表面相关联（至少对于小位移）：其中是相关表面的小位移。

(1 + u_{x}^{2}) u_{y y} - 2 u_{x} u_{y} u_{x y} + (1 + u_{y}^{2}) u_{x x} = 0

$(1+u_x^2)u_{yy} - 2u_x u_y u_{xy} + (1+u_y^2)u_{xx} = 0$

u (x, y)

$u(x,y)$

您可以通过将边界点向下投影到平面并将其值用作边界上的使用边界点的坐标作为域的边界，您可以在内部创建一个三角形网格，然后使用上述边界条件求解上述方程。 $(x,y)$ $z$ $u(x,y)$ $(x,y)$

deal.II 的示例#15做到了这一点，尽管我认为你需要一个四边形网格而不是三角形。我无法从图中看出，但如果您的边界名义上是一个圆形，那么为其生成四边形网格应该很简单。如果它是任意的，那么你将不得不做更多的工作。

如果您更喜欢 Python 而不是 C++，则可能还有一个FEniCS示例，但我无法通过网络上的快速搜索来判断。如果您可以访问，还有一个使用 PDE 工具箱的快速MATLAB 示例。它甚至会为您生成网格。

您可以使用 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。）

其它你可能感兴趣的问题

上一篇“伞式抽样”的词源下一篇FEniCS：自定义正交规则