之前的两篇文章,讲解了使用 Gmsh 如何构建几何模型与进行网格生成:
但在网格生成后,我们通常还需要利用网格的信息去进行计算,因此本文主要讲解如何得到 Gmsh 网格的信息。
首先生成一个网格作为示例
import gmsh
import numpy as np
from fealpy.mesh import TriangleMesh
import matplotlib.pyplot as plt
gmsh.initialize()
gmsh.model.add("t1")
lc=0.1
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(1, 0, 0, lc, 2)
gmsh.model.geo.addPoint(1, 1, 0, lc, 3)
gmsh.model.geo.addPoint(0, 1, 0, lc, 4)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
for i in range(1,5):
gmsh.model.geo.mesh.setTransfiniteCurve(i,3)
gmsh.model.geo.mesh.setTransfiniteSurface(1,"Right")
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
生成网格如下,其中节点的 tag 和颜色可以在左上角的 Tools/Options/Mesh
中的 Visibility
和 Color
中进行调整。
获得节点的命令如下
vtags,vxyz, vp = gmsh.model.mesh.getNodes()
其中 vtags
和 vxyz
是一维数组,分别返回所有点的 tag 和点的坐标,vp
返回点的参数坐标,一般情况下返回的都是空数组。除此以外,还可以输入参数 dim
,tag
,includeBoundary
,例如
vtags,vxyz, vp = gmsh.model.mesh.getNodes(dim=1,tag=-1,includeBoundary=True)
此时 dim=1
,vtags
返回值为 [5 1 2 6 2 3 7 3 4 8 4 1]
,按照边返回了所有的边界点。默认情况下 includeBoundary
为 False
,此时返回 [5 6 7 8]
,只返回边界的内部点。默认 tag=-1
,返回该维数下所有实体的点,若指定 tag,则只返回对应实体上的点。
通过以下命令,可以获得单元信息
tritype = gmsh.model.mesh.getElementType("Triangle",1)
etags,evtags = gmsh.model.mesh.getElementsByType(tritype)
在 Gmsh 中,不同的单元用不同的整数表示,使用函数 getElementType
输入单元的类型与单元阶数,则可以返回表示该单元的数字。将数字传入函数 getELementsByType
中,则可以返回该单元的 tag etags
与单元中点的 tag evtags
,它们均是一维数组。一阶三角形用整数 2
来表示,因此输入参数 2
也能达到同样的效果。
由于使用 getNodes
和 getElementsByType
可以获得所有节点与单元,因此可将网格传入到 FEALPy 中。
ntags, vxyz, _ = gmsh.model.mesh.getNodes()
node = vxyz.reshape((-1,3))
node = node[:,:2]
vmap = dict({j:i for i,j in enumerate(ntags)})
tris_tags,evtags = gmsh.model.mesh.getElementsByType(2)
evid = np.array([vmap[j] for j in evtags])
cell = evid.reshape((tris_tags.shape[-1],-1))
同样的,四边形、四面体、六面体网格也可以通过 getElementType
获得其表示的值,再通过 getNodes
和 getElementsByType
获得节点和单元信息,传入 FEALPy 中。
Gmsh 中通过节点和单元就可以描述网格,边与面的信息是不必要的,因此网格中并不存储边与面的实体。不过 Gmsh 中仍然可以获得边与面的信息。
通过如下代码可以获得所有边的信息
gmsh.model.mesh.createEdges()
edgeTags,edgeNodes = gmsh.model.mesh.getAllEdges()
edge = edgeNodes.reshape((-1,2))
由于 Gmsh 网格默认没有存储边的实体,因此需要使用 createEdges
生成边实体,并通过 getAllEdges
得到所有边的 tag edgeTags
与对应 node 的 tag edgeNodes
,edgeTags
与 edgeNodes
均是一维数组。
除了获取所有边的信息,Gmsh 还可以获得边界处的边与单元所对应的边。
linetype = gmsh.model.mesh.getElementType("Line",1)
tritype = gmsh.model.mesh.getElementType("Triangle",1)
LedgeNodes = gmsh.model.mesh.getElementEdgeNodes(linetype)
EedgeNodes = gmsh.model.mesh.getElementEdgeNodes(tritype)
LedgeTags, _ = gmsh.model.mesh.getEdges(LedgeNodes)
EedgeTags, _ = gmsh.model.mesh.getEdges(EedgeNodes)
上述代码通过 getElementType
获得表示边单元和三角形单元的数字,并将其传入 getElementEdgeNodes
中,得到组成实体的边的节点 tag,并将节点 tag 传入 getEdges
中,得到边tag。
getElementEdgeNodes
中,可以输入实体的 tag,从而返回指定实体上的边节点 tag,再结合函数 gmsh.model.getEntities
获得实体的 tag,则可以达到我们的目的。仍以文章开头处的网格为例,结合如下代码gmsh.model.mesh.generate(2)
gmsh.model.mesh.createEdges()
linetype = gmsh.model.mesh.getElementType("Line",1)
dim1tag = gmsh.model.getEntities(1)
dim1tag=np.asarray(dim1tag)
Linetag = dim1tag[:,1]
for i in Linetag:
LEdgeNodes = gmsh.model.mesh.getElementEdgeNodes(linetype,i)
print("Line%d上的edge:"%i)
print(LEdgeNodes.reshape(-1,2))
gmsh.fltk.run()
gmsh.finalize()
结果如下
三维下情况类似,这里仅给出代码,不作说明
import gmsh
import numpy as np
gmsh.initialize()
gmsh.model.add("t1")
gmsh.model.occ.addBox(0,0,0,1,1,1)
gmsh.model.occ.synchronize() # 同步到模型
gmsh.option.setNumber("Mesh.MeshSizeMin", 2.)
gmsh.model.mesh.generate(3) # 生成网格
gmsh.model.mesh.createEdges()
gmsh.model.mesh.createFaces()
tritype = gmsh.model.mesh.getElementType("Triangle",1)
dim2tag = gmsh.model.getEntities(2)
dim2tag=np.asarray(dim2tag)
Surfacetag = dim2tag[:,1]
for i in Surfacetag:
StriEdgeNodes = gmsh.model.mesh.getElementEdgeNodes(tritype,i)
print("Surface%d上的三角单元:"%i)
print(StriEdgeNodes.reshape(-1,3,2))
gmsh.fltk.run()
gmsh.finalize()
通过这三篇文章,虽然称不上面面俱到,但也能对 Gmsh 有一个较为具体的理解,能使用 Gmsh 进行一些工作。不过 Gmsh 的功能还不止文章里提到的这些,Gmsh 还具有生成自适应网格、有限元计算等功能,如果想进一步了解 Gmsh,可以去 Gmsh 的官方网站 https://gmsh.info/ 查找更多资料。