我正在使用行进立方体从体积中提取二维表面。在这种情况下是陀螺仪。
import numpy as np
from numpy import sin, cos, pi
from skimage import measure
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def gyroid(x, y, z, t):
return cos(x)*sin(y) + cos(y)*sin(z) + cos(z)*sin(x) - t
lattice_param = 1.0
strut_param = 0.0
resolution = 31j
x, y, z = pi*np.mgrid[-1:1:resolution, -1:1:resolution, -1:1:resolution] * lattice_param
vol = gyroid(x, y, z, strut_param)
verts, faces = measure.marching_cubes(vol, 0, spacing=(0.1, 0.1, 0.1)) # , normals, values
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='ocean', lw=1)
这一切都工作正常,但网格质量在很多地方都令人震惊。我无法在网格上运行任何 FEA,因为许多元素/面的面积接近于零或高度扭曲。
是否有一种方法可以在给定顶点的情况下重新划分网格并确保特定的元素面/面指标(例如纵横比)或强制行进立方体执行此类操作?
只要网格是一个公平的近似值,我就不会担心移动顶点。
一种选择可能是使用表面网格划分包“重新网格化”行进立方体输出。本质上,这意味着行进立方体三角剖分将作为要重新三角剖分的初始表面定义。
可以采用多种技术来做到这一点。一些可能有用的选项(所有
C++ / C
实现):
JIGSAW
:一种restricted、frontal-delaunay算法^^,通常会构建非常高质量的表面Delaunay三角剖分。对于所示的对象类型,我希望它能很好地工作。在包含的演示中(MATLAB
中提供),几个示例解决了行进立方体输出的重新网格化问题。 CGAL
:一种restricted、delaunay-refinement方法,也可以构建表面Delaunay三角剖分,但使用与JIGSAW
略有不同的算法,并且还包括CVT
型网格优化方案。MMG
:重新网格划分/优化策略的集合(据我所知),可用于通过局部修改的迭代应用来转换(从而改进)初始网格。^^
我是JIGSAW
的作者,所以,这里基本上是无耻的推广。
我建议使用带有平滑功能的 3D 自适应网格重新划分。 Geogram 编程库 (http://alice.loria.fr/software/geogram/doc/html/index.html) 提供了一个很好的实现。请参阅:https://twitter.com/brunolevy01/status/1132343120690122752?lang=en
也可以使用pygalmesh进行表面重新网格化。 (这是一个易于使用的 CGAL 界面。)
import pygalmesh
# create verts, faces
meshio.write_points_cells("in.vtu", verts, [("triangle", faces)])
mesh = pygalmesh.remesh_surface(
"in.vtu",
edge_size=0.025,
facet_angle=25,
facet_size=0.1,
facet_distance=0.001,
verbose=False,
)
# mesh.points, mesh.cells
奇怪的三角形来自奇怪的数据,而不是来自三角测量所使用的方法。
我可以说德劳内三角剖分实现了最佳的三角形面积/三角形周长比(最佳的理论比例是等边三角形)。但您不能将 if 用于网格,因为 Delaunay 三角剖分会输出凸网格。
您面临着一项艰巨的任务。一些想法:
z
都可以遍历网格,直到找到包含网格中点的 x,y
坐标的三角形,然后对 z
进行插值。如果您有一些额外的信息,例如邻居或层次结构,则搜索可以比检查每个三角形更快地完成。对网格进行三角测量非常简单。我没有足够的声誉来发表评论,所以我将此作为答案提交......
添加到 Nico 的答案中,似乎
pygalmesh.remesh_surface()
的关键字参数已更改。链接中提供的示例是:
mesh = pygalmesh.remesh_surface(
"lion-head.off",
max_edge_size_at_feature_edges=0.025,
min_facet_angle=25,
max_radius_surface_delaunay_ball=0.1,
max_facet_distance=0.001,
verbose=False,
)
此外,
meshio.write_points_cells
的第三个参数定义了面,对于三角形元素,faces
必须具有形状(n_faces, 3)
。就我而言,我使用的网格与另一个包 (PyVista) 的格式不同,并且必须重新调整数组的形状。