From bf800fcaebb472f15b055e95bf6a0fc8f51a2651 Mon Sep 17 00:00:00 2001 From: jasminium Date: Wed, 4 Sep 2019 17:13:21 +0100 Subject: [PATCH] multiple grid geometry outputs --- gprMax/geometry_outputs.py | 293 +++++++++++++++++++++++++++++++++++++ 1 file changed, 293 insertions(+) diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index ee0139bb..78ddf4e5 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -388,3 +388,296 @@ class GeometryObjects(object): dispersionstr += '{:g} {:g} '.format(material.tau[pole], material.alpha[pole]) dispersionstr += material.ID fmaterials.write(dispersionstr + '\n') + +class GeometryViewFineMultiGrid: + """Geometry view for all grids in the simulation.""" + + """"Slicing is not supported by this class :( - only the full extent of the grids + are output. The subgrids are output without the non-working regions If you + require domainslicing GeometryView seperately for each grid you require and + view them at once in Paraview.""" + + if sys.byteorder == 'little': + byteorder = 'LittleEndian' + else: + byteorder = 'BigEndian' + + def __init__(self, xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, G): + """ + Args: + xs, xf, ys, yf, zs, zf (int): Extent of the volume in cells. + dx, dy, dz (int): Spatial discretisation in cells. + filename (str): Filename to save to. + fileext (str): File extension of VTK file - either '.vti' for a per cell + geometry view, or '.vtp' for a per cell edge geometry view. + """ + self.G = G + self.nx = G.nx + self.ny = G.ny + self.nz = G.nz + self.basefilename = filename + self.fileext = '.vtp' + self.sg_views = [] + + self.additional_lines = 0 + self.additional_points = 0 + + for sg in G.subgrids: + # create an object to contain data relevant to the geometry processing + sg_gv = SubgridGeometryView(sg) + self.sg_views.append(sg_gv) + # total additional lines required for subgrid + self.additional_lines += sg_gv.n_total_lines + # total additional points required for subgrid + + self.additional_points += sg_gv.n_total_points + + self.vtk_numpoints = self.additional_points + (self.nx + 1) * (self.ny + 1) * (self.nz + 1) + self.vtk_numpoint_components = 3 + self.vtk_numlines = self.additional_lines + 2 * self.nx * self.ny + 2 * self.ny * self.nz + 2 * self.nx * self.nz + 3 * self.nx * self.ny * self.nz + self.nx + self.ny + self.nz + self.vtk_numline_components = 2 + self.vtk_connectivity_offset = round_value(int((self.vtk_numpoints * self.vtk_numpoint_components * np.dtype(np.float32).itemsize) + np.dtype(np.uint32).itemsize)) + self.vtk_offsets_offset = round_value(int(self.vtk_connectivity_offset + (self.vtk_numlines * self.vtk_numline_components * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)) + self.vtk_materials_offset = round_value(int(self.vtk_offsets_offset + (self.vtk_numlines * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)) + self.datawritesize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components + np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components + np.dtype(np.uint32).itemsize * self.vtk_numlines + np.dtype(np.uint32).itemsize * self.vtk_numlines + + def set_filename(self, appendmodelnumber): + """ + Construct filename from user-supplied name and model run number. + + Args: + appendmodelnumber (str): Text to append to filename. + """ + + self.filename = os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(config.general['inputfilepath'])), self.basefilename + appendmodelnumber)) + self.filename += self.fileext + + def write_vtk(self, *args): + """ + Writes the geometry information to a VTK file. Either ImageData (.vti) for a + per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view. + N.B. No Python 3 support for VTK at time of writing (03/2015) + Args: + G (class): Grid class instance - holds essential parameters describing the model. + """ + + G = self.G + + with open(self.filename, 'wb') as f: + + # refine parameters for subgrid + + f.write('\n'.encode('utf-8')) + f.write('\n'.format(GeometryViewFineMultiGrid.byteorder).encode('utf-8')) + f.write('\n\n'.format(self.vtk_numpoints, self.vtk_numlines).encode('utf-8')) + + f.write('\n\n\n'.encode('utf-8')) + f.write('\n\n'.format(self.vtk_connectivity_offset).encode('utf-8')) + f.write('\n\n'.format(self.vtk_offsets_offset).encode('utf-8')) + f.write('\n'.encode('utf-8')) + f.write('\n'.format(self.vtk_materials_offset).encode('utf-8')) + f.write('\n'.encode('utf-8')) + + f.write('\n\n\n_'.encode('utf-8')) + + # Write points + print('writing points main grid') + datasize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components + f.write(pack('I', datasize)) + for i in range(0, G.nx + 1): + for j in range(0, G.ny + 1): + for k in range(0, self.G.nz + 1): + f.write(pack('fff', i * G.dx, j * G.dy, k * G.dz)) + + for sg_v in self.sg_views: + print('writing points subgrid') + sg_v.write_points(f, G) + + n_x_lines = self.nx * (self.ny + 1) * (self.nz + 1) + x_lines = np.zeros((n_x_lines, 2), dtype=np.uint32) + 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) + + print('calculate connectivity main grid') + label = 0 + counter_x = 0 + counter_y = 0 + counter_z = 0 + for i in range(self.nx + 1): + for j in range(self.ny + 1): + for k in range(self.nz + 1): + + if i < self.nx: + # x connectivity + label_x = label + (self.ny + 1) * (self.nz + 1) + x_lines[counter_x][0] = label + x_lines[counter_x][1] = label_x + # material for the line + x_materials[counter_x] = G.ID[0, i, j, k] + counter_x += 1 + if j < self.ny: + label_y = label + self.nz + 1 + y_lines[counter_y][0] = label + y_lines[counter_y][1] = label_y + y_materials[counter_y] = G.ID[1, i, j, k] + counter_y += 1 + if k < self.nz: + label_z = label + 1 + z_lines[counter_z][0] = label + z_lines[counter_z][1] = label_z + z_materials[counter_z] = G.ID[2, i, j, k] + counter_z += 1 + + label = label + 1 + + print('calculate connectivity subgrids') + for sg_v in self.sg_views: + sg_v.populate_connectivity_and_materials(label) + # use the last subgrids label for the next view + label = sg_v.label + + datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components + f.write(pack('I', datasize)) + + f.write(x_lines.tostring()) + for sg_v in self.sg_views: + f.write(sg_v.x_s_lines.tostring()) + f.write(y_lines.tostring()) + for sg_v in self.sg_views: + f.write(sg_v.y_s_lines.tostring()) + f.write(z_lines.tostring()) + for sg_v in self.sg_views: + f.write(sg_v.z_s_lines.tostring()) + + # Write cell type (line) offsets + vtk_cell_pts = 2 + datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines + f.write(pack('I', datasize)) + for vtk_offsets in range(vtk_cell_pts, (self.vtk_numline_components * self.vtk_numlines) + vtk_cell_pts, vtk_cell_pts): + f.write(pack('I', vtk_offsets)) + + datasize = np.dtype(np.uint32).itemsize * self.vtk_numlines + f.write(pack('I', datasize)) + + f.write(x_materials.tostring()) + for sg_v in self.sg_views: + f.write(sg_v.x_s_materials.tostring()) + f.write(y_materials.tostring()) + for sg_v in self.sg_views: + f.write(sg_v.y_s_materials.tostring()) + f.write(z_materials.tostring()) + for sg_v in self.sg_views: + f.write(sg_v.z_s_materials.tostring()) + + f.write('\n\n'.encode('utf-8')) + + #self.write_gprmax_info(f, G, materialsonly=True) + + def write_gprmax_info(self, f, G, materialsonly=False): + """ + Writes gprMax specific information relating material, source, + and receiver names to numeric identifiers. + Args: + f (filehandle): VTK file. + G (class): Grid class instance - holds essential parameters describing the model. + materialsonly (boolean): Only write information on materials + """ + + f.write('\n\n\n'.encode('utf-8')) + for material in G.materials: + f.write('{}\n'.format(material.ID, material.numID).encode('utf-8')) + if not materialsonly: + f.write('1\n'.encode('utf-8')) + for index, src in enumerate(G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines): + f.write('{}\n'.format(src.ID, index + 2).encode('utf-8')) + for index, rx in enumerate(G.rxs): + f.write('{}\n'.format(rx.ID, index + 1).encode('utf-8')) + f.write('\n'.encode('utf-8')) + return None + + +class SubgridGeometryView: + def __init__(self, sg): + + self.sg = sg + # n component lines in each direction required for subgrid in the working region + n_sx_lines = sg.nwx * (sg.nwy + 1) * (sg.nwz + 1) + n_sy_lines = sg.nwy * (sg.nwx + 1) * (sg.nwz + 1) + n_sz_lines = sg.nwz * (sg.nwx + 1) * (sg.nwy + 1) + + n_total_lines = n_sx_lines + n_sy_lines + n_sz_lines + self.n_total_lines = n_total_lines.astype(np.int32) + + # n points in the the working region + n_total_points = (sg.nwx + 1) * (sg.nwy + 1) * (sg.nwz + 1) + self.n_total_points = n_total_points.astype(np.int32) + + # connectivity array. 2 labels form an x component connection + self.x_s_lines = np.zeros((n_sx_lines, 2), dtype=np.uint32) + # material array. Each index contains a material index + self.x_s_materials = np.zeros((n_sx_lines), dtype=np.uint32) + + self.y_s_lines = np.zeros((n_sy_lines, 2), dtype=np.uint32) + self.y_s_materials = np.zeros((n_sy_lines), dtype=np.uint32) + + self.z_s_lines = np.zeros((n_sz_lines, 2), dtype=np.uint32) + self.z_s_materials = np.zeros((n_sz_lines), dtype=np.uint32) + + self.label = 0 + + def write_points(self, f, G): + sg = self.sg + for i in range(sg.i0, sg.i0 + sg.nwx + 1): + for j in range(sg.j0, sg.j0 + sg.nwy + 1): + for k in range(sg.k0, sg.k0 + sg.nwz + 1): + p_x = (sg.i0 * G.dx) + ((i - sg.i0) * sg.dx) + p_y = (sg.j0 * G.dy) + ((j - sg.j0) * sg.dy) + p_z = (sg.k0 * G.dz) + ((k - sg.k0) * sg.dz) + f.write(pack('fff', p_x, p_y, p_z)) + + def populate_connectivity_and_materials(self, label): + """Label is the starting label. 0 if no other grids are present but +1 the last label used + for a multigrid view""" + sg = self.sg + self.label = label + + # counters to to index numpy arrays + counter_x = 0 + counter_y = 0 + counter_z = 0 + + for i in range(sg.nwx + 1): + for j in range(sg.nwy + 1): + for k in range(sg.nwz + 1): + i_s = i + sg.n_boundary_cells_x + j_s = j + sg.n_boundary_cells_y + k_s = k + sg.n_boundary_cells_z + if i < sg.nwx: + # x connectivity + label_x = self.label + (sg.nwy + 1) * (sg.nwz + 1) + self.x_s_lines[counter_x][0] = self.label + self.x_s_lines[counter_x][1] = label_x + # material for the line + self.x_s_materials[counter_x] = sg.ID[0, i_s, j_s, k_s] + counter_x += 1 + if j < sg.nwy: + label_y = self.label + sg.nwz + 1 + self.y_s_lines[counter_y][0] = self.label + self.y_s_lines[counter_y][1] = label_y + self.y_s_materials[counter_y] = sg.ID[1, i_s, j_s, k_s] + counter_y += 1 + if k < sg.nwz: + label_z = self.label + 1 + self.z_s_lines[counter_z][0] = self.label + self.z_s_lines[counter_z][1] = label_z + self.z_s_materials[counter_z] = sg.ID[2, i_s, j_s, k_s] + counter_z += 1 + + self.label = self.label + 1