geometry output speed up

gprMax/geometry_out.pyx

@@ -0,0 +1,63 @@
+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

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_out import define_fine_geometry
+
 
 class GeometryView(object):
     """Views of the geometry of the model."""
@@ -181,80 +183,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
+                # 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 = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components
                 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.tostring())
+                pbar.update(n=points.nbytes)
 
-                # 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.tostring())
+                pbar.update(n=x_lines.nbytes)
+                f.write(y_lines.tostring())
+                pbar.update(n=y_lines.nbytes)
+                f.write(z_lines.tostring())
+                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.tostring())
+                pbar.update(n=x_materials.nbytes)
+                f.write(y_materials.tostring())
+                pbar.update(n=y_materials.nbytes)
+                f.write(z_materials.tostring())
+                pbar.update(n=z_materials.nbytes)
 
                 f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))