diff --git a/gprMax/geometry_views.py b/gprMax/geometry_views.py index 12fe2574..12a34cde 100644 --- a/gprMax/geometry_views.py +++ b/gprMax/geometry_views.py @@ -58,7 +58,7 @@ class GeometryView(object): def write_xdmf(self, modelrun, numbermodelruns, G): filename = self.filename[:-4] - write_output_file(filename, G) + write_output_file(filename, G, self.type) def write_vtk(self, modelrun, numbermodelruns, G): """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. diff --git a/gprMax/xdmf.py b/gprMax/xdmf.py index e6e1eee0..ab1878d9 100644 --- a/gprMax/xdmf.py +++ b/gprMax/xdmf.py @@ -39,41 +39,6 @@ class Coordinates: self.coordinate_count += 1 -def hexCellPicker(grid, i, j, k): - """ - This is the ordering of nodes in the hexahedron cell. - - 7 --------- 6 - / /| - 4 --------- 5 2 - | 3 | / - | / |/ - 0 --------- 1 - - 0 1 2 3 4 5 6 7 - """ - - cell = [ - grid[i][j][k], - # 1 - grid[i + 1][j][k], - # 2 - grid[i + 1][j + 1][k], - # 3 - grid[i][j + 1][k], - # 4 - grid[i][j][k + 1], - # 5 - grid[i + 1][j][k + 1], - # 6 - grid[i + 1][j + 1][k + 1], - # 7 - grid[i][j + 1][k + 1] - ] - - return cell - - class Solids: def __init__(self, fdtd_grid): @@ -96,9 +61,44 @@ class SolidLabels(): self.total_solids = label_grid.n_cells() self.solid_labels = np.zeros((self.total_solids, 8), np.float32) + def hexCellPicker(self, grid, i, j, k): + + """ + This is the ordering of nodes in the hexahedron cell. + + 7 --------- 6 + / /| + 4 --------- 5 2 + | 3 | / + | / |/ + 0 --------- 1 + + 0 1 2 3 4 5 6 7 + """ + + cell = [ + grid[i][j][k], + # 1 + grid[i + 1][j][k], + # 2 + grid[i + 1][j + 1][k], + # 3 + grid[i][j + 1][k], + # 4 + grid[i][j][k + 1], + # 5 + grid[i + 1][j][k + 1], + # 6 + grid[i + 1][j + 1][k + 1], + # 7 + grid[i][j + 1][k + 1] + ] + + return cell + def add(self, i, j, k): - solid_labels = hexCellPicker(self.label_grid.grid, i, j, k) + solid_labels = self.hexCellPicker(self.label_grid.grid, i, j, k) self.solid_labels[self.count] = solid_labels self.count += 1 @@ -119,25 +119,23 @@ class Materials: self.material_count += 1 - -def process_grid(fdtd_grid): +def process_grid(fdtd_grid, res): # Dimensions of the problem domain. nx = fdtd_grid.nx ny = fdtd_grid.ny nz = fdtd_grid.nz + # useful indices + i_max = nx - 1 + j_max = ny - 1 + k_max = nz - 1 + # label each node in the space labels = np.arange(nx * ny * nz).reshape(nx, ny, nz) label_grid = Grid(labels) - # Edges define the connectivity of the grid. - edges = Edges(label_grid) - - # Material for each edge - edge_materials = Materials(fdtd_grid) - # Define coordinates for each node coordinates = Coordinates(fdtd_grid) @@ -147,9 +145,12 @@ def process_grid(fdtd_grid): # Connectivity for hexhahedron grid solid_labels = SolidLabels(label_grid) - i_max = nx - 1 - j_max = ny - 1 - k_max = nz - 1 + if res == 'f': + # Edges define the connectivity of the grid. + edges = Edges(label_grid) + + # Material for each edge + edge_materials = Materials(fdtd_grid) for i, ix in enumerate(labels): for j, jx in enumerate(ix): @@ -157,80 +158,86 @@ def process_grid(fdtd_grid): label = labels[i][j][k] - # Each vertex can have varying numbers of edges - - # Type 1 vertex if i < i_max and j < j_max and k < k_max: - edges.add_edge(label, i + 1, j, k) - edges.add_edge(label, i, j + 1, k) - edges.add_edge(label, i, j, k + 1) - - edge_materials.add_material(i, j, k, 0) - edge_materials.add_material(i, j, k, 1) - edge_materials.add_material(i, j, k, 2) - - # Only this node can support a cell solids.add_solid(i, j, k) solid_labels.add(i, j, k) - # Type 2 vertex - elif i < i_max and j == j_max and k == k_max: - edges.add_edge(label, i + 1, j, k) - edge_materials.add_material(i, j, k, 0) + if res == 'f': + # Each vertex can have varying numbers of edges + # Type 1 vertex + if i < i_max and j < j_max and k < k_max: + edges.add_edge(label, i + 1, j, k) + edges.add_edge(label, i, j + 1, k) + edges.add_edge(label, i, j, k + 1) - # Type 7 vertex - elif i < i_max and j == j_max and k < k_max: - edges.add_edge(label, i + 1, j, k) - edges.add_edge(label, i, j, k + 1) - edge_materials.add_material(i, j, k, 0) - edge_materials.add_material(i, j, k, 2) + edge_materials.add_material(i, j, k, 0) + edge_materials.add_material(i, j, k, 1) + edge_materials.add_material(i, j, k, 2) - # Type 6 vertex - elif i == i_max and j == j_max and k < k_max: - edges.add_edge(label, i, j, k + 1) - edge_materials.add_material(i, j, k, 2) + # Only this node can support a cell - # Type 5 vertex - elif i == i_max and j < j_max and k < k_max: - edges.add_edge(label, i, j, k + 1) - edges.add_edge(label, i, j + 1, k) - edge_materials.add_material(i, j, k, 2) - edge_materials.add_material(i, j, k, 1) + # Type 2 vertex + elif i < i_max and j == j_max and k == k_max: + edges.add_edge(label, i + 1, j, k) + edge_materials.add_material(i, j, k, 0) - # Type 4 vertex - elif i == i_max and j < j_max and k == k_max: - edges.add_edge(label, i, j + 1, k) - edge_materials.add_material(i, j, k, 1) + # Type 7 vertex + elif i < i_max and j == j_max and k < k_max: + edges.add_edge(label, i + 1, j, k) + edges.add_edge(label, i, j, k + 1) + edge_materials.add_material(i, j, k, 0) + edge_materials.add_material(i, j, k, 2) - # Type 8 vertex - elif i < i_max and j < j_max and k == k_max: - edges.add_edge(label, i, j + 1, k) - edges.add_edge(label, i + 1, j, k) - edge_materials.add_material(i, j, k, 1) - edge_materials.add_material(i, j, k, 0) + # Type 6 vertex + elif i == i_max and j == j_max and k < k_max: + edges.add_edge(label, i, j, k + 1) + edge_materials.add_material(i, j, k, 2) - # Type 3 vertex - # Has no new connectivity - elif i == i_max and j == j_max and k == k_max: - pass - else: - print('oh no') + # Type 5 vertex + elif i == i_max and j < j_max and k < k_max: + edges.add_edge(label, i, j, k + 1) + edges.add_edge(label, i, j + 1, k) + edge_materials.add_material(i, j, k, 2) + edge_materials.add_material(i, j, k, 1) + + # Type 4 vertex + elif i == i_max and j < j_max and k == k_max: + edges.add_edge(label, i, j + 1, k) + edge_materials.add_material(i, j, k, 1) + + # Type 8 vertex + elif i < i_max and j < j_max and k == k_max: + edges.add_edge(label, i, j + 1, k) + edges.add_edge(label, i + 1, j, k) + edge_materials.add_material(i, j, k, 1) + edge_materials.add_material(i, j, k, 0) + + # Type 3 vertex + # Has no new connectivity + elif i == i_max and j == j_max and k == k_max: + pass + else: + print('oh no') # Add the coordinates coordinates.add_coordinate(i, j, k) - return { + data = { 'coordinates': coordinates, 'solids': solids, 'solid_labels': solid_labels, - 'edges': edges, - 'edge_materials': edge_materials, } + if res == 'f': + data['edges'] = edges + data['edge_materials'] = edge_materials -def write_output_file(filename, grid): + return data - data = process_grid(grid) + +def write_output_file(filename, grid, res): + + data = process_grid(grid, res) data['filename'] = filename data['xml_doc'] = create_xdmf_markup(data) @@ -246,14 +253,17 @@ def write_xml_doc(options): def write_H5file(options): - f = h5py.File(options['filename'] + '.h5', "w") - coords = f.create_group("mesh") - coords.create_dataset('coordinates', data=options['coordinates'].coordinates) - coords.create_dataset('connectivity', data=options['edges'].edges) - coords.create_dataset('solid_connectivity', data=options['solid_labels'].solid_labels) - data = f.create_group("data") + f = h5py.File(options['filename'] + '.h5', "w") + coords = f.create_group("mesh") + data = f.create_group("data") + + coords.create_dataset('coordinates', data=options['coordinates'].coordinates) + coords.create_dataset('solid_connectivity', data=options['solid_labels'].solid_labels) + data.create_dataset('solids', data=options['solids'].solids) + + if 'edges' in options: data.create_dataset('materials', data=options['edge_materials'].materials) - data.create_dataset('solids', data=options['solids'].solids) + coords.create_dataset('connectivity', data=options['edges'].edges) def create_xdmf_markup(options): @@ -264,36 +274,40 @@ def create_xdmf_markup(options): domain_el = etree.Element("Domain") xdmf_el.append(domain_el) - grid_el = etree.Element("Grid", Name="Edges", GridType="Uniform") - domain_el.append(grid_el) - - # Create the grid node - topology_el = etree.Element("Topology", TopologyType="Polyline", NumberOfElements=str(options['edges'].total_edges)) - grid_el.append(topology_el) - - topology_dimensions = "{} 2".format(options['edges'].total_edges) - top_data_el = etree.Element("DataItem", Dimensions=topology_dimensions, NumberType="Float", Precision="8", Format="HDF") - top_data_el.text = "{}:/mesh/connectivity".format(options['filename'] + '.h5') - topology_el.append(top_data_el) - - # Create the Geometry node geometry_el = etree.Element("Geometry", GeometryType="XYZ") - grid_el.append(geometry_el) - - # Create the origin coordinates coordinates_dimensions = "{} 3".format(options['coordinates'].total_coordinates) origin_el = etree.Element("DataItem", Dimensions=coordinates_dimensions, NumberType="Float", Precision="8", Format="HDF") origin_el.text = "{}:/mesh/coordinates".format(options['filename'] + '.h5') geometry_el.append(origin_el) - # Create the materials attribute - attr_el = etree.Element("Attribute", Center="Cell", Name="Edge_Materials") - grid_el.append(attr_el) + # Check if there are edges to write + if 'edges' in options: - materials_dimensions = "{} 1".format(options['edge_materials'].material_count) - materials_el = etree.Element("DataItem", Dimensions=materials_dimensions, NumberType="Float", Precision="8", Format="HDF") - materials_el.text = "{}:/data/materials".format(options['filename'] + '.h5') - attr_el.append(materials_el) + grid_el = etree.Element("Grid", Name="Edges", GridType="Uniform") + domain_el.append(grid_el) + + # Create the grid node + topology_el = etree.Element("Topology", TopologyType="Polyline", NumberOfElements=str(options['edges'].total_edges)) + grid_el.append(topology_el) + + topology_dimensions = "{} 2".format(options['edges'].total_edges) + top_data_el = etree.Element("DataItem", Dimensions=topology_dimensions, NumberType="Float", Precision="8", Format="HDF") + top_data_el.text = "{}:/mesh/connectivity".format(options['filename'] + '.h5') + topology_el.append(top_data_el) + + # Create the Geometry node + grid_el.append(copy.deepcopy(geometry_el)) + + # Create the origin coordinates + + # Create the materials attribute + attr_el = etree.Element("Attribute", Center="Cell", Name="Edge_Materials") + grid_el.append(attr_el) + + materials_dimensions = "{} 1".format(options['edge_materials'].material_count) + materials_el = etree.Element("DataItem", Dimensions=materials_dimensions, NumberType="Float", Precision="8", Format="HDF") + materials_el.text = "{}:/data/materials".format(options['filename'] + '.h5') + attr_el.append(materials_el) v_grid_el = etree.Element("Grid", Name="Voxel", GridType="Uniform") domain_el.append(v_grid_el)