基于python的三维射线追踪库-ttcrpy详解
Posted mica fish
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了基于python的三维射线追踪库-ttcrpy详解相关的知识，希望对你有一定的参考价值。
基于python的三维射线追踪库-ttcrpy详解（5）
ttcrpy是加拿大学者伯纳德·吉鲁（Bernard Giroux）于2021年发布的开源python库，详见（https://github.com/groupeLIAMG），参考文献（Giroux B. 2021. ttcrpy: A Python package for traveltime computation and raytracing.
SoftwareX, vol. 16, 100834. doi:10.1016/j.softx.2021.100834
）。
ttcrpy库包含了三种射线追踪方法：快速扫描算法（FSM）、最短路径法（SPM）、动节点最短路径法（DSPM）。包含其二维与三维的实现。
ttcrpy库中给出了2D矩形网格和三角形网格、3D正六面体与四面体网格等网格剖分形式，对于非规则网格，要利用python中的vtk库和pygmsh库生成，本博文研究ttcrpy中三角网格射线追踪。
文章目录
1、vtk库
Vtk，（visualization toolkit）是一个开源的免费软件系统，主要用于三维计算机图形学、图像处理和可视化。Vtk是在面向对象原理的基础上设计和实现的，它的内核是用C++构建的，包含有大约250,000行代码，2000多个类，还包含有几个转换界面，因此也可以自由的通过Java，Tcl/Tk和Python各种语言使用vtk。
1.1、vtk库的安装
安装：直接在cmd中pip即可。（正常情况下会报错，可以在vtk官网下载支持python的*.whl文件安装）
pip install vtk
csdn有不少博客介绍python vtk的安装方法，本人安装的版本是9.1.0版本（目前最新版本）。
1.2、vtk库的使用
在python vtk官网（地址：https://vtk.org/download/），有一个简单示例，代码如下：
import vtk
cone_a=vtk.vtkConeSource()
coneMapper = vtk.vtkPolyDataMapper()
coneMapper.SetInputConnection(cone_a.GetOutputPort())
coneActor = vtk.vtkActor()
coneActor.SetMapper(coneMapper)
ren1= vtk.vtkRenderer()
ren1.AddActor( coneActor )
ren1.SetBackground( 0.1, 0.2, 0.4 )
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer( ren1 )
renWin.SetSize( 600, 600 )
renWin.Render()
iren=vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.Initialize()
iren.Start()
在spyder下运行结果：
vtk库包含的方法单词较长，ctrl+V比较方便，vtk说明书里面语法讲解的很详细，较易上手，但是熟练掌握要花费很多时间，ttcrpy库只需要部分vtk功能。
2、pygmsh库
Gmsh是一个强大的网格生成工具，使用脚本语言 这是出了名的难写。
pygmsh的目标是将gmsh的强大功能与python的多功能性结合起来，并 从gmsh脚本语言中提供有用的抽象，以便您可以创建复杂的 几何图形更容易。
2.1、pygmsh库安装
直接pip安装。
pip install pygmsh
我安装的版本是7.1.17版本。
2.2、pygmsh库使用
pygmsh官方文档给出了一个简单的示例（见官网https://pypi.org/project/pygmsh/），示例代码如下：
import pygmsh
with pygmsh.geo.Geometry() as geom:
geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)
mesh = geom.generate_mesh()
# 将网格文件保存为 vtk 格式
mesh.write("test.vtk")
将生成的vtk文件用paraview打开，显示模型：
如果每次生成vtk模型，都要用paraview打开，将非常麻烦，利用python代码可以将vtk模型显示出来。
2.2.1、二维矩形网格
利用python vtk生成二维矩形网格如下：
对应的python代码如下：
# -*- coding: utf-8 -*-
"""
Created on Thu May 5 14:48:05 2022
@author: 86159
"""
'''
此代码是用vtk绘制一个平面网格
'''
# 导入vtk库
import vtk
# 设置点与网格信息
points = vtk.vtkPoints()
cells = vtk.vtkCellArray()
polydata = vtk.vtkPolyData()
mapper = vtk.vtkPolyDataMapper()
# 图的范围
rangeX = [-10,10]
rangeY = [-10,10]
intervalX = 2
intervalY = 2
for gridY in range(rangeY[0], rangeY[1] + intervalY, intervalY):
lineStart = [rangeX[0], gridY, 1.0]
lineEnd = [rangeX[1], gridY, 1.0]
pointIdStart = points.InsertNextPoint(lineStart)
pointIdEnd = points.InsertNextPoint(lineEnd)
singleLineCell = [pointIdStart, pointIdEnd]
cells.InsertNextCell(2,singleLineCell)
print(lineStart)
print(lineEnd)
print(pointIdStart)
print(pointIdEnd)
print(singleLineCell)
for gridX in range(rangeX[0], rangeX[1] + intervalX, intervalX):
lineStart = [gridX, rangeY[0], 1.0]
lineEnd = [gridX, rangeY[1], 1.0]
pointIdStart = points.InsertNextPoint(lineStart)
pointIdEnd = points.InsertNextPoint(lineEnd)
singleLineCell = [pointIdStart, pointIdEnd]
cells.InsertNextCell(2,singleLineCell)
polydata.SetLines(cells)
polydata.SetPoints(points)
mapper.SetInputData(polydata)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
ren1 = vtk.vtkRenderer()
ren1.AddActor(actor)
ren1.SetBackground(0.1,0.2,0.4)
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren1);
renWin.SetSize(300, 300);
#enderWindowInteractor
iren = vtk.vtkRenderWindowInteractor();
iren.SetRenderWindow(renWin);
style = vtk.vtkInteractorStyleTrackballCamera();
iren.SetInteractorStyle(style);
renWin.SetSize(800, 800);
renWin.SetWindowName('test ttcrpy vtk')
renWin.Render();
iren.Initialize()
iren.Start()
2.2.2、二维三角形网格
利用python vtk生成二维三角形网格如下：
对应的python代码：
# -*- coding: utf-8 -*-
"""
Created on Fri May 6 17:03:16 2022
@author: 86159
"""
import pygmsh
import vtk
with pygmsh.geo.Geometry() as geom:
geom.add_polygon(
[
[0.0, 0.0],
[1.0, 0.0],
[1.0, 1.0],
[0.0, 1.0],
],
mesh_size=0.1,
)
mesh = geom.generate_mesh()
a = mesh.points
b = mesh.cells_dict['triangle']
points = vtk.vtkPoints()
for ai in a:
points.InsertNextPoint(ai)
cells = vtk.vtkCellArray()
for bi in b:
cells.InsertNextCell(3, bi)
cellColor = vtk.vtkUnsignedCharArray()
cellColor.SetNumberOfComponents(3)
for tmp in b:
cellColor.InsertNextTuple(tmp)
cube = vtk.vtkPolyData()
cube.SetPoints(points)
cube.SetPolys(cells)
# # 2.创建Mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetColorModeToDefault()
mapper.SetInputData(cube)
# 3.创建Actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetEdgeColor(1, 0, 1)
actor.GetProperty().SetEdgeVisibility(1) # 显示边
# 4.创建Renderer
renderer = vtk.vtkRenderer()
renderer.SetBackground(0.1, 0.2, 0.4) # 背景白色
renderer.AddActor(actor) # 将actor加入
renderer.ResetCamera() # 调整显示
# 5.渲染窗口
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(renderer)
renWin.SetSize(800,800)
renWin.SetWindowName('test ttcrpy vtk')
renWin.Render()
# 6.交互
renWinInteractor = vtk.vtkRenderWindowInteractor()
renWinInteractor.SetRenderWindow(renWin)
renWinInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
renWinInteractor.Start()
2.2.3、一些其他图形
以下图形来自pygmsh官网上给出的示例，将三角形其填充颜色后显示为下列图片，具体的代码可参照pygmsh7.1.17官网，利用python vtk显示图片。
with pygmsh.geo.Geometry() as geom:
geom.add_polygon(
[
[0.0, 0.0],
[1.0, -0.2],
[1.1, 1.2],
[0.1, 0.7],
],
mesh_size=0.1,
)
mesh = geom.generate_mesh()
# ellpsoid with holes
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_max = 0.1
ellipsoid = geom.add_ellipsoid([0.0, 0.0, 0.0], [1.0, 0.7, 0.5])
cylinders = [
geom.add_cylinder([-1.0, 0.0, 0.0], [2.0, 0.0, 0.0], 0.3),
geom.add_cylinder([0.0, -1.0, 0.0], [0.0, 2.0, 0.0], 0.3),
geom.add_cylinder([0.0, 0.0, -1.0], [0.0, 0.0, 2.0], 0.3),
]
geom.boolean_difference(ellipsoid, geom.boolean_union(cylinders))
mesh = geom.generate_mesh()
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_min = 0.1
geom.characteristic_length_max = 0.1
rectangle = geom.add_rectangle([-1.0, -1.0, 0.0], 2.0, 2.0)
disk1 = geom.add_disk([-1.2, 0.0, 0.0], 0.5)
disk2 = geom.add_disk([+1.2, 0.0, 0.0], 0.5)
disk3 = geom.add_disk([0.0, -0.9, 0.0], 0.5)
disk4 = geom.add_disk([0.0, +0.9, 0.0], 0.5)
flat = geom.boolean_difference(
geom.boolean_union([rectangle, disk1, disk2]),
geom.boolean_union([disk3, disk4]),
)
geom.extrude(flat, [0, 0, 0.3])
mesh = geom.generate_mesh()
with pygmsh.geo.Geometry() as geom:
poly = geom.add_polygon(
[
[0.0, 0.0],
[2.0, 0.0],
[3.0, 1.0],
[1.0, 2.0],
[0.0, 1.0],
],
mesh_size=0.3,
)
field0 = geom.add_boundary_layer(
edges_list=[poly.curves[0]],
lcmin=0.05,
lcmax=0.2,
distmin=0.0,
distmax=0.2,
)
field1 = geom.add_boundary_layer(
nodes_list=[poly.points[2]],
lcmin=0.05,
lcmax=0.2,
distmin=0.1,
distmax=0.4,
)
geom.set_background_mesh([field0, field1], operator="Min")
mesh = geom.generate_mesh()
from math import pi
import pygmsh
with pygmsh.geo.Geometry() as geom:
poly = geom.add_polygon(
[
[+0.0, +0.5],
[-0.1, +0.1],
[-0.5, +0.0],
[-0.1, -0.1],
[+0.0, -0.5],
[+0.1, -0.1],
[+0.5, +0.0],
[+0.1, +0.1],
],
mesh_size=0.05,
)
geom.twist(
poly,
translation_axis=[0, 0, 1],
rotation_axis=[0, 0, 1],
point_on_axis=[0, 0, 0],
angle=pi / 3,
)
mesh = geom.generate_mesh()
3、ttrpy中tmesh三角网射线追踪
3.1、导入对应的库
需要导入的库有time、numpy、matplotlib.pyplot、pygmsh、ttcrpy中的tmesh、numpy、vtk等。
import time
import numpy as np
import matplotlib.pyplot as plt
import pygmsh
from ttcrpy import tmesh
from numpy.random import default_rng
import vtk
3.2、定义vtk中显示网格的函数
此部分代码如下，基本上复制上述代码即可。
def draw_vtk(mesh):
'''
在vtk中显示网格模型
'''
a = mesh.points
b = mesh.cells_dict['triangle']
points = vtk.vtkPoints()
for ai in a:
points.InsertNextPoint(ai)
cells = vtk.vtkCellArray()
for bi in b:
cells.InsertNextCell(3, bi)
cellColor = vtk.vtkUnsignedCharArray()
cellColor.SetNumberOfComponents(3)
for tmp in b:
cellColor.InsertNextTuple(tmp)
cube = vtk.vtkPolyData()
cube.SetPoints(points)
cube.SetPolys(cells)
# # 2.创建Mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetColorModeToDefault()
mapper.SetInputData(cube)
# 3.创建Actor
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetEdgeColor(1, 0, 1)
actor.GetProperty().SetEdgeVisibility(1) # 显示边
# 4.创建Renderer
renderer = vtk.vtkRenderer()
renderer.SetBackground(1, 0, 0) # 背景白色
renderer.AddActor(actor) # 将actor加入
renderer.ResetCamera() # 调整显示
# 5.渲染窗口
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(renderer)
renWin.SetSize(600,600)
renWin.SetWindowName('meshio example')
renWin.Render()
# 6.交互
renWinInteractor = vtk.vtkRenderWindowInteractor()
renWinInteractor.SetRenderWindow(renWin)
renWinInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
renWinInteractor.Start()
3.3、设置网格
设置一个大小为1×1×1的三维立体网格，三角形单元的边长为0.042，显示图形
with pygmsh.geo.Geometry() as geom:
geom.add_box(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
mesh_size = 0.042)
msh = geom.generate_mesh()
# 显示模型
draw_vtk(msh)
模型显示结果如下所示：
3.4、设置慢度模型
根据ttcrpy官方文档中示例，设置模型参数如下：
def f(z):
return 1.0 / (1.5 + 4.5*z)
slowness = f(msh.points[:,2])
rng = default_rng(1966)
src = rng.uniform(0.01, 0.99, (10, 3)) # 10 src points
rcv = rng.uniform(0.05, 0.95, (20, 3)) # 20 receivers
# src & rcv must have the same number of rows, with rows corresponding to src-rcv pairs
src = np.kron(np.ones((20, 1)), src)
rcv = np.kron(rcv, np.ones((10, 1)))
3.5、射线追踪并计算走时
运行以下代码，输出计算得到的走时。
tm = tmesh.Mesh3d(msh.points,msh.cells_dict['tetra'],
n_threads = 1, cell_slowness=0, method='DSPM',
gradient_method = 1, tt_from_rp = 1,
process_vel = False,
maxit = 20, min_dist = 1e-5,
n_secondary = 2, n_tertiary = 2,
radius_factor_tertiary=3, translate_grid = False)
tm.set_slowness(slowness)
tstart = time.time()
tt1 = tm.raytrace(src, rcv)
tend1 = time.time() - tstart
print('time serial = ', tend1)
计算结果如下
计算耗时：time serial = 13.381624937057495
走时曲线
以上就是完整过程了，详见官方文档。
4、二维tmesh
创建二维三角网格的代码2D_maker如下：
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 26 13:54:30 2021
@author: giroux
"""
import numpy as np
import vtk
from math import floor
import gmsh
N = 20
dx = 1.0
# Rcv
with open('rcv2d.dat', 'w') as f:
xx = np.arange(0, N+1)
print('0:d'.format(xx.size*xx.size), file=f)
for x in xx:
for z in xx:
print('0:d 1:d'.format(x, z), file=f)
with open('rcv2d_in.dat', 'w') as f:
xx = np.arange(1, N)
print('0:d'.format(xx.size*xx.size), file=f)
for x in xx:
for z in xx:
print('0:d 1:d'.format(x, z), file=f)
a = 1.0
V20 = 3.0
b = (V20-a)/20.0
#
# Coarse grid
#
frac = 2
N *= frac
dx /= frac
x = np.arange(0.0, 20.0+0.01, dx)
z = np.arange(0.0, 20.0+0.01, dx)
# %% gradient
ugrid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()
slowness = vtk.vtkDoubleArray()
slowness.SetName('Slowness')
for x in np.arange(0.0, 20.0+0.01, dx):
for z in np.arange(0.0, 20.0+0.01, dx):
pts.InsertNextPoint(x, 0.0, z)
s = 1.0 / (a+b*z)
slowness.InsertNextValue(s)
ugrid.SetPoints(pts)
ugrid.GetPointData().SetScalars(slowness)
tri = vtk.vtkTriangle()
Np = N + 1
for i in np.arange(N):
for j in np.arange(N):
tri.GetPointIds().SetId(0, i*Np+j)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, i*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
tri.GetPointIds().SetId(0, i*Np+j+1)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, (i+1)*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName('gradient_coarse2d.vtu')
writer.SetInputData( ugrid )
writer.SetDataModeToBinary();
writer.Update()
xCoords = vtk.vtkDoubleArray()
for n in range(0, N+1):
xCoords.InsertNextValue( n*dx )
yCoords = vtk.vtkDoubleArray()
yCoords.InsertNextValue(0.0)
zCoords = xCoords
rgrid = vtk.vtkRectilinearGrid()
rgrid.SetDimensions( N+1, 1, N+1 )
rgrid.SetXCoordinates(xCoords)
rgrid.SetYCoordinates(yCoords)
rgrid.SetZCoordinates(zCoords)
slowness = vtk.vtkDoubleArray()
slowness.SetName("Slowness")
slowness.SetNumberOfComponents(1)
slowness.SetNumberOfTuples( rgrid.GetNumberOfPoints() )
for n in range(0, rgrid.GetNumberOfPoints()):
x = rgrid.GetPoint(n)
s = 1.0 / (a+b*x[2])
slowness.SetTuple1(n, s)
rgrid.GetPointData().SetScalars( slowness );
writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetFileName( "gradient_coarse2d.vtr" )
writer.SetInputData( rgrid )
writer.SetDataModeToBinary()
writer.Update()
# %% layers
ugrid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()
for x in np.arange(0.0, 20.0+0.01, dx):
for z in np.arange(0.0, 20.0+0.01, dx):
pts.InsertNextPoint(x, 0.0, z)
ugrid.SetPoints(pts)
tri = vtk.vtkTriangle()
slowness = vtk.vtkDoubleArray()
slowness.SetName('Slowness')
Np = N + 1
for i in np.arange(N):
for j in np.arange(N):
tri.GetPointIds().SetId(0, i*Np+j)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, i*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
z = floor(j*dx + 0.001) + 0.5
s = 1.0 / (a+b*z)
slowness.InsertNextValue(s)
tri.GetPointIds().SetId(0, i*Np+j+1)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, (i+1)*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
slowness.InsertNextValue(s)
ugrid.GetCellData().SetScalars(slowness)
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName('layers_coarse2d.vtu')
writer.SetInputData( ugrid )
writer.SetDataModeToBinary();
writer.Update()
xCoords = vtk.vtkDoubleArray()
for n in range(0, N+1):
xCoords.InsertNextValue( n*dx )
yCoords = vtk.vtkDoubleArray()
yCoords.InsertNextValue(0.0)
zCoords = xCoords
rgrid = vtk.vtkRectilinearGrid()
rgrid.SetDimensions( N+1, 1, N+1 )
rgrid.SetXCoordinates(xCoords)
rgrid.SetYCoordinates(yCoords)
rgrid.SetZCoordinates(zCoords)
slowness = vtk.vtkDoubleArray()
slowness.SetName("Slowness")
slowness.SetNumberOfComponents(1)
slowness.SetNumberOfTuples( rgrid.GetNumberOfCells() )
for n in range(0, rgrid.GetNumberOfCells()):
bo = rgrid.GetCell(n).GetBounds()
z = floor(bo[4]) + 0.5
s = 1.0 / (a+b*z)
slowness.SetTuple1(n, s)
rgrid.GetCellData().SetScalars( slowness );
writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetFileName( "layers_coarse2d.vtr" )
writer.SetInputData( rgrid )
writer.SetDataModeToBinary()
writer.Update()
# %%
# Finer grid
#
frac = 2.5
N = int(frac*N + 0.000000001)
dx /= frac
x = np.arange(0.0, 20.0+0.01, dx)
z = np.arange(0.0, 20.0+0.01, dx)
# %% gradient
ugrid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()
slowness = vtk.vtkDoubleArray()
slowness.SetName('Slowness')
for x in np.arange(0.0, 20.0+0.01, dx):
for z in np.arange(0.0, 20.0+0.01, dx):
pts.InsertNextPoint(x, 0.0, z)
s = 1.0 / (a+b*z)
slowness.InsertNextValue(s)
ugrid.SetPoints(pts)
ugrid.GetPointData().SetScalars(slowness)
tri = vtk.vtkTriangle()
Np = N + 1
for i in np.arange(N):
for j in np.arange(N):
tri.GetPointIds().SetId(0, i*Np+j)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, i*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
tri.GetPointIds().SetId(0, i*Np+j+1)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, (i+1)*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName('gradient_fine2d.vtu')
writer.SetInputData( ugrid )
writer.SetDataModeToBinary();
writer.Update()
xCoords = vtk.vtkDoubleArray()
for n in range(0, N+1):
xCoords.InsertNextValue( n*dx )
yCoords = vtk.vtkDoubleArray()
yCoords.InsertNextValue(0.0)
zCoords = xCoords
rgrid = vtk.vtkRectilinearGrid()
rgrid.SetDimensions( N+1, 1, N+1 )
rgrid.SetXCoordinates(xCoords)
rgrid.SetYCoordinates(yCoords)
rgrid.SetZCoordinates(zCoords)
slowness = vtk.vtkDoubleArray()
slowness.SetName("Slowness")
slowness.SetNumberOfComponents(1)
slowness.SetNumberOfTuples( rgrid.GetNumberOfPoints() )
for n in range(0, rgrid.GetNumberOfPoints()):
x = rgrid.GetPoint(n)
s = 1.0 / (a+b*x[2])
slowness.SetTuple1(n, s)
rgrid.GetPointData().SetScalars( slowness );
writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetFileName( "gradient_fine2d.vtr" )
writer.SetInputData( rgrid )
writer.SetDataModeToBinary()
writer.Update()
# %% layers
ugrid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()
for x in np.arange(0.0, 20.0+0.01, dx):
for z in np.arange(0.0, 20.0+0.01, dx):
pts.InsertNextPoint(x, 0.0, z)
ugrid.SetPoints(pts)
tri = vtk.vtkTriangle()
slowness = vtk.vtkDoubleArray()
slowness.SetName('Slowness')
Np = N + 1
for i in np.arange(N):
for j in np.arange(N):
tri.GetPointIds().SetId(0, i*Np+j)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, i*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
z = floor(j*dx + 0.001) + 0.5
s = 1.0 / (a+b*z)
slowness.InsertNextValue(s)
tri.GetPointIds().SetId(0, i*Np+j+1)
tri.GetPointIds().SetId(1, (i+1)*Np+j)
tri.GetPointIds().SetId(2, (i+1)*Np+j+1)
ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
slowness.InsertNextValue(s)
ugrid.GetCellData().SetScalars(slowness)
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName('layers_fine2d.vtu')
writer.SetInputData( ugrid )
writer.SetDataModeToBinary();
writer.Update()
xCoords = vtk.vtkDoubleArray()
for n in range(0, N+1):
xCoords.InsertNextValue( n*dx )
yCoords = vtk.vtkDoubleArray()
yCoords.InsertNextValue(0.0)
zCoords = xCoords
rgrid = vtk.vtkRectilinearGrid()
rgrid.SetDimensions( N+1, 1, N+1 )
rgrid.SetXCoordinates(xCoords)
rgrid.SetYCoordinates(yCoords)
rgrid.SetZCoordinates(zCoords)
slowness = vtk.vtkDoubleArray()
slowness.SetName("Slowness")
slowness.SetNumberOfComponents(1)
slowness.SetNumberOfTuples( rgrid.GetNumberOfCells() )
for n in range(0, rgrid.GetNumberOfCells()):
bo = rgrid.GetCell(n).GetBounds()
z = floor(bo[4]) + 0.5
s = 1.0 / (a+b*z)
slowness.SetTuple1(n, s)
rgrid.GetCellData().SetScalars( slowness );
writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetFileName( "layers_fine2d.vtr" )
writer.SetInputData( rgrid )
writer.SetDataModeToBinary()
writer.Update()
# %%
gmsh.initialize()
# %%
meshSize = 0.43
gmsh.clear()
tag = 1
# for x in range(21):
# for z in range(21):
# gmsh.model.geo.addPoint(x, 0, z, meshSize=meshSize, tag=tag)
# tag += 1
# pts = []
# for k in range(21):
# pts.append(k)
# for i in range(1, 21):
# pts.append(i*21+20)
# for k in range(19, -1, -1):
# pts.append(20*21 + k)
# for i in range(19, -1, -1):
# pts.append(i*21)
gmsh.model.geo.addPoint(0, 0, 0, meshSize=meshSize, tag=1)
gmsh.model.geo.addPoint(0, 0, 20, meshSize=meshSize, tag=2)
gmsh.model.geo.addPoint(20, 0, 20, meshSize=meshSize, tag=3)
gmsh.model.geo.addPoint(20, 0, 0, meshSize=meshSize, tag=4)
# add a few points "randomly" for have a more irregular mesh
gmsh.model.geo.addPoint(13.43, 0, 3.3312, meshSize=meshSize, tag=5)
gmsh.model.geo.addPoint(3.43, 0, 13.3312, meshSize=meshSize, tag=6)
gmsh.model.geo.addPoint(8.43, 0, 8.3312, meshSize=meshSize, tag=7)
gmsh.model.geo.addPoint(17.43, 0, 7.3312, meshSize=meshSize, tag=8)
gmsh.model.geo.addPoint(1.43, 0, 17.12, meshSize=meshSize, tag=9)
gmsh.model.geo.addPoint(1.43, 0, 1.3312, meshSize=meshSize, tag=10)
gmsh.model.geo.addPoint(9.43, 0, 15.3312, meshSize=meshSize, tag=11)
pts = (0, 1, 2, 3, 0)
lines = []
for n in range(len(pts)-1):
lines.append(gmsh.model.geo.addLine(pts[n]+1, pts[n+1]+1))
loop = gmsh.model.geo.addCurveLoop(lines)
srf = gmsh.model.geo.addPlaneSurface((loop,))
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(0, [5, 6, 7, 8, 9, 10, 11], 2, srf)
gmsh.model.geo.synchronize()
gmsh.model.geo.addPhysicalGroup(2, (srf,))
gmsh.model.mesh.generate(2)
gmsh.write('s_mesh.vtk')
# %%
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName('s_mesh.vtk')
reader.Update()
slowness = vtk.vtkDoubleArray()
slowness.SetName('Slowness')
ugrid = reader.GetOutput()
for n in range(ugrid.GetNumberOfPoints()):
x, y, z = ugrid.GetPoint(n)
s = 1.0 / (a+b*z)
slowness.InsertNextValue(s)
ugrid.GetPointData().SetScalars(slowness)
ugrid2 = vtk.vtkUnstructuredGrid()
ugrid2.SetPoints(ugrid.GetPoints())
ugrid2.GetPointData().SetScalars(slowness)
for n in range(ugrid.GetNumberOfCells()):
cell = ugrid.GetCell(n)
if type(cell) is vtk.vtkTriangle:
ugrid2.InsertNextCell(cell.GetCellType(), cell.GetPointIds())
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetFileName('gradient_gmsh2d.vtu')
writer.SetInputData( ugrid2 )
writer.SetDataModeToBinary();
writer.Update()
# %%
# nodeTags, nodeCoords, nodeParams = gmsh.model.mesh.getNodes()
# elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(2)
# elemTags = elemTags[0]
# elemNodeTags = elemNodeTags[0]
# nodeCoords = np.reshape(nodeCoords, (-1, 3))
# elemNodeTags = np.reshape(elemNodeTags, (-1, 3))
# ugrid = vtk.vtkUnstructuredGrid()
# pts = vtk.vtkPoints()
# slowness = vtk.vtkDoubleArray()
# slowness.SetName('Slowness')
# for n in range(nodeCoords.shape[0]):
# x = nodeCoords[n, 0]
# z = nodeCoords[n, 2]
# pts.InsertNextPoint(x, 0.0, z)
# s = 1.0 / (a+b*z)
# slowness.InsertNextValue(s)
# ugrid.SetPoints(pts)
# ugrid.GetPointData().SetScalars(slowness)
# tri = vtk.vtkTriangle()
# for n in np.arange(elemNodeTags.shape[0]):
# tri.GetPointIds().SetId(0, int(elemNodeTags[n, 0]-1))
# tri.GetPointIds().SetId(1, int(elemNodeTags[n, 1]-1))
# tri.GetPointIds().SetId(2, int(elemNodeTags[n, 2]-1))
# ugrid.InsertNextCell( tri.GetCellType(), tri.GetPointIds() )
# writer = vtk.vtkXMLUnstructuredGridWriter()
# writer.SetFileName('gradient_gmsh2d.vtu')
# writer.SetInputData( ugrid )
# writer.SetDataModeToBinary();
# writer.Update()
2D模型的测试代码如下：
# -*- coding: utf-8 -*-
"""Tests for verifying python wrappers, module tmesh in 2D"""
import unittest
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import ttcrpy.tmesh as tm
def get_tt(filename):
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
data = reader.GetOutput()
names = ('Travel Time', 'Travel time', 'travel time')
for name in names:
if data.GetPointData().HasArray(name):
break
tt = vtk_to_numpy(data.GetPointData().GetArray(name))
return tt
class TestMesh2dc(unittest.TestCase):
def setUp(self):
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName('layers_fine2d.vtu')
reader.Update()
self.nodes = np.empty((reader.GetOutput().GetNumberOfPoints(), 2 ))
for n in range(reader.GetOutput().GetNumberOfPoints()):
x = reader.GetOutput().GetPoint(n)
self.nodes[n, 0] = x[0]
self.nodes[n, 1] = x[2]
self.tri = np.empty((reader.GetOutput().GetNumberOfCells(), 3 ),
dtype=int)
ind = vtk.vtkIdList()
for n in range(reader.GetOutput().GetNumberOfCells()):
reader.GetOutput().GetCellPoints(n, ind)
self.tri[n, 0] = ind.GetId(0)
self.tri[n, 1] = ind.GetId(1)
self.tri[n, 2] = ind.GetId(2)
data = reader.GetOutput()
self.slowness = vtk_to_numpy(data.GetCellData().GetArray('Slowness'))
self.src = np.loadtxt('src2d.dat',skiprows=1)
self.src = self.src.reshape((1, 3))
self.rcv = np.loadtxt('rcv2d.dat',skiprows=1)
def test_Mesh2Dfs(self):
g = tm.Mesh2d(self.nodes, self.tri.astype('int64'), method='FSM')
tt = g.raytrace(self.src, self.rcv, slowness=self.slowness)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('Grid2Ducfs_tt_grid.vtu')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'FSM accuracy failed (slowness in cells)')
def test_Mesh2Dsp(self):
g = tm.Mesh2d(self.nodes, self.tri.astype('int64'), method='SPM', n_secondary=10)
tt = g.raytrace(self.src, self.rcv, slowness=self.slowness)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('Grid2Ducsp_tt_grid.vtu')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'SPM accuracy failed (slowness in cells)')
def test_Mesh2Ddsp(self):
g = tm.Mesh2d(self.nodes, self.tri.astype('int64'), method='DSPM', n_secondary=3,
n_tertiary=3, radius_factor_tertiary=3.0)
tt = g.raytrace(self.src, self.rcv, slowness=self.slowness)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('Grid2Ducdsp_tt_grid.vtu')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'SPM accuracy failed (slowness in cells)')
class TestMesh2dn(unittest.TestCase):
def setUp(self):
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName('gradient_fine2d.vtu')
reader.Update()
self.nodes = np.empty((reader.GetOutput().GetNumberOfPoints(), 2 ))
for n in range(reader.GetOutput().GetNumberOfPoints()):
x = reader.GetOutput().GetPoint(n)
self.nodes[n, 0] = x[0]
self.nodes[n, 1] = x[2]
self.tri = np.empty((reader.GetOutput().GetNumberOfCells(), 3 ),
dtype=int)
ind = vtk.vtkIdList()
for n in range(reader.GetOutput().GetNumberOfCells()):
reader.GetOutput().GetCellPoints(n, ind)
self.tri[n, 0] = ind.GetId(0)
self.tri[n, 1] = ind.GetId(1)
self.tri[n, 2] = ind.GetId(2)
data = reader.GetOutput()
self.slowness = vtk_to_numpy(data.GetPointData().GetArray('Slowness'))
self.src = np.loadtxt('src2d.dat',skiprows=1)
self.src = self.src.reshape((1, 3))
self.rcv = np.loadtxt('rcv2d.dat',skiprows=1)
def test_Mesh2Dfs(self):
g = tm.Mesh2d(self.nodes, self.tri.astype('int64'), method='FSM', cell_slowness=0)
tt = g.raytrace(self.src, self.rcv, slowness=self.slowness)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('Grid2Dunfs_tt_grid.vtu')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'FSM accuracy failed (slowness at nodes)')
def test_Mesh2Dsp(self):
g = tm.Mesh2d(self.nodes, self.tri.astype('int64'), method='SPM', n_secondary=10,
cell_slowness=0)
tt = g.raytrace(self.src, self.rcv, slowness=self.slowness)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('Grid2Dunsp_tt_grid.vtu')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'SPM accuracy failed (slowness at nodes)')
def test_Mesh2Ddsp(self):
g = tm.Mesh2d(self.nodes, self.tri.astype('int64'), method='DSPM', n_secondary=3,
n_tertiary=3, radius_factor_tertiary=3.0, cell_slowness=0)
tt = g.raytrace(self.src, self.rcv, slowness=self.slowness)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('Grid2Dundsp_tt_grid.vtu')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'SPM accuracy failed (slowness at nodes)')
class Data_kernel(unittest.TestCase):
def test_2d(self):
V = np.ones((11, 13))
V[:, 7:] = 2
slowness = 1. / V.flatten()
grx = np.arange(12.)
grz = np.arange(14.)
z = 0.5 + np.arange(13.)
Tx = np.vstack((0.5+np.zeros((13,)),
z)).T
Rx = np.vstack((10.5+np.zeros((13,)),
z)).T
nTx = Tx.shape[0]
nRx = Rx.shape[0]
Tx = np.kron(Tx, np.ones((nRx,1)))
Rx = np.kron(np.ones((nTx,1)), Rx)
L = tm.Mesh2d.data_kernel_straight_rays(Tx, Rx, grx, grz)
tt = L.dot(slowness)
tt2 = np.zeros(tt.shape)
d = np.sqrt(np.sum((Tx-Rx)**2, axis=1))
ind = np.logical_and(Tx[:,1]>7, Rx[:,1]>7)
tt2[ind] = d[ind]/2
ind2 = np.logical_and(Tx[:,1]<7, Rx[:,1]<7)
tt2[ind2] = d[ind2]
ind3 = np.logical_and(np.logical_not(ind), np.logical_not(ind2))
f = (7-Tx[ind3,1]) / (Rx[ind3,1]-Tx[ind3,1])
ind = (Rx[ind3,1]-Tx[ind3,1]) < 0
f[ind] = 1-f[ind]
tt2[ind3] = d[ind3]*f + d[ind3]*(1-f)/2
self.assertAlmostEqual(np.sum(np.abs(tt-tt2)), 0.0 )
if __name__ == '__main__':
unittest.main()
4、参考资料
Giroux B. 2021. ttcrpy: A Python package for traveltime computation and raytracing.
SoftwareX, vol. 16, 100834. doi: 10.1016/j.softx.2021.100834
https://github.com/groupeLIAMG/ttcr
搬砖不易，走过路过，点个赞可好
以上是关于基于python的三维射线追踪库-ttcrpy详解的主要内容，如果未能解决你的问题，请参考以下文章