default voxel and optional edge writing

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.

gprMax/xdmf.py

@@ -39,7 +39,30 @@ class Coordinates:
         self.coordinate_count += 1
 
 
-def hexCellPicker(grid, i, j, k):
+class Solids:
+
+    def __init__(self, fdtd_grid):
+        self.count = 0
+        self.fdtd_grid = fdtd_grid
+        self.total_solids = fdtd_grid.n_cells()
+        self.solids = np.zeros((self.total_solids), np.float32)
+
+    def add_solid(self, i, j, k):
+
+        self.solids[self.count] = self.fdtd_grid.solid[i][j][k]
+        self.count += 1
+
+
+class SolidLabels():
+
+    def __init__(self, label_grid):
+        self.count = 0
+        self.label_grid = label_grid
+        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.
 
@@ -73,32 +96,9 @@ def hexCellPicker(grid, i, j, k):
 
         return cell
 
-
-class Solids:
-
-    def __init__(self, fdtd_grid):
-        self.count = 0
-        self.fdtd_grid = fdtd_grid
-        self.total_solids = fdtd_grid.n_cells()
-        self.solids = np.zeros((self.total_solids), np.float32)
-
-    def add_solid(self, i, j, k):
-
-        self.solids[self.count] = self.fdtd_grid.solid[i][j][k]
-        self.count += 1
-
-
-class SolidLabels():
-
-    def __init__(self, label_grid):
-        self.count = 0
-        self.label_grid = label_grid
-        self.total_solids = label_grid.n_cells()
-        self.solid_labels = np.zeros((self.total_solids, 8), np.float32)
-
     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,8 +158,12 @@ def process_grid(fdtd_grid):
 
                 label = labels[i][j][k]
 
-                # Each vertex can have varying numbers of edges
+                if i < i_max and j < j_max and k < k_max:
+                    solids.add_solid(i, j, k)
+                    solid_labels.add(i, j, k)
 
+                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)
@@ -170,8 +175,6 @@ def process_grid(fdtd_grid):
                         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:
@@ -219,18 +222,22 @@ def process_grid(fdtd_grid):
                 # 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)
 
@@ -248,13 +255,16 @@ 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")
-        data.create_dataset('materials', data=options['edge_materials'].materials)
+
+    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)
+        coords.create_dataset('connectivity', data=options['edges'].edges)
+
 
 def create_xdmf_markup(options):
 
@@ -264,6 +274,15 @@ def create_xdmf_markup(options):
     domain_el = etree.Element("Domain")
     xdmf_el.append(domain_el)
 
+    geometry_el = etree.Element("Geometry", GeometryType="XYZ")
+    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)
+
+    # Check if there are edges to write
+    if 'edges' in options:
+
         grid_el = etree.Element("Grid", Name="Edges", GridType="Uniform")
         domain_el.append(grid_el)
 
@@ -277,14 +296,9 @@ def create_xdmf_markup(options):
         topology_el.append(top_data_el)
 
         # Create the Geometry node
-    geometry_el = etree.Element("Geometry", GeometryType="XYZ")
-    grid_el.append(geometry_el)
+        grid_el.append(copy.deepcopy(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")