From 2f7c2edbfda0995448fbb187f49cc401358c1ce4 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 30 Sep 2016 14:23:36 +0100 Subject: [PATCH] Added averaging (dielectric smoothing) option when reading geometry objects from file. Added checking for presence of rigid arrays when reading geometry objects from file. Change material 'average' parameter name to 'averrable' for clarity. Added averaging (dielectric smoothing) option for #fractal_box, default is False (off). --- gprMax/input_cmds_geometry.py | 80 ++++++++++++++++++++++++++++------- 1 file changed, 65 insertions(+), 15 deletions(-) diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index 282efeae..7dd61b4b 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -50,9 +50,22 @@ def process_geometrycmds(geometry, G): for object in tqdm(geometry, desc='Processing geometry related cmds', unit='cmds', ncols=get_terminal_width() - 1, file=sys.stdout, disable=tqdmdisable): tmp = object.split() - if tmp[0] == '#geometry_objects_file:': - if len(tmp) != 6: - raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly five parameters') + if tmp[0] == '#geometry_objects_read:': + averaging = G.averagevolumeobjects + + if len(tmp) < 6: + raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least five parameters') + + elif len(tmp) == 7: + if tmp[6].lower() == 'y': + averaging = True + elif tmp[6].lower() == 'n': + averaging = False + else: + raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n') + + elif len(tmp) > 7: + raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') xs = round_value(float(tmp[1]) / G.dx) ys = round_value(float(tmp[2]) / G.dy) @@ -77,6 +90,11 @@ def process_geometrycmds(geometry, G): # Process parameters for commands that can occur multiple times in the model process_multicmds(multicmds, G) + # Update material type + for material in G.materials: + if material.numID >= numexistmaterials: + material.type += ', imported' + # See if geometry object file exists at specified path and if not try input file directory if not os.path.isfile(geofile): geofile = os.path.abspath(os.path.join(G.inputdirectory, geofile)) @@ -93,10 +111,23 @@ def process_geometrycmds(geometry, G): if data.dtype != 'int16': data = data.astype('int16') - build_voxels_from_array(xs, ys, zs, numexistmaterials, data, G.solid, G.rigidE, G.rigidH, G.ID) + # Look rigid arrays if present (these should be present if the original geometry objects were written from gprMax) + try: + rigidE = f['/rigidE'][:] + rigidH = f['/rigidH'][:] + G.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE + G.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH + except KeyError: + averaging = False + + build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, data, G.solid, G.rigidE, G.rigidH, G.ID) if G.messages: - tqdm.write('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}, dielectric smoothing is off for these objects.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile)) + if averaging: + dielectricsmoothing = 'on' + else: + dielectricsmoothing = 'off' + tqdm.write('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}, dielectric smoothing is {} for these objects.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile, dielectricsmoothing)) elif tmp[0] == '#edge:': if len(tmp) != 8: @@ -334,7 +365,7 @@ def process_geometrycmds(geometry, G): if thickness > 0: # Isotropic case if len(materials) == 1: - averaging = materials[0].average and averagetriangularprism + averaging = materials[0].averagable and averagetriangularprism numID = numIDx = numIDy = numIDz = materials[0].numID # Uniaxial anisotropic case @@ -442,7 +473,7 @@ def process_geometrycmds(geometry, G): # Isotropic case if len(materials) == 1: - averaging = materials[0].average and averagebox + averaging = materials[0].averagable and averagebox numID = numIDx = numIDy = numIDz = materials[0].numID # Uniaxial anisotropic case @@ -523,7 +554,7 @@ def process_geometrycmds(geometry, G): # Isotropic case if len(materials) == 1: - averaging = materials[0].average and averagecylinder + averaging = materials[0].averagable and averagecylinder numID = numIDx = numIDy = numIDz = materials[0].numID # Uniaxial anisotropic case @@ -619,7 +650,7 @@ def process_geometrycmds(geometry, G): if thickness > 0: # Isotropic case if len(materials) == 1: - averaging = materials[0].average and averagecylindricalsector + averaging = materials[0].averagable and averagecylindricalsector numID = numIDx = numIDy = numIDz = materials[0].numID # Uniaxial anisotropic case @@ -729,7 +760,7 @@ def process_geometrycmds(geometry, G): # Isotropic case if len(materials) == 1: - averaging = materials[0].average and averagesphere + averaging = materials[0].averagable and averagesphere numID = numIDx = numIDy = numIDz = materials[0].numID # Uniaxial anisotropic case @@ -765,12 +796,23 @@ def process_geometrycmds(geometry, G): tqdm.write('Sphere with centre {:g}m, {:g}m, {:g}m, radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(xc * G.dx, yc * G.dy, zc * G.dz, r, ', '.join(materialsrequested), dielectricsmoothing)) elif tmp[0] == '#fractal_box:': + # Default is no dielectric smoothing for a fractal box + averagefractalbox = False + if len(tmp) < 14: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least thirteen parameters') elif len(tmp) == 14: seed = None elif len(tmp) == 15: seed = int(tmp[14]) + elif len(tmp) == 16: + seed = int(tmp[14]) + if tmp[15].lower() == 'y': + averagefractalbox = True + elif tmp[15].lower() == 'n': + averagefractalbox = False + else: + raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires averaging to be either y or n') else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') @@ -825,9 +867,17 @@ def process_geometrycmds(geometry, G): volume.nbins = nbins volume.seed = seed volume.weighting = (float(tmp[8]), float(tmp[9]), float(tmp[10])) + try: + volume.averaging = averagefractalbox + except: + pass if G.messages: - tqdm.write('Fractal box {} from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with {}, fractal dimension {:g}, fractal weightings {:g}, {:g}, {:g}, fractal seeding {} and with {} values created.'.format(volume.ID, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, volume.operatingonID, volume.dimension, volume.weighting[0], volume.weighting[1], volume.weighting[2], volume.seed, volume.nbins)) + if volume.averaging: + dielectricsmoothing = 'on' + else: + dielectricsmoothing = 'off' + tqdm.write('Fractal box {} from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with {}, fractal dimension {:g}, fractal weightings {:g}, {:g}, {:g}, fractal seeding {}, with {} material(s) created, dielectric smoothing is {}.'.format(volume.ID, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, volume.operatingonID, volume.dimension, volume.weighting[0], volume.weighting[1], volume.weighting[2], volume.seed, volume.nbins, dielectricsmoothing)) G.fractalvolumes.append(volume) @@ -1038,7 +1088,7 @@ def process_geometrycmds(geometry, G): # Check to see if water has been already defined as a material if not any(x.ID == 'water' for x in G.materials): m = Material(len(G.materials), 'water') - m.average = False + m.averagable = False m.type = 'builtin, debye' m.er = Material.watereri m.deltaer.append(Material.waterdeltaer) @@ -1186,7 +1236,7 @@ def process_geometrycmds(geometry, G): # Check to see if grass has been already defined as a material if not any(x.ID == 'grass' for x in G.materials): m = Material(len(G.materials), 'grass') - m.average = False + m.averagable = False m.type = 'builtin, debye' m.er = Material.grasseri m.deltaer.append(Material.grassdeltaer) @@ -1455,7 +1505,7 @@ def process_geometrycmds(geometry, G): grassnumID = next((x.numID for x in G.materials if x.ID == 'grass'), 0) data = volume.fractalvolume.astype('int16', order='C') mask = volume.mask.copy(order='C') - build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs, waternumID, grassnumID, mask, data, G.solid, G.rigidE, G.rigidH, G.ID) + build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs, waternumID, grassnumID, volume.averaging, mask, data, G.solid, G.rigidE, G.rigidH, G.ID) else: if volume.nbins == 1: @@ -1465,4 +1515,4 @@ def process_geometrycmds(geometry, G): volume.fractalvolume += mixingmodel.startmaterialnum data = volume.fractalvolume.astype('int16', order='C') - build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0, data, G.solid, G.rigidE, G.rigidH, G.ID) + build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0, volume.averaging, data, G.solid, G.rigidE, G.rigidH, G.ID)