From c2dd4ba3480555b2a0c1580a1feb205174397c1e Mon Sep 17 00:00:00 2001 From: jasminium Date: Fri, 15 Dec 2017 19:54:46 +0000 Subject: [PATCH] geometry output speed up --- gprMax/geometry_out.pyx | 63 ++++++++++++++++++++ gprMax/geometry_outputs.py | 117 ++++++++++++++++++++----------------- 2 files changed, 125 insertions(+), 55 deletions(-) create mode 100644 gprMax/geometry_out.pyx diff --git a/gprMax/geometry_out.pyx b/gprMax/geometry_out.pyx new file mode 100644 index 00000000..3c4f3ffa --- /dev/null +++ b/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 diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index e2357a87..d72330f0 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_out import define_fine_geometry + class GeometryView(object): """Views of the geometry of the model.""" @@ -181,80 +183,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 + # 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\n'.encode('utf-8'))