Gmsh教程

网卡了

于 2024-07-17 06:00:00 发布

阅读量104

点赞数

分类专栏： Gmsh 文章标签： Gmsh

本文链接：https://blog.csdn.net/weixin_44897632/article/details/140450916

版权

Gmsh 专栏收录该内容

5 篇文章 0 订阅

订阅专栏

1、几何和网格数据

# 导入Gmsh模块和sys模块用于处理命令行参数  
import gmsh  
import sys  
  
# 初始化Gmsh  
gmsh.initialize()  
  
# 检查命令行参数，如果提供了参数且不是以'-'开头，则尝试打开该文件  
if len(sys.argv) > 1 and sys.argv[1][0] != '-':  
    gmsh.open(sys.argv[1])  
else:  
    # 否则，创建一个简单的圆锥几何体并生成网格  
    gmsh.model.occ.addCone(1, 0, 0, 1, 0, 0, 0.5, 0.1)  # 创建一个圆锥，参数依次为：标签、底部中心点x,y,z、顶部中心点x,y,z、半径、高度  
    gmsh.model.occ.synchronize()  # 同步OCC内核与Gmsh模型  
    gmsh.model.mesh.generate()  # 生成网格  
  
# 打印当前模型的名称和维度  
print('Model ' + gmsh.model.getCurrent() + ' (' + str(gmsh.model.getDimension()) + 'D)')  
  
# 获取模型中的所有基本实体（点、曲线、表面、体积）  
entities = gmsh.model.getEntities()  
  
for e in entities:  
    dim, tag = e  # 分解出实体的维度和标签  
  
    # 获取实体的网格节点信息  
    nodeTags, nodeCoords, nodeParams = gmsh.model.mesh.getNodes(dim, tag)  
  
    # 获取实体的网格元素信息  
    elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)  
  
    # 获取实体的类型和名称  
    type = gmsh.model.getType(dim, tag)  
    name = gmsh.model.getEntityName(dim, tag)  
    if len(name): name += ' '  
  
    # 打印实体信息  
    print("Entity " + name + str(e) + " of type " + type)  
  
    # 计算并打印网格中的节点和元素数量  
    numElem = sum(len(i) for i in elemTags)  
    print(" - Mesh has " + str(len(nodeTags)) + " nodes and " + str(numElem) + " elements")  
  
    # 获取并打印实体的上下相邻实体信息  
    up, down = gmsh.model.getAdjacencies(dim, tag)  
    if len(up):  
        print(" - Upward adjacencies: " + str(up))  
    if len(down):  
        print(" - Downward adjacencies: " + str(down))  
  
    # 检查实体是否属于物理组，并打印物理组信息  
    physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)  
    if len(physicalTags):  
        s = ''  
        for p in physicalTags:  
            n = gmsh.model.getPhysicalName(dim, p)  
            if n: n += ' '  
            s += n + '(' + str(dim) + ', ' + str(p) + ') '  
        print(" - Physical groups: " + s)  
  
    # 检查实体是否是分区实体，并打印分区和父实体信息  
    partitions = gmsh.model.getPartitions(dim, tag)  
    if len(partitions):  
        print(" - Partition tags: " + str(partitions) + " - parent entity " + str(gmsh.model.getParent(dim, tag)))  
  
    # 打印构成实体网格的所有元素类型信息  
    for t in elemTypes:  
        name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(t)  
        print(" - Element type: " + name + ", order " + str(order) + " (" + str(numv) + " nodes in param coord: " + str(parv) + ")")  
  
# 如果命令行参数中没有'-nopopup'，则启动Gmsh GUI以显示模型  
if '-nopopup' not in sys.argv:  
    gmsh.fltk.run()  
  
# 清除所有模型数据  
gmsh.clear()  
  
# 结束Gmsh会话  
gmsh.finalize()

2、网格导入、离散实体、混合模型、地形网格

import gmsh  # 导入Gmsh库  
import sys  # 导入系统库，用于访问命令行参数  
import math  # 导入数学库  
  
gmsh.initialize()  # 初始化Gmsh库  
  
# 创建一个新的模型，命名为"x2"  
gmsh.model.add("x2")  
  
N = 100  # 设置网格的细分程度  
  
# 定义一个辅助函数，用于计算节点的标签  
def tag(i, j):  
    return (N + 1) * i + j + 1  
  
coords = []  # 用于存储节点坐标的列表  
nodes = []  # 用于存储节点标签的列表  
tris = []  # 用于存储三角形元素标签的列表  
lin = [[], [], [], []]  # 用于存储线性元素（边）标签的列表  
  
# 定义四个角点的标签  
pnt = [tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)]  
  
# 遍历每个网格点，计算其坐标，并存储节点标签和坐标  
for i in range(N + 1):  
    for j in range(N + 1):  
        nodes.append(tag(i, j))  
        # 计算节点的x, y, z坐标，z坐标根据x和y的某种函数关系计算  
        coords.extend([float(i) / N, float(j) / N, 0.05 * math.sin(10 * float(i + j) / N)])  
        # 添加三角形元素  
        if i > 0 and j > 0:  
            tris.extend([tag(i - 1, j - 1), tag(i, j - 1), tag(i - 1, j)])  
            tris.extend([tag(i, j - 1), tag(i, j), tag(i - 1, j)])  
        # 添加边界线性元素  
        if (i == 0 or i == N) and j > 0:  
            lin[3 if i == 0 else 1].extend([tag(i, j - 1), tag(i, j)])  
        if (j == 0 or j == N) and i > 0:  
            lin[0 if j == 0 else 2].extend([tag(i - 1, j), tag(i, j)])  
  
# 添加四个离散实体（点）  
for i in range(4):  
    gmsh.model.addDiscreteEntity(0, i + 1)  
# 设置四个角点的坐标  
gmsh.model.setCoordinates(1, 0, 0, coords[3 * tag(0, 0) - 1])  
gmsh.model.setCoordinates(2, 1, 0, coords[3 * tag(N, 0) - 1])  
gmsh.model.setCoordinates(3, 1, 1, coords[3 * tag(N, N) - 1])  
gmsh.model.setCoordinates(4, 0, 1, coords[3 * tag(0, N) - 1])  
  
# 添加四条线性实体（边）  
for i in range(4):  
    gmsh.model.addDiscreteEntity(1, i + 1, [i + 1, i + 2 if i < 3 else 1])  
  
# 添加一个曲面实体  
gmsh.model.addDiscreteEntity(2, 1, [1, 2, -3, -4])  
  
# 添加节点和坐标到网格  
gmsh.model.mesh.addNodes(2, 1, nodes, coords)  
  
# 添加不同类型的元素到网格  
for i in range(4):  
    gmsh.model.mesh.addElementsByType(i + 1, 15, [], [pnt[i]])  # 添加点元素  
    gmsh.model.mesh.addElementsByType(i + 1, 1, [], lin[i])     # 添加线性元素  
  
# 添加三角形元素到网格  
gmsh.model.mesh.addElementsByType(1, 2, [], tris)  
  
# 重新分类节点和创建几何体  
gmsh.model.mesh.reclassifyNodes()  
gmsh.model.mesh.createGeometry()  
  
# 添加一个额外的几何体，用于演示如何添加更复杂的几何形状  
# ...（此处省略了具体的几何体添加代码，因为主要是重复的添加点和线，然后形成面和体的过程）  
  
# 同步几何体，准备网格化  
gmsh.model.geo.synchronize()  
  
# 设置网格化参数  
transfinite = False  # 是否使用Transfinite网格化方法  
transfiniteAuto = False  # 是否使用自动Transfinite网格化方法  
  
# 根据不同的网格化参数设置，进行网格化  
if transfinite:  
    # 使用Transfinite网格化方法  
    NN = 30  
    for c in gmsh.model.getEntities(1):  
        gmsh.model.mesh.setTransfiniteCurve(c[1], NN)  
    for s in gmsh.model.getEntities(2):  
        gmsh.model.mesh.setTransfiniteSurface(s[1])  
        gmsh.model.mesh.setRecombine(s[0], s[1])  
        gmsh.model.mesh.setSmoothing(s[0], s[1], 100)  
    gmsh.model.mesh.setTransfiniteVolume(v1)  
elif transfiniteAuto:  
    # 使用自动Transfinite网格化方法  
    gmsh.option.setNumber('Mesh.MeshSizeMin', 0.5)  
    gmsh.option.setNumber('Mesh.MeshSizeMax', 0.5)  
    gmsh.model.mesh.setTransfiniteAutomatic()  
else:  
    # 使用默认的网格化方法  
    gmsh.option.setNumber('Mesh.MeshSizeMin', 0.05)  
    gmsh.option.setNumber('Mesh.MeshSizeMax', 0.05)  
  
# 生成网格  
gmsh.model.mesh.generate(3)  
# 将网格写入文件  
gmsh.write('x2.msh')  
  
# 如果命令行参数中没有'-nopopup'，则运行Gmsh的图形用户界面  
if '-nopopup' not in sys.argv:  
    gmsh.fltk.run()  
  
# 清理资源，结束Gmsh的使用  
gmsh.finalize()

3、后处理数据导入：基于列表

import gmsh
import sys

gmsh.initialize(sys.argv)

# 下面的循环尝试动态修改三角形坐标，
t1 = gmsh.view.add("A list-based view")
	triangle1 = [0., 1., 1., # x coordinates of the 3 triangle nodes
	0., 0., 1., # y coordinates of the 3 triangle nodes
	0., 0., 0.] # z coordinates of the 3 triangle nodes
	triangle2 = [0., 1., 0., 0., 1., 1., 0., 0., 0.]
	for step in range(0, 10) :
		triangle1.extend([10., 11. - step, 12.])  
		triangle2.extend([11., 12., 13. + step])
gmsh.view.addListData(t1, "ST", 2, triangle1 + triangle2)


line = [
	0., 1., # x coordinate of the 2 line nodes
		1.2, 1.2, # y coordinate of the 2 line nodes
		0., 0.    # z coordinate of the 2 line nodes
]
for step in range(0, 10) :

	line.extend([10. + step, 0., 0.,
		10. + step, 0., 0.])
	gmsh.view.addListData(t1, "VL", 1, line)

			gmsh.view.addListDataString(t1, [20., -20.], ["Created with Gmsh"])
			gmsh.view.addListDataString(t1, [0.5, 1.5, 0.],
				["A multi-step list-based view"],
				["Align", "Center", "Font", "Helvetica"])


			gmsh.view.option.setNumber(t1, "TimeStep", 5)
			gmsh.view.option.setNumber(t1, "IntervalsType", 3)
			ns = gmsh.view.option.getNumber(t1, "NbTimeStep")
			print("View " + str(t1) + " has " + str(ns) + " time steps")


				print("Value at (0.9, 0.1, 0)", gmsh.view.probe(t1, 0.9, 0.1, 0))


				gmsh.view.write(t1, "x3.pos")

t2 = gmsh.view.add("Second order quad")


quad = [0., 1., 1., 0., # x coordinates of the 4 quadrangle nodes
- 1.2, -1.2, -0.2, -0.2, # y coordinates of the 4 quadrangle nodes
0., 0., 0., 0.] # z coordinates of the 4 quadrangle nodes


quad.extend([1., 1., 1., 1., 3., 3., 3., 3., -3.])


gmsh.view.setInterpolationMatrices(t2, "Quadrangle", 9,
	[0, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0.25,
	0, 0, 0.25, 0, 0, -0.25, 0.25, 0, -0.25,
	0, 0, 0.25, 0, 0, 0.25, 0.25, 0, 0.25,
	0, 0, 0.25, 0, 0, 0.25, -0.25, 0, -0.25,
	0, 0, -0.5, 0.5, 0, 0.5, 0, -0.5, 0,
	0, 0.5, -0.5, 0, 0.5, 0, -0.5, 0, 0,
	0, 0, -0.5, 0.5, 0, -0.5, 0, 0.5, 0,
	0, 0.5, -0.5, 0, -0.5, 0, 0.5, 0, 0,
	1, -1, 1, -1, 0, 0, 0, 0, 0],
	[0, 0, 0,
	2, 0, 0,
	2, 2, 0,
	0, 2, 0,
	1, 0, 0,
	2, 1, 0,
	1, 2, 0,
	0, 1, 0,
	1, 1, 0])

gmsh.view.addListData(t2, "SQ", 1, quad)
	gmsh.view.option.setNumber(t2, "AdaptVisualizationGrid", 1)
	gmsh.view.option.setNumber(t2, "TargetError", 1e-2)
	gmsh.view.option.setNumber(t2, "MaxRecursionLevel", 5)


if '-nopopup' not in sys.argv :
	gmsh.fltk.run()

	gmsh.finalize()

4、后处理数据导入：基于模型

import gmsh
import sys

gmsh.initialize(sys.argv)
创建第一个模型：
gmsh.model.add("simple model")
surf = gmsh.model.addDiscreteEntity(2)

gmsh.model.mesh.addNodes(2, surf, [1, 2, 3, 4],
                         [0., 0., 0., 1., 0., 0., 1., 1., 0., 0., 1., 0.])
gmsh.model.mesh.addElementsByType(surf, 2, [1, 2], [1, 2, 3, 1, 3, 4])

添加和可视化节点数据：

t1 = gmsh.view.add("Continuous")
for step in range(0, 10):
    gmsh.view.addHomogeneousModelData(
        t1, step, "simple model", "NodeData",
        [1, 2, 3, 4],  # tags of nodes
        [10., 10., 12. + step, 13. + step])  # data, per node

t2 = gmsh.view.add("Discontinuous")
for step in range(0, 10):
    gmsh.view.addHomogeneousModelData(
        t2, step, "simple model", "ElementNodeData",
        [1, 2],  
        [10., 10., 12. + step, 14., 15., 13. + step]) 
创建第二个模型：
gmsh.model.add("another model")
gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)
获取和可视化第二个模型的节点数据：

nodes, coord, _ = gmsh.model.mesh.getNodes()
for step in range(11, 20):
    gmsh.view.addHomogeneousModelData(
        t1, step, "another model", "NodeData", nodes,
        [step * coord[i] for i in range(0, len(coord), 3)])

gmsh.view.write(t1, "x4_t1.msh")
gmsh.view.write(t2, "x4_t2.msh")

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()

5、附加几何数据：参数化、法线、曲率

import gmsh
import sys
import math

gmsh.initialize(sys.argv)

gmsh.model.add("x5")
s = gmsh.model.occ.addSphere(0, 0, 0, 1)
b = gmsh.model.occ.addBox(0.5, 0, 0, 1.3, 2, 3)
gmsh.model.occ.fuse([(3, s)], [(3, b)])
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)

normals = []
curvatures = []
计算法线和曲率：
for e in gmsh.model.getEntities(2):

    s = e[1]

    tags, coord, param = gmsh.model.mesh.getNodes(2, s, True)

    norm = gmsh.model.getNormal(s, param)

    curv = gmsh.model.getCurvature(2, s, param)

    for i in range(0, len(coord), 3):
        normals.append(coord[i])
        normals.append(coord[i + 1])
        normals.append(coord[i + 2])
        normals.append(norm[i])
        normals.append(norm[i + 1])
        normals.append(norm[i + 2])
        curvatures.append(coord[i])
        curvatures.append(coord[i + 1])
        curvatures.append(coord[i + 2])
        curvatures.append(curv[i // 3])

可视化法线和曲率
vn = gmsh.view.add("normals")
gmsh.view.addListData(vn, "VP", len(normals) // 6, normals)
gmsh.view.option.setNumber(vn, 'ShowScale', 0)
gmsh.view.option.setNumber(vn, 'ArrowSizeMax', 30)
gmsh.view.option.setNumber(vn, 'ColormapNumber', 19)
vc = gmsh.view.add("curvatures")
gmsh.view.addListData(vc, "SP", len(curvatures) // 4, curvatures)
gmsh.view.option.setNumber(vc, 'ShowScale', 0)
参数化和坐标转换：
bounds = gmsh.model.getParametrizationBounds(1, 5)
N = 20
t = [bounds[0][0] + i * (bounds[1][0] - bounds[0][0]) / N for i in range(N)]
xyz1 = gmsh.model.getValue(1, 5, t)

uv = gmsh.model.reparametrizeOnSurface(1, 5, t, 1)
xyz2 = gmsh.model.getValue(2, 1, uv)

if max([abs(a - b) for (a, b) in zip(xyz1, xyz2)]) < 1e-12:
    gmsh.logger.write('Evaluation on curve and surface match!')
else:
    gmsh.logger.write('Evaluation on curve and surface do not match!', 'error')

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()
注释代码

6、附加网格数据：积分点、雅可比矩阵和基函数

import gmsh
import sys

gmsh.initialize(sys.argv)

gmsh.model.add("x6")

添加几何体并生成网格：
gmsh.model.occ.addRectangle(0, 0, 0, 1, 0.1)：在原点添加一个矩形，长度为1，高度为0.1。
gmsh.model.occ.synchronize()：同步几何体，确保所有更改都被考虑。
gmsh.model.mesh.setTransfiniteAutomatic()：设置网格生成为自动超限插值。
gmsh.model.mesh.generate(2)：生成2维网格。
设置网格元素阶数：
elementOrder = 1
interpolationOrder = 2
gmsh.model.mesh.setOrder(elementOrder)

获取并打印网格元素信息：
def pp(label, v, mult):
    print(" * " + str(len(v) / mult) + " " + label + ": " + str(v))

elementTypes = gmsh.model.mesh.getElementTypes()
遍历所有元素类型，获取并打印每种元素类型的名称、维度、阶数、节点数、局部节点坐标、基本节点数等信息。
对于每种元素类型，获取并打印积分点、基函数及其梯度、雅可比矩阵的行列式等信息。
结束：
for t in elementTypes:
    elementName, dim, order, numNodes, localNodeCoord, numPrimNodes =\
    gmsh.model.mesh.getElementProperties(t)
    print("\n** " + elementName + " **\n")

    localCoords, weights =\
    gmsh.model.mesh.getIntegrationPoints(t, "Gauss" + str(interpolationOrder))
    pp("integration points to integrate order " +
       str(interpolationOrder) + " polynomials", localCoords, 3)

    numComponents, basisFunctions, numOrientations =\
    gmsh.model.mesh.getBasisFunctions(t, localCoords, "Lagrange")
    pp("basis functions at integration points", basisFunctions, 1)
    numComponents, basisFunctions, numOrientations =\
    gmsh.model.mesh.getBasisFunctions(t, localCoords, "GradLagrange")
    pp("basis function gradients at integration points", basisFunctions, 3)

    jacobians, determinants, coords =\
    gmsh.model.mesh.getJacobians(t, localCoords)
    pp("Jacobian determinants at integration points", determinants, 1)

gmsh.finalize()

7、其他网格数据：内部边缘和面

import sys  
import gmsh  # 导入gmsh库，用于几何建模和网格生成  
  
gmsh.initialize(sys.argv)  # 初始化gmsh环境  
  
gmsh.model.add("x7")  # 添加一个新的模型，命名为"x7"  
  
# 在模型中添加一个立方体，左下角位于原点，边长为1  
gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)  
gmsh.model.occ.synchronize()  # 同步几何实体，准备进行网格划分  
  
# 设置网格的最小尺寸为2  
gmsh.option.setNumber("Mesh.MeshSizeMin", 2.)  
gmsh.model.mesh.generate(3)  # 生成三维网格  
  
# 获取四面体元素类型  
elementType = gmsh.model.mesh.getElementType("tetrahedron", 1)  
edgeNodes = gmsh.model.mesh.getElementEdgeNodes(elementType)  # 获取四面体的边节点  
faceNodes = gmsh.model.mesh.getElementFaceNodes(elementType, 3)  # 获取四面体的面节点  
  
# 创建网格的边和面  
gmsh.model.mesh.createEdges()  
gmsh.model.mesh.createFaces()  
  
# 获取所有边的标签和方向  
edgeTags, edgeOrientations = gmsh.model.mesh.getEdges(edgeNodes)  
# 获取所有面的标签和方向，只考虑3维的面  
faceTags, faceOrientations = gmsh.model.mesh.getFaces(3, faceNodes)  
  
# 获取所有四面体的标签和节点标签  
elementTags, elementNodeTags = gmsh.model.mesh.getElementsByType(elementType)  
edges2Elements = {}  
faces2Elements = {}  
# 遍历所有边和面，记录每个边和面所属的四面体  
for i in range(len(edgeTags)): # 每个四面体有6条边  
    if not edgeTags[i] in edges2Elements:  
        edges2Elements[edgeTags[i]] = [elementTags[i // 6]]  
    else:  
        edges2Elements[edgeTags[i]].append(elementTags[i // 6])  
for i in range(len(faceTags)): # 每个四面体有4个面  
    if not faceTags[i] in faces2Elements:  
        faces2Elements[faceTags[i]] = [elementTags[i // 4]]  
    else:  
        faces2Elements[faceTags[i]].append(elementTags[i // 4])  
  
# 添加一个新的离散实体，维度为2  
s = gmsh.model.addDiscreteEntity(2)  
  
maxElementTag = gmsh.model.mesh.getMaxElementTag()  # 获取当前最大的元素标签  
uniqueFaceTags = set()  
tagsForTriangles = []  
faceNodesForTriangles = []  
# 遍历所有面的标签，为每一个面创建一个三角形元素  
for i in range(len(faceTags)):  
    if faceTags[i] not in uniqueFaceTags:  
        uniqueFaceTags.add(faceTags[i])  
        tagsForTriangles.append(faceTags[i] + maxElementTag)  
        faceNodesForTriangles.append(faceNodes[3 * i])  
        faceNodesForTriangles.append(faceNodes[3 * i + 1])  
        faceNodesForTriangles.append(faceNodes[3 * i + 2])  
# 获取三角形元素类型  
elementType2D = gmsh.model.mesh.getElementType("triangle", 1)  
# 添加三角形元素到模型中  
gmsh.model.mesh.addElementsByType(s, elementType2D, tagsForTriangles,  
                                  faceNodesForTriangles)  
  
# 打印每个三角形连接的四面体  
for t in tagsForTriangles:  
    print("triangle " + str(int(t)) + " is connected to tetrahedra " +  
          str(faces2Elements[t - maxElementTag]))  
  
# 获取所有边的标签和节点  
edgeTags, edgeNodes = gmsh.model.mesh.getAllEdges()  
# 获取所有面的标签和节点，只考虑3维的面  
faceTags, faceNodes = gmsh.model.mesh.getAllFaces(3)  
  
# 如果命令行参数中没有'-nopopup'，则启动gmsh的图形界面  
if '-nopopup' not in sys.argv:  
    gmsh.fltk.run()  
  
gmsh.finalize()  # 清理并退出gmsh环境