当与高阶元素(地质力学模拟所需)进行网格划分时,以stl文件形式给出的相当中等形状的地形会导致网格不可用(jac.<0 )。我在官方 Gmsh 演示模型terrain_stl.py中演示了这个问题。
使用 (tetra mesh)将行添加gmsh.model.mesh.setOrder(2)
到此脚本transfinite=False
会导致警告
Warning : Surface mesh: worst distortion = -0.148472 (avg = 0.990486, 2 elements with jac. < 0); worst gamma = 0.532537
Warning : Volume mesh: worst distortion = -0.124232 (avg = 0.997014, 1 elements with jac. < 0)
并且生成的网格无法用于后续的 FEM 模拟。
设置transfinite=True
(十六进制网格)会发出警告Failed to compute equidistant parameters
,但生成的网格有效。然而,这种超限网格的制备通常需要更多的努力,甚至并不总是可行的。
这些扭曲元素背后的原因是什么?如何规避它们?是否有助于引入中间层、所有面的 STL 网格,或者是否有特别适合的网格划分算法或选项?
为了完整起见,这里是完整的代码(Python)
import gmsh
import math
import os
import sys
gmsh.initialize(sys.argv)
path = os.path.dirname(os.path.abspath(__file__))
# load an STL surface
gmsh.merge(os.path.join(path, 'terrain_stl_data.stl'))
# classify the surface mesh according to given angle, and create discrete model
# entities (surfaces, curves and points) accordingly; curveAngle forces bounding
# curves to be split on sharp corners
gmsh.model.mesh.classifySurfaces(math.pi, curveAngle=math.pi / 3)
# create a geometry for the discrete curves and surfaces
gmsh.model.mesh.createGeometry()
# retrieve the surface, its boundary curves and corner points
s = gmsh.model.getEntities(2)
c = gmsh.model.getBoundary(s)
if (len(c) != 4):
gmsh.logger.write('Should have 4 boundary curves!', level='error')
p = []
xyz = []
for i in range(len(c)):
pt = gmsh.model.getBoundary([c[i]], combined=False)
p.extend([pt[0][1]])
xyz.extend(gmsh.model.getValue(0, pt[0][1], []))
# create other CAD entities to form one volume below the terrain surface; beware
# that only built-in CAD entities can be hybrid, i.e. have discrete entities on
# their boundary: OpenCASCADE does not support this feature
z = -1000
p1 = gmsh.model.geo.addPoint(xyz[0], xyz[1], z)
p2 = gmsh.model.geo.addPoint(xyz[3], xyz[4], z)
p3 = gmsh.model.geo.addPoint(xyz[6], xyz[7], z)
p4 = gmsh.model.geo.addPoint(xyz[9], xyz[10], z)
c1 = gmsh.model.geo.addLine(p1, p2)
c2 = gmsh.model.geo.addLine(p2, p3)
c3 = gmsh.model.geo.addLine(p3, p4)
c4 = gmsh.model.geo.addLine(p4, p1)
c10 = gmsh.model.geo.addLine(p1, p[0])
c11 = gmsh.model.geo.addLine(p2, p[1])
c12 = gmsh.model.geo.addLine(p3, p[2])
c13 = gmsh.model.geo.addLine(p4, p[3])
ll1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
s1 = gmsh.model.geo.addPlaneSurface([ll1])
ll3 = gmsh.model.geo.addCurveLoop([c1, c11, -c[0][1], -c10])
s3 = gmsh.model.geo.addPlaneSurface([ll3])
ll4 = gmsh.model.geo.addCurveLoop([c2, c12, -c[1][1], -c11])
s4 = gmsh.model.geo.addPlaneSurface([ll4])
ll5 = gmsh.model.geo.addCurveLoop([c3, c13, -c[2][1], -c12])
s5 = gmsh.model.geo.addPlaneSurface([ll5])
ll6 = gmsh.model.geo.addCurveLoop([c4, c10, -c[3][1], -c13])
s6 = gmsh.model.geo.addPlaneSurface([ll6])
sl1 = gmsh.model.geo.addSurfaceLoop([s1, s3, s4, s5, s6, s[0][1]])
v1 = gmsh.model.geo.addVolume([sl1])
gmsh.model.geo.synchronize()
# set this to True to build a fully hex mesh
transfinite = False
if 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)
else:
gmsh.option.setNumber('Mesh.MeshSizeMin', 100)
gmsh.option.setNumber('Mesh.MeshSizeMax', 100)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.setOrder(2) # <-- THAT IS THE ADDITIONAL LINE
gmsh.write('terrain.msh')
#gmsh.fltk.run()
gmsh.finalize() # to remove all objects from memory