diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py
index eefcf989..b0a12e9b 100644
--- a/gprMax/geometry_outputs.py
+++ b/gprMax/geometry_outputs.py
@@ -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('\n'.format(vtk_rxs_offset).encode('utf-8'))
f.write('\n'.encode('utf-8'))
f.write('\n\n\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\n'.encode('utf-8'))
@@ -181,80 +197,85 @@ class GeometryView(object):
f.write('\n'.encode('utf-8'))
f.write('\n\n\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\n'.encode('utf-8'))
diff --git a/gprMax/geometry_outputs_ext.pyx b/gprMax/geometry_outputs_ext.pyx
new file mode 100644
index 00000000..c1182a81
--- /dev/null
+++ b/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
+
+
+