diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py
index a1791473..7804f0c1 100644
--- a/gprMax/input_cmds_geometry.py
+++ b/gprMax/input_cmds_geometry.py
@@ -16,7 +16,10 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
+import os
+import h5py
 import numpy as np
+import xml.etree.ElementTree as ET
 
 from gprMax.constants import floattype
 from gprMax.exceptions import CmdInputError
@@ -34,10 +37,66 @@ def process_geometrycmds(geometry, G):
         
     """
 
-    for x, object in enumerate(geometry):
+    for object in geometry:
         tmp = object.split()
         
-        if tmp[0] == '#edge:':
+        if tmp[0] == '#xdmf_geometry_file:':
+            if len(tmp) != 2:
+                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly one parameter')
+            xdmffile = tmp[1]
+
+            # See if file exists at specified path and if not try input file directory
+            if not os.path.isfile(xdmffile):
+                xdmffile = os.path.join(G.inputdirectory, xdmffile)
+
+            # Read Xdmf file for information on dataset
+            tree = ET.parse(xdmffile)
+            root = tree.getroot()
+            print(root.tag)
+            data = root.findall('./UnstructuredGrid/Piece/CellData/DataArray')
+            data = data[0].text
+            data = data.split()
+            print(data)
+
+            # Find Geometry element and get information on origin and spatial sampling of data
+#            geoinfo = root.findall('./Domain/Grid/Geometry/DataItem')    
+#            origin = geoinfo[0].text.split()
+#            dxdydz = geoinfo[1].text.split()
+#            origin = [float(item) for item in origin]
+#            dxdydz = [float(item) for item in dxdydz]
+
+            # Get name of dataset stored in HDF5 file
+#            datainfo = root.find('./Domain/Grid/Attribute/DataItem')
+#            datainfo = datainfo.text
+#            datasetname = datainfo.split(':')[1]
+
+            # Load HDF5 file and dataset
+#            xdmfh5file = os.path.splitext(xdmffile)[0] + '_old.h5'
+#            f = h5py.File(xdmfh5file, 'r')
+            nx = 40
+            ny = 70
+            nz = 233
+#            data = f['Data2']
+#
+            for index in range(len(data)):
+                subscripts = np.unravel_index(index, (nz, ny, nx))
+                i = subscripts[0]
+                j = subscripts[1]
+                k = subscripts[2]
+                numIDlocal = int(data[index])
+                for material in G.materials:
+                    try:
+                        material.ID.split('_')[1]
+                        if int(material.ID.split('_')[1]) == numIDlocal:
+                            numID = numIDx = numIDy = numIDz = material.numID
+                            build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, False, G.solid, G.rigidE, G.rigidH, G.ID)
+                        else:
+                            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' cannot find material ID {}'.format(numIDlocal))
+                    except:
+                        pass
+
+        
+        elif tmp[0] == '#edge:':
             if len(tmp) != 8:
                 raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters')
             
@@ -1181,7 +1240,7 @@ def process_geometrycmds(geometry, G):
                 if volume.nbins == 1:
                     volume.fractalvolume = np.ones((volume.nx + 1, volume.ny + 1, volume.nz + 1), dtype=floattype)
                     materialnumID = next(x.numID for x in G.materials if x.ID == volume.operatingonID)
-                    volume.fractalvolume = volume.fractalvolume * materialnumID
+                    volume.fractalvolume *= materialnumID
                 else:
                     volume.generate_fractal_volume(G)
                     volume.fractalvolume += mixingmodel.startmaterialnum