Initial investigation of Xdmf-based file reading.

这个提交包含在:
craig-warren
2016-03-13 12:37:17 +00:00
父节点 c548e715f8
当前提交 0c90153a21

查看文件

@@ -16,7 +16,10 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with gprMax. If not, see <http://www.gnu.org/licenses/>. # along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import os
import h5py
import numpy as np import numpy as np
import xml.etree.ElementTree as ET
from gprMax.constants import floattype from gprMax.constants import floattype
from gprMax.exceptions import CmdInputError 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() 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: if len(tmp) != 8:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters') raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters')
@@ -1181,7 +1240,7 @@ def process_geometrycmds(geometry, G):
if volume.nbins == 1: if volume.nbins == 1:
volume.fractalvolume = np.ones((volume.nx + 1, volume.ny + 1, volume.nz + 1), dtype=floattype) 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) materialnumID = next(x.numID for x in G.materials if x.ID == volume.operatingonID)
volume.fractalvolume = volume.fractalvolume * materialnumID volume.fractalvolume *= materialnumID
else: else:
volume.generate_fractal_volume(G) volume.generate_fractal_volume(G)
volume.fractalvolume += mixingmodel.startmaterialnum volume.fractalvolume += mixingmodel.startmaterialnum