0

当与高阶元素(地质力学模拟所需)进行网格划分时,以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
4

0 回答 0