你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-06 12:36:51 +08:00
Craig Warren
@@ -26,6 +26,8 @@ from struct import pack
|
||||
from gprMax._version import __version__
|
||||
from gprMax.utilities import round_value
|
||||
|
||||
from .geometry_outputs_ext import define_fine_geometry, define_normal_geometry
|
||||
|
||||
|
||||
class GeometryView(object):
|
||||
"""Views of the geometry of the model."""
|
||||
@@ -62,16 +64,17 @@ class GeometryView(object):
|
||||
|
||||
if self.fileext == '.vti':
|
||||
# Calculate number of cells according to requested sampling for geometry view
|
||||
self.vtk_xscells = round_value(self.xs / self.dx)
|
||||
self.vtk_xfcells = round_value(self.xf / self.dx)
|
||||
self.vtk_yscells = round_value(self.ys / self.dy)
|
||||
self.vtk_yfcells = round_value(self.yf / self.dy)
|
||||
self.vtk_zscells = round_value(self.zs / self.dz)
|
||||
self.vtk_zfcells = round_value(self.zf / self.dz)
|
||||
self.vtk_nxcells = self.vtk_xfcells - self.vtk_xscells
|
||||
self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells
|
||||
self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells
|
||||
self.datawritesize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + 2 * (int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells))
|
||||
self.vtk_xscells = self.xs // self.dx
|
||||
self.vtk_xfcells = self.xf // self.dx
|
||||
self.vtk_yscells = self.ys // self.dy
|
||||
self.vtk_yfcells = self.yf // self.dy
|
||||
self.vtk_zscells = self.zs // self.dz
|
||||
self.vtk_zfcells = self.zf // self.dz
|
||||
self.vtk_nxcells = self.nx // self.dx
|
||||
self.vtk_nycells = self.ny // self.dy
|
||||
self.vtk_nzcells = self.nz // self.dz
|
||||
self.n_vtk_cells = self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells
|
||||
self.datawritesize = int(np.dtype(np.uint32).itemsize * self.n_vtk_cells) + 2 * (int(np.dtype(np.int8).itemsize * self.n_vtk_cells))
|
||||
|
||||
elif self.fileext == '.vtp':
|
||||
self.vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
|
||||
@@ -133,34 +136,47 @@ class GeometryView(object):
|
||||
f.write('<DataArray type="Int8" Name="Receivers" format="appended" offset="{}" />\n'.format(vtk_rxs_offset).encode('utf-8'))
|
||||
f.write('</CellData>\n'.encode('utf-8'))
|
||||
f.write('</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
||||
|
||||
solid_geometry = np.zeros((self.n_vtk_cells), dtype=np.uint32)
|
||||
srcs_pml_geometry = np.zeros((self.n_vtk_cells), dtype=np.int8)
|
||||
rxs_geometry = np.zeros((self.n_vtk_cells), dtype=np.int8)
|
||||
|
||||
define_normal_geometry(
|
||||
self.xs,
|
||||
self.xf,
|
||||
self.ys,
|
||||
self.yf,
|
||||
self.zs,
|
||||
self.zf,
|
||||
self.dx,
|
||||
self.dy,
|
||||
self.dz,
|
||||
G.solid,
|
||||
self.srcs_pml,
|
||||
self.rxs,
|
||||
solid_geometry,
|
||||
srcs_pml_geometry,
|
||||
rxs_geometry)
|
||||
|
||||
# Write material IDs
|
||||
datasize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||
datasize = solid_geometry.nbytes
|
||||
# Write number of bytes of appended data as UInt32
|
||||
f.write(pack('I', datasize))
|
||||
for k in range(self.zs, self.zf, self.dz):
|
||||
for j in range(self.ys, self.yf, self.dy):
|
||||
for i in range(self.xs, self.xf, self.dx):
|
||||
pbar.update(n=4)
|
||||
f.write(pack('I', G.solid[i, j, k]))
|
||||
pbar.update(n=4)
|
||||
f.write(solid_geometry)
|
||||
pbar.update(n=datasize)
|
||||
|
||||
# Write source/PML IDs
|
||||
datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||
datasize = srcs_pml_geometry.nbytes
|
||||
f.write(pack('I', datasize))
|
||||
for k in range(self.zs, self.zf, self.dz):
|
||||
for j in range(self.ys, self.yf, self.dy):
|
||||
for i in range(self.xs, self.xf, self.dx):
|
||||
pbar.update()
|
||||
f.write(pack('b', self.srcs_pml[i, j, k]))
|
||||
pbar.update(n=4)
|
||||
f.write(srcs_pml_geometry)
|
||||
pbar.update(n=datasize)
|
||||
|
||||
# Write receiver IDs
|
||||
datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||
datasize = rxs_geometry.nbytes
|
||||
f.write(pack('I', datasize))
|
||||
for k in range(self.zs, self.zf, self.dz):
|
||||
for j in range(self.ys, self.yf, self.dy):
|
||||
for i in range(self.xs, self.xf, self.dx):
|
||||
pbar.update()
|
||||
f.write(pack('b', self.rxs[i, j, k]))
|
||||
pbar.update(n=4)
|
||||
f.write(rxs_geometry)
|
||||
pbar.update(n=datasize)
|
||||
|
||||
f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
|
||||
|
||||
@@ -181,80 +197,85 @@ class GeometryView(object):
|
||||
f.write('</CellData>\n'.encode('utf-8'))
|
||||
|
||||
f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
||||
|
||||
# Coordinates of each point
|
||||
n_points = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
|
||||
points = np.zeros((n_points, 3), dtype=np.float32)
|
||||
|
||||
# Write points
|
||||
datasize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components
|
||||
# Number of x components
|
||||
n_x_lines = self.nx * (self.ny + 1) * (self.nz + 1)
|
||||
# Node connectivity. Each index contains a pair of connected x nodes
|
||||
x_lines = np.zeros((n_x_lines, 2), dtype=np.uint32)
|
||||
# Material at Ex location in Yee cell.
|
||||
x_materials = np.zeros((n_x_lines), dtype=np.uint32)
|
||||
|
||||
n_y_lines = self.ny * (self.nx + 1) * (self.nz + 1)
|
||||
y_lines = np.zeros((n_y_lines, 2), dtype=np.uint32)
|
||||
y_materials = np.zeros((n_y_lines), dtype=np.uint32)
|
||||
|
||||
n_z_lines = self.nz * (self.nx + 1) * (self.ny + 1)
|
||||
z_lines = np.zeros((n_z_lines, 2), dtype=np.uint32)
|
||||
z_materials = np.zeros((n_z_lines), dtype=np.uint32)
|
||||
|
||||
define_fine_geometry(self.nx,
|
||||
self.ny,
|
||||
self.nz,
|
||||
self.xs,
|
||||
self.xf,
|
||||
self.ys,
|
||||
self.yf,
|
||||
self.zs,
|
||||
self.zf,
|
||||
G.dx,
|
||||
G.dy,
|
||||
G.dz,
|
||||
G.ID,
|
||||
points,
|
||||
x_lines,
|
||||
x_materials,
|
||||
y_lines,
|
||||
y_materials,
|
||||
z_lines,
|
||||
z_materials)
|
||||
|
||||
# Write point data
|
||||
datasize = points.nbytes
|
||||
f.write(pack('I', datasize))
|
||||
for i in range(self.xs, self.xf + 1):
|
||||
for j in range(self.ys, self.yf + 1):
|
||||
for k in range(self.zs, self.zf + 1):
|
||||
pbar.update(n=12)
|
||||
f.write(pack('fff', i * G.dx, j * G.dy, k * G.dz))
|
||||
f.write(points)
|
||||
pbar.update(n=datasize)
|
||||
|
||||
# Write cell type (line) connectivity for x components
|
||||
# Write connectivity data
|
||||
datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components
|
||||
f.write(pack('I', datasize))
|
||||
vtk_x2 = (self.ny + 1) * (self.nz + 1)
|
||||
for vtk_x1 in range(self.nx * (self.ny + 1) * (self.nz + 1)):
|
||||
pbar.update(n=8)
|
||||
f.write(pack('II', vtk_x1, vtk_x2))
|
||||
# print('x {} {}'.format(vtk_x1, vtk_x2))
|
||||
vtk_x2 += 1
|
||||
pbar.update(n=4)
|
||||
|
||||
# Write cell type (line) connectivity for y components
|
||||
vtk_ycnt1 = 1
|
||||
vtk_ycnt2 = 0
|
||||
for vtk_y1 in range((self.nx + 1) * (self.ny + 1) * (self.nz + 1)):
|
||||
if vtk_y1 >= (vtk_ycnt1 * (self.ny + 1) * (self.nz + 1)) - (self.nz + 1) and vtk_y1 < vtk_ycnt1 * (self.ny + 1) * (self.nz + 1):
|
||||
vtk_ycnt2 += 1
|
||||
else:
|
||||
vtk_y2 = vtk_y1 + self.nz + 1
|
||||
pbar.update(n=8)
|
||||
f.write(pack('II', vtk_y1, vtk_y2))
|
||||
# print('y {} {}'.format(vtk_y1, vtk_y2))
|
||||
if vtk_ycnt2 == self.nz + 1:
|
||||
vtk_ycnt1 += 1
|
||||
vtk_ycnt2 = 0
|
||||
|
||||
# Write cell type (line) connectivity for z components
|
||||
vtk_zcnt = self.nz
|
||||
for vtk_z1 in range((self.nx + 1) * (self.ny + 1) * self.nz + (self.nx + 1) * (self.ny + 1)):
|
||||
if vtk_z1 != vtk_zcnt:
|
||||
vtk_z2 = vtk_z1 + 1
|
||||
pbar.update(n=8)
|
||||
f.write(pack('II', vtk_z1, vtk_z2))
|
||||
# print('z {} {}'.format(vtk_z1, vtk_z2))
|
||||
else:
|
||||
vtk_zcnt += self.nz + 1
|
||||
f.write(x_lines)
|
||||
pbar.update(n=x_lines.nbytes)
|
||||
f.write(y_lines)
|
||||
pbar.update(n=y_lines.nbytes)
|
||||
f.write(z_lines)
|
||||
pbar.update(n=z_lines.nbytes)
|
||||
|
||||
# Write cell type (line) offsets
|
||||
vtk_cell_pts = 2
|
||||
datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines
|
||||
f.write(pack('I', datasize))
|
||||
pbar.update(n=4)
|
||||
for vtk_offsets in range(vtk_cell_pts, (self.vtk_numline_components * self.vtk_numlines) + vtk_cell_pts, vtk_cell_pts):
|
||||
pbar.update(n=4)
|
||||
f.write(pack('I', vtk_offsets))
|
||||
pbar.update(n=4)
|
||||
|
||||
# Write material IDs per-cell-edge, i.e. from ID array
|
||||
datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines
|
||||
f.write(pack('I', datasize))
|
||||
for i in range(self.xs, self.xf):
|
||||
for j in range(self.ys, self.yf + 1):
|
||||
for k in range(self.zs, self.zf + 1):
|
||||
pbar.update(n=4)
|
||||
f.write(pack('I', G.ID[0, i, j, k]))
|
||||
pbar.update(n=4)
|
||||
|
||||
for i in range(self.xs, self.xf + 1):
|
||||
for j in range(self.ys, self.yf):
|
||||
for k in range(self.zs, self.zf + 1):
|
||||
pbar.update(n=4)
|
||||
f.write(pack('I', G.ID[1, i, j, k]))
|
||||
|
||||
for i in range(self.xs, self.xf + 1):
|
||||
for j in range(self.ys, self.yf + 1):
|
||||
for k in range(self.zs, self.zf):
|
||||
pbar.update(n=4)
|
||||
f.write(pack('I', G.ID[2, i, j, k]))
|
||||
f.write(x_materials)
|
||||
pbar.update(n=x_materials.nbytes)
|
||||
f.write(y_materials)
|
||||
pbar.update(n=y_materials.nbytes)
|
||||
f.write(z_materials)
|
||||
pbar.update(n=z_materials.nbytes)
|
||||
|
||||
f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
|
||||
|
||||
|
||||
96
gprMax/geometry_outputs_ext.pyx
普通文件
96
gprMax/geometry_outputs_ext.pyx
普通文件
@@ -0,0 +1,96 @@
|
||||
cimport numpy as np
|
||||
|
||||
from gprMax.constants cimport floattype_t
|
||||
|
||||
cpdef void define_fine_geometry(
|
||||
int nx,
|
||||
int ny,
|
||||
int nz,
|
||||
int xs,
|
||||
int xf,
|
||||
int ys,
|
||||
int yf,
|
||||
int zs,
|
||||
int zf,
|
||||
float dx,
|
||||
float dy,
|
||||
float dz,
|
||||
np.uint32_t[:, :, :, :] ID,
|
||||
floattype_t[:, :] points,
|
||||
np.uint32_t[:, :] x_lines,
|
||||
np.uint32_t[:] x_materials,
|
||||
np.uint32_t[:, :] y_lines,
|
||||
np.uint32_t[:] y_materials,
|
||||
np.uint32_t[:, :] z_lines,
|
||||
np.uint32_t[:] z_materials
|
||||
):
|
||||
|
||||
cdef Py_ssize_t i, j, k
|
||||
cdef Py_ssize_t label = 0
|
||||
cdef Py_ssize_t counter_x = 0
|
||||
cdef Py_ssize_t counter_y = 0
|
||||
cdef Py_ssize_t counter_z = 0
|
||||
|
||||
cdef int label_x, label_y, label_z
|
||||
|
||||
for i in range(xs, xf + 1):
|
||||
for j in range(ys, yf + 1):
|
||||
for k in range(zs, zf + 1):
|
||||
points[label][0] = i * dx
|
||||
points[label][1] = j * dy
|
||||
points[label][2] = k * dz
|
||||
if i < xf:
|
||||
# x connectivity
|
||||
label_x = label + (ny + 1) * (nz + 1)
|
||||
x_lines[counter_x][0] = label
|
||||
x_lines[counter_x][1] = label_x
|
||||
# material for the line
|
||||
x_materials[counter_x] = ID[0, i, j, k]
|
||||
counter_x += 1
|
||||
if j < yf:
|
||||
label_y = label + nz + 1
|
||||
y_lines[counter_y][0] = label
|
||||
y_lines[counter_y][1] = label_y
|
||||
y_materials[counter_y] = ID[1, i, j, k]
|
||||
counter_y += 1
|
||||
if k < zf:
|
||||
label_z = label + 1
|
||||
z_lines[counter_z][0] = label
|
||||
z_lines[counter_z][1] = label_z
|
||||
z_materials[counter_z] = ID[2, i, j, k]
|
||||
counter_z += 1
|
||||
|
||||
label = label + 1
|
||||
|
||||
cpdef void define_normal_geometry(
|
||||
int xs,
|
||||
int xf,
|
||||
int ys,
|
||||
int yf,
|
||||
int zs,
|
||||
int zf,
|
||||
int dx,
|
||||
int dy,
|
||||
int dz,
|
||||
np.uint32_t[:, :, :] solid,
|
||||
np.int8_t[:, :, :] srcs_pml,
|
||||
np.int8_t[:, :, :] rxs,
|
||||
np.uint32_t[:] solid_geometry,
|
||||
np.int8_t[:] srcs_pml_geometry,
|
||||
np.int8_t[:] rxs_geometry,
|
||||
|
||||
):
|
||||
|
||||
cdef Py_ssize_t i, j, k
|
||||
cdef Py_ssize_t counter = 0
|
||||
|
||||
for k in range(zs, zf, dz):
|
||||
for j in range(ys, yf, dy):
|
||||
for i in range(xs, xf, dx):
|
||||
solid_geometry[counter] = solid[i, j, k]
|
||||
srcs_pml_geometry[counter] = srcs_pml[i, j, k]
|
||||
rxs_geometry[counter] = rxs[i, j, k]
|
||||
counter = counter + 1
|
||||
|
||||
|
||||
|
||||
在新工单中引用
屏蔽一个用户