4. Gmsh 示例 Gmsh Examples#
Overview
Generating meshes with the Gmsh Python API: a Möbius strip, extrusion for PML layers, and notes on Gmsh node ordering for high-order elements. Requires the gmsh package (this notebook is not executed during the book build).
4.1. 莫比乌斯带#
Firedrake 没有内置的莫比乌斯带网格. 它是嵌在三维空间里的二维曲面, 只能借助外部
网格生成器按参数方程一段段拼出来. 本节先用 Gmsh 的 Python API 把这条带子造出来,
再把生成的 .msh 文件交给 Firedrake, 在它上面解一个曲面 Poisson 问题.
4.1.1. 用 matplotlib 绘图#
先按参数方程把曲面画出来, 对下一段脚本里的坐标公式有个直观印象. 取 \(u \in [0, 2\pi]\) 沿中心圆周绕行, \(v \in [-1/2, 1/2]\) 沿带宽方向, 令 \(r = R - W v \sin(u/2)\), 则曲面上的点为 \((r\cos u,\ r\sin u,\ W v \cos(u/2))\). \(u\) 每绕一整圈, 截面就转过半圈, 首尾 接上时带子的两条边正好互换, 这就是莫比乌斯带.
本单元只是画图, 与后面的网格生成互不影响, 所以这里的 R、W 取值和下一段脚本
里的也不相同.
import gmsh
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
R = 10
W = 5
_u = np.linspace(0, 2*np.pi, 201, endpoint=True)
_v = np.linspace(-1/2, 1/2, 10)
u, v = np.meshgrid(_u, _v)
r = R - W*v*np.sin(u/2)
x = r*np.cos(u)
y = r*np.sin(u)
z = W*v*np.cos(u/2)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_surface(x, y, z)
zlim = ax.set_zlim([-2*W, 2*W])
4.1.2. Python 脚本#
下面用 Gmsh 的 Python API 把这条带子”缝”出来: 沿 \(u\) 方向分成 M 段, 每段在带宽
两侧各放一个点, 相邻两段的四个点围成一个曲线环 (add_curve_loop), 再填成一个
平面片 (add_plane_surface). 关键在收尾的两条线: 最后一段与第 0 段是交叉相连的
(points[0][-1] 接 points[1][0], points[1][-1] 接 points[0][0]), 半圈的
扭转就是这样闭合成莫比乌斯带的.
Mesh.MeshSizeMin/Mesh.MeshSizeMax 把网格尺寸控制在 h 附近,
gmsh.model.mesh.generate() 才真正开始剖分. 剖分前把 General.Verbosity 调低
以免刷屏; 剖分失败时 Gmsh 会抛出异常, 因此用 try/finally 保证输出等级一定
被还原.
写文件的 gmsh.write 被注释掉了: 仓库里已经存有生成好的 gmsh/mobius.msh,
取消注释会把它覆盖掉. 末尾的 gmsh.fltk.run() 只在命令行带 -popup 参数时才
打开 Gmsh 的图形界面, 在 Notebook 里不会触发.
import sys
import gmsh
import numpy as np
gmsh.initialize()
fac = gmsh.model.geo
R = 5
W = 2
M = 50
_u = np.linspace(0, 2*np.pi, M+1, endpoint=True)
_v = np.linspace(-1/2, 1/2, 2)
u, v = np.meshgrid(_u, _v)
r = R - W*v*np.sin(u/2)
x = r*np.cos(u)
y = r*np.sin(u)
z = W*v*np.cos(u/2)
points = [[], []]
line_current = None
ss = []
for i in range(M):
points[0].append(fac.add_point(x[0, i], y[0, i], z[0, i]))
points[1].append(fac.add_point(x[1, i], y[1, i], z[1, i]))
line_pre = line_current
line_current = fac.add_line(points[0][-1], points[1][-1])
if i == 0:
line0 = line_current
else:
j = i - 1
line3 = fac.add_line(points[0][i-1], points[0][i])
line4 = fac.add_line(points[1][i-1], points[1][i])
cl = fac.add_curve_loop([line_pre, line4, -line_current, -line3])
ss.append(fac.add_plane_surface([cl]))
line3 = fac.add_line(points[0][-1], points[1][0])
line4 = fac.add_line(points[1][-1], points[0][0])
cl = fac.add_curve_loop([line_current, line4, line0, -line3])
ss.append(fac.add_plane_surface([cl]))
fac.add_surface_loop(ss)
fac.synchronize()
h = 0.4
gmsh.option.setNumber("Mesh.MeshSizeMin", h)
gmsh.option.setNumber("Mesh.MeshSizeMax", h)
# 剖分失败时 gmsh 会抛出异常, 用 try/finally 保证输出等级被还原
old_verbosity = gmsh.option.getNumber("General.Verbosity")
gmsh.option.setNumber("General.Verbosity", 1)
try:
gmsh.model.mesh.generate()
finally:
gmsh.option.setNumber("General.Verbosity", old_verbosity)
# gmsh.write('gmsh/mobius.msh')
# gui or not
if '-popup' in sys.argv:
gmsh.fltk.run()
gmsh.finalize()
4.1.3. 在莫比乌斯带上求解 Poisson 问题#
Mesh 读入 .msh 时用 dim=3 指明几何维数是 3: 这张网格的单元是三角形
(拓扑维数为 2), 但顶点坐标是三维的. triplot 画出读进来的网格, 可以和上面
matplotlib 画的曲面对照着看.
from firedrake import *
msh = Mesh('gmsh/mobius.msh', dim=3)
triplot(msh)
曲面上的 Poisson 问题写法和平面上完全一样: dx 在这里是曲面上的测度, grad 是
切向梯度, 因此下面的变分形式离散的其实是曲面上的 Laplace-Beltrami 算子. 右端项
\(f = x^2 + y^2 + z^2\) 由 SpatialCoordinate 给出,
DirichletBC(V, 0, 'on_boundary') 在带子的外边界上取零.
这里把双线性形式和右端项写在同一个表达式 a 里, 再用 lhs/rhs 把它拆开.
注意 u 先后被绑定了两次: 先是 TrialFunction, 后是存放解的 Function; a
在重新绑定之前就已经建好了, 所以不受影响, 但初学时容易看错.
x, y, z = SpatialCoordinate(msh)
f = x**2 + y**2 + z**2
V = FunctionSpace(msh, 'CG', 1)
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v))*dx - inner(f, v)*dx
u = Function(V)
bc = DirichletBC(V, 0, 'on_boundary')
solve(lhs(a) == rhs(a), u, bcs=bc)
c = trisurf(u)
4.2. 为 PML 拉伸生成网格#
完美匹配层 (PML, perfectly matched layer) 是波动问题里常用的人工吸收层: 在物理 区域外面再套一圈薄层, 让外行波在层内衰减而几乎不反射回来. 这对网格提出两个 要求: 物理区域外要多出一层贴着边界的单元, 并且物理区域和吸收层必须能分开标记.
Gmsh 的 extrude_boundary_layer 正好做这件事: 把一条边界曲线 (二维) 或一张
边界曲面 (三维) 沿法向拉伸出若干层单元. 下面两小节的结构完全一样: 先造几何体,
取出它的边界, 沿边界拉伸出 PML 层, 再给物理区域、PML 层和内外边界分别打上
physical group 标记.
这些标记读进 Firedrake 之后, 物理区域和 PML 层的标记就是 dx(1)、dx(2) 里的
子域号. 边界标记则要当心: 物理区域与 PML 层的交界面拉伸之后是内部面, 只能用
dS(...) 积分, 用 ds(...) 会得到 0. 本地实测 gmsh/disk_pml.msh 的 tag 1
(inner bdy) ds 积分为 0、dS 积分为 3.13, tag 2 (outer bdy) 才是真正的外边界
(ds 积分 3.76); gmsh/sphere_pml.msh 同理, tag 2 是内部面, tag 1、3、4 才是
外边界. 弄错时 Firedrake 只给一条 Subdomain is empty 的 warning, 积分会悄悄
算成 0.
4.2.1. 三维#
gmsh.model.occ.add_sphere(0, 0, 0, 0.5, angle1=0) 用 OCC 内核造一个半径 0.5 的
半球体 (angle1=0 把极角的起点从 \(-\pi/2\) 抬到 \(0\), 只保留上半部分), 因此它的
边界既有球面又有底面圆盘, 需要用 gmsh.model.get_type 按类型区分开. 随后只对
球面 (inner_bdy) 拉伸出边界层作为 PML, 并给半球体、PML 层以及内外边界分别
编号. 仓库里存好的 gmsh/sphere_pml.msh 实测外半径为 0.6, 即球半径 0.5 加上
0.1 厚的一层.
import sys
import gmsh
import numpy as np
gmsh.initialize()
sphere = gmsh.model.occ.add_sphere(0, 0, 0, 0.5, angle1=0)
gmsh.model.occ.synchronize()
bdy = gmsh.model.get_boundary([[3, sphere]])
plane = []
for dim, tag in bdy:
name = gmsh.model.get_type(dim, abs(tag))
if name == 'Sphere':
inner_bdy = tag
elif name == 'Plane':
plane.append(abs(tag))
# heights = [0.1 for _ in range(8)]
# change recombine to True to get mix-cell mesh
top = gmsh.model.geo.extrude_boundary_layer([[2, inner_bdy]], numElements=[2], heights=[0.1], recombine=False)
gmsh.model.geo.synchronize()
gmsh.model.add_physical_group(2, plane, tag=1)
gmsh.model.set_physical_name(2, 1, "plane")
gmsh.model.add_physical_group(2, [inner_bdy], tag=2)
gmsh.model.set_physical_name(2, 2, "inner bdy")
gmsh.model.add_physical_group(2, [top[3][1]], tag=3)
gmsh.model.set_physical_name(2, 3, "plane2")
gmsh.model.add_physical_group(2, [top[0][1]], tag=4)
gmsh.model.set_physical_name(2, 4, "outer bdy")
gmsh.model.add_physical_group(3, [sphere], tag=1)
gmsh.model.set_physical_name(3, 1, "Domain")
gmsh.model.add_physical_group(3, [top[1][1]], tag=2)
gmsh.model.set_physical_name(3, 2, "PML")
gmsh.model.geo.synchronize()
# for dim, tag in top:
# name = gmsh.model.get_type(dim, abs(tag))
# print(name, dim, tag)
# gmsh.fltk.run()
gmsh.model.mesh.generate()
# gmsh.write('gmsh/sphere_pml.msh')
gmsh.finalize()
Info : Skipping boundary layer extrusion of degenerate curve 1
Info : Meshing 1D...
Info : [ 20%] Meshing curve 2 (Circle)
Info : [ 30%] Meshing curve 3 (Circle)
Info : Done meshing 1D (Wall 0.0001125s, CPU 0.000119s)
Info : Meshing 2D...
Info : 0 dependencies in mesh of source faces
Info : Meshing surface 1 (Sphere, Frontal-Delaunay)
Info : Meshing curve 5
Info : Meshing curve 5 (Extruded)
Info : Meshing curve 6
Info : Meshing curve 6 (Extruded)
Info : Meshing curve 9
Info : Meshing curve 9 (Extruded)
Info : Meshing curve 10
Info : Meshing curve 10 (Extruded)
Info : Meshing surface 2 (Plane, Frontal-Delaunay)
Info : Meshing surface 11 (Extruded)
Info : Meshing surface 15 (Extruded)
Info : Meshing surface 20 (Extruded)
Info : Meshing surface 11 (Boundary layer surface)
Info : Meshing surface 11 (Extruded)
Info : Meshing surface 15 (Boundary layer surface)
Info : Meshing surface 15 (Extruded)
Info : Meshing surface 20 (Boundary layer surface)
Info : Meshing surface 20 (Extruded)
Info : Done meshing 2D (Wall 0.00580604s, CPU 0.005009s)
Info : Meshing 3D...
Info : Meshing volume 2 (Extruded)
Info : Subdividing extruded mesh
Info : Swapping 6
Info : Swapping 0
Info : Remeshing surface 11
Info : Meshing surface 11 (Extruded)
Info : Remeshing surface 15
Info : Meshing surface 15 (Extruded)
Info : 3D Meshing 1 volume with 1 connected component
Info : Tetrahedrizing 149 nodes...
Info : Done tetrahedrizing 157 nodes (Wall 0.000950333s, CPU 0.000659s)
Info : Reconstructing mesh...
Info : - Creating surface mesh
Info : - Identifying boundary edges
Info : - Recovering boundary
Info : Done reconstructing mesh (Wall 0.00207496s, CPU 0.001475s)
Info : Found volume 1
Info : It. 0 - 0 nodes created - worst tet radius 1.64893 (nodes removed 0 0)
Info : 3D refinement terminated (303 nodes total):
Info : - 0 Delaunay cavities modified for star shapeness
Info : - 0 nodes could not be inserted
Info : - 495 tetrahedra created in 0.00060025 sec. (824656 tets/s)
Info : 0 node relocations
Info : Done meshing 3D (Wall 0.00691087s, CPU 0.005343s)
Info : Optimizing mesh...
Info : Optimizing volume 1
Info : Optimization starts (volume = 0.25424) with worst = 0.0255552 / average = 0.716624:
Info : 0.00 < quality < 0.10 : 3 elements
Info : 0.10 < quality < 0.20 : 4 elements
Info : 0.20 < quality < 0.30 : 12 elements
Info : 0.30 < quality < 0.40 : 12 elements
Info : 0.40 < quality < 0.50 : 17 elements
Info : 0.50 < quality < 0.60 : 36 elements
Info : 0.60 < quality < 0.70 : 94 elements
Info : 0.70 < quality < 0.80 : 159 elements
Info : 0.80 < quality < 0.90 : 110 elements
Info : 0.90 < quality < 1.00 : 48 elements
Info : 19 edge swaps, 0 node relocations (volume = 0.25424): worst = 0.313387 / average = 0.738245 (Wall 0.000205958s, CPU 0.000167s)
Info : No ill-shaped tets in the mesh :-)
Info : 0.00 < quality < 0.10 : 0 elements
Info : 0.10 < quality < 0.20 : 0 elements
Info : 0.20 < quality < 0.30 : 0 elements
Info : 0.30 < quality < 0.40 : 11 elements
Info : 0.40 < quality < 0.50 : 17 elements
Info : 0.50 < quality < 0.60 : 37 elements
Info : 0.60 < quality < 0.70 : 95 elements
Info : 0.70 < quality < 0.80 : 159 elements
Info : 0.80 < quality < 0.90 : 112 elements
Info : 0.90 < quality < 1.00 : 48 elements
Info : Done optimizing mesh (Wall 0.000420666s, CPU 0.000389s)
Info : 388 nodes 2324 elements
4.2.2. 二维#
二维版本与三维完全平行: add_disk 造一个半径 0.5 的圆盘, 它的边界是一条曲线,
对这条曲线拉伸出边界层. numElements=[2]、heights=[0.1] 表示沿法向分 2 层、
总厚度 0.1, 生成的 gmsh/disk_pml.msh 外半径实测正是 0.6. 把 recombine 改成
True 可以得到含四边形的混合网格.
import sys
import gmsh
import numpy as np
gmsh.initialize()
disk = gmsh.model.occ.add_disk(0, 0, 0, 0.5, 0.5)
gmsh.model.occ.synchronize()
bdy = gmsh.model.get_boundary([[2, disk]])
inner_bdy = bdy[0][1]
# heights = [0.1 for _ in range(8)]
# change recombine to True to get mix-cell mesh
top = gmsh.model.geo.extrude_boundary_layer([[1, inner_bdy]], numElements=[2], heights=[0.1], recombine=False)
gmsh.model.geo.synchronize()
gmsh.model.add_physical_group(1, [inner_bdy], tag=1)
gmsh.model.set_physical_name(1, 1, "inner bdy")
gmsh.model.add_physical_group(1, [top[0][1]], tag=2)
gmsh.model.set_physical_name(1, 2, "outer bdy")
gmsh.model.add_physical_group(2, [disk], tag=1)
gmsh.model.set_physical_name(2, 1, "Domain")
gmsh.model.add_physical_group(2, [top[1][1]], tag=2)
gmsh.model.set_physical_name(2, 2, "PML")
gmsh.model.geo.synchronize()
# for dim, tag in top:
# name = gmsh.model.get_type(dim, abs(tag))
# print(name, dim, tag)
# gmsh.fltk.run()
gmsh.model.mesh.generate()
# gmsh.write('gmsh/disk_pml.msh')
gmsh.finalize()
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Ellipse)
Info : Done meshing 1D (Wall 5.58328e-05s, CPU 7.6e-05s)
Info : Meshing 2D...
Info : 0 dependencies in mesh of source faces
Info : Meshing curve 2
Info : Meshing curve 2 (Extruded)
Info : Meshing curve 3
Info : Meshing curve 3 (Extruded)
Info : Meshing surface 1 (Plane, Frontal-Delaunay)
Info : Meshing surface 5 (Extruded)
Info : Meshing surface 5 (Boundary layer surface)
Info : Meshing surface 5 (Extruded)
Info : Done meshing 2D (Wall 0.00300558s, CPU 0.003076s)
Info : Meshing 3D...
Info : Done meshing 3D (Wall 1.0829e-06s, CPU 1e-06s)
Info : 119 nodes 262 elements
4.2.2.1. 在 HDivT 空间中绘制区域标记#
网格读进来之后, 常常想确认标记有没有落在预期的单元和面上. 这里用一个自由度长在
面上的空间 HDivT (H(div) 的迹空间; 一阶时每条边上是一次多项式, 即每条边两个
自由度) 把子域标记”画”出来.
做法是用 par_loop 在子域 1 (dx(1)) 上给每个单元的面自由度累加 1/2 (INC
表示累加): 两侧单元都属于子域 1 的边被加了两次, 得到 1; 只有一侧属于子域 1 的边
只加一次, 得到 1/2; 完全不挨着子域 1 的边保持 0. 再把小于 3/4 的值清零, 剩下
取值为 1 的就正好是子域 1 的内部边. 本例实测 656 个自由度中, 334 个为 1,
46 个为 1/2, 276 个为 0.
后面的绘图部分没有用 Firedrake 的画图函数, 而是直接从 mesh.topology_dm 取出
每条边的两个端点坐标画线段, 再通过 HDivT 空间的 Section 找到这条边对应的
自由度, 把标记值写在边的中点上. 输出文件名带上进程数和进程号, 并行运行时每个
进程各画自己那一份.
from firedrake import *
import numpy as np
mesh = Mesh('gmsh/disk_pml.msh')
V = FunctionSpace(mesh, 'HDivT', 1)
marker = Function(V, name='f')
par_loop(('{[i] : 0 <= i < f.dofs}', 'f[i, 0] = 1/2'), dx(1), {'f': (marker, INC)})
index = marker.dat.data_with_halos < 3/4
marker.dat.data_with_halos[index] = 0
plex = mesh.topology_dm
s, e = plex.getHeightStratum(1)
coords = mesh.coordinates
csec = coords.function_space().dm.getSection()
sec = V.dm.getSection()
import matplotlib.pyplot as plt
plt.figure(figsize=[8, 8])
for i in range(s, e):
a, b = plex.getCone(i)
off_a = csec.getOffset(a)
off_b = csec.getOffset(b)
x1, y1 = coords.dat.data_ro_with_halos[off_a].real
x2, y2 = coords.dat.data_ro_with_halos[off_b].real
plt.plot([x1, x2], [y1, y2])
off_i = sec.getOffset(i)
v = marker.dat.data_ro_with_halos[off_i].real
plt.text((x1 + x2)/2, (y1 + y2)/2, round(v), ha='center', va='center')
plt.axis('equal')
rank, size = mesh.comm.rank, mesh.comm.size
plt.savefig(f'figures/hdivt-marker-{size}-{rank}.pdf')
4.3. Gmsh 单纯形节点顺序#
高阶网格的一个单元上有多个节点, 而不同软件对这些节点的排列约定并不一致: Gmsh 先列顶点, 再列各条边上的节点, 然后才是面内、体内的节点; 有限元库里则更常用按 格点坐标字典序排列的编号. 要把 Gmsh 的高阶网格读进来, 就得在两套编号之间建立 对应关系.
下面这组函数就是干这件事的. GmshLexOrder_SEG、GmshLexOrder_TRI、
GmshLexOrder_TET 返回一个二元组 (lex, node): lex 是数组, lex[k] 表示
字典序第 k 个格点在 Gmsh 编号里排第几; node 是递归时用的下一个可用编号.
辅助函数分三类: SN1/SN2/SN3 给出 \(p\) 阶线段、三角形、四面体上的节点总数;
SI1/SI2/SI3 把格点坐标映射成字典序下标; SL1/SL2/SL3 枚举内部格点.
三个 GmshLexOrder_* 按 Gmsh 的约定先编顶点、再编各条边, 最后把面内和体内的
节点递归交给降阶后的同名函数处理.
例如 GmshLexOrder_TRI(2)[0] 是 [0., 3., 1., 5., 4., 2.]: 二阶三角形上字典序
的 6 个格点, 依次对应 Gmsh 的第 0、3、1、5、4、2 号节点. 注意取的是返回值的
第一个分量, 直接写 lex = GmshLexOrder_TRI(2) 拿到的是整个二元组.
import numpy as np
def SN1(p):
return p + 1
def SN2(p):
return SN1(p) * SN1(p + 1) // 2
def SN3(p):
return SN2(p) * SN1(p + 2) // 3
def SI1(p, i):
return i
def SI2(p, i, j):
return i + (SN2(p) - SN2(p - j))
def SI3(p, i, j, k):
return SI2(p - k, i, j) + SN3(p) - SN3(p - k)
def SL1(p):
for i in range(1, p):
yield i
def SL2(p):
for i in range(1, p - 1):
for j in range(1, p - i):
yield i, j
def SL3(p):
for i in range(1, p - 2):
for j in range(1, p - i):
for k in range(1, p - i - j):
yield i, j, k
def GmshLexOrder_SEG(p, node=0):
index = lambda i: SI1(p, i)
lex = - np.ones(SN1(p))
if p == 0:
lex[0] = node; node += 1
return lex, node
lex[index(0)] = node; node += 1
lex[index(p)] = node; node += 1
if p == 1:
return lex, node
for i in SL1(p):
lex[index(i)] = node; node += 1
return lex, node
def GmshLexOrder_TRI(p, node=0):
index = lambda i, j: SI2(p, i, j)
lex = - np.ones(SN2(p))
if p == 0:
lex[0] = node; node += 1
return lex, node
lex[index(0, 0)] = node; node += 1
lex[index(p, 0)] = node; node += 1
lex[index(0, p)] = node; node += 1
if p == 1:
return lex, node
for i in SL1(p):
lex[index(i, 0)] = node; node += 1
for j in SL1(p):
lex[index(p - j, j)] = node; node += 1
for j in SL1(p):
lex[index(0, p - j)] = node; node += 1
if p == 2:
return lex, node
sub, node = GmshLexOrder_TRI(p - 3, node);
for _, (j, i) in enumerate(SL2(p)):
lex[index(i, j)] = sub[_];
return lex, node
def GmshLexOrder_TET(p, node=0):
index = lambda i, j, k: SI3(p, i, j, k)
lex = - np.ones(SN3(p))
if p == 0:
lex[0] = node; node += 1
return lex, node
lex[index(0, 0, 0)] = node; node += 1
lex[index(p, 0, 0)] = node; node += 1
lex[index(0, p, 0)] = node; node += 1
lex[index(0, 0, p)] = node; node += 1
if p == 1:
return lex, node
# internal edge nodes
for i in SL1(p): lex[index( i, 0, 0)] = node; node += 1
for j in SL1(p): lex[index(p - j, j, 0)] = node; node += 1
for j in SL1(p): lex[index( 0, p - j, 0)] = node; node += 1
for k in SL1(p): lex[index( 0, 0, p - k)] = node; node += 1
for j in SL1(p): lex[index( 0, j, p - j)] = node; node += 1
for i in SL1(p): lex[index( i, 0, p - i)] = node; node += 1
if p == 2:
return lex, node
# /* internal face nodes */
sub, node = GmshLexOrder_TRI(p - 3, node)
for _, (i, j) in enumerate(SL2(p)):
lex[index(i, j, 0)] = sub[_]
sub, node = GmshLexOrder_TRI(p - 3, node);
for _, (k, i) in enumerate(SL2(p)):
lex[index(i, 0, k)] = sub[_]
sub, node = GmshLexOrder_TRI(p - 3, node);
for _, (j, k) in enumerate(SL2(p)):
lex[index(0, j, k)] = sub[_]
sub, node = GmshLexOrder_TRI(p - 3, node);
for _, (j, i) in enumerate(SL2(p)):
lex[index(i, j, p - i - j)] = sub[_]
if p == 3:
return lex, node
# internal cell nodes */
sub, node = GmshLexOrder_TET(p - 4, node);
for _, (k, j, i) in enumerate(SL3(p)):
lex[index(i, j, k)] = sub[_];
return lex, node
4.3.1. 高阶单元测试#
Firedrake 读 .msh 文件时只会给网格配一个一阶的坐标单元 (见 firedrake/mesh.py
的 make_mesh_from_mesh_topology, 那里至今留着一条
“meshfile might indicates higher-order coordinate element” 的 TODO). 因此直接读
一个二阶网格会失败, 在 Firedrake 2026.4.1 下报的是坐标数组形状对不上:
>>> Mesh('gmsh/cube_2rd.msh')
ValueError: cannot reshape array of size 720 into shape (14,3)
下面几个单元演示如何绕过这一步: 先从 DMPlex 上读出坐标单元的阶数, 再自己拼一个
相应阶数的坐标单元, 最后把 Firedrake 的 make_mesh_from_mesh_topology 换掉.
已知问题
这段演示只解决了”读得进来”的问题, 读进来的几何是不对的. 用
-dm_plex_gmsh_project 让 PETSc 把逐单元存放的高阶坐标投影成连续坐标之后, 坐标
向量里的数值并不是各节点的坐标. 在 Firedrake 2026.4.1 + PETSc 3.25 下实测, 读
单位立方体网格 gmsh/cube_2rd.msh 得到的坐标三个分量最大值分别是 9、7 和 8.5
(本应都是 1), 由此建出来的网格体积约为 75.41 (本应是 1); 具体输出见下面打印坐标
范围的那个单元. 因此本节只是”问题出在哪里”的记录, 还不是一个可用的高阶网格读取
方案.
from firedrake import *
import firedrake as fd
import numpy as np
try:
from petsctools.options import OptionsManager
except ImportError:
from firedrake.petsc import OptionsManager
下面两个辅助函数各管一件事: getCoordinateFESpaceOrder 从 DMPlex 的坐标 DM
上取出 PETSc 给坐标场配的有限元名字 (形如 P2), 从中解析出阶数, 取不到就按
一阶处理; coordinates_from_plex 把 DMPlex 坐标向量里的值按 Firedrake 的自由度
编号重排, 装进一个 CoordinatelessFunction.
coordinates_from_plex 里的 assert 是有意加的: 它每个网格点只搬运一个节点的
坐标, 因此只适用于单纯形上最高二阶的情形. 三阶网格每条边有 2 个节点, 若少了这条
assert, 多出来的那些行会悄悄留成零坐标 (实测读 gmsh/cube_3rd.msh 时
172 行里有 49 行全零), 却不报任何错.
这一段照着 firedrake/mesh.py 的 make_mesh_from_mesh_topology 改写, 与原版有
两处区别:
坐标单元的阶数改成从网格文件里读出来的, 而不是固定为一阶;
原版里处理周期网格的那个分支被丢掉了: 原版会先看
dm.getCoordinatesLocalized(), 对周期网格调用_fully_localize_coordinates并改用逐单元的DG/DQ一阶坐标单元 (variant="equispaced"), 这里没有照搬.
改写后的函数被塞回 firedrake.mesh 模块. Mesh() 内部是按模块级名字调用它的,
因此替换掉模块属性就能生效 —— 但这也意味着替换之后, 本进程里所有的网格构造都会
走这段代码, 包括与本节毫无关系的那些.
这段代码是实验性的, 而且会改动 Firedrake 的内部状态
它有三条硬边界, 越界时大多不会报错, 只会悄悄给出错误的网格:
不支持周期网格. 上面第 2 条丢掉的正是周期网格分支. 本地实测, 装上这个 替换之后
PeriodicUnitSquareMesh(4, 4)的面积从 1.0 变成 2.25,PeriodicIntervalMesh(8, 1.0)的长度从 1.0 变成 1.75, 全程不报任何错; 非周期网格不受影响 (UnitSquareMesh(4, 4)仍是 1.0).只支持每个网格点带一个坐标节点的情形, 即单纯形上最高二阶. 这一条由
coordinates_from_plex里的assert兜住, 会明确报错.必须配合
-dm_plex_gmsh_project和petscdualspace_lagrange_continuity一起使用, 否则报ValueError: could not broadcast input array from shape (0,) into shape (3,).
因此用完一定要把原函数还原 (见下面 triplot 之后的那个还原单元), 否则同一次会话
里后面所有的网格都会受影响.
-dm_plex_gmsh_project 触发 PETSc 在读入网格后重新投影坐标, 但它默认的目标
空间仍然是不连续的: 实测只给这一个选项时, 坐标 Section 的存储量还是 720、
自由度仍旧挂在单元上, 与不给选项没有区别. 必须再加上
petscdualspace_lagrange_continuity, 坐标才会变成挂在顶点和边上的连续 P2
(存储量 189, 即 (14 个顶点 + 49 条边) \(\times\) 3).
这两个都是 PETSc 的命令行选项, 用 OptionsManager 的 inserted_options() 临时
塞进 PETSc 的选项数据库, 退出 with 块后自动移除. 对上面那段替换过的
make_mesh_from_mesh_topology 来说这两个选项是必需的, 少了就会报
ValueError: could not broadcast input array from shape (0,) into shape (3,).
project_parameter = {
"dm_plex_gmsh_project": True,
"dm_plex_gmsh_project_": {
# "fe_view": "ascii",
"petscdualspace_lagrange_continuity": True,
}
}
om = OptionsManager(project_parameter, options_prefix="")
with om.inserted_options():
mesh = Mesh('gmsh/cube_2rd.msh')
下面把读进来的网格画出来, 同时打印坐标范围和网格体积. gmsh/cube_2rd.msh 是
单位立方体, 坐标本应落在 \([0, 1]^3\) 内、体积本应等于 1, 实际输出与此相差很远,
这正是上面”已知问题”里说的那种情况.
xs = mesh.coordinates.dat.data_ro.real
print("coordinate range:", xs.min(axis=0), xs.max(axis=0))
print("mesh volume:", assemble(Constant(1.0)*dx(domain=mesh)))
c = triplot(mesh)
coordinate range: [-3.73461584e-17 -1.11457691e-16 -3.69532952e-17] [9. 7. 8.5]
mesh volume: 75.41047563906
实验做完, 立刻把 Firedrake 原本的 make_mesh_from_mesh_topology 还原回去.
若不还原, 同一次会话里后面建的每一张网格都会走上面那段实验代码, 周期网格会被
静默算错.
fd.mesh.make_mesh_from_mesh_topology = _orig_make_mesh_from_mesh_topology
4.3.1.1. Gmsh 网格读取的问题#
这一小节绕开 Firedrake, 直接用 PETSc.DMPlex 读同一个文件, 对比”投影前”和
“投影后”两种坐标存放方式, 看清上面那个问题出在哪一步.
show_plex_coordinates 打印投影前的坐标: 此时高阶坐标是逐单元存放的, 坐标
Section 的自由度挂在单元上, 每个四面体存 10 个节点的坐标, 数值都落在 \([0, 1]\)
内, 是对的.
show_plex_coordinates_continue 打印投影后的坐标: 此时自由度改挂到了顶点和边上
(二阶四面体每个顶点、每条边各一个节点), 所以要沿单元的传递闭包
(getTransitiveClosure) 把这些点找出来再取值. 打印出的数值明显超出了 \([0, 1]\),
可见问题出在这一步投影上, 而不是出在前面替换 make_mesh_from_mesh_topology 的
那段代码里.
from firedrake.petsc import PETSc
import numpy as np
try:
from petsctools.options import OptionsManager
except ImportError:
from firedrake.petsc import OptionsManager
def show_plex_coordinates_continue(plex):
dim = plex.getCoordinateDim()
dm_sec = plex.getCoordinateSection()
dm_coords = plex.getCoordinatesLocal().array_r.real.reshape([-1, dim])
cs, ce = plex.getHeightStratum(0)
print(f"plex {plex.getName()}")
for i in range(cs, ce):
fs = plex.getCone(i)
print(" cell:", i, tuple(fs))
cl, o = plex.getTransitiveClosure(i)
print(" coordinates:")
for p in cl:
if dm_sec.getDof(p) == 0: # 只有带坐标自由度的点才需要打印
continue
offset = dm_sec.getOffset(p)
print(" ", p, "[%.2f, %.2f, %.2f]"%tuple(dm_coords[offset//dim, :]))
def show_plex_coordinates(plex):
dim = plex.getCoordinateDim()
dm_sec = plex.getCoordinateSection()
dm_coords = plex.getCoordinatesLocal().array_r.real
cs, ce = plex.getHeightStratum(0)
print(f"plex {plex.getName()}")
for i in range(cs, ce):
offset = dm_sec.getOffset(i)
dof = dm_sec.getDof(i)
fs = plex.getCone(i)
print(" cell:", i, tuple(fs))
print(" coordinates:")
coords = dm_coords[offset:offset+dof].reshape([-1, dim])
for k in coords:
print(" ", "[%.2f, %.2f, %.2f]"%tuple(k))
plex = PETSc.DMPlex().createFromFile('gmsh/cube_2rd.msh', 'gmsh_2rd')
show_plex_coordinates(plex)
project_parameter = {
"dm_plex_gmsh_project": True,
"dm_plex_gmsh_project_": {
# "fe_view": "ascii",
"petscdualspace_lagrange_continuity": True,
}
}
om = OptionsManager(project_parameter, options_prefix="")
with om.inserted_options():
plex_proj = PETSc.DMPlex().createFromFile('gmsh/cube_2rd.msh', 'gmsh_2rd_with_projection')
show_plex_coordinates_continue(plex_proj)
plex gmsh_2rd
cell: 0 (np.int32(38), np.int32(39), np.int32(40), np.int32(41))
coordinates:
[0.50, 1.00, 0.50]
[0.25, 0.75, 0.50]
[0.00, 0.50, 0.50]
[0.50, 0.50, 0.50]
[0.25, 0.25, 0.50]
[0.50, 0.00, 0.50]
[0.50, 0.75, 0.75]
[0.25, 0.50, 0.75]
[0.50, 0.25, 0.75]
[0.50, 0.50, 1.00]
cell: 1 (np.int32(42), np.int32(43), np.int32(44), np.int32(38))
coordinates:
[0.00, 0.50, 0.50]
[0.25, 0.50, 0.25]
[0.50, 0.50, 0.00]
[0.25, 0.75, 0.50]
[0.50, 0.75, 0.25]
[0.50, 1.00, 0.50]
[0.25, 0.25, 0.50]
[0.50, 0.25, 0.25]
[0.50, 0.50, 0.50]
[0.50, 0.00, 0.50]
cell: 2 (np.int32(44), np.int32(45), np.int32(46), np.int32(47))
coordinates:
[0.50, 0.50, 0.00]
[0.50, 0.75, 0.25]
[0.50, 1.00, 0.50]
[0.50, 0.25, 0.25]
[0.50, 0.50, 0.50]
[0.50, 0.00, 0.50]
[0.75, 0.50, 0.25]
[0.75, 0.75, 0.50]
[0.75, 0.25, 0.50]
[1.00, 0.50, 0.50]
cell: 3 (np.int32(41), np.int32(48), np.int32(46), np.int32(49))
coordinates:
[0.50, 0.50, 1.00]
[0.50, 0.25, 0.75]
[0.50, 0.00, 0.50]
[0.50, 0.75, 0.75]
[0.50, 0.50, 0.50]
[0.50, 1.00, 0.50]
[0.75, 0.50, 0.75]
[0.75, 0.25, 0.50]
[0.75, 0.75, 0.50]
[1.00, 0.50, 0.50]
cell: 4 (np.int32(50), np.int32(51), np.int32(52), np.int32(53))
coordinates:
[0.50, 0.50, 0.00]
[0.25, 0.75, 0.00]
[0.00, 1.00, 0.00]
[0.50, 0.75, 0.25]
[0.25, 1.00, 0.25]
[0.50, 1.00, 0.50]
[0.75, 0.75, 0.00]
[0.50, 1.00, 0.00]
[0.75, 1.00, 0.25]
[1.00, 1.00, 0.00]
cell: 5 (np.int32(54), np.int32(55), np.int32(56), np.int32(57))
coordinates:
[0.50, 0.50, 0.00]
[0.25, 0.25, 0.00]
[0.00, 0.00, 0.00]
[0.25, 0.50, 0.25]
[0.00, 0.25, 0.25]
[0.00, 0.50, 0.50]
[0.25, 0.75, 0.00]
[0.00, 0.50, 0.00]
[0.00, 0.75, 0.25]
[0.00, 1.00, 0.00]
cell: 6 (np.int32(58), np.int32(59), np.int32(60), np.int32(61))
coordinates:
[0.50, 1.00, 0.50]
[0.25, 1.00, 0.25]
[0.00, 1.00, 0.00]
[0.25, 0.75, 0.50]
[0.00, 0.75, 0.25]
[0.00, 0.50, 0.50]
[0.25, 1.00, 0.75]
[0.00, 1.00, 0.50]
[0.00, 0.75, 0.75]
[0.00, 1.00, 1.00]
cell: 7 (np.int32(62), np.int32(63), np.int32(64), np.int32(65))
coordinates:
[0.00, 0.00, 1.00]
[0.00, 0.25, 0.75]
[0.00, 0.50, 0.50]
[0.00, 0.50, 1.00]
[0.00, 0.75, 0.75]
[0.00, 1.00, 1.00]
[0.25, 0.25, 1.00]
[0.25, 0.50, 0.75]
[0.25, 0.75, 1.00]
[0.50, 0.50, 1.00]
cell: 8 (np.int32(66), np.int32(67), np.int32(68), np.int32(69))
coordinates:
[0.00, 0.50, 0.50]
[0.00, 0.25, 0.25]
[0.00, 0.00, 0.00]
[0.25, 0.25, 0.50]
[0.25, 0.00, 0.25]
[0.50, 0.00, 0.50]
[0.00, 0.25, 0.75]
[0.00, 0.00, 0.50]
[0.25, 0.00, 0.75]
[0.00, 0.00, 1.00]
cell: 9 (np.int32(70), np.int32(71), np.int32(72), np.int32(73))
coordinates:
[0.00, 0.00, 0.00]
[0.25, 0.25, 0.00]
[0.50, 0.50, 0.00]
[0.25, 0.00, 0.25]
[0.50, 0.25, 0.25]
[0.50, 0.00, 0.50]
[0.50, 0.00, 0.00]
[0.75, 0.25, 0.00]
[0.75, 0.00, 0.25]
[1.00, 0.00, 0.00]
cell: 10 (np.int32(74), np.int32(75), np.int32(76), np.int32(77))
coordinates:
[0.50, 0.50, 1.00]
[0.25, 0.25, 1.00]
[0.00, 0.00, 1.00]
[0.50, 0.25, 0.75]
[0.25, 0.00, 0.75]
[0.50, 0.00, 0.50]
[0.75, 0.25, 1.00]
[0.50, 0.00, 1.00]
[0.75, 0.00, 0.75]
[1.00, 0.00, 1.00]
cell: 11 (np.int32(78), np.int32(79), np.int32(80), np.int32(81))
coordinates:
[0.50, 1.00, 0.50]
[0.25, 1.00, 0.75]
[0.00, 1.00, 1.00]
[0.50, 0.75, 0.75]
[0.25, 0.75, 1.00]
[0.50, 0.50, 1.00]
[0.75, 1.00, 0.75]
[0.50, 1.00, 1.00]
[0.75, 0.75, 1.00]
[1.00, 1.00, 1.00]
cell: 12 (np.int32(82), np.int32(83), np.int32(84), np.int32(85))
coordinates:
[1.00, 1.00, 1.00]
[1.00, 0.50, 1.00]
[1.00, 0.00, 1.00]
[0.75, 0.75, 1.00]
[0.75, 0.25, 1.00]
[0.50, 0.50, 1.00]
[1.00, 0.75, 0.75]
[1.00, 0.25, 0.75]
[0.75, 0.50, 0.75]
[1.00, 0.50, 0.50]
cell: 13 (np.int32(86), np.int32(87), np.int32(88), np.int32(89))
coordinates:
[1.00, 1.00, 0.00]
[1.00, 1.00, 0.50]
[1.00, 1.00, 1.00]
[0.75, 1.00, 0.25]
[0.75, 1.00, 0.75]
[0.50, 1.00, 0.50]
[1.00, 0.75, 0.25]
[1.00, 0.75, 0.75]
[0.75, 0.75, 0.50]
[1.00, 0.50, 0.50]
cell: 14 (np.int32(90), np.int32(91), np.int32(92), np.int32(93))
coordinates:
[1.00, 0.50, 0.50]
[1.00, 0.25, 0.75]
[1.00, 0.00, 1.00]
[0.75, 0.25, 0.50]
[0.75, 0.00, 0.75]
[0.50, 0.00, 0.50]
[1.00, 0.25, 0.25]
[1.00, 0.00, 0.50]
[0.75, 0.00, 0.25]
[1.00, 0.00, 0.00]
cell: 15 (np.int32(94), np.int32(95), np.int32(96), np.int32(97))
coordinates:
[1.00, 0.50, 0.50]
[0.75, 0.50, 0.25]
[0.50, 0.50, 0.00]
[1.00, 0.75, 0.25]
[0.75, 0.75, 0.00]
[1.00, 1.00, 0.00]
[1.00, 0.25, 0.25]
[0.75, 0.25, 0.00]
[1.00, 0.50, 0.00]
[1.00, 0.00, 0.00]
cell: 16 (np.int32(69), np.int32(74), np.int32(40), np.int32(63))
coordinates:
[0.00, 0.00, 1.00]
[0.25, 0.00, 0.75]
[0.50, 0.00, 0.50]
[0.00, 0.25, 0.75]
[0.25, 0.25, 0.50]
[0.00, 0.50, 0.50]
[0.25, 0.25, 1.00]
[0.50, 0.25, 0.75]
[0.25, 0.50, 0.75]
[0.50, 0.50, 1.00]
cell: 17 (np.int32(61), np.int32(39), np.int32(78), np.int32(64))
coordinates:
[0.00, 0.50, 0.50]
[0.25, 0.75, 0.50]
[0.50, 1.00, 0.50]
[0.00, 0.75, 0.75]
[0.25, 1.00, 0.75]
[0.00, 1.00, 1.00]
[0.25, 0.50, 0.75]
[0.50, 0.75, 0.75]
[0.25, 0.75, 1.00]
[0.50, 0.50, 1.00]
cell: 18 (np.int32(57), np.int32(42), np.int32(50), np.int32(58))
coordinates:
[0.00, 0.50, 0.50]
[0.25, 0.50, 0.25]
[0.50, 0.50, 0.00]
[0.00, 0.75, 0.25]
[0.25, 0.75, 0.00]
[0.00, 1.00, 0.00]
[0.25, 0.75, 0.50]
[0.50, 0.75, 0.25]
[0.25, 1.00, 0.25]
[0.50, 1.00, 0.50]
cell: 19 (np.int32(54), np.int32(66), np.int32(70), np.int32(43))
coordinates:
[0.00, 0.50, 0.50]
[0.00, 0.25, 0.25]
[0.00, 0.00, 0.00]
[0.25, 0.50, 0.25]
[0.25, 0.25, 0.00]
[0.50, 0.50, 0.00]
[0.25, 0.25, 0.50]
[0.25, 0.00, 0.25]
[0.50, 0.25, 0.25]
[0.50, 0.00, 0.50]
cell: 20 (np.int32(77), np.int32(48), np.int32(84), np.int32(90))
coordinates:
[0.50, 0.00, 0.50]
[0.50, 0.25, 0.75]
[0.50, 0.50, 1.00]
[0.75, 0.00, 0.75]
[0.75, 0.25, 1.00]
[1.00, 0.00, 1.00]
[0.75, 0.25, 0.50]
[0.75, 0.50, 0.75]
[1.00, 0.25, 0.75]
[1.00, 0.50, 0.50]
cell: 21 (np.int32(81), np.int32(85), np.int32(49), np.int32(88))
coordinates:
[1.00, 1.00, 1.00]
[0.75, 0.75, 1.00]
[0.50, 0.50, 1.00]
[0.75, 1.00, 0.75]
[0.50, 0.75, 0.75]
[0.50, 1.00, 0.50]
[1.00, 0.75, 0.75]
[0.75, 0.50, 0.75]
[0.75, 0.75, 0.50]
[1.00, 0.50, 0.50]
cell: 22 (np.int32(53), np.int32(45), np.int32(94), np.int32(89))
coordinates:
[0.50, 1.00, 0.50]
[0.50, 0.75, 0.25]
[0.50, 0.50, 0.00]
[0.75, 1.00, 0.25]
[0.75, 0.75, 0.00]
[1.00, 1.00, 0.00]
[0.75, 0.75, 0.50]
[0.75, 0.50, 0.25]
[1.00, 0.75, 0.25]
[1.00, 0.50, 0.50]
cell: 23 (np.int32(47), np.int32(93), np.int32(72), np.int32(95))
coordinates:
[1.00, 0.50, 0.50]
[0.75, 0.25, 0.50]
[0.50, 0.00, 0.50]
[0.75, 0.50, 0.25]
[0.50, 0.25, 0.25]
[0.50, 0.50, 0.00]
[1.00, 0.25, 0.25]
[0.75, 0.00, 0.25]
[0.75, 0.25, 0.00]
[1.00, 0.00, 0.00]
plex gmsh_2rd_with_projection
cell: 0 (np.int32(38), np.int32(39), np.int32(40), np.int32(41))
coordinates:
98 [1.25, 3.75, 2.00]
99 [1.75, 1.50, 1.75]
100 [1.75, 1.50, 2.25]
101 [1.50, 2.00, 3.50]
102 [2.25, 2.75, 4.00]
103 [2.00, 1.50, 4.00]
32 [2.50, 5.50, 4.00]
35 [5.50, 7.00, 5.50]
34 [5.50, 3.00, 7.50]
37 [5.50, 5.00, 8.50]
cell: 1 (np.int32(42), np.int32(43), np.int32(44), np.int32(38))
coordinates:
104 [0.50, 2.00, 1.00]
98 [1.25, 3.75, 2.00]
105 [2.00, 3.50, 1.50]
106 [2.00, 1.25, 1.75]
100 [1.75, 1.50, 2.25]
99 [1.75, 1.50, 1.75]
36 [3.00, 4.00, 4.00]
32 [2.50, 5.50, 4.00]
35 [5.50, 7.00, 5.50]
34 [5.50, 3.00, 7.50]
cell: 2 (np.int32(44), np.int32(45), np.int32(46), np.int32(47))
coordinates:
105 [2.00, 3.50, 1.50]
106 [2.00, 1.25, 1.75]
99 [1.75, 1.50, 1.75]
107 [4.25, 2.75, 2.50]
108 [3.00, 2.00, 1.50]
109 [4.00, 2.00, 3.00]
35 [5.50, 7.00, 5.50]
36 [3.00, 4.00, 4.00]
34 [5.50, 3.00, 7.50]
33 [9.00, 5.00, 4.00]
cell: 3 (np.int32(41), np.int32(48), np.int32(46), np.int32(49))
coordinates:
103 [2.00, 1.50, 4.00]
102 [2.25, 2.75, 4.00]
99 [1.75, 1.50, 1.75]
109 [4.00, 2.00, 3.00]
110 [3.50, 2.00, 2.50]
107 [4.25, 2.75, 2.50]
34 [5.50, 3.00, 7.50]
37 [5.50, 5.00, 8.50]
35 [5.50, 7.00, 5.50]
33 [9.00, 5.00, 4.00]
cell: 4 (np.int32(50), np.int32(51), np.int32(52), np.int32(53))
coordinates:
111 [0.75, 2.00, 0.75]
105 [2.00, 3.50, 1.50]
112 [1.00, 2.25, 0.75]
113 [0.75, 0.75, -0.00]
114 [2.00, 2.25, 0.75]
115 [2.00, 2.75, 0.75]
27 [1.50, 3.00, 0.50]
36 [3.00, 4.00, 4.00]
35 [5.50, 7.00, 5.50]
31 [2.50, 3.00, 0.50]
cell: 5 (np.int32(54), np.int32(55), np.int32(56), np.int32(57))
coordinates:
116 [0.50, 0.75, 0.75]
104 [0.50, 2.00, 1.00]
117 [0.75, 1.00, 0.75]
118 [0.25, 0.75, -0.00]
111 [0.75, 2.00, 0.75]
119 [0.50, 2.25, 0.25]
25 [1.00, 1.50, 1.50]
36 [3.00, 4.00, 4.00]
32 [2.50, 5.50, 4.00]
27 [1.50, 3.00, 0.50]
cell: 6 (np.int32(58), np.int32(59), np.int32(60), np.int32(61))
coordinates:
112 [1.00, 2.25, 0.75]
98 [1.25, 3.75, 2.00]
119 [0.50, 2.25, 0.25]
120 [0.25, 1.00, 0.75]
121 [0.75, 2.25, 2.00]
122 [0.25, 2.25, 2.00]
27 [1.50, 3.00, 0.50]
35 [5.50, 7.00, 5.50]
32 [2.50, 5.50, 4.00]
26 [1.00, 3.50, 2.50]
cell: 7 (np.int32(62), np.int32(63), np.int32(64), np.int32(65))
coordinates:
123 [0.50, 0.75, 2.25]
124 [-0.00, 0.75, 0.75]
122 [0.25, 2.25, 2.00]
101 [1.50, 2.00, 3.50]
125 [1.00, 1.50, 2.50]
126 [1.00, 2.25, 2.25]
32 [2.50, 5.50, 4.00]
24 [0.50, 2.00, 3.50]
26 [1.00, 3.50, 2.50]
37 [5.50, 5.00, 8.50]
cell: 8 (np.int32(66), np.int32(67), np.int32(68), np.int32(69))
coordinates:
117 [0.75, 1.00, 0.75]
100 [1.75, 1.50, 2.25]
127 [0.75, 0.75, 1.00]
128 [-0.00, 0.25, 0.75]
123 [0.50, 0.75, 2.25]
129 [0.25, 0.50, 2.25]
25 [1.00, 1.50, 1.50]
32 [2.50, 5.50, 4.00]
34 [5.50, 3.00, 7.50]
24 [0.50, 2.00, 3.50]
cell: 9 (np.int32(70), np.int32(71), np.int32(72), np.int32(73))
coordinates:
116 [0.50, 0.75, 0.75]
127 [0.75, 0.75, 1.00]
106 [2.00, 1.25, 1.75]
130 [2.25, 0.25, 0.50]
131 [0.75, 0.00, 0.25]
132 [2.75, 0.50, 0.75]
36 [3.00, 4.00, 4.00]
25 [1.00, 1.50, 1.50]
34 [5.50, 3.00, 7.50]
29 [4.00, -0.00, 0.00]
cell: 10 (np.int32(74), np.int32(75), np.int32(76), np.int32(77))
coordinates:
125 [1.00, 1.50, 2.50]
103 [2.00, 1.50, 4.00]
129 [0.25, 0.50, 2.25]
133 [0.75, 0.25, 1.00]
134 [2.25, 0.75, 2.50]
135 [2.25, 0.50, 2.75]
24 [0.50, 2.00, 3.50]
37 [5.50, 5.00, 8.50]
34 [5.50, 3.00, 7.50]
28 [3.50, 2.00, 3.50]
cell: 11 (np.int32(78), np.int32(79), np.int32(80), np.int32(81))
coordinates:
121 [0.75, 2.25, 2.00]
102 [2.25, 2.75, 4.00]
126 [1.00, 2.25, 2.25]
136 [0.75, 1.00, 0.75]
137 [2.25, 2.50, 2.25]
138 [2.00, 2.25, 2.75]
26 [1.00, 3.50, 2.50]
35 [5.50, 7.00, 5.50]
37 [5.50, 5.00, 8.50]
30 [3.00, 3.50, 2.50]
cell: 12 (np.int32(82), np.int32(83), np.int32(84), np.int32(85))
coordinates:
139 [0.75, 0.75, 1.00]
138 [2.00, 2.25, 2.75]
134 [2.25, 0.75, 2.50]
140 [2.50, 1.50, 2.00]
141 [2.50, 2.00, 1.50]
110 [3.50, 2.00, 2.50]
28 [3.50, 2.00, 3.50]
30 [3.00, 3.50, 2.50]
37 [5.50, 5.00, 8.50]
33 [9.00, 5.00, 4.00]
cell: 13 (np.int32(86), np.int32(87), np.int32(88), np.int32(89))
coordinates:
142 [0.75, 1.00, 0.25]
115 [2.00, 2.75, 0.75]
137 [2.25, 2.50, 2.25]
141 [2.50, 2.00, 1.50]
143 [2.25, 2.00, 0.75]
107 [4.25, 2.75, 2.50]
30 [3.00, 3.50, 2.50]
31 [2.50, 3.00, 0.50]
35 [5.50, 7.00, 5.50]
33 [9.00, 5.00, 4.00]
cell: 14 (np.int32(90), np.int32(91), np.int32(92), np.int32(93))
coordinates:
140 [2.50, 1.50, 2.00]
109 [4.00, 2.00, 3.00]
135 [2.25, 0.50, 2.75]
144 [1.00, 0.25, 0.25]
145 [2.50, 0.75, 0.25]
132 [2.75, 0.50, 0.75]
28 [3.50, 2.00, 3.50]
33 [9.00, 5.00, 4.00]
34 [5.50, 3.00, 7.50]
29 [4.00, -0.00, 0.00]
cell: 15 (np.int32(94), np.int32(95), np.int32(96), np.int32(97))
coordinates:
108 [3.00, 2.00, 1.50]
143 [2.25, 2.00, 0.75]
114 [2.00, 2.25, 0.75]
130 [2.25, 0.25, 0.50]
145 [2.50, 0.75, 0.25]
146 [0.75, 0.25, 0.00]
36 [3.00, 4.00, 4.00]
33 [9.00, 5.00, 4.00]
31 [2.50, 3.00, 0.50]
29 [4.00, -0.00, 0.00]
cell: 16 (np.int32(69), np.int32(74), np.int32(40), np.int32(63))
coordinates:
129 [0.25, 0.50, 2.25]
123 [0.50, 0.75, 2.25]
100 [1.75, 1.50, 2.25]
103 [2.00, 1.50, 4.00]
125 [1.00, 1.50, 2.50]
101 [1.50, 2.00, 3.50]
34 [5.50, 3.00, 7.50]
24 [0.50, 2.00, 3.50]
32 [2.50, 5.50, 4.00]
37 [5.50, 5.00, 8.50]
cell: 17 (np.int32(61), np.int32(39), np.int32(78), np.int32(64))
coordinates:
98 [1.25, 3.75, 2.00]
122 [0.25, 2.25, 2.00]
121 [0.75, 2.25, 2.00]
102 [2.25, 2.75, 4.00]
101 [1.50, 2.00, 3.50]
126 [1.00, 2.25, 2.25]
35 [5.50, 7.00, 5.50]
32 [2.50, 5.50, 4.00]
26 [1.00, 3.50, 2.50]
37 [5.50, 5.00, 8.50]
cell: 18 (np.int32(57), np.int32(42), np.int32(50), np.int32(58))
coordinates:
104 [0.50, 2.00, 1.00]
119 [0.50, 2.25, 0.25]
111 [0.75, 2.00, 0.75]
105 [2.00, 3.50, 1.50]
98 [1.25, 3.75, 2.00]
112 [1.00, 2.25, 0.75]
36 [3.00, 4.00, 4.00]
32 [2.50, 5.50, 4.00]
27 [1.50, 3.00, 0.50]
35 [5.50, 7.00, 5.50]
cell: 19 (np.int32(54), np.int32(66), np.int32(70), np.int32(43))
coordinates:
117 [0.75, 1.00, 0.75]
104 [0.50, 2.00, 1.00]
116 [0.50, 0.75, 0.75]
127 [0.75, 0.75, 1.00]
100 [1.75, 1.50, 2.25]
106 [2.00, 1.25, 1.75]
25 [1.00, 1.50, 1.50]
32 [2.50, 5.50, 4.00]
36 [3.00, 4.00, 4.00]
34 [5.50, 3.00, 7.50]
cell: 20 (np.int32(77), np.int32(48), np.int32(84), np.int32(90))
coordinates:
103 [2.00, 1.50, 4.00]
135 [2.25, 0.50, 2.75]
134 [2.25, 0.75, 2.50]
110 [3.50, 2.00, 2.50]
109 [4.00, 2.00, 3.00]
140 [2.50, 1.50, 2.00]
37 [5.50, 5.00, 8.50]
34 [5.50, 3.00, 7.50]
28 [3.50, 2.00, 3.50]
33 [9.00, 5.00, 4.00]
cell: 21 (np.int32(81), np.int32(85), np.int32(49), np.int32(88))
coordinates:
138 [2.00, 2.25, 2.75]
137 [2.25, 2.50, 2.25]
102 [2.25, 2.75, 4.00]
110 [3.50, 2.00, 2.50]
141 [2.50, 2.00, 1.50]
107 [4.25, 2.75, 2.50]
37 [5.50, 5.00, 8.50]
30 [3.00, 3.50, 2.50]
35 [5.50, 7.00, 5.50]
33 [9.00, 5.00, 4.00]
cell: 22 (np.int32(53), np.int32(45), np.int32(94), np.int32(89))
coordinates:
105 [2.00, 3.50, 1.50]
115 [2.00, 2.75, 0.75]
114 [2.00, 2.25, 0.75]
108 [3.00, 2.00, 1.50]
107 [4.25, 2.75, 2.50]
143 [2.25, 2.00, 0.75]
36 [3.00, 4.00, 4.00]
35 [5.50, 7.00, 5.50]
31 [2.50, 3.00, 0.50]
33 [9.00, 5.00, 4.00]
cell: 23 (np.int32(47), np.int32(93), np.int32(72), np.int32(95))
coordinates:
109 [4.00, 2.00, 3.00]
108 [3.00, 2.00, 1.50]
106 [2.00, 1.25, 1.75]
132 [2.75, 0.50, 0.75]
145 [2.50, 0.75, 0.25]
130 [2.25, 0.25, 0.50]
34 [5.50, 3.00, 7.50]
33 [9.00, 5.00, 4.00]
36 [3.00, 4.00, 4.00]
29 [4.00, -0.00, 0.00]