4. Gmsh 示例 Gmsh Examples#
概要
用 Gmsh 的 Python API 生成网格: 莫比乌斯带、为 PML 层做的拉伸, 以及高阶单元的
Gmsh 节点顺序. 需要 gmsh 包 (本笔记本在构建时不执行).
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. 莫比乌斯带 Möbius strip#
Firedrake 没有内置的莫比乌斯带网格. 它是嵌在三维空间里的二维曲面, 只能借助外部网格生成器按参数方程一段段拼出来. 本节先用 Gmsh 的 Python API 把这条带子造出来,
再把生成的 .msh 文件交给 Firedrake, 在它上面解一个曲面 Poisson 问题.
Firedrake has no built-in mesh for the Möbius strip. It is a two-dimensional
surface embedded in three-dimensional space and can only be stitched together
segment by segment from its parametric equations with an external mesh
generator. This section first builds the strip with the Python API of Gmsh, and
then hands the resulting .msh file to Firedrake to solve a surface Poisson
problem on it.
4.1.1. 用 matplotlib 绘图 Plotting with 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 取值和下一段脚本里的也不相同.
The surface is first drawn from its parametric equations, to give an idea of the coordinate formulas used in the script of the next section. Take \(u \in [0, 2\pi]\) going around the central circle and \(v \in [-1/2, 1/2]\) across the width of the strip, and let \(r = R - W v \sin(u/2)\); the points on the surface are then \((r\cos u,\ r\sin u,\ W v \cos(u/2))\). Each full turn in \(u\) rotates the cross section by half a turn, so that the two edges of the strip are exchanged when the ends are joined: this is the Möbius strip.
This cell only draws a picture and does not interact with the mesh generation
below, so the values of R and W used here differ from those in the script of
the next section.
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 脚本 Python script#
下面用 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 里不会触发.
The Python API of Gmsh is used below to “stitch” the strip together: it is cut
into M segments along \(u\), each segment carries one point on either side across
the width, and the four points of two neighboring segments bound a curve loop
(add_curve_loop) which is then filled into a plane surface
(add_plane_surface). The key lies in the two closing lines: the last segment is
connected crosswise to segment 0 (points[0][-1] to points[1][0],
points[1][-1] to points[0][0]), and this is how the half turn closes up into
a Möbius strip.
Mesh.MeshSizeMin/Mesh.MeshSizeMax keep the mesh size around h, and
gmsh.model.mesh.generate() is what actually starts the meshing.
General.Verbosity is lowered before meshing to avoid flooding the output; Gmsh
raises an exception when meshing fails, so try/finally makes sure the
verbosity is restored in any case.
The call to gmsh.write is commented out: the generated gmsh/mobius.msh is
already stored in the repository, and uncommenting the call would overwrite it.
The gmsh.fltk.run() at the end only opens the graphical interface of Gmsh when
the command line carries a -popup argument, which never happens in a
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 问题 Solving a Poisson problem on the Möbius strip#
Mesh 读入 .msh 时用 dim=3 指明几何维数是 3: 这张网格的单元是三角形
(拓扑维数为 2), 但顶点坐标是三维的. triplot 画出读进来的网格, 可以和上面
matplotlib 画的曲面对照着看.
Mesh is given dim=3 when reading the .msh file to state that the geometric
dimension is 3: the cells of this mesh are triangles (the topological dimension
is 2), but the vertex coordinates are three-dimensional. triplot draws the mesh
that has been read in, which can be compared with the surface drawn with
matplotlib above.
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
在重新绑定之前就已经建好了, 所以不受影响, 但初学时容易看错.
A Poisson problem on a surface is written exactly as in the plane: dx is the
measure on the surface here and grad is the tangential gradient, so the
variational form below in fact discretizes the Laplace-Beltrami operator on the
surface. The right-hand side \(f = x^2 + y^2 + z^2\) is given by
SpatialCoordinate, and DirichletBC(V, 0, 'on_boundary') sets zero on the
outer boundary of the strip.
The bilinear form and the right-hand side are written in a single expression a
here and then split apart with lhs/rhs. Note that u is bound twice: first
as a TrialFunction and then as the Function holding the solution; a has
already been built before the rebinding and is therefore unaffected, but this is
easy to misread at first.
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 拉伸生成网格 Generating a mesh with extrusion for a 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.
A perfectly matched layer (PML) is an artificial absorbing layer commonly used for wave problems: a thin layer is wrapped around the physical region so that outgoing waves decay inside it and are hardly reflected back. This puts two requirements on the mesh: there has to be an extra layer of cells outside the physical region hugging its boundary, and the physical region and the absorbing layer must be markable separately.
The extrude_boundary_layer of Gmsh does exactly this: it extrudes a boundary
curve (in two dimensions) or a boundary surface (in three dimensions) along the
normal into several layers of cells. The two subsections below have exactly the
same structure: build a geometry, take its boundary, extrude the PML layer along
that boundary, and then mark the physical region, the PML layer and the inner and
outer boundaries with physical groups.
Once these markers have been read into Firedrake, the markers of the physical
region and of the PML layer are the subdomain numbers appearing in dx(1) and
dx(2). The boundary markers need care: the interface between the physical
region and the PML layer becomes an interior facet after the extrusion, so it
can only be integrated with dS(...); ds(...) gives 0. Measured locally on
gmsh/disk_pml.msh, the ds integral over tag 1 (inner bdy) is 0 while the dS
integral is 3.13, and tag 2 (outer bdy) is the real outer boundary (ds integral
3.76); the same holds for gmsh/sphere_pml.msh, where tag 2 is an interior facet
and tags 1, 3 and 4 are the outer boundaries. When these are confused, Firedrake
only emits a Subdomain is empty warning and the integral is silently computed
as 0.
4.2.1. 三维 Three dimensions#
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 厚的一层.
gmsh.model.occ.add_sphere(0, 0, 0, 0.5, angle1=0) builds a half ball of radius
0.5 with the OCC kernel (angle1=0 raises the start of the polar angle from
\(-\pi/2\) to \(0\), keeping only the upper half), so its boundary consists of both a
spherical surface and a disk at the bottom, which have to be told apart by type
with gmsh.model.get_type. Only the spherical surface (inner_bdy) is then
extruded into a boundary layer serving as the PML, and the half ball, the PML
layer and the inner and outer boundaries are numbered separately. The
gmsh/sphere_pml.msh stored in the repository has a measured outer radius of
0.6 — the ball radius 0.5 plus a layer of thickness 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. 二维 Two dimensions#
二维版本与三维完全平行: add_disk 造一个半径 0.5 的圆盘, 它的边界是一条曲线,
对这条曲线拉伸出边界层. numElements=[2]、heights=[0.1] 表示沿法向分 2 层、
总厚度 0.1, 生成的 gmsh/disk_pml.msh 外半径实测正是 0.6. 把 recombine 改成
True 可以得到含四边形的混合网格.
The two-dimensional version is exactly parallel to the three-dimensional one:
add_disk builds a disk of radius 0.5 whose boundary is a single curve, and a
boundary layer is extruded from that curve. numElements=[2] and
heights=[0.1] mean 2 layers along the normal with a total thickness of 0.1, and
the outer radius of the generated gmsh/disk_pml.msh is measured to be exactly
0.6. Setting recombine to True gives a mixed mesh containing
quadrilaterals.
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 空间中绘制区域标记 Plotting subdomain markers in an HDivT space#
网格读进来之后, 常常想确认标记有没有落在预期的单元和面上. 这里用一个自由度长在面上的空间 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 找到这条边对应的自由度, 把标记值写在边的中点上. 输出文件名带上进程数和进程号, 并行运行时每个进程各画自己那一份.
Once a mesh has been read in, one often wants to check that the markers landed on
the expected cells and facets. Here the subdomain markers are “drawn” with
HDivT, a space whose degrees of freedom live on the facets (the trace space of
H(div); at first order it is a linear polynomial on each edge, so two degrees of
freedom per edge).
The trick is to use par_loop on subdomain 1 (dx(1)) to add 1/2 to the facet
degrees of freedom of every cell (INC means accumulate): an edge whose two
neighboring cells both belong to subdomain 1 receives the increment twice and gets 1; an
edge with only one side in subdomain 1 is added to once and gets 1/2; an edge not
touching subdomain 1 at all stays 0. Values below 3/4 are then zeroed out, and
what is left with value 1 is exactly the set of interior edges of subdomain 1. In
this example, out of 656 degrees of freedom, 334 are measured to be 1, 46 are 1/2
and 276 are 0.
The plotting part below does not use the plotting functions of Firedrake; it
takes the coordinates of the two endpoints of every edge directly from
mesh.topology_dm to draw the segment, then finds the degree of freedom
belonging to that edge through the Section of the HDivT space and writes the
marker value at the midpoint of the edge. The output file names carry the number
of processes and the process number, so that in a parallel run every process
draws its own share.
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 simplex node ordering#
高阶网格的一个单元上有多个节点, 而不同软件对这些节点的排列约定并不一致: 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) 拿到的是整个二元组.
A cell of a higher-order mesh carries several nodes, and different pieces of software do not agree on how these nodes are arranged: Gmsh lists the vertices first, then the nodes on each edge, and only then the nodes interior to the faces and to the volume; finite element libraries more commonly number them in lexicographic order of the lattice coordinates. Reading a higher-order Gmsh mesh therefore requires a correspondence between the two numberings.
The following group of functions does exactly that. GmshLexOrder_SEG,
GmshLexOrder_TRI and GmshLexOrder_TET return a pair (lex, node): lex is
an array in which lex[k] is the position in the Gmsh numbering of the k-th
lattice point in lexicographic order, and node is the next available number
used during the recursion. The auxiliary functions come in three groups:
SN1/SN2/SN3 give the total number of nodes on a segment, a triangle and a
tetrahedron of degree \(p\); SI1/SI2/SI3 map lattice coordinates to the
lexicographic index; SL1/SL2/SL3 enumerate the interior lattice points. The
three GmshLexOrder_* follow the convention of Gmsh by numbering the vertices
first and then each edge, and finally hand the nodes interior to the faces and to
the volume over recursively to the same functions at reduced degree.
For instance GmshLexOrder_TRI(2)[0] is [0., 3., 1., 5., 4., 2.]: the 6
lattice points of a second-order triangle in lexicographic order correspond in
turn to the Gmsh nodes numbered 0, 3, 1, 5, 4 and 2. Note that it is the first
component of the return value that is taken; writing
lex = GmshLexOrder_TRI(2) directly gives the whole pair.
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. 高阶单元测试 Higher-order element test#
Firedrake 读 .msh 文件时只会给网格配一个一阶的坐标单元 (见 firedrake/mesh.py
的 make_mesh_from_mesh_topology, 那里至今留着一条
“meshfile might indicates higher-order coordinate element” 的 TODO). 因此直接读一个二阶网格会失败, 在 Firedrake 2026.4.1 下报的是坐标数组形状对不上:
When Firedrake reads a .msh file it only equips the mesh with a first-order
coordinate element (see make_mesh_from_mesh_topology in firedrake/mesh.py,
where a TODO reading “meshfile might indicates higher-order coordinate element”
is still left). Reading a second-order mesh directly therefore fails; under
Firedrake 2026.4.1 what is reported is a mismatch in the shape of the coordinate
array:
>>> Mesh('gmsh/cube_2rd.msh')
ValueError: cannot reshape array of size 720 into shape (14,3)
下面几个单元演示如何绕过这一步: 先从 DMPlex 上读出坐标单元的阶数, 再自己拼一个相应阶数的坐标单元, 最后把 Firedrake 的 make_mesh_from_mesh_topology 换掉.
The next few cells show how to work around this step: the degree of the
coordinate element is first read off the DMPlex, a coordinate element of that
degree is then built by hand, and finally the make_mesh_from_mesh_topology of
Firedrake is replaced.
已知问题
这段演示只解决了”读得进来”的问题, 读进来的几何是不对的. 用
-dm_plex_gmsh_project 让 PETSc 把逐单元存放的高阶坐标投影成连续坐标之后, 坐标向量里的数值并不是各节点的坐标. 在 Firedrake 2026.4.1 + PETSc 3.25 下实测, 读单位立方体网格 gmsh/cube_2rd.msh 得到的坐标三个分量最大值分别是 9、7 和 8.5
(本应都是 1), 由此建出来的网格体积约为 75.41 (本应是 1); 具体输出见下面打印坐标范围的那个单元. 因此本节只是”问题出在哪里”的记录, 还不是一个可用的高阶网格读取方案.
Known issue
This demonstration only solves the problem of getting the mesh read in; the
geometry that is read in is wrong. After -dm_plex_gmsh_project makes PETSc
project the cell-wise stored higher-order coordinates into continuous
coordinates, the numbers in the coordinate vector are not the coordinates of the
nodes. Measured under Firedrake 2026.4.1 with PETSc 3.25, reading the unit cube
mesh gmsh/cube_2rd.msh gives maxima of the three coordinate components of 9, 7
and 8.5 (all of which should be 1), and the mesh built from them has a volume of
about 75.41 (which should be 1); the actual output is shown by the cell below
that prints the coordinate range. This section is therefore only a record of
where the problem lies, not yet a usable scheme for reading higher-order
meshes.
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 行全零), 却不报任何错.
The two helper functions below each do one thing: getCoordinateFESpaceOrder
takes the name of the finite element that PETSc assigned to the coordinate field
(of the form P2) from the coordinate DM of the DMPlex and parses the degree
out of it, falling back to first order when it cannot be obtained;
coordinates_from_plex reorders the values of the DMPlex coordinate vector
according to the degree-of-freedom numbering of Firedrake and packs them into a
CoordinatelessFunction.
The assert inside coordinates_from_plex is there on purpose: it moves the
coordinates of only one node per mesh point, and is therefore only applicable to
at most second order on simplices. A third-order mesh has 2 nodes on each edge,
and without this assert the extra rows would silently be left as zero
coordinates (49 of the 172 rows are measured to be all zeros when reading
gmsh/cube_3rd.msh) without any error being raised.
这一段照着 firedrake/mesh.py 的 make_mesh_from_mesh_topology 改写, 与原版有
两处区别:
坐标单元的阶数改成从网格文件里读出来的, 而不是固定为一阶;
原版里处理周期网格的那个分支被丢掉了: 原版会先看
dm.getCoordinatesLocalized(), 对周期网格调用_fully_localize_coordinates并改用逐单元的DG/DQ一阶坐标单元 (variant="equispaced"), 这里没有照搬.
改写后的函数被塞回 firedrake.mesh 模块. Mesh() 内部是按模块级名字调用它的,
因此替换掉模块属性就能生效 —— 但这也意味着替换之后, 本进程里所有的网格构造都会走这段代码, 包括与本节毫无关系的那些.
This part is adapted from make_mesh_from_mesh_topology in firedrake/mesh.py,
with two differences from the original:
the degree of the coordinate element is read from the mesh file instead of being fixed at one;
the branch of the original that handles periodic meshes has been dropped: the original first looks at
dm.getCoordinatesLocalized(), calls_fully_localize_coordinatesfor a periodic mesh and switches to cell-wise first-orderDG/DQcoordinate elements (variant="equispaced"), which is not carried over here.
The rewritten function is put back into the firedrake.mesh module. Mesh()
calls it internally by its module-level name, so replacing the module attribute
is enough to take effect. But this also means that after the replacement every
mesh construction in this process goes through this code, including those that
have nothing to do with this section.
这段代码是实验性的, 而且会改动 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 之后的那个还原单元), 否则同一次会话里后面所有的网格都会受影响.
This code is experimental, and it modifies the internal state of Firedrake
It has three hard limits, and crossing them mostly raises no error at all; it just silently gives a wrong mesh:
Periodic meshes are not supported. What point 2 above drops is exactly the periodic-mesh branch. Measured locally, with this replacement installed the area of
PeriodicUnitSquareMesh(4, 4)changes from 1.0 to 2.25 and the length ofPeriodicIntervalMesh(8, 1.0)from 1.0 to 1.75, without any error being raised at any point; non-periodic meshes are unaffected (UnitSquareMesh(4, 4)is still 1.0).Only the case of one coordinate node per mesh point is supported, which means at most second order on simplices. This one is caught by the
assertincoordinates_from_plex, which does report a clear error.It must be used together with
-dm_plex_gmsh_projectandpetscdualspace_lagrange_continuity, otherwise it reportsValueError: could not broadcast input array from shape (0,) into shape (3,).
The original function must therefore be restored after use (see the restoring
cell after triplot below), otherwise every mesh built later in the same session
is affected.
-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,).
-dm_plex_gmsh_project makes PETSc reproject the coordinates after the mesh has
been read, but its default target space is still discontinuous: measured with
only this option given, the storage size of the coordinate Section is still 720
and the degrees of freedom still hang on the cells, no different from not giving
the option at all. petscdualspace_lagrange_continuity has to be added as well
before the coordinates become a continuous P2 hanging on the vertices and edges
(storage size 189, which is (14 vertices + 49 edges) \( imes\) 3).
Both of these are PETSc command-line options, and inserted_options() of
OptionsManager is used to insert them temporarily into the options database of
PETSc; they are removed automatically on leaving the with block. For the
replaced make_mesh_from_mesh_topology above these two options are
required; without them it reports
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, 实际输出与此相差很远,
这正是上面”已知问题”里说的那种情况.
The mesh that has been read in is drawn below, and the coordinate range and the
mesh volume are printed at the same time. gmsh/cube_2rd.msh is the unit cube,
so the coordinates should lie in \([0, 1]^3\) and the volume should equal 1; the
actual output is far from that, which is exactly the situation described in the
“Known issue” above.
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 还原回去.
若不还原, 同一次会话里后面建的每一张网格都会走上面那段实验代码, 周期网格会被静默算错.
Now that the experiment is done, the original make_mesh_from_mesh_topology of
Firedrake is restored immediately. Without restoring it, every mesh built later
in the same session would go through the experimental code above, and periodic
meshes would be silently computed wrong.
fd.mesh.make_mesh_from_mesh_topology = _orig_make_mesh_from_mesh_topology
4.3.1.1. Gmsh 网格读取的问题 Problems in reading Gmsh meshes#
这一小节绕开 Firedrake, 直接用 PETSc.DMPlex 读同一个文件, 对比”投影前”和
“投影后”两种坐标存放方式, 看清上面那个问题出在哪一步.
show_plex_coordinates 打印投影前的坐标: 此时高阶坐标是逐单元存放的, 坐标
Section 的自由度挂在单元上, 每个四面体存 10 个节点的坐标, 数值都落在 \([0, 1]\)
内, 是对的.
show_plex_coordinates_continue 打印投影后的坐标: 此时自由度改挂到了顶点和边上
(二阶四面体每个顶点、每条边各一个节点), 所以要沿单元的传递闭包
(getTransitiveClosure) 把这些点找出来再取值. 打印出的数值明显超出了 \([0, 1]\),
可见问题出在这一步投影上, 而不是出在前面替换 make_mesh_from_mesh_topology 的那段代码里.
This subsection bypasses Firedrake and reads the same file directly with
PETSc.DMPlex, comparing the two ways of storing the coordinates, “before
projection” and “after projection”, in order to see at which step the problem
above arises.
show_plex_coordinates prints the coordinates before the projection: the
higher-order coordinates are stored cell by cell at this point, the degrees of
freedom of the coordinate Section hang on the cells, every tetrahedron stores
the coordinates of 10 nodes, and the values all lie within \([0, 1]\), which is
correct.
show_plex_coordinates_continue prints the coordinates after the projection: the
degrees of freedom have moved to the vertices and the edges (one node on each
vertex and each edge of a second-order tetrahedron), so the transitive closure of
the cell (getTransitiveClosure) has to be walked to find these points and read
their values. The printed values clearly exceed \([0, 1]\), which shows that the
problem lies in this projection step, and not in the code replacing
make_mesh_from_mesh_topology earlier.
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]