From c2dd4ba3480555b2a0c1580a1feb205174397c1e Mon Sep 17 00:00:00 2001 From: jasminium Date: Fri, 15 Dec 2017 19:54:46 +0000 Subject: [PATCH 1/3] 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')) From 6ab0eafe4a6c9c5ecf1864ffab6052094797b690 Mon Sep 17 00:00:00 2001 From: jasminium Date: Thu, 21 Dec 2017 11:39:42 +0000 Subject: [PATCH 2/3] voxel building --- gprMax/geometry_outputs.py | 94 ++++++++++++++++++-------------- gprMax/geometry_outputs_ext.pyx | 96 +++++++++++++++++++++++++++++++++ 2 files changed, 150 insertions(+), 40 deletions(-) create mode 100644 gprMax/geometry_outputs_ext.pyx diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index d72330f0..ed72df8f 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -26,7 +26,7 @@ from struct import pack from gprMax._version import __version__ from gprMax.utilities import round_value -from .geometry_out import define_fine_geometry +from .geometry_outputs_ext import define_fine_geometry, define_normal_geometry class GeometryView(object): @@ -64,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) @@ -135,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')) @@ -225,21 +239,21 @@ class GeometryView(object): z_materials) # Write point data - datasize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components + datasize = points.nbytes f.write(pack('I', datasize)) - f.write(points.tostring()) - pbar.update(n=points.nbytes) + f.write(points) + pbar.update(n=datasize) # Write connectivity data datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components f.write(pack('I', datasize)) pbar.update(n=4) - f.write(x_lines.tostring()) + f.write(x_lines) pbar.update(n=x_lines.nbytes) - f.write(y_lines.tostring()) + f.write(y_lines) pbar.update(n=y_lines.nbytes) - f.write(z_lines.tostring()) + f.write(z_lines) pbar.update(n=z_lines.nbytes) # Write cell type (line) offsets @@ -256,11 +270,11 @@ class GeometryView(object): f.write(pack('I', datasize)) pbar.update(n=4) - f.write(x_materials.tostring()) + f.write(x_materials) pbar.update(n=x_materials.nbytes) - f.write(y_materials.tostring()) + f.write(y_materials) pbar.update(n=y_materials.nbytes) - f.write(z_materials.tostring()) + 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 + + + From d6c718588252d96e8b0e188ac53c1d4ac6589b30 Mon Sep 17 00:00:00 2001 From: jasminium Date: Thu, 21 Dec 2017 11:57:50 +0000 Subject: [PATCH 3/3] removed file --- gprMax/geometry_out.pyx | 63 ----------------------------------------- 1 file changed, 63 deletions(-) delete mode 100644 gprMax/geometry_out.pyx diff --git a/gprMax/geometry_out.pyx b/gprMax/geometry_out.pyx deleted file mode 100644 index 3c4f3ffa..00000000 --- a/gprMax/geometry_out.pyx +++ /dev/null @@ -1,63 +0,0 @@ -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