首页/文章/ 详情

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

1年前浏览3417

  之前的两篇文章,讲解了使用 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 中的 VisibilityColor 中进行调整。

       获得节点的命令如下

vtags,vxyz, vp = gmsh.model.mesh.getNodes()

      其中 vtagsvxyz 是一维数组,分别返回所有点的 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=1vtags 返回值为 [5 1 2 6 2 3 7 3 4 8 4 1],按照边返回了所有的边界点。默认情况下 includeBoundaryFalse,此时返回 [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 也能达到同样的效果。

       由于使用 getNodesgetElementsByType 可以获得所有节点与单元,因此可将网格传入到 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 获得其表示的值,再通过 getNodesgetElementsByType 获得节点和单元信息,传入 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 edgeNodesedgeTagsedgeNodes 均是一维数组。

  除了获取所有边的信息,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/ 查找更多资料。



来源:易木木响叮当
科普网格处理
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2023-06-14
最近编辑:1年前
易木木响叮当
硕士 有限元爱好者
获赞 220粉丝 261文章 349课程 2
点赞
收藏
未登录
1条评论
@hdil
签名征集中
1年前
你好,最后的代码运行之后,打印出来的面上的三角单元是个空列表,是什么原因呢,我一直遇到这个问题,想请教一下,谢谢了
回复
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈