Gmsh 使用入门：获取网格信息

  之前的两篇文章，讲解了使用 Gmsh 如何构建几何模型与进行网格生成:

       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，返回所有点
• 若 dim=0，则返回角点
• 若 dim=1，则返回边上的点
• 若 dim=2，则返回面上的点
• 三维情况下，若 dim=3 返回几何体上的点

      此时 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。  

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/ 查找更多资料。