From b92ce556f6f747d59d375f7149828f74efa4355d Mon Sep 17 00:00:00 2001 From: jasminium Date: Tue, 10 May 2016 16:12:40 +0100 Subject: [PATCH] hdf5 xdmf support --- gprMax/geometry_views.py | 47 ++--- gprMax/gprMax.py | 73 ++++---- gprMax/grid.py | 67 +++++--- gprMax/xdmf.py | 361 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 469 insertions(+), 79 deletions(-) create mode 100644 gprMax/xdmf.py diff --git a/gprMax/geometry_views.py b/gprMax/geometry_views.py index eade70b9..c3ad4094 100644 --- a/gprMax/geometry_views.py +++ b/gprMax/geometry_views.py @@ -19,13 +19,14 @@ import sys import numpy as np from struct import pack +from gprMax.xdmf import write_output_file from gprMax.utilities import round_value class GeometryView: """Views of the geometry of the model.""" - + if sys.byteorder == 'little': byteorder = 'LittleEndian' else: @@ -39,7 +40,7 @@ class GeometryView: filename (str): Filename to save to. type (str): Either 'n' for a per cell geometry view, or 'f' for a per cell edge geometry view. """ - + self.xs = xs self.ys = ys self.zs = zs @@ -55,26 +56,30 @@ class GeometryView: self.filename = filename self.type = type + def write_xdmf(self, modelrun, numbermodelruns, G): + filename = self.filename[:-4] + write_output_file(filename, G) + 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. - + N.B. No Python 3 support for VTK at time of writing (03/2015) - + Args: modelrun (int): Current model run number. numbermodelruns (int): Total number of model runs. G (class): Grid class instance - holds essential parameters describing the model. """ - + # Construct filename from user-supplied name and model run number if numbermodelruns == 1: self.filename = G.inputdirectory + self.filename else: self.filename = G.inputdirectory + self.filename + str(modelrun) - + if self.type == 'n': self.filename += '.vti' - + # Calculate number of cells according to requested sampling for geometry view self.vtk_xscells = round_value(self.xs / self.dx) self.vtk_xfcells = round_value(self.xf / self.dx) @@ -85,16 +90,16 @@ class GeometryView: self.vtk_nxcells = self.vtk_xfcells - self.vtk_xscells self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells - + # Create an array and add numeric IDs for PML, sources and receivers self.srcs_rxs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8) for pml in G.pmls: self.srcs_rxs_pml[pml.xs:pml.xf, pml.ys:pml.yf, pml.zs:pml.zf] = 1 for index, srcrx in enumerate(G.rxs + G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines): self.srcs_rxs_pml[srcrx.xcoord, srcrx.ycoord, srcrx.zcoord] = index + 2 - + vtk_srcs_rxs_pml_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize) - + with open(self.filename, 'wb') as f: f.write('\n'.encode('utf-8')) f.write('\n'.format(GeometryView.byteorder).encode('utf-8')) @@ -105,7 +110,7 @@ class GeometryView: f.write('\n'.format(vtk_srcs_rxs_pml_offset).encode('utf-8')) f.write('\n'.encode('utf-8')) f.write('\n\n\n_'.encode('utf-8')) - + # Write material IDs datasize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) # Write number of bytes of appended data as UInt32 @@ -114,7 +119,7 @@ class GeometryView: for j in range(self.ys, self.yf, self.dy): for i in range(self.xs, self.xf, self.dx): f.write(pack('I', G.solid[i, j, k])) - + # Write source/receiver IDs datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) f.write(pack('I', datasize)) @@ -124,34 +129,34 @@ class GeometryView: f.write(pack('b', self.srcs_rxs_pml[i, j, k])) f.write('\n\n'.encode('utf-8')) - + self.write_gprmax_info(f, G) elif self.type == 'f': self.filename += '.vtp' - + vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1) vtk_numpoint_components = 3 vtk_numlines = 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 vtk_numline_components = 2; - + vtk_connectivity_offset = round_value((vtk_numpoints * vtk_numpoint_components * np.dtype(np.float32).itemsize) + np.dtype(np.uint32).itemsize) vtk_offsets_offset = round_value(vtk_connectivity_offset + (vtk_numlines * vtk_numline_components * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize) vtk_materials_offset = round_value(vtk_offsets_offset + (vtk_numlines * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize) - + with open(self.filename, 'wb') as f: f.write('\n'.encode('utf-8')) f.write('\n'.format(GeometryView.byteorder).encode('utf-8')) f.write('\n\n'.format(vtk_numpoints, vtk_numlines).encode('utf-8')) - + f.write('\n\n\n'.encode('utf-8')) f.write('\n\n'.format(vtk_connectivity_offset).encode('utf-8')) f.write('\n\n'.format(vtk_offsets_offset).encode('utf-8')) - + f.write('\n'.encode('utf-8')) f.write('\n'.format(vtk_materials_offset).encode('utf-8')) f.write('\n'.encode('utf-8')) - + f.write('\n\n\n_'.encode('utf-8')) # Write points @@ -209,7 +214,7 @@ class GeometryView: for j in range(self.ys, self.yf + 1): for k in range(self.zs, self.zf + 1): f.write(pack('I', G.ID[0, i, j, k])) - + for i in range(self.xs, self.xf + 1): for j in range(self.ys, self.yf): for k in range(self.zs, self.zf + 1): @@ -227,7 +232,7 @@ class GeometryView: 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. diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index bf126edf..1cb32c53 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -42,7 +42,7 @@ from gprMax.yee_cell_build import build_electric_components, build_magnetic_comp def main(): """This is the main function for gprMax.""" - + # Print gprMax logo, version, and licencing/copyright information logo(gprMax.__version__ + ' (Bowmore)') @@ -59,7 +59,7 @@ def main(): numbermodelruns = args.n inputdirectory = os.path.dirname(os.path.abspath(args.inputfile)) + os.sep inputfile = inputdirectory + os.path.basename(args.inputfile) - + # Create a separate namespace that users can access in any Python code blocks in the input file usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'inputdirectory': inputdirectory} @@ -92,7 +92,7 @@ def main(): def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None): """Run standard simulation - models are run one after another and each model is parallelised with OpenMP - + Args: args (dict): Namespace with command line arguments numbermodelruns (int): Total number of model runs. @@ -100,7 +100,7 @@ def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. optparams (dict): Optional argument. For Taguchi optimisation it provides the parameters to optimise and their values. """ - + tsimstart = perf_counter() for modelrun in range(1, numbermodelruns + 1): if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace @@ -117,20 +117,20 @@ def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) def run_benchmark_sim(args, inputfile, usernamespace): """Run standard simulation in benchmarking mode - models are run one after another and each model is parallelised with OpenMP - + Args: args (dict): Namespace with command line arguments inputfile (str): Name of the input file to open. usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. """ - + # Number of threads to test - start from max physical CPU cores and divide in half until 1 thread = psutil.cpu_count(logical=False) threads = [thread] while not thread%2: thread /= 2 threads.append(int(thread)) - + benchtimes = np.zeros(len(threads)) numbermodelruns = len(threads) @@ -150,7 +150,7 @@ def run_benchmark_sim(args, inputfile, usernamespace): def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None): """Run mixed mode MPI/OpenMP simulation - MPI task farm for models with each model parallelised with OpenMP - + Args: args (dict): Namespace with command line arguments numbermodelruns (int): Total number of model runs. @@ -180,7 +180,7 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) data = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=status) source = status.Get_source() tag = status.Get_tag() - + if tag == tags.READY.value: # Worker is ready, so send it a task if modelrun < numbermodelruns + 1: comm.send(modelrun, dest=source, tag=tags.START.value) @@ -188,10 +188,10 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) modelrun += 1 else: comm.send(None, dest=source, tag=tags.EXIT.value) - + elif tag == tags.DONE.value: print('Worker {}: completed.'.format(source)) - + elif tag == tags.EXIT.value: print('Worker {}: exited.'.format(source)) closedworkers += 1 @@ -202,7 +202,7 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) comm.send(None, dest=0, tag=tags.READY.value) modelrun = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) # Receive a model number to run from the master tag = status.Get_tag() - + # Run a model if tag == tags.START.value: if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace @@ -212,10 +212,10 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) modelusernamespace.update({'optparams': tmp}) else: modelusernamespace = usernamespace - + run_model(args, modelrun, numbermodelruns, inputfile, modelusernamespace) comm.send(None, dest=0, tag=tags.DONE.value) - + elif tag == tags.EXIT.value: break @@ -224,34 +224,34 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None) def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): """Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop. - + Args: args (dict): Namespace with command line arguments modelrun (int): Current model run number. numbermodelruns (int): Total number of model runs. inputfile (str): Name of the input file to open. usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file. - + Returns: tsolve (int): Length of time (seconds) of main FDTD calculations """ - + # Monitor memory usage p = psutil.Process() - + print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile)) - + # Add the current model run to namespace that can be accessed by user in any Python code blocks in input file usernamespace['current_model_run'] = modelrun print('Constants/variables available for Python scripting: {}\n'.format(usernamespace)) - + # Process any user input Python commands processedlines = process_python_include_code(inputfile, usernamespace) - + # Write a file containing the input commands after Python blocks have been processed if args.write_processed: write_processed_file(inputfile, modelrun, numbermodelruns, processedlines) - + # Check validity of command names & that essential commands are present singlecmds, multicmds, geometry = check_cmd_names(processedlines) @@ -284,7 +284,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): # Build the PML and calculate initial coefficients build_pmls(G) - + # Build the model, i.e. set the material properties (ID) for every edge of every Yee cell tbuildstart = perf_counter() build_electric_components(G.solid, G.rigidE, G.ID, G) @@ -295,7 +295,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): # Process any voltage sources (that have resistance) to create a new material at the source location for voltagesource in G.voltagesources: voltagesource.create_material(G) - + # Initialise arrays of update coefficients to pass to update functions G.initialise_std_updatecoeff_arrays() @@ -309,22 +309,22 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): print('ID\tName\t\tProperties') print('{}'.format('-'*50)) for material in G.materials: - + # Calculate update coefficients for material material.calculate_update_coeffsE(G) material.calculate_update_coeffsH(G) - + # Store all update coefficients together G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm - + # Store coefficients for any dispersive materials if Material.maxpoles != 0: z = 0 for pole in range(Material.maxpoles): G.updatecoeffsdispersive[material.numID, z:z+3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole] z += 3 - + if G.messages: if material.deltaer and material.tau: tmp = 'delta_epsr={}, tau={} secs; '.format(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer), ', '.join('{:g}'.format(tau) for tau in material.tau)) @@ -340,7 +340,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): resolution = dispersion_check(G) if resolution != 0 and max((G.dx, G.dy, G.dz)) > resolution: print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution)) - + # Write files for any geometry views if not G.geometryviews and args.geometry_only: raise GeneralError('No geometry views found.') @@ -348,16 +348,17 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): tgeostart = perf_counter() for geometryview in G.geometryviews: geometryview.write_vtk(modelrun, numbermodelruns, G) + geometryview.write_xdmf(modelrun, numbermodelruns, G) tgeoend = perf_counter() print('\nGeometry file(s) written in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tgeoend - tgeostart)))) # Run simulation if not doing only geometry if not args.geometry_only: - + # Prepare any snapshot files for snapshot in G.snapshots: snapshot.prepare_vtk_imagedata(modelrun, numbermodelruns, G) - + # Adjust position of sources and receivers if required if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0: for source in itertools.chain(G.hertziandipoles, G.magneticdipoles, G.voltagesources, G.transmissionlines): @@ -396,10 +397,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): for timestep in range(G.iterations): if timestep == 0: tstepstart = perf_counter() - + # Write field outputs to file write_hdf5(f, timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G) - + # Write any snapshots to file for snapshot in G.snapshots: if snapshot.time == timestep + 1: @@ -432,7 +433,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): # Increment absolute time value abstime += 0.5 * G.dt - + # Update magnetic field components update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz) @@ -447,7 +448,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): # Increment absolute time value abstime += 0.5 * G.dt - + # Calculate time for two iterations, used to estimate overall runtime if timestep == 1: tstepend = perf_counter() @@ -457,7 +458,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): sys.stdout.flush() elif timestep > 1: update_progress((timestep + 1) / G.iterations) - + # Close output file f.close() diff --git a/gprMax/grid.py b/gprMax/grid.py index 365b2d08..6e970695 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -22,10 +22,33 @@ import matplotlib.pyplot as plt from gprMax.constants import c, floattype, complextype from gprMax.materials import Material +class Grid(): -class FDTDGrid: + def __init__(self, grid): + self.nx = grid.shape[0] + self.ny = grid.shape[1] + self.nz = grid.shape[2] + self.grid = grid + + def n_edges(self): + l = self.nx + m = self.ny + n = self.nz + e = (l * m * (n - 1)) + (m * n * (l - 1)) + (l * n * (m - 1)) + return e + + def n_nodes(self): + return self.nx * self.ny * self.nz + + def n_cells(self): + return (self.nx - 1) * (self.ny - 1) * (self.nz - 1) + + def get(self, i, j, k): + return self.grid[i, j, k] + +class FDTDGrid(Grid): """Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters.""" - + def __init__(self): self.inputdirectory = '' self.title = '' @@ -61,7 +84,7 @@ class FDTDGrid: self.rxstepz = 0 self.rxs = [] self.snapshots = [] - + def initialise_std_arrays(self): """Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid); an array for cell edge IDs (ID); and arrays for the electric and magnetic field components. Solid and ID arrays are initialised to free_space (one); rigid arrays @@ -78,7 +101,7 @@ class FDTDGrid: self.Hx = np.zeros((self.nx + 1, self.ny, self.nz), dtype=floattype) self.Hy = np.zeros((self.nx, self.ny + 1, self.nz), dtype=floattype) self.Hz = np.zeros((self.nx, self.ny, self.nz + 1), dtype=floattype) - + def initialise_std_updatecoeff_arrays(self): """Initialise arrays for storing update coefficients.""" self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=floattype) @@ -94,25 +117,25 @@ class FDTDGrid: def dispersion_check(G): """Check for potential numerical dispersion. Is the smallest wavelength present in the simulation discretised by at least a factor of 10 - + Args: G (class): Grid class instance - holds essential parameters describing the model. - + Returns: resolution (float): Potential numerical dispersion """ - + # Minimum number of spatial steps to resolve smallest wavelength resolvedsteps = 10 - + # Find maximum frequency maxfreqs = [] for waveform in G.waveforms: - + # User-defined waveform if waveform.uservalues is not None: waveformvalues = waveform.uservalues - + # Built-in waveform else: time = np.linspace(0, 1, G.iterations) @@ -130,24 +153,24 @@ def dispersion_check(G): # Shift powers so that frequency with maximum power is at zero decibels power -= np.amax(power) - + # Set maximum frequency to -60dB from maximum power freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1 maxfreqs.append(freqs[freq]) if maxfreqs: maxfreq = max(maxfreqs) - + # Find minimum wavelength ers = [material.er for material in G.materials] maxer = max(ers) # Minimum velocity minvelocity = c / np.sqrt(maxer) - + # Minimum wavelength minwavelength = minvelocity / maxfreq - + # Resolution of minimum wavelength resolution = minwavelength / resolvedsteps @@ -159,10 +182,10 @@ def dispersion_check(G): def get_other_directions(direction): """Return the two other directions from x, y, z given a single direction - + Args: direction (str): Component x, y or z - + Returns: (tuple): Two directions from x, y, z """ @@ -174,7 +197,7 @@ def get_other_directions(direction): def Ix(x, y, z, Hy, Hz, G): """Calculates the x-component of current at a grid position. - + Args: x, y, z (float): Coordinates of position in grid. Hy, Hz (memory view): numpy array of magnetic field values. @@ -184,14 +207,14 @@ def Ix(x, y, z, Hy, Hz, G): if y == 0 or z == 0: Ix = 0 return Ix - + else: Ix = G.dy * (Hy[x, y, z - 1] - Hy[x, y, z]) + G.dz * (Hz[x, y, z] - Hz[x, y - 1, z]) return Ix def Iy(x, y, z, Hx, Hz, G): """Calculates the y-component of current at a grid position. - + Args: x, y, z (float): Coordinates of position in grid. Hx, Hz (memory view): numpy array of magnetic field values. @@ -201,14 +224,14 @@ def Iy(x, y, z, Hx, Hz, G): if x == 0 or z == 0: Iy = 0 return Iy - + else: Iy = G.dx * (Hx[x, y, z] - Hx[x, y, z - 1]) + G.dz * (Hz[x - 1, y, z] - Hz[x, y, z]) return Iy def Iz(x, y, z, Hx, Hy, G): """Calculates the z-component of current at a grid position. - + Args: x, y, z (float): Coordinates of position in grid. Hx, Hy (memory view): numpy array of magnetic field values. @@ -218,7 +241,7 @@ def Iz(x, y, z, Hx, Hy, G): if x == 0 or y == 0: Iz = 0 return Iz - + else: Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z]) return Iz diff --git a/gprMax/xdmf.py b/gprMax/xdmf.py new file mode 100644 index 00000000..326ff6ab --- /dev/null +++ b/gprMax/xdmf.py @@ -0,0 +1,361 @@ +import numpy as np +import h5py +from lxml import etree +from gprMax.grid import Grid +import copy + +class Edges: + + def __init__(self, grid): + + """ + Class to define some connectivity of for an n x l x m + grid + """ + self.total_edges = grid.n_edges() + self.edges = np.zeros((self.total_edges, 2), np.float32) + self.edge_count = 0 + self.grid = grid + """ + Adds the the edge specified by in_node and the i,j,k position of the outnode + """ + def add_edge(self, in_label, i, j, k): + + out_label = self.grid.get(i, j, k) + self.edges[self.edge_count] = np.array([in_label, out_label]) + self.edge_count += 1 + +class Coordinates: + + def __init__(self, grid): + self.total_coordinates = grid.nx * grid.ny * grid.nz + self.coordinate_count = 0 + self.coordinates = np.zeros((self.total_coordinates, 3), np.float32) + + def add_coordinate(self, x, y, z): + self.coordinates[self.coordinate_count] = np.array([x, y, z]) + 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): + 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) + self.solid_labels[self.count] = solid_labels + self.count += 1 + + +class Materials: + + def __init__(self, fdtd_grid): + self.fdtd_grid = fdtd_grid + self.n_edges = fdtd_grid.n_edges() + self.materials = np.zeros((self.n_edges), np.float32) + self.material_count = 0 + + # direction x->0 y->1 z->2 + def add_material(self, i, j, k, direction): + + material = self.fdtd_grid.ID[direction, i, j, k] + self.materials[self.material_count] = material + + self.material_count += 1 + +def process_grid(fdtd_grid): + + # Dimensions of the problem domain. + nx = fdtd_grid.nx + ny = fdtd_grid.ny + nz = fdtd_grid.nz + + # 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) + + # Material for each solid + solids = Solids(fdtd_grid) + + # Connectivity for hexhahedron grid + solid_labels = SolidLabels(label_grid) + + i_max = nx - 1 + j_max = ny - 1 + k_max = nz - 1 + + for i, ix in enumerate(labels): + for j, jx in enumerate(ix): + for k, kx in enumerate(jx): + + 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) + + # 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 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 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) + + #x = np.arange(fdtd_grid.nx) + #y = np.arange(fdtd_grid.ny) + #z = np.arange(fdtd_grid.nz) + + return { + 'coordinates': coordinates, + 'solids': solids, + 'solid_labels': solid_labels, + 'edges': edges, + 'edge_materials': edge_materials, + #'x': x, + #'y': y, + #'z': z + } + +def write_output_file(filename, grid): + + data = process_grid(grid) + data['filename'] = filename + data['xml_doc'] = create_xdmf_markup(data) + + write_H5file(data) + write_xml_doc(data) + +def write_xml_doc(options): + #write xml to file + with open(options['filename'] + '.xdmf', 'wb') as xdmf_f: + xdmf_f.write(options['xml_doc']) + +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) + #coords.create_dataset('x', data=options['x']) + #coords.create_dataset('y', data=options['y']) + #coords.create_dataset('z', data=options['z']) + data = f.create_group("data") + data.create_dataset('materials', data=options['edge_materials'].materials) + data.create_dataset('solids', data=options['solids'].solids) + +def create_xdmf_markup(options): + + # Write the XDMF markup for edge style grid + xdmf_el = etree.Element("Xdmf", Version="2.0") + + 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) + + 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) + + """ + # VOXEL style markup + v_grid_el = etree.Element("Grid", Name="Voxel", GridType="Uniform") + domain_el.append(v_grid_el) + + noe = "{} {} {}".format(options['x'].size, options['y'].size, options['y'].size) + v_topology_el = etree.Element("Topology", TopologyType="3DRectMesh", NumberOfElements=noe) + v_grid_el.append(v_topology_el) + + v_geometry = etree.Element("Geometry", GeometryType="VXVYVZ") + v_grid_el.append(v_geometry) + + d1 = etree.Element("DataItem", Dimensions=str(options['x'].size), NumberType="Float", Precision="4", Format="HDF") + d1.text = "{}:/mesh/x".format(options['filename'] + '.h5') + v_geometry.append(d1) + + d2 = etree.Element("DataItem", Dimensions=str(options['y'].size), NumberType="Float", Precision="4", Format="HDF") + d2.text = "{}:/mesh/y".format(options['filename'] + '.h5') + v_geometry.append(d2) + + d3 = etree.Element("DataItem", Dimensions=str(options['z'].size), NumberType="Float", Precision="4", Format="HDF") + d3.text = "{}:/mesh/z".format(options['filename'] + '.h5') + v_geometry.append(d3) + + v_attr = etree.Element("Attribute", Name="material-blocks", Center="Cell") + v_grid_el.append(v_attr) + + d4 = etree.Element("DataItem", Format="HDF", NumberType="Float", Precision="4", Dimensions=str(options['solids'].solids.size)) + d4.text = "{}:/data/solids".format(options['filename'] + '.h5') + v_attr.append(d4) + + """ + + v_grid_el = etree.Element("Grid", Name="Voxel", GridType="Uniform") + domain_el.append(v_grid_el) + + n_solids = str(options['solids'].solids.size) + v_topology_el = etree.Element("Topology", TopologyType="Hexahedron", NumberOfElements=str(options['solids'].solids.size)) + v_grid_el.append(v_topology_el) + + solid_label_d = "{} {}".format(n_solids, 8) + solid_labels_el = etree.Element("DataItem", Dimensions=solid_label_d, Format="HDF") + solid_labels_el.text = "{}:/mesh/solid_connectivity".format(options['filename'] + '.h5') + v_topology_el.append(solid_labels_el) + + # Same geometry as edges + v_grid_el.append(copy.deepcopy(geometry_el)) + + v_attr = etree.Element("Attribute", Name="Voxel_Materials", Center="Cell") + v_grid_el.append(v_attr) + + d4 = etree.Element("DataItem", Format="HDF", NumberType="Float", Precision="4", Dimensions=str(options['solids'].solids.size)) + d4.text = "{}:/data/solids".format(options['filename'] + '.h5') + v_attr.append(d4) + + # Define a doctype - useful for parsers + doc_type = '' + + # Serialize elements + xml_doc = etree.tostring(xdmf_el, xml_declaration=True, + encoding="utf-8", doctype=doc_type, pretty_print=True) + + return xml_doc +