From e8ebd1498f016eebf00b65ad4b0663b153395aa6 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 21 Jul 2016 16:17:39 +0100 Subject: [PATCH] Whitespace cleanups. --- gprMax/exceptions.py | 2 + gprMax/fractals.py | 20 ++-- gprMax/geometry_views.py | 3 +- gprMax/grid.py | 5 +- gprMax/input_cmd_funcs.py | 3 +- gprMax/input_cmds_geometry.py | 163 +++++++++++++++------------------ gprMax/input_cmds_multiuse.py | 18 +--- gprMax/input_cmds_singleuse.py | 33 +++---- gprMax/materials.py | 18 ++-- gprMax/optimisation_taguchi.py | 9 +- gprMax/pml.py | 26 +++--- gprMax/sources.py | 6 +- gprMax/utilities.py | 8 +- gprMax/writer_hdf5.py | 2 - gprMax/xdmf.py | 6 +- 15 files changed, 147 insertions(+), 175 deletions(-) diff --git a/gprMax/exceptions.py b/gprMax/exceptions.py index e6165032..e4d84d58 100644 --- a/gprMax/exceptions.py +++ b/gprMax/exceptions.py @@ -16,10 +16,12 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . + class GeneralError(ValueError): """Handles general errors. Subclasses the ValueError class.""" pass + class CmdInputError(ValueError): """Handles errors in user specified commands. Subclasses the ValueError class.""" pass diff --git a/gprMax/fractals.py b/gprMax/fractals.py index 9f36ce5e..bf9a5bcb 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -76,7 +76,7 @@ class FractalSurface(object): self.fractalsurface = np.zeros(surfacedims, dtype=complextype) # Positional vector at centre of array, scaled by weighting - v1 = np.array([self.weighting[0]*(surfacedims[0])/2, self.weighting[1]*(surfacedims[1])/2]) + v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1] * (surfacedims[1]) / 2]) # 2D array of random numbers to be convolved with the fractal function R = np.random.RandomState(self.seed) @@ -88,13 +88,13 @@ class FractalSurface(object): for i in range(surfacedims[0]): for j in range(surfacedims[1]): # Positional vector for current position - v2 = np.array([self.weighting[0]*i, self.weighting[1]*j]) + v2 = np.array([self.weighting[0] * i, self.weighting[1] * j]) rr = np.linalg.norm(v2 - v1) try: - self.fractalsurface[i, j] = A[i, j] * 1/(rr**self.b) + self.fractalsurface[i, j] = A[i, j] * 1 / (rr**self.b) except FloatingPointError: rr = 0.9 - self.fractalsurface[i, j] = A[i, j] * 1/(rr**self.b) + self.fractalsurface[i, j] = A[i, j] * 1 / (rr**self.b) # Shift the zero frequency component to the centre of the spectrum self.fractalsurface = np.fft.ifftshift(self.fractalsurface) @@ -104,7 +104,7 @@ class FractalSurface(object): fractalmin = np.amin(self.fractalsurface) fractalmax = np.amax(self.fractalsurface) fractalrange = fractalmax - fractalmin - self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0])/fractalrange) + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0])/fractalrange) * fractalmin + self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin class FractalVolume(object): @@ -118,7 +118,7 @@ class FractalVolume(object): """ self.ID = None - self.operatingonID= None + self.operatingonID = None self.xs = xs self.xf = xf self.ys = ys @@ -146,7 +146,7 @@ class FractalVolume(object): self.fractalvolume = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=complextype) # Positional vector at centre of array, scaled by weighting - v1 = np.array([self.weighting[0]*(self.nx + 1)/2, self.weighting[1]*(self.ny + 1)/2, self.weighting[2]*(self.nz + 1)/2]) + v1 = np.array([self.weighting[0] * (self.nx + 1) / 2, self.weighting[1] * (self.ny + 1) / 2, self.weighting[2] * (self.nz + 1) / 2]) # 3D array of random numbers to be convolved with the fractal function R = np.random.RandomState(self.seed) @@ -159,13 +159,13 @@ class FractalVolume(object): for j in range(self.ny + 1): for k in range(self.nz + 1): # Positional vector for current position - v2 = np.array([self.weighting[0]*i, self.weighting[1]*j, self.weighting[2]*k]) + v2 = np.array([self.weighting[0] * i, self.weighting[1] * j, self.weighting[2] * k]) rr = np.linalg.norm(v2 - v1) try: - self.fractalvolume[i, j, k] = A[i, j, k] * 1/(rr**self.b) + self.fractalvolume[i, j, k] = A[i, j, k] * 1 / (rr**self.b) except FloatingPointError: rr = 0.9 - self.fractalvolume[i, j, k] = A[i, j, k] * 1/(rr**self.b) + self.fractalvolume[i, j, k] = A[i, j, k] * 1 / (rr**self.b) # Shift the zero frequency component to the centre of the spectrum self.fractalvolume = np.fft.ifftshift(self.fractalvolume) diff --git a/gprMax/geometry_views.py b/gprMax/geometry_views.py index 3b137ff2..4c06e4d3 100644 --- a/gprMax/geometry_views.py +++ b/gprMax/geometry_views.py @@ -153,7 +153,7 @@ class GeometryView(object): 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_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) @@ -244,7 +244,6 @@ class GeometryView(object): self.write_gprmax_info(f, G, materialsonly=True) - def write_gprmax_info(self, f, G, materialsonly=False): """Writes gprMax specific information relating material, source, and receiver names to numeric identifiers. diff --git a/gprMax/grid.py b/gprMax/grid.py index 670ae354..a9172c38 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -63,6 +63,7 @@ class Grid(object): co = round_value(float(val) / getattr(self, 'd' + coord)) return co + class FDTDGrid(Grid): """Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters.""" @@ -155,7 +156,7 @@ def dispersion_check(G): if waveform.type == 'sine' or waveform.type == 'contsine': maxfreqs.append(4 * waveform.freq) - elif waveform.type =='impulse': + elif waveform.type == 'impulse': pass else: @@ -239,6 +240,7 @@ def Ix(x, y, z, Hy, Hz, G): 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. @@ -256,6 +258,7 @@ def Iy(x, y, z, Hx, Hz, G): 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. diff --git a/gprMax/input_cmd_funcs.py b/gprMax/input_cmd_funcs.py index 019e67b7..c4f51f8b 100644 --- a/gprMax/input_cmd_funcs.py +++ b/gprMax/input_cmd_funcs.py @@ -35,6 +35,7 @@ Coordinate(x=0.1, y=0.2, z=0.3) from collections import namedtuple Coordinate_tuple = namedtuple('Coordinate', ['x', 'y', 'z']) + class Coordinate(Coordinate_tuple): """Subclass of a namedtuple where __str__ outputs 'x y z'""" def __str__(self): @@ -62,7 +63,7 @@ def command(cmd, *parameters): except TypeError as e: # append info about cmd and parameters to the exception: if not e.args: - e.args=('', ) + e.args = ('', ) additional_info = "Creating cmd = #{} with parameters {} -> {} failed".format(cmd, parameters, filtered_list) e.args = e.args + (additional_info,) raise e diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index 2cbe5247..be7ebec3 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -46,9 +46,9 @@ def process_geometrycmds(geometry, G): if len(tmp) != 6: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly five parameters') - xs = round_value(float(tmp[1])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - zs = round_value(float(tmp[3])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) geofile = tmp[4] matfile = tmp[5] @@ -90,17 +90,16 @@ def process_geometrycmds(geometry, G): if G.messages: print('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile)) - elif tmp[0] == '#edge:': if len(tmp) != 8: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters') - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) if xs < 0 or xs > G.nx: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx)) @@ -150,7 +149,6 @@ def process_geometrycmds(geometry, G): if G.messages: print('Edge from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, tmp[7])) - elif tmp[0] == '#plate:': if len(tmp) < 8: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least seven parameters') @@ -166,12 +164,12 @@ def process_geometrycmds(geometry, G): else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) if xs < 0 or xs > G.nx: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx)) @@ -259,7 +257,6 @@ def process_geometrycmds(geometry, G): if G.messages: print('Plate from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested))) - elif tmp[0] == '#triangle:': if len(tmp) < 12: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eleven parameters') @@ -286,15 +283,15 @@ def process_geometrycmds(geometry, G): else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') - x1 = round_value(float(tmp[1])/G.dx) * G.dx - y1 = round_value(float(tmp[2])/G.dy) * G.dy - z1 = round_value(float(tmp[3])/G.dz) * G.dz - x2 = round_value(float(tmp[4])/G.dx) * G.dx - y2 = round_value(float(tmp[5])/G.dy) * G.dy - z2 = round_value(float(tmp[6])/G.dz) * G.dz - x3 = round_value(float(tmp[7])/G.dx) * G.dx - y3 = round_value(float(tmp[8])/G.dy) * G.dy - z3 = round_value(float(tmp[9])/G.dz) * G.dz + x1 = round_value(float(tmp[1]) / G.dx) * G.dx + y1 = round_value(float(tmp[2]) / G.dy) * G.dy + z1 = round_value(float(tmp[3]) / G.dz) * G.dz + x2 = round_value(float(tmp[4]) / G.dx) * G.dx + y2 = round_value(float(tmp[5]) / G.dy) * G.dy + z2 = round_value(float(tmp[6]) / G.dz) * G.dz + x3 = round_value(float(tmp[7]) / G.dx) * G.dx + y3 = round_value(float(tmp[8]) / G.dy) * G.dy + z3 = round_value(float(tmp[9]) / G.dz) * G.dz thickness = float(tmp[10]) if x1 < 0 or x2 < 0 or x3 < 0 or x1 > G.nx or x2 > G.nx or x3 > G.nx: @@ -379,7 +376,6 @@ def process_geometrycmds(geometry, G): else: print('Triangle with coordinates {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m of material(s) {} created.'.format(x1, y1, z1, x2, y2, z2, x3, y3, z3, ', '.join(materialsrequested))) - elif tmp[0] == '#box:': if len(tmp) < 8: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least seven parameters') @@ -406,12 +402,12 @@ def process_geometrycmds(geometry, G): else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) if xs < 0 or xs > G.nx: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx)) @@ -471,7 +467,6 @@ def process_geometrycmds(geometry, G): dielectricsmoothing = 'off' print('Box from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created, dielectric smoothing is {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested), dielectricsmoothing)) - elif tmp[0] == '#cylinder:': if len(tmp) < 9: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eight parameters') @@ -498,12 +493,12 @@ def process_geometrycmds(geometry, G): else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') - x1 = round_value(float(tmp[1])/G.dx) * G.dx - y1 = round_value(float(tmp[2])/G.dy) * G.dy - z1 = round_value(float(tmp[3])/G.dz) * G.dz - x2 = round_value(float(tmp[4])/G.dx) * G.dx - y2 = round_value(float(tmp[5])/G.dy) * G.dy - z2 = round_value(float(tmp[6])/G.dz) * G.dz + x1 = round_value(float(tmp[1]) / G.dx) * G.dx + y1 = round_value(float(tmp[2]) / G.dy) * G.dy + z1 = round_value(float(tmp[3]) / G.dz) * G.dz + x2 = round_value(float(tmp[4]) / G.dx) * G.dx + y2 = round_value(float(tmp[5]) / G.dy) * G.dy + z2 = round_value(float(tmp[6]) / G.dz) * G.dz r = float(tmp[7]) if r <= 0: @@ -552,7 +547,6 @@ def process_geometrycmds(geometry, G): dielectricsmoothing = 'off' print('Cylinder with face centres {:g}m, {:g}m, {:g}m and {:g}m, {:g}m, {:g}m, with radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(x1, y1, z1, x2, y2, z2, r, ', '.join(materialsrequested), dielectricsmoothing)) - elif tmp[0] == '#cylindrical_sector:': if len(tmp) < 10: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least nine parameters') @@ -586,8 +580,8 @@ def process_geometrycmds(geometry, G): extent2 = float(tmp[5]) thickness = extent2 - extent1 r = float(tmp[6]) - sectorstartangle = 2*np.pi*(float(tmp[7])/360) - sectorangle = 2*np.pi*(float(tmp[8])/360) + sectorstartangle = 2 * np.pi * (float(tmp[7]) / 360) + sectorangle = 2 * np.pi * (float(tmp[8]) / 360) if normal != 'x' and normal != 'y' and normal != 'z': raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the normal direction must be either x, y or z.') @@ -648,21 +642,21 @@ def process_geometrycmds(geometry, G): # yz-plane cylindrical sector if normal == 'x': - ctr1 = round_value(ctr1/G.dy) * G.dy - ctr2 = round_value(ctr2/G.dz) * G.dz - level = round_value(extent1/G.dx) + ctr1 = round_value(ctr1 / G.dy) * G.dy + ctr2 = round_value(ctr2 / G.dz) * G.dz + level = round_value(extent1 / G.dx) # xz-plane cylindrical sector elif normal == 'y': - ctr1 = round_value(ctr1/G.dx) * G.dx - ctr2 = round_value(ctr2/G.dz) * G.dz - level = round_value(extent1/G.dy) + ctr1 = round_value(ctr1 / G.dx) * G.dx + ctr2 = round_value(ctr2 / G.dz) * G.dz + level = round_value(extent1 / G.dy) # xy-plane cylindrical sector elif normal == 'z': - ctr1 = round_value(ctr1/G.dx) * G.dx - ctr2 = round_value(ctr2/G.dy) * G.dy - level = round_value(extent1/G.dz) + ctr1 = round_value(ctr1 / G.dx) * G.dx + ctr2 = round_value(ctr2 / G.dy) * G.dy + level = round_value(extent1 / G.dz) build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID) @@ -672,10 +666,9 @@ def process_geometrycmds(geometry, G): dielectricsmoothing = 'on' else: dielectricsmoothing = 'off' - print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, thickness {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(ctr1, ctr2, r, (sectorstartangle/(2*np.pi))*360, (sectorangle/(2*np.pi))*360, thickness, ', '.join(materialsrequested), dielectricsmoothing)) + print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, thickness {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, thickness, ', '.join(materialsrequested), dielectricsmoothing)) else: - print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, of material(s) {} created.'.format(ctr1, ctr2, r, (sectorstartangle/(2*np.pi))*360, (sectorangle/(2*np.pi))*360, ', '.join(materialsrequested))) - + print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, of material(s) {} created.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, ', '.join(materialsrequested))) elif tmp[0] == '#sphere:': if len(tmp) < 6: @@ -704,9 +697,9 @@ def process_geometrycmds(geometry, G): raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') # Centre of sphere - xc = round_value(float(tmp[1])/G.dx) - yc = round_value(float(tmp[2])/G.dy) - zc = round_value(float(tmp[3])/G.dz) + xc = round_value(float(tmp[1]) / G.dx) + yc = round_value(float(tmp[2]) / G.dy) + zc = round_value(float(tmp[3]) / G.dz) r = float(tmp[4]) # Look up requested materials in existing list of material instances @@ -752,7 +745,6 @@ def process_geometrycmds(geometry, G): dielectricsmoothing = 'off' print('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:': if len(tmp) < 14: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least thirteen parameters') @@ -763,12 +755,12 @@ def process_geometrycmds(geometry, G): else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) if xs < 0 or xs > G.nx: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx)) @@ -836,12 +828,12 @@ def process_geometrycmds(geometry, G): # Only process rough surfaces for this fractal volume if tmp[12] == volume.ID: - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) if xs < 0 or xs > G.nx: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx)) @@ -943,12 +935,12 @@ def process_geometrycmds(geometry, G): # Only process surfaces for this fractal volume if tmp[8] == volume.ID: - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) depth = float(tmp[7]) if xs < 0 or xs > G.nx: @@ -1044,7 +1036,6 @@ def process_geometrycmds(geometry, G): if G.messages: print('Water on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with depth {:g}m, added to {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, filldepth, surface.operatingonID)) - if tmp[0] == '#add_grass:': if len(tmp) < 12: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eleven parameters') @@ -1057,12 +1048,12 @@ def process_geometrycmds(geometry, G): # Only process grass for this fractal volume if tmp[11] == volume.ID: - xs = round_value(float(tmp[1])/G.dx) - xf = round_value(float(tmp[4])/G.dx) - ys = round_value(float(tmp[2])/G.dy) - yf = round_value(float(tmp[5])/G.dy) - zs = round_value(float(tmp[3])/G.dz) - zf = round_value(float(tmp[6])/G.dz) + xs = round_value(float(tmp[1]) / G.dx) + xf = round_value(float(tmp[4]) / G.dx) + ys = round_value(float(tmp[2]) / G.dy) + yf = round_value(float(tmp[5]) / G.dy) + zs = round_value(float(tmp[3]) / G.dz) + zf = round_value(float(tmp[6]) / G.dz) numblades = int(tmp[10]) if xs < 0 or xs > G.nx: @@ -1194,7 +1185,6 @@ def process_geometrycmds(geometry, G): if G.messages: print('{} blades of grass on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(numblades, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, surface.dimension, surface.seed, float(tmp[8]), float(tmp[9]), surface.operatingonID)) - # Process any modifications to the original fractal box then generate it if volume.fractalsurfaces: volume.originalxs = volume.xs @@ -1439,7 +1429,6 @@ def process_geometrycmds(geometry, G): depth += 1 k -= 1 - # Build voxels from any true values of the 3D mask array waternumID = next((x.numID for x in G.materials if x.ID == 'water'), 0) grassnumID = next((x.numID for x in G.materials if x.ID == 'grass'), 0) diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 20495e33..7feb5288 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -71,7 +71,6 @@ def process_multicmds(multicmds, G): G.waveforms.append(w) - # Voltage source cmdname = '#voltage_source' if multicmds[cmdname] != 'None': @@ -98,7 +97,7 @@ def process_multicmds(multicmds, G): raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' there is no waveform with the identifier {}'.format(tmp[5])) v = VoltageSource() - v.polarisation= tmp[0] + v.polarisation = tmp[0] v.xcoord = xcoord v.ycoord = ycoord v.zcoord = zcoord @@ -132,7 +131,6 @@ def process_multicmds(multicmds, G): G.voltagesources.append(v) - # Hertzian dipole cmdname = '#hertzian_dipole' if multicmds[cmdname] != 'None': @@ -192,7 +190,6 @@ def process_multicmds(multicmds, G): G.hertziandipoles.append(h) - # Magnetic dipole cmdname = '#magnetic_dipole' if multicmds[cmdname] != 'None': @@ -325,9 +322,9 @@ def process_multicmds(multicmds, G): raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' has an incorrect number of parameters') # Check position parameters - xcoord = round_value(float(tmp[0])/G.dx) - ycoord = round_value(float(tmp[1])/G.dy) - zcoord = round_value(float(tmp[2])/G.dz) + xcoord = round_value(float(tmp[0]) / G.dx) + ycoord = round_value(float(tmp[1]) / G.dy) + zcoord = round_value(float(tmp[2]) / G.dz) check_coordinates(xcoord, ycoord, zcoord) if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]: print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n') @@ -359,7 +356,6 @@ def process_multicmds(multicmds, G): G.rxs.append(r) - # Receiver array cmdname = '#rx_array' if multicmds[cmdname] != 'None': @@ -425,7 +421,6 @@ def process_multicmds(multicmds, G): print(' Receiver at {:g}m, {:g}m, {:g}m with output(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs))) G.rxs.append(r) - # Snapshot cmdname = '#snapshot' if multicmds[cmdname] != 'None': @@ -511,7 +506,6 @@ def process_multicmds(multicmds, G): # Append the new material object to the materials list G.materials.append(m) - cmdname = '#add_dispersion_debye' if multicmds[cmdname] != 'None': for cmdinstance in multicmds[cmdname]: @@ -583,7 +577,6 @@ def process_multicmds(multicmds, G): if G.messages: print('Lorentz disperion added to {} with delta_epsr={}, omega={} secs, and gamma={} created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha))) - cmdname = '#add_dispersion_drude' if multicmds[cmdname] != 'None': for cmdinstance in multicmds[cmdname]: @@ -619,7 +612,6 @@ def process_multicmds(multicmds, G): if G.messages: print('Drude disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha))) - cmdname = '#soil_peplinski' if multicmds[cmdname] != 'None': for cmdinstance in multicmds[cmdname]: @@ -650,7 +642,6 @@ def process_multicmds(multicmds, G): # Append the new material object to the materials list G.mixingmodels.append(s) - # Geometry views (creates VTK-based geometry files) cmdname = '#geometry_view' if multicmds[cmdname] != 'None': @@ -695,7 +686,6 @@ def process_multicmds(multicmds, G): # Append the new GeometryView object to the geometry views list G.geometryviews.append(g) - # Complex frequency shifted (CFS) PML parameter cmdname = '#pml_cfs' if multicmds[cmdname] != 'None': diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/input_cmds_singleuse.py index 57a5f058..9466a8e1 100644 --- a/gprMax/input_cmds_singleuse.py +++ b/gprMax/input_cmds_singleuse.py @@ -50,7 +50,6 @@ def process_singlecmds(singlecmds, G): else: raise CmdInputError(cmd + ' requires input values of either y or n') - # Title cmd = '#title' if singlecmds[cmd] != 'None': @@ -58,7 +57,6 @@ def process_singlecmds(singlecmds, G): if G.messages: print('Model title: {}'.format(G.title)) - # Number of processors to run on (OpenMP) cmd = '#num_threads' os.environ['OMP_WAIT_POLICY'] = 'active' @@ -85,7 +83,6 @@ def process_singlecmds(singlecmds, G): if G.nthreads > psutil.cpu_count(logical=False): print('\nWARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, psutil.cpu_count(logical=False))) - # Spatial discretisation cmd = '#dx_dy_dz' tmp = [float(x) for x in singlecmds[cmd].split()] @@ -103,15 +100,14 @@ def process_singlecmds(singlecmds, G): if G.messages: print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz)) - # Domain cmd = '#domain' tmp = [float(x) for x in singlecmds[cmd].split()] if len(tmp) != 3: raise CmdInputError(cmd + ' requires exactly three parameters') - G.nx = round_value(tmp[0]/G.dx) - G.ny = round_value(tmp[1]/G.dy) - G.nz = round_value(tmp[2]/G.dz) + G.nx = round_value(tmp[0] / G.dx) + G.ny = round_value(tmp[1] / G.dy) + G.nz = round_value(tmp[2] / G.dz) if G.nx == 0 or G.ny == 0 or G.nz == 0: raise CmdInputError(cmd + ' requires at least one cell in every dimension') if G.messages: @@ -120,7 +116,6 @@ def process_singlecmds(singlecmds, G): mem = (((G.nx + 1) * (G.ny + 1) * (G.nz + 1) * 13 * np.dtype(floattype).itemsize + (G.nx + 1) * (G.ny + 1) * (G.nz + 1) * 18) * 1.1) + 30e6 print('Memory (RAM) usage: ~{} required, {} available'.format(human_size(mem), human_size(psutil.virtual_memory().total))) - # Time step CFL limit (use either 2D or 3D) and default PML thickness if G.nx == 1: G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz))) @@ -145,7 +140,6 @@ def process_singlecmds(singlecmds, G): if G.messages: print('Time step (at {} CFL limit): {:g} secs'.format(G.dtlimit, G.dt)) - # Time step stability factor cmd = '#time_step_stability_factor' if singlecmds[cmd] != 'None': @@ -158,7 +152,6 @@ def process_singlecmds(singlecmds, G): if G.messages: print('Time step (modified): {:g} secs'.format(G.dt)) - # Time window cmd = '#time_window' tmp = singlecmds[cmd].split() @@ -182,7 +175,6 @@ def process_singlecmds(singlecmds, G): if G.messages: print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations)) - # PML cmd = '#pml_cells' if singlecmds[cmd] != 'None': @@ -193,36 +185,33 @@ def process_singlecmds(singlecmds, G): G.pmlthickness = (int(tmp[0]), int(tmp[0]), int(tmp[0]), int(tmp[0]), int(tmp[0]), int(tmp[0])) else: G.pmlthickness = (int(tmp[0]), int(tmp[1]), int(tmp[2]), int(tmp[3]), int(tmp[4]), int(tmp[5])) - if 2*G.pmlthickness[0] >= G.nx or 2*G.pmlthickness[1] >= G.ny or 2*G.pmlthickness[2] >= G.nz or 2*G.pmlthickness[3] >= G.nx or 2*G.pmlthickness[4] >= G.ny or 2*G.pmlthickness[5] >= G.nz: + if 2 * G.pmlthickness[0] >= G.nx or 2 * G.pmlthickness[1] >= G.ny or 2 * G.pmlthickness[2] >= G.nz or 2 * G.pmlthickness[3] >= G.nx or 2 * G.pmlthickness[4] >= G.ny or 2 * G.pmlthickness[5] >= G.nz: raise CmdInputError(cmd + ' has too many cells for the domain size') - # src_steps cmd = '#src_steps' if singlecmds[cmd] != 'None': tmp = singlecmds[cmd].split() if len(tmp) != 3: raise CmdInputError(cmd + ' requires exactly three parameters') - G.srcstepx = round_value(float(tmp[0])/G.dx) - G.srcstepy = round_value(float(tmp[1])/G.dy) - G.srcstepz = round_value(float(tmp[2])/G.dz) + G.srcstepx = round_value(float(tmp[0]) / G.dx) + G.srcstepy = round_value(float(tmp[1]) / G.dy) + G.srcstepz = round_value(float(tmp[2]) / G.dz) if G.messages: print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcstepx * G.dx, G.srcstepy * G.dy, G.srcstepz * G.dz)) - # rx_steps cmd = '#rx_steps' if singlecmds[cmd] != 'None': tmp = singlecmds[cmd].split() if len(tmp) != 3: raise CmdInputError(cmd + ' requires exactly three parameters') - G.rxstepx = round_value(float(tmp[0])/G.dx) - G.rxstepy = round_value(float(tmp[1])/G.dy) - G.rxstepz = round_value(float(tmp[2])/G.dz) + G.rxstepx = round_value(float(tmp[0]) / G.dx) + G.rxstepy = round_value(float(tmp[1]) / G.dy) + G.rxstepz = round_value(float(tmp[2]) / G.dz) if G.messages: print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxstepx * G.dx, G.rxstepy * G.dy, G.rxstepz * G.dz)) - # Excitation file for user-defined source waveforms cmd = '#excitation_file' if singlecmds[cmd] != 'None': @@ -251,7 +240,7 @@ def process_singlecmds(singlecmds, G): if len(waveformvalues.shape) == 1: w.uservalues = waveformvalues[:] else: - w.uservalues = waveformvalues[:,waveform] + w.uservalues = waveformvalues[:, waveform] if G.messages: print('User waveform {} created.'.format(w.ID)) diff --git a/gprMax/materials.py b/gprMax/materials.py index fa6efc68..45cf6bfa 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -74,8 +74,8 @@ class Material(object): G (class): Grid class instance - holds essential parameters describing the model. """ - HA = (m0*self.mr / G.dt) + 0.5*self.sm - HB = (m0*self.mr / G.dt) - 0.5*self.sm + HA = (m0 * self.mr / G.dt) + 0.5 * self.sm + HB = (m0 * self.mr / G.dt) - 0.5 * self.sm self.DA = HB / HA self.DBx = (1 / G.dx) * 1 / HA self.DBy = (1 / G.dy) * 1 / HA @@ -122,12 +122,12 @@ class Material(object): self.zt[x] = (self.w[x] / self.q[x]) * (1 - self.eqt[x]) / G.dt self.zt2[x] = (self.w[x] / self.q[x]) * (1 - self.eqt2[x]) - EA = (e0*self.er / G.dt) + 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real) - EB = (e0*self.er / G.dt) - 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real) + EA = (e0 * self.er / G.dt) + 0.5 * self.se - (e0 / G.dt) * np.sum(self.zt2.real) + EB = (e0 * self.er / G.dt) - 0.5 * self.se - (e0 / G.dt) * np.sum(self.zt2.real) else: - EA = (e0*self.er / G.dt) + 0.5*self.se - EB = (e0*self.er / G.dt) - 0.5*self.se + EA = (e0 * self.er / G.dt) + 0.5 * self.se + EB = (e0 * self.er / G.dt) - 0.5 * self.se if self.ID == 'pec' or self.se == float('inf'): self.CA = 0 @@ -198,15 +198,15 @@ class PeplinskiSoil(object): muiter = np.nditer(mumaterials, flags=['c_index']) while not muiter.finished: # Real part for frequencies in the range 1.4GHz to 18GHz - er1 = (1 + (self.rb/self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1/a) + er1 = (1 + (self.rb / self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1 / a) # Real part for frequencies in the range 0.3GHz to 1.3GHz er2 = 1.15 * er1 - 0.68 # Imaginary part for frequencies in the range 0.3GHz to 1.3GHz - eri = er2 - (muiter[0]**(b2/a) * Material.waterdeltaer) + eri = er2 - (muiter[0]**(b2 / a) * Material.waterdeltaer) # Effective conductivity - sig = muiter[0]**(b2/a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0])) + sig = muiter[0]**(b2 / a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0])) # Check to see if the material already exists before creating a new one requiredID = '|{:.4f}|'.format(float(muiter[0])) diff --git a/gprMax/optimisation_taguchi.py b/gprMax/optimisation_taguchi.py index 783f2e1e..102eb7fc 100644 --- a/gprMax/optimisation_taguchi.py +++ b/gprMax/optimisation_taguchi.py @@ -68,7 +68,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): # Select OA OA, N, cols, k, s, t = construct_OA(optparams) - print('\n{}\nTaguchi optimisation...\n'.format(68*'*')) + print('\n{}\nTaguchi optimisation...\n'.format(68 * '*')) print('\tOrthogonal array: {:g} experiments per iteration, {:g} parameters ({:g} will be used), {:g} levels, and strength {:g}'.format(N, cols, k, s, t)) tmp = [(k, v) for k, v in optparams.items()] print('\tParameters to optimise with ranges: {}'.format(str(tmp).strip('[]'))) @@ -131,7 +131,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): # Rename confirmation experiment output file so that it is retained for each iteraction os.rename(outputfile, os.path.splitext(outputfile)[0] + '_final' + str(iteration + 1) + '.out') - print('\nTaguchi optimisation, iteration {} completed. History of optimal parameter values {} and of fitness values {}'.format(iteration + 1, dict(optparamshist), fitnessvalueshist, 68*'*')) + print('\nTaguchi optimisation, iteration {} completed. History of optimal parameter values {} and of fitness values {}'.format(iteration + 1, dict(optparamshist), fitnessvalueshist, 68 * '*')) iteration += 1 # Stop optimisation if stopping criterion has been reached @@ -154,7 +154,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace): pickle.dump(fitnessvalueshist, f) pickle.dump(optparamsinit, f) - print('\n{}\nTaguchi optimisation completed after {} iteration(s).\nHistory of optimal parameter values {} and of fitness values {}\n{}\n'.format(68*'*', iteration, dict(optparamshist), fitnessvalueshist, 68*'*')) + print('\n{}\nTaguchi optimisation completed after {} iteration(s).\nHistory of optimal parameter values {} and of fitness values {}\n{}\n'.format(68 * '*', iteration, dict(optparamshist), fitnessvalueshist, 68 * '*')) def taguchi_code_blocks(inputfile, taguchinamespace): @@ -313,7 +313,7 @@ def calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, le # Gaussian reduction function used for calculating levels T = 18 # Usually values between 15 - 20 - RR = np.exp(-(i/T)**2) + RR = np.exp(-(i / T)**2) # Calculate levels for each parameter for p in range(k): @@ -341,7 +341,6 @@ def calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, le levels[0, p] = levels[1, p] - levelsdiff[p] levels[2, p] = levels[1, p] + levelsdiff[p] - # Update dictionary of parameters to optimise with lists of new values; clear dictionary first optparams = OrderedDict((key, list()) for key in optparams) p = 0 diff --git a/gprMax/pml.py b/gprMax/pml.py index 52ea5b36..4c94e476 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -98,7 +98,7 @@ class CFS(object): Hvalues (float): numpy array holding scaling profile values for magnetic PML update. """ - tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2*len(Evalues)) / (len(Evalues) - 1)) ** order + tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2 * len(Evalues)) / (len(Evalues) - 1)) ** order Evalues = tmp[0:-1:2] Hvalues = tmp[1::2] return Evalues, Hvalues @@ -115,7 +115,7 @@ class CFS(object): Hvalues (float): numpy array holding profile value for magnetic PML update. """ - Evalues= np.zeros(thickness + 1, dtype=floattype) + Evalues = np.zeros(thickness + 1, dtype=floattype) Hvalues = np.zeros(thickness + 1, dtype=floattype) if parameter.scalingprofile == 'constant': @@ -218,18 +218,18 @@ class PML(object): Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma) # Electric PML update coefficients - tmp = (2*e0*Ekappa) + G.dt * (Ealpha * Ekappa + Esigma) - self.ERA[x, :] = (2*e0 + G.dt*Ealpha) / tmp - self.ERB[x, :] = (2*e0*Ekappa) / tmp - self.ERE[x, :] = ((2*e0*Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp - self.ERF[x, :] = (2*Esigma*G.dt) / (Ekappa * tmp) + tmp = (2 * e0 * Ekappa) + G.dt * (Ealpha * Ekappa + Esigma) + self.ERA[x, :] = (2 * e0 + G.dt * Ealpha) / tmp + self.ERB[x, :] = (2 * e0 * Ekappa) / tmp + self.ERE[x, :] = ((2 * e0 * Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp + self.ERF[x, :] = (2 * Esigma * G.dt) / (Ekappa * tmp) # Magnetic PML update coefficients - tmp = (2*e0*Hkappa) + G.dt * (Halpha * Hkappa + Hsigma) - self.HRA[x, :] = (2*e0 + G.dt*Halpha) / tmp - self.HRB[x, :] = (2*e0*Hkappa) / tmp - self.HRE[x, :] = ((2*e0*Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp - self.HRF[x, :] = (2*Hsigma*G.dt) / (Hkappa * tmp) + tmp = (2 * e0 * Hkappa) + G.dt * (Halpha * Hkappa + Hsigma) + self.HRA[x, :] = (2 * e0 + G.dt * Halpha) / tmp + self.HRB[x, :] = (2 * e0 * Hkappa) / tmp + self.HRE[x, :] = ((2 * e0 * Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp + self.HRF[x, :] = (2 * Hsigma * G.dt) / (Hkappa * tmp) def build_pmls(G): @@ -244,7 +244,7 @@ def build_pmls(G): if pmlthickness > 0: sumer = 0 # Sum of relative permittivities in PML slab summr = 0 # Sum of relative permeabilities in PML slab - pmldirection=PML.directions[index] + pmldirection = PML.directions[index] if pmldirection[0] == 'x': if pmldirection == 'xminus': diff --git a/gprMax/sources.py b/gprMax/sources.py index 66b3e63d..7e3fabbe 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -185,13 +185,13 @@ class MagneticDipole(Source): waveform = next(x for x in G.waveforms if x.ID == self.waveformID) if self.polarisation is 'x': - Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) + Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) elif self.polarisation is 'y': - Hy[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) + Hy[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) elif self.polarisation is 'z': - Hz[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) + Hz[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz)) class TransmissionLine(Source): diff --git a/gprMax/utilities.py b/gprMax/utilities.py index 5ba3f364..d47e431d 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -57,9 +57,9 @@ along with gprMax. If not, see .""" url = 'www.gprmax.com' print('\nElectromagnetic modelling software based on the Finite-Difference \nTime-Domain (FDTD) method') - print('\n{} {} {}'.format('*'*round((width - len(url))/2), url, '*'*round((width - len(url))/2))) + print('\n{} {} {}'.format('*' * round((width - len(url)) / 2), url, '*' * round((width - len(url)) / 2))) print('{}'.format(gprMaxlogo)) - print('{} v{} {}'.format('*'*round((width - len(version))/2), (version), '*'*round((width - len(version))/2))) + print('{} v{} {}'.format('*' * round((width - len(version)) / 2), (version), '*' * round((width - len(version)) / 2))) print(licenseinfo) @@ -91,12 +91,12 @@ def round_value(value, decimalplaces=0): # Rounds to nearest integer (half values are rounded downwards) if decimalplaces == 0: - rounded = int(d.Decimal(value).quantize(d.Decimal('1'),rounding=d.ROUND_HALF_DOWN)) + rounded = int(d.Decimal(value).quantize(d.Decimal('1'), rounding=d.ROUND_HALF_DOWN)) # Rounds down to nearest float represented by number of decimal places else: precision = '1.{places}'.format(places='0' * decimalplaces) - rounded = float(d.Decimal(value).quantize(d.Decimal(precision),rounding=d.ROUND_FLOOR)) + rounded = float(d.Decimal(value).quantize(d.Decimal(precision), rounding=d.ROUND_FLOOR)) return rounded diff --git a/gprMax/writer_hdf5.py b/gprMax/writer_hdf5.py index 2ce2ffa0..26206c22 100644 --- a/gprMax/writer_hdf5.py +++ b/gprMax/writer_hdf5.py @@ -43,7 +43,6 @@ def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G): f.attrs['srcsteps'] = (G.srcstepx, G.srcstepy, G.srcstepz) f.attrs['rxsteps'] = (G.rxstepx, G.rxstepy, G.rxstepz) - # Create group for sources (except transmission lines); add type and positional data attributes srclist = G.voltagesources + G.hertziandipoles + G.magneticdipoles for srcindex, src in enumerate(srclist): @@ -65,7 +64,6 @@ def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G): f['/tls/tl' + str(tlindex + 1) + '/Vtotal'] = tl.Vtotal f['/tls/tl' + str(tlindex + 1) + '/Itotal'] = tl.Itotal - # Create group, add positional data and write field component arrays for receivers for rxindex, rx in enumerate(G.rxs): grp = f.create_group('/rxs/rx' + str(rxindex + 1)) diff --git a/gprMax/xdmf.py b/gprMax/xdmf.py index 92841401..07ef71b2 100644 --- a/gprMax/xdmf.py +++ b/gprMax/xdmf.py @@ -6,6 +6,7 @@ import numpy as np from gprMax.grid import Grid + class ListCounter(): def __init__(self, npArray): @@ -16,10 +17,10 @@ class ListCounter(): self.array[self.count] = item self.count += 1 + class EdgeLabels: def __init__(self, grid): - """ Class to define some connectivity of for an n x l x m grid @@ -37,6 +38,7 @@ class EdgeLabels: edge = np.array([in_label, out_label]) self.edge_counter.add(edge) + class EdgeMaterials: def __init__(self, fdtd_grid): @@ -85,7 +87,6 @@ class SolidLabels(): self.label_counter = ListCounter(self.solid_labels) def hexCellPicker(self, grid, i, j, k): - """ This is the ordering of nodes in the hexahedron cell. @@ -124,6 +125,7 @@ class SolidLabels(): solid_labels = self.hexCellPicker(self.label_grid.grid, i, j, k) self.label_counter.add(solid_labels) + class SolidManager(Grid): def __init__(self, label_grid, fdtd_grid):