working

gprMax/cmds_geometry/add_grass.py

@@ -160,5 +160,5 @@ class AddGrass(UserObjectGeometry):
 
         volume.fractalsurfaces.append(surface)
 
-        if grid.messages:
+        if config.is_messages():
             tqdm.write('{} 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(n_blades, xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, surface.dimension, surface.seed, limits[0], limits[1], surface.operatingonID))

gprMax/cmds_geometry/add_surface_roughness.py

@@ -123,5 +123,5 @@ class AddSurfaceRoughness(UserObjectGeometry):
         surface.generate_fractal_surface(grid)
         volume.fractalsurfaces.append(surface)
 
-        if config.general['messages']:
+        if config.is_messages():
             tqdm.write('Fractal surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal weightings {:g}, {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, surface.dimension, surface.weighting[0], surface.weighting[1], surface.seed, limits[0], limits[1], surface.operatingonID))

gprMax/cmds_geometry/add_surface_water.py

@@ -114,5 +114,5 @@ class AddSurfaceWater(UserObjectGeometry):
         if testwater:
             raise CmdInputError(self.__str__() + ' requires the time step for the model to be less than the relaxation time required to model water.')
 
-        if grid.messages:
+        if config.is_messages():
             tqdm.write('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 * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, filldepth, surface.operatingonID))

gprMax/cmds_geometry/box.py

@@ -3,6 +3,7 @@ from .cmds_geometry import UserObjectGeometry
 from ..exceptions import CmdInputError
 from ..materials import Material
 from ..cython.geometry_primitives import build_box
+import gprMax.config as config
 
 from tqdm import tqdm
 import numpy as np
@@ -84,7 +85,7 @@ class Box(UserObjectGeometry):
 
         build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             if averaging:
                 dielectricsmoothing = 'on'
             else:

gprMax/cmds_geometry/cylinder.py

@@ -6,6 +6,7 @@ from ..cython.geometry_primitives import build_cylinder
 
 from tqdm import tqdm
 import numpy as np
+import gprMax.config as config
 
 
 class Cylinder(UserObjectGeometry):
@@ -89,7 +90,7 @@ class Cylinder(UserObjectGeometry):
 
         build_cylinder(x1, y1, z1, x2, y2, z2, r, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             if averaging:
                 dielectricsmoothing = 'on'
             else:

gprMax/cmds_geometry/cylindrical_sector.py

@@ -126,7 +126,7 @@ class CylindricalSector(UserObjectGeometry):
 
         build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             if thickness > 0:
                 if averaging:
                     dielectricsmoothing = 'on'

gprMax/cmds_geometry/edge.py

@@ -60,5 +60,5 @@ class Edge(UserObjectGeometry):
                 for k in range(zs, zf):
                     build_edge_z(xs, ys, k, material.numID, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             tqdm.write('Edge from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, material_id))

gprMax/cmds_geometry/fractal_box.py

@@ -25,7 +25,7 @@ class FractalBox(UserObjectGeometry):
             weighting = np.array(self.kwargs['weighting'])
             n_materials = self.kwargs['n_materials']
             mixing_model_id = self.kwargs['mixing_model_id']
-            ID = self.kwargs['ID']
+            ID = self.kwargs['id']
 
         except KeyError:
             raise CmdInputError(self.__str__() + ' Incorrect parameters')
@@ -83,7 +83,7 @@ class FractalBox(UserObjectGeometry):
         volume.averaging = averagefractalbox
         volume.mixingmodel = mixingmodel
 
-        if config.general['messages']:
+        if config.is_messages():
             if volume.averaging:
                 dielectricsmoothing = 'on'
             else:

gprMax/cmds_geometry/plate.py

@@ -108,5 +108,5 @@ class Plate(UserObjectGeometry):
                 for j in range(ys, yf):
                     build_face_xy(i, j, zs, numIDx, numIDy, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             tqdm.write('Plate from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, ', '.join(materialsrequested)))

gprMax/cmds_geometry/sphere.py

@@ -83,7 +83,7 @@ class Sphere(UserObjectGeometry):
 
         build_sphere(xc, yc, zc, r, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             if averaging:
                 dielectricsmoothing = 'on'
             else:

gprMax/cmds_geometry/triangle.py

@@ -121,7 +121,7 @@ class Triangle(UserObjectGeometry):
 
         build_triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, normal, thickness, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
 
-        if grid.messages:
+        if config.is_messages():
             if thickness > 0:
                 if averaging:
                     dielectricsmoothing = 'on'

gprMax/cmds_multiple.py

@@ -106,7 +106,7 @@ class Waveform(UserObjectMulti):
         w.amp = amp
         w.freq = freq
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Waveform {} of type {} with maximum amplitude scaling {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq))
 
         grid.waveforms.append(w)
@@ -183,7 +183,7 @@ class VoltageSource(UserObjectMulti):
 
         v.calculate_waveform_values(grid)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * grid.dx, v.ycoord * grid.dy, v.zcoord * grid.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID))
 
         grid.voltagesources.append(v)
@@ -268,7 +268,7 @@ class HertzianDipole(UserObjectMulti):
 
         h.calculate_waveform_values(grid)
 
-        if config.general['messages']:
+        if config.is_messages():
             if grid.mode == '2D':
                 print('Hertzian dipole is a line source in 2D with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * grid.dx, h.ycoord * grid.dy, h.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
             else:
@@ -345,7 +345,7 @@ class MagneticDipole(UserObjectMulti):
 
         m.calculate_waveform_values(grid)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Magnetic dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(m.polarisation, m.xcoord * grid.dx, m.ycoord * grid.dy, m.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(m.waveformID))
 
         grid.magneticdipoles.append(m)
@@ -427,7 +427,7 @@ class TransmissionLine(UserObjectMulti):
         t.calculate_waveform_values(grid)
         t.calculate_incident_V_I(grid)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Transmission line with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(t.polarisation, t.xcoord * grid.dx, t.ycoord * grid.dy, t.zcoord * grid.dz, t.resistance) + startstop + 'using waveform {} created.'.format(t.waveformID))
 
         grid.transmissionlines.append(t)
@@ -456,7 +456,7 @@ class Rx(UserObjectMulti):
         r.xcoordorigin, r.ycoordorigin, r.zcoordorigin = p
 
         try:
-            r.ID = self.kwargs['ID']
+            r.ID = self.kwargs['id']
             outputs = self.kwargs['outputs']
             # Get allowable outputs
             if grid.gpu is not None:
@@ -475,8 +475,7 @@ class Rx(UserObjectMulti):
             r.ID = r.__class__.__name__ + '(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')'
             for key in RxUser.defaultoutputs:
                 r.outputs[key] = np.zeros(grid.iterations, dtype=floattype)
-
+        if config.is_messages():
-        if config.general['messages']:
             print('Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs)))
 
         grid.rxs.append(r)
@@ -526,7 +525,7 @@ class RxArray(UserObjectMulti):
             else:
                 raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Receiver array {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g}m'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz))
 
         for x in range(xs, xf + 1, dx):
@@ -542,7 +541,7 @@ class RxArray(UserObjectMulti):
                     r.ID = r.__class__.__name__ + '(' + str(x) + ',' + str(y) + ',' + str(z) + ')'
                     for key in RxUser.defaultoutputs:
                         r.outputs[key] = np.zeros(grid.iterations, dtype=floattype)
-                    if config.general['messages']:
+                    if config.is_messages():
                         print('  Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs)))
                     grid.rxs.append(r)
 
@@ -601,7 +600,7 @@ class Snapshot(UserObjectMulti):
         else:
             s = SnapshotUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, filename)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Snapshot from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, at {:g} secs with filename {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, s.time * grid.dt, s.basefilename))
 
         grid.snapshots.append(s)
@@ -651,7 +650,7 @@ class Material(UserObjectMulti):
         if m.se == float('inf'):
             m.averagable = False
 
-        if config.general['messages']:
+        if config.is_messages():
             tqdm.write('Material {} with eps_r={:g}, sigma={:g} S/m; mu_r={:g}, sigma*={:g} Ohm/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm))
 
         # Append the new material object to the materials list
@@ -702,7 +701,7 @@ class AddDebyeDispersion(UserObjectMulti):
             if material.poles > MaterialUser.maxpoles:
                 MaterialUser.maxpoles = material.poles
 
-            if config.general['messages']:
+            if config.is_messages():
                 tqdm.write('Debye disperion added to {} with delta_eps_r={}, and tau={} secs created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau)))
 
 
@@ -750,7 +749,7 @@ class AddLorentzDispersion(UserObjectMulti):
             if material.poles > MaterialUser.maxpoles:
                 MaterialUser.maxpoles = material.poles
 
-            if config.general['messages']:
+            if config.is_messages():
                 tqdm.write('Lorentz disperion added to {} with delta_eps_r={}, 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)))
 
 
@@ -795,7 +794,7 @@ class AddDrudeDispersion(UserObjectMulti):
             if material.poles > MaterialUser.maxpoles:
                 MaterialUser.maxpoles = material.poles
 
-            if config.general['messages']:
+            if config.is_messages():
                 tqdm.write('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)))
 
 
@@ -817,7 +816,7 @@ class SoilPeplinski(UserObjectMulti):
             sand_density = self.kwargs['sand_density']
             water_fraction_lower = self.kwargs['water_fraction_lower']
             water_fraction_upper = self.kwargs['water_fraction_upper']
-            ID = self.kwargs['ID']
+            ID = self.kwargs['id']
 
         except KeyError:
             raise CmdInputError(self.__str__() + ' requires at exactly seven parameters')
@@ -840,7 +839,7 @@ class SoilPeplinski(UserObjectMulti):
         # Create a new instance of the Material class material (start index after pec & free_space)
         s = PeplinskiSoilUser(ID, sand_fraction, clay_fraction, bulk_density, sand_density, (water_fraction_lower, water_fraction_upper))
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Mixing model (Peplinski) used to create {} with sand fraction {:g}, clay fraction {:g}, bulk density {:g}g/cm3, sand particle density {:g}g/cm3, and water volumetric fraction {:g} to {:g} created.'.format(s.ID, s.S, s.C, s.rb, s.rs, s.mu[0], s.mu[1]))
 
         # Append the new material object to the materials list
@@ -912,7 +911,7 @@ class GeometryView(UserObjectMulti):
 
         g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, grid)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, multi_grid {}, grid={}, with filename base {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, self.multi_grid, grid.name, g.basefilename))
 
         # Append the new GeometryView object to the geometry views list
@@ -942,7 +941,7 @@ class GeometryObjectsWrite(UserObjectMulti):
 
         g = GeometryObjectsUser(x0, y0, z0, x1, y1, z1, filename)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Geometry objects in the volume from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, will be written to {}, with materials written to {}'.format(p1[0] * grid.dx, p1[1] * grid.dy, p1[2] * grid.dz, p2[0] * grid.dx, p2[1] * grid.dy, p2[2] * grid.dz, g.filename, g.materialsfilename))
 
         # Append the new GeometryView object to the geometry objects to write list
@@ -1017,7 +1016,7 @@ class PMLCFS(UserObjectMulti):
         cfs.kappa = cfskappa
         cfs.sigma = cfssigma
 
-        if config.general['messages']:
+        if config.is_messages():
             print('PML CFS parameters: alpha (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), kappa (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), sigma (scaling: {}, scaling direction: {}, min: {:g}, max: {}) created.'.format(cfsalpha.scalingprofile, cfsalpha.scalingdirection, cfsalpha.min, cfsalpha.max, cfskappa.scalingprofile, cfskappa.scalingdirection, cfskappa.min, cfskappa.max, cfssigma.scalingprofile, cfssigma.scalingdirection, cfssigma.min, cfssigma.max))
 
         grid.cfs.append(cfs)

gprMax/cmds_single_use.py

@@ -86,7 +86,7 @@ class DomainSingle(UserObjectSingle):
 
         if G.nx == 0 or G.ny == 0 or G.nz == 0:
             raise CmdInputError(self.__str__ + ' requires at least one cell in every dimension')
-        if config.general['messages']:
+        if config.is_messages():
             print('Domain size: {:g} x {:g} x {:g}m ({:d} x {:d} x {:d} = {:g} cells)'.format(self.kwargs['p1'][0], self.kwargs['p1'][1], self.kwargs['p1'][2], G.nx, G.ny, G.nz, (G.nx * G.ny * G.nz)))
 
         # Time step CFL limit (either 2D or 3D); switch off appropriate PMLs for 2D
@@ -114,7 +114,7 @@ class DomainSingle(UserObjectSingle):
         # Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
         G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Mode: {}'.format(G.mode))
             print('Time step (at CFL limit): {:g} secs'.format(G.dt))
 
@@ -139,7 +139,7 @@ class DomainSingle(UserObjectSingle):
             G.nthreads = hostinfo['physicalcores']
             os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
         if G.nthreads > hostinfo['physicalcores']:
             print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
@@ -194,7 +194,7 @@ class Discretisation(UserObjectSingle):
         if G.dl[2] <= 0:
             raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero')
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl))
 
 
@@ -243,7 +243,7 @@ class TimeWindow(UserObjectSingle):
         if not G.timewindow:
             raise CmdInputError('TimeWindow: Specify a time or number of iterations')
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
 
 
@@ -291,7 +291,7 @@ class Title(UserObjectSingle):
         except KeyError:
             pass
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Model title: {}'.format(G.title))
 
 class NumThreads(UserObjectSingle):
@@ -321,13 +321,13 @@ class NumThreads(UserObjectSingle):
         G.nthreads = n
         os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
 
-        if config.general['messages']:
+        if config.is_messages():
             print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
         if G.nthreads > hostinfo['physicalcores']:
             print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
 
         # Print information about any GPU in use
-        if config.general['messages']:
+        if config.is_messages():
             if G.gpu is not None:
                 print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))
 
@@ -358,7 +358,7 @@ class TimeStepStabilityFactor(UserObjectSingle):
         if f <= 0 or f > 1:
             raise CmdInputError(self.__str__() + ' requires the value of the time step stability factor to be between zero and one')
         G.dt = G.dt * f
-        if config.general['messages']:
+        if config.is_messages():
             print('Time step (modified): {:g} secs'.format(G.dt))
 
 
@@ -412,7 +412,7 @@ class SrcSteps(UserObjectSingle):
         except KeyError:
             raise CmdInputError('#src_steps: requires exactly three parameters')
         # src_steps
-        if config.general['messages']:
+        if config.is_messages():
             print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcsteps[0] * G.dx, G.srcsteps[1] * G.dy, G.srcsteps[2] * G.dz))
 
 
@@ -429,7 +429,7 @@ class RxSteps(UserObjectSingle):
             G.rxsteps = uip.discretise_point(self.kwargs['p1'])
         except KeyError:
             raise CmdInputError('#rx_steps: requires exactly three parameters')
-        if config.general['messages']:
+        if config.is_messages():
             print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxsteps[0] * G.dx, G.rxsteps[1] * G.dy, G.rxsteps[2] * G.dz))
 
 
@@ -455,7 +455,7 @@ class ExcitationFile(UserObjectSingle):
             if not os.path.isfile(excitationfile):
                 excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
 
-            if config.general['messages']:
+            if config.is_messages():
                 print('\nExcitation file: {}'.format(excitationfile))
 
             # Get waveform names
@@ -495,7 +495,7 @@ class ExcitationFile(UserObjectSingle):
                 # Interpolate waveform values
                 w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)
 
-                if config.general['messages']:
+                if config.is_messages():
                     print('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))
 
                 G.waveforms.append(w)

gprMax/config.py

@@ -51,6 +51,12 @@ z0 = np.sqrt(m0 / e0)
 general = {'inputfilepath': 'gprMax', 'outputfilepath': 'gprMax', 'messages': True,
            'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False, 'autotranslate': False}
 
+
+def is_messages():
+    """Function to return messages."""
+    return general['messages']
+
+
 # Store information about host machine
 hostinfo = get_host_info()
 

gprMax/cython/fields_updates_hsg.pyx

@@ -1,8 +1,24 @@
+# Copyright (C) 2015-2017: The University of Edinburgh
+#                 Authors: Craig Warren and Antonis Giannopoulos
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 import numpy as np
 cimport numpy as np
 from cython.parallel import prange
 
-
 cpdef void cython_update_electric_os(
         np.float64_t[:, :] updatecoeffsE,
         np.uint32_t[:, :, :, :] ID,
@@ -208,7 +224,7 @@ cpdef void cython_update_is(
         int nwx,
         int nwy,
         int nwz,
-        double[:, :] updatecoeffsE,
+        np.float64_t[:, :] updatecoeffsE,
         np.uint32_t[:, :, :, :] ID,
         int n,
         int offset,
@@ -216,9 +232,9 @@ cpdef void cython_update_is(
         int nwm,
         int nwn,
         int face,
-        double[:, :, :] field,
+        np.float64_t[:, :, :] field,
-        double[:, :] inc_field_l,
+        np.float64_t[:, :] inc_field_l,
-        double[:, :] inc_field_u,
+        np.float64_t[:, :] inc_field_u,
         Py_ssize_t lookup_id,
         int sign_l,
         int sign_u,

gprMax/fields_outputs.py

@@ -17,6 +17,7 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 from string import Template
+from pathlib import Path
 
 import h5py
 
@@ -93,6 +94,53 @@ __global__ void store_outputs(int NRX, int iteration, const int* __restrict__ rx
 
 
 def write_hdf5_outputfile(outputfile, G):
+    write_hdf5_main_grid_outputfile(outputfile, G)
+    write_hdf5_sub_grid_outputfile(outputfile, G)
+
+
+def write_hdf5_main_grid_outputfile(outputfile, G):
+    """Write an output file in HDF5 format.
+
+    Args:
+        outputfile (str): Name of the output file.
+        Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
+        G (class): Grid class instance - holds essential parameters describing the model.
+    """
+
+    write_data(outputfile, G)
+
+
+def write_hdf5_sub_grid_outputfile(outputfile, G):
+    """Write an output file in HDF5 format.
+
+    Args:
+        outputfile (str): Name of the output file.
+        Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
+        G (class): Grid class instance - holds essential parameters describing the model.
+    """
+
+    stem = outputfile.stem
+    suffix = outputfile.suffix
+    parent = outputfile.parent
+
+    for sg in G.subgrids:
+
+        # create an outputfile for each subgrid
+        fp = stem + '_' + sg.name + suffix
+        fp = parent / Path(fp)
+
+        f = write_data(fp, sg)
+
+        # write some additional meta data about the subgrid
+        f.attrs['is_os_sep'] = sg.is_os_sep
+        f.attrs['pml_separation'] = sg.pml_separation
+        f.attrs['subgrid_pml_thickness'] = sg.pmlthickness['x0']
+        f.attrs['filter'] = sg.filter
+        f.attrs['ratio'] = sg.ratio
+        f.attrs['interpolation'] = sg.interpolation
+
+
+def write_data(outputfile, G):
     """Write an output file in HDF5 format.
 
     Args:
@@ -146,3 +194,5 @@ def write_hdf5_outputfile(outputfile, G):
 
         for output in rx.outputs:
             f['/rxs/rx' + str(rxindex + 1) + '/' + output] = rx.outputs[output]
+
+    return f

gprMax/input_cmds_geometry.py

@@ -246,13 +246,13 @@ def process_geometrycmds(geometry):
             ID = tmp[13]
             # without seed
             if len(tmp) == 14:
-                fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, ID=ID, n_materials=n_materials)
+                fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials)
             # with seed
             elif len(tmp) == 15:
-                fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, ID=ID, n_materials=n_materials, seed=tmp[14])
+                fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials, seed=tmp[14])
             # user specified averaging
             elif len(tmp) == 16:
-                fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, ID=ID, n_materials=n_materials, seed=tmp[14], averaging=tmp[15].lower())
+                fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials, seed=tmp[14], averaging=tmp[15].lower())
             else:
                 raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
 

gprMax/input_cmds_multiuse.py

@@ -143,11 +143,10 @@ def process_multicmds(multicmds):
             tmp = cmdinstance.split()
             if len(tmp) != 3 and len(tmp) < 5:
                 raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' has an incorrect number of parameters')
-
             if len(tmp) == 3:
                 rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])))
             else:
-                rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])), ID=tmp[3], outputs=tmp[4:])
+                rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])), id=tmp[3], outputs=tmp[4:])
             scene_objects.append(rx)
 
     # Receiver array
@@ -196,7 +195,7 @@ def process_multicmds(multicmds):
             if len(tmp) != 5:
                 raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly five parameters')
 
-            material = Material(er=float(tmp[0]), se=float(tmp[1]), mr=float(tmp[2]), sm=float(tmp[3]), ID=tmp[4])
+            material = Material(er=float(tmp[0]), se=float(tmp[1]), mr=float(tmp[2]), sm=float(tmp[3]), id=tmp[4])
             scene_objects.append(material)
 
     cmdname = '#add_dispersion_debye'
@@ -274,7 +273,7 @@ def process_multicmds(multicmds):
                                  sand_density=float(tmp[3]),
                                  water_fraction_lower=float(tmp[4]),
                                  water_fraction_upper=float(tmp[5]),
-                                 ID=tmp[6])
+                                 id=tmp[6])
             scene_objects.append(soil)
 
     # Geometry views (creates VTK-based geometry files)

gprMax/model_build_run.py

@@ -81,7 +81,7 @@ class ModelBuildRun:
         self.G = G
         self.sim_config = sim_config
         self.model_config = model_config
-        self.printer = Printer(sim_config)
+        self.printer = Printer(config)
         # Monitor memory usage
         self.p = None
 
@@ -130,9 +130,9 @@ class ModelBuildRun:
                 receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2]
 
         # Write files for any geometry views and geometry object outputs
-        if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.general['messages']:
+        if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.is_messages():
             print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL)
-        if config.general['messages']: print()
+        if config.is_messages(): print()
         for i, geometryview in enumerate(G.geometryviews):
             geometryview.set_filename(self.model_config.appendmodelnumber)
             pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry view file {}/{}, {}'.format(i + 1, len(G.geometryviews), os.path.split(geometryview.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
@@ -153,7 +153,7 @@ class ModelBuildRun:
 
         G = self.G
 
-        printer = Printer(sim_config)
+        printer = Printer(config)
         printer.print(model_config.next_model)
 
         scene = self.build_scene()
@@ -300,7 +300,7 @@ class ModelBuildRun:
         """
         G = self.G
 
-        if config.general['messages']:
+        if config.is_messages():
             iterator = tqdm(range(G.iterations), desc='Running simulation, model ' + str(self.model_config
                             .i + 1) + '/' + str(self.sim_config.model_end), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
         else:

gprMax/solvers.py

@@ -21,7 +21,7 @@ from gprMax.utilities import timer
 from .grid import FDTDGrid
 from .grid import GPUGrid
 import gprMax.config as config
-from .subgrids.updates import create_subgrid_updates
+from .subgrids.solver import create_solver as create_subgrid_solver
 
 
 def create_G(sim_config):
@@ -40,20 +40,20 @@ def create_solver(G, sim_config):
     """Returns the configured solver."""
     if sim_config.gpu:
         updates = GPUUpdates(G)
+        solver = Solver(updates)
     elif sim_config.subgrid:
-        updates = create_subgrid_updates(G)
+        solver = create_subgrid_solver(G)
     else:
         updates = CPUUpdates(G)
+        solver = Solver(updates)
 
-    solver = Solver(updates)
+    return solver
 
     # a large range of function exist to advance the time step for dispersive
     # materials. The correct function is set here  based on the
     # the required numerical precision and dispersive material type.
-    props = updates.adapt_dispersive_config(config)
+    #props = updates.adapt_dispersive_config(config)
-    updates.set_dispersive_updates(props)
+    #updates.set_dispersive_updates(props)
-
-    return solver
 
 
 class Solver:
@@ -76,7 +76,6 @@ class Solver:
     def solve(self, iterator):
         """Time step the FDTD model."""
         tsolvestart = timer()
-
         for iteration in iterator:
             self.updates.store_outputs(iteration)
             self.updates.store_snapshots(iteration)

gprMax/sources.py

@@ -169,7 +169,6 @@ class HertzianDipole(Source):
             j = self.ycoord
             k = self.zcoord
             componentID = 'E' + self.polarisation
-
             if self.polarisation == 'x':
                 Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz))
 

gprMax/subgrids/base.py

@@ -1,24 +1,18 @@
-# Copyright (C) 2015-2019: The University of Edinburgh
-#                 Authors: Craig Warren and Antonis Giannopoulos
-#
-# This file is part of gprMax.
-#
-# gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# gprMax is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
-
 from ..grid import FDTDGrid
+
+from ..materials import Material
+
+from scipy.constants import mu_0
+from scipy.constants import epsilon_0
+from scipy.constants import c
+
 import numpy as np
 
+from colorama import init
+from colorama import Fore
+from colorama import Style
+init()
+
 
 class SubGridBase(FDTDGrid):
 
@@ -26,8 +20,6 @@ class SubGridBase(FDTDGrid):
         super().__init__()
 
         self.mode = '3D'
-
-        # subgridding ratio
         self.ratio = kwargs['ratio']
 
         if self.ratio % 2 == 0:
@@ -60,17 +52,21 @@ class SubGridBase(FDTDGrid):
         self.n_boundary_cells_y = d_to_pml + self.pmlthickness['y0']
         self.n_boundary_cells_z = d_to_pml + self.pmlthickness['z0']
 
-        # vectorise coordinates
-        #self.p0 = np.array(self.i0, self.j0, self.k0)
-        self.n_boundary_cells_p = np.array([self.n_boundary_cells_x, self.n_boundary_cells_y, self.n_boundary_cells_z])
-
-        # interpolation scheme
         self.interpolation = kwargs['interpolation']
 
-    def main_grid_index_to_subgrid_index(self, p):
+    def calculate_dt(self):
-        """
+        self.dt = (1 / (c * np.sqrt(
-        Return the equivalent spatial index of the global (i, j, k) in the subgrid.
+                    (1 / self.dx) * (1 / self.dx) +
-        Args:
+                    (1 / self.dy) * (1 / self.dy) +
-              p (numpy.array): i, j, k indices of a point in the main grid.
+                    (1 / self.dz) * (1 / self.dz))))
-        """
+
-        return self.n_boundary_cells_p + (p - self.p0) * self.ratio
+
+    def main_grid_index_to_subgrid_index(self, i, j, k):
+        i_s = self.n_boundary_cells_x + (i - self.i0) * self.ratio
+        j_s = self.n_boundary_cells_y + (j - self.j0) * self.ratio
+        k_s = self.n_boundary_cells_z + (k - self.k0) * self.ratio
+
+        return (i_s, j_s, k_s)
+
+    def initialise_geometry_arrays(self):
+        super().initialise_geometry_arrays()

gprMax/subgrids/precursor_nodes.py

@@ -1,23 +1,6 @@
-# Copyright (C) 2015-2019: The University of Edinburgh
-#                 Authors: Craig Warren and Antonis Giannopoulos
-#
-# This file is part of gprMax.
-#
-# gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# gprMax is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
-
 import numpy as np
 from scipy import interpolate
+import sys
 
 
 def calculate_weighting_coefficients(x1, x):
@@ -26,14 +9,13 @@ def calculate_weighting_coefficients(x1, x):
     return (c1, c2)
 
 
-class PrecusorNodesBase:
+class PrecusorNodes2dBase(object):
 
     def __init__(self, fdtd_grid, sub_grid):
         self.G = fdtd_grid
         self.ratio = sub_grid.ratio
         self.nwx = sub_grid.nwx
         self.nwy = sub_grid.nwy
-        self.nwz = sub_grid.nwz
         self.sub_grid = sub_grid
         self.interpolation = sub_grid.interpolation
 
@@ -60,6 +42,300 @@ class PrecusorNodesBase:
 
         self.l_weight = self.ratio // 2
         self.r_weight = self.ratio - self.l_weight
+        #self.l_weight = 1
+        #self.r_weight = 2
+
+    def _initialize_fields(self):
+        print('dont call me')
+        sys.exit()
+        pass
+
+    def _initialize_field_names(self):
+        print('dont call me')
+        sys.exit()
+        pass
+
+    def interpolate_magnetic_in_time(self, m):
+        self.weight_pre_and_current_fields(m, self.fn_m)
+
+    def interpolate_electric_in_time(self, m):
+        self.weight_pre_and_current_fields(m, self.fn_e)
+
+    def weight_pre_and_current_fields(self, m, field_names):
+        c1, c2 = calculate_weighting_coefficients(m, self.ratio)
+
+        for f in field_names:
+            try:
+                val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
+            except ValueError:
+                print(self.ex_front_0.shape)
+                print(self.ex_front_1.shape)
+                raise Exception(f)
+            setattr(self, f, val)
+
+    def calc_exact_field(self, field_names):
+        """
+            Function to set the fields used in update calculations to the
+            values at the current main time step.
+            i.e. ey_left = copy.ey_left_1
+        """
+        for f in field_names:
+            val = np.copy(getattr(self, f + '_1'))
+            setattr(self, f, val)
+
+    def calc_exact_magnetic_in_time(self):
+        self.calc_exact_field(self.fn_m)
+
+    def calc_exact_electric_in_time(self):
+        self.calc_exact_field(self.fn_e)
+
+
+class PrecursorNodes2dTM(PrecusorNodes2dBase):
+
+    def __init__(self, fdtd_grid, sub_grid):
+        super().__init__(fdtd_grid, sub_grid)
+
+    def _initialize_fields(self):
+
+        # H precursors
+        self.hx_front_0 = np.zeros(self.nwx + 1)
+        self.hx_front_1 = np.zeros(self.nwx + 1)
+
+        self.hx_back_0 = np.zeros(self.nwx + 1)
+        self.hx_back_1 = np.zeros(self.nwx + 1)
+
+        self.hy_left_0 = np.zeros(self.nwy + 1)
+        self.hy_left_1 = np.zeros(self.nwy + 1)
+
+        self.hy_right_0 = np.zeros(self.nwy + 1)
+        self.hy_right_1 = np.zeros(self.nwy + 1)
+
+        # E precursors
+        self.ez_front_0 = np.zeros(self.nwx + 1)
+        self.ez_front_1 = np.zeros(self.nwx + 1)
+
+        self.ez_back_0 = np.zeros(self.nwx + 1)
+        self.ez_back_1 = np.zeros(self.nwx + 1)
+
+        self.ez_left_0 = np.zeros(self.nwy + 1)
+        self.ez_left_1 = np.zeros(self.nwy + 1)
+
+        self.ez_right_0 = np.zeros(self.nwy + 1)
+        self.ez_right_1 = np.zeros(self.nwy + 1)
+
+    def _initialize_field_names(self):
+        self.fn_m = ['hy_left', 'hy_right', 'hx_back', 'hx_front']
+        self.fn_e = ['ez_left', 'ez_right', 'ez_back', 'ez_front']
+
+    def interpolate_across_sub_cells(self, y, n0, n1):
+        n = n1 - n0
+        x = np.arange(0, y.shape[0], 1.0)
+        x_new = np.linspace(1, y.shape[0] - 2, n * self.ratio + 1)
+        f = interpolate.interp1d(x, y, kind='linear')
+        y_new = f(x_new)
+        """plt.plot(x, y, 'x', x_new, y_new, '.')
+        plt.show()
+        sys.exit()"""
+        return y_new
+
+    def update_electric(self):
+        # f = np.sin(np.arange(0, 13))
+
+        # line of Ez nodes along the left edge of the IS
+        self.ez_left_0 = np.copy(self.ez_left_1)
+        # interpolate nodes between two Ez nodes 1 cell beyond/infront the IS
+        f = self.Ez[self.i0, self.j0 - 1:self.j1 + 2, 0]
+        self.ez_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        self.ez_right_0 = np.copy(self.ez_right_1)
+        f = self.Ez[self.i1, self.j0 - 1:self.j1 + 2, 0]
+        self.ez_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        self.ez_front_0 = np.copy(self.ez_front_1)
+        f = self.Ez[self.i0 - 1:self.i1 + 2, self.j0, 0]
+        self.ez_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+        self.ez_back_0 = np.copy(self.ez_back_1)
+        f = self.Ez[self.i0 - 1:self.i1 + 2, self.j1, 0]
+        self.ez_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+    def update_hy_left_1(self):
+
+        self.hy_left_0 = np.copy(self.hy_left_1)
+        c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+        l = self.Hy[self.i0 - 1, self.j0 - 1:self.j1 + 2, 0]
+        r = self.Hy[self.i0, self.j0 - 1:self.j1 + 2, 0]
+
+        # Tranverse interpolated hz
+        hz_t_i = c1 * l + c2 * r
+        # Tangential interpolated hz
+        self.hy_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hy_right_1(self):
+
+        self.hy_right_0 = np.copy(self.hy_right_1)
+        c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+        l = self.Hy[self.i1 - 1, self.j0 - 1:self.j1 + 2, 0]
+        r = self.Hy[self.i1, self.j0 - 1:self.j1 + 2, 0]
+        hz_t_i = c1 * l + c2 * r
+        self.hy_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hx_back_1(self):
+        self.hx_back_0 = np.copy(self.hx_back_1)
+        c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+        l = self.Hx[self.i0 - 1:self.i1 + 2, self.j1 - 1, 0]
+        r = self.Hx[self.i0 - 1:self.i1 + 2, self.j1, 0]
+        hz_t_i = c1 * l + c2 * r
+        self.hx_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_hx_front_1(self):
+
+        self.hx_front_0 = np.copy(self.hx_front_1)
+        c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+        l = self.Hx[self.i0 - 1:self.i1 + 2, self.j0 - 1, 0]
+        r = self.Hx[self.i0 - 1:self.i1 + 2, self.j0, 0]
+        hz_t_i = c1 * l + c2 * r
+        self.hx_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_magnetic(self):
+        self.update_hy_left_1()
+        self.update_hy_right_1()
+        self.update_hx_back_1()
+        self.update_hx_front_1()
+
+
+class PrecursorNodes2d(PrecusorNodes2dBase):
+
+    def __init__(self, fdtd_grid, sub_grid):
+        super().__init__(fdtd_grid, sub_grid)
+
+    def _initialize_fields(self):
+        # E precursors
+        self.ex_front_0 = np.zeros(self.nwx)
+        self.ex_front_1 = np.zeros(self.nwx)
+
+        self.ex_back_0 = np.zeros(self.nwx)
+        self.ex_back_1 = np.zeros(self.nwx)
+
+        self.ey_left_0 = np.zeros(self.nwy)
+        self.ey_left_1 = np.zeros(self.nwy)
+
+        self.ey_right_0 = np.zeros(self.nwy)
+        self.ey_right_1 = np.zeros(self.nwy)
+
+        # H precursors
+        self.hz_front_0 = np.zeros(self.nwx)
+        self.hz_front_1 = np.zeros(self.nwx)
+
+        self.hz_back_0 = np.zeros(self.nwx)
+        self.hz_back_1 = np.zeros(self.nwx)
+
+        self.hz_left_0 = np.zeros(self.nwy)
+        self.hz_left_1 = np.zeros(self.nwy)
+
+        self.hz_right_0 = np.zeros(self.nwy)
+        self.hz_right_1 = np.zeros(self.nwy)
+
+    def _initialize_field_names(self):
+        # Field names
+        self.fn_m = ['hz_left', 'hz_right', 'hz_back', 'hz_front']
+        self.fn_e = ['ey_left', 'ey_right', 'ex_back', 'ex_front']
+
+    def update_hz_left_1(self):
+        self.hz_left_0 = np.copy(self.hz_left_1)
+        c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+        l = self.Hz[self.i0 - 1, self.j0 - 1:self.j1 + 1, 1]
+        r = self.Hz[self.i0, self.j0 - 1:self.j1 + 1, 1]
+
+        # Tranverse interpolated hz
+        hz_t_i = c1 * l + c2 * r
+        # Tangential interpolated hz
+        self.hz_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hz_right_1(self):
+
+        self.hz_right_0 = np.copy(self.hz_right_1)
+        c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+        l = self.Hz[self.i1 - 1, self.j0 - 1:self.j1 + 1, 1]
+        r = self.Hz[self.i1, self.j0 - 1:self.j1 + 1, 1]
+        hz_t_i = c1 * l + c2 * r
+        self.hz_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hz_back_1(self):
+        self.hz_back_0 = np.copy(self.hz_back_1)
+        c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+        l = self.Hz[self.i0 - 1:self.i1 + 1, self.j1 - 1, 1]
+        r = self.Hz[self.i0 - 1:self.i1 + 1, self.j1, 1]
+        hz_t_i = c1 * l + c2 * r
+        self.hz_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_hz_front_1(self):
+        self.hz_front_0 = np.copy(self.hz_front_1)
+        c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+        l = self.Hz[self.i0 - 1:self.i1 + 1, self.j0 - 1, 1]
+        r = self.Hz[self.i0 - 1:self.i1 + 1, self.j0, 1]
+        hz_t_i = c1 * l + c2 * r
+        self.hz_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_magnetic(self):
+        self.update_hz_left_1()
+        self.update_hz_right_1()
+        self.update_hz_back_1()
+        self.update_hz_front_1()
+
+    def update_electric(self):
+
+        """Function to calculate the precursor nodes at the next main
+            grid time-step
+        """
+
+        # LEFT BOUNDARY
+        # Save the last precursor node values - these will be used interpolate
+        # field values at sub grid time step values.
+        self.ey_left_0 = np.copy(self.ey_left_1)
+        # line of Ex nodes along the left edge of the IS
+        f = self.Ey[self.i0, self.j0 - 1:self.j1 + 1, 1]
+        self.ey_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        # RIGHT BOUNDARY
+        self.ey_right_0 = np.copy(self.ey_right_1)
+        f = self.Ey[self.i1, self.j0 - 1:self.j1 + 1, 1]
+        self.ey_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        # BACK BOUNDARY
+        self.ex_back_0 = np.copy(self.ex_back_1)
+        f = self.Ex[self.i0 - 1:self.i1 + 1, self.j1, 1]
+        self.ex_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+        # FRONT BOUNDARY
+        self.ex_front_0 = np.copy(self.ex_front_1)
+        f = self.Ex[self.i0 - 1:self.i1 + 1, self.j0, 1]
+        self.ex_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+    def interpolate_across_sub_cells(self, y, n0, n1):
+        n = n1 - n0
+        x = np.arange(0.5, y.shape[0], 1.0)
+        x_new = np.linspace(1 + self.d, (1 + n) - self.d, n * self.ratio)
+        f = interpolate.interp1d(x, y, kind='linear')
+        y_new = f(x_new)
+
+        return y_new
+
+
+class PrecursorNodes(PrecusorNodes2dBase):
+
+    def __init__(self, fdtd_grid, sub_grid):
+
+        self.nwz = sub_grid.nwz
+        super().__init__(fdtd_grid, sub_grid)
+        self.initialize_magnetic_slices_array()
+        self.initialize_electric_slices_array()
+
+    def build_coefficient_matrix(self, update_coefficients, lookup_id, inc_field, name):
+        """Function builds a 2d matrix of update coefficients for each face and field.
+            This allows the is nodes to be updated via slicing rather than iterating which is much faster
+        """
+        pass
 
     def _initialize_fields(self):
 
@@ -157,124 +433,12 @@ class PrecusorNodesBase:
             'ex_bottom', 'ey_bottom'
         ]
 
-    def interpolate_magnetic_in_time(self, m):
-        self.weight_pre_and_current_fields(m, self.fn_m)
-
-    def interpolate_electric_in_time(self, m):
-        self.weight_pre_and_current_fields(m, self.fn_e)
-
-    def weight_pre_and_current_fields(self, m, field_names):
-        c1, c2 = calculate_weighting_coefficients(m, self.ratio)
-
-        for f in field_names:
-            try:
-                val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
-            except ValueError:
-                print(self.ex_front_0.shape)
-                print(self.ex_front_1.shape)
-                raise Exception(f)
-            setattr(self, f, val)
-
-    def calc_exact_field(self, field_names):
-        """
-            Function to set the fields used in update calculations to the
-            values at the current main time step.
-            i.e. ey_left = copy.ey_left_1
-        """
-        for f in field_names:
-            val = np.copy(getattr(self, f + '_1'))
-            setattr(self, f, val)
-
     def update_previous_timestep_fields(self, field_names):
         for fn in field_names:
             val = getattr(self, fn + '_1')
             val_c = np.copy(val)
             setattr(self, fn + '_0', val_c)
 
-    def calc_exact_magnetic_in_time(self):
-        self.calc_exact_field(self.fn_m)
-
-    def calc_exact_electric_in_time(self):
-        self.calc_exact_field(self.fn_e)
-
-    def create_interpolated_coords(self, mid, field):
-
-        n_x = field.shape[0]
-        n_y = field.shape[1]
-
-        if mid:
-            x = np.arange(0.5, n_x, 1.0)
-            z = np.arange(0, n_y, 1.0)
-
-            # Coordinates that require interpolated values
-            x_sg = np.linspace(self.d, n_x - self.d, n_x * self.ratio)
-            z_sg = np.linspace(0, n_y - 1, (n_y - 1) * self.ratio + 1)
-
-        else:
-            x = np.arange(0, n_x, 1.0)
-            z = np.arange(0.5, n_y, 1.0)
-
-            # Coordinates that require interpolated values
-            x_sg = np.linspace(0, n_x - 1, (n_x - 1) * self.ratio + 1)
-            z_sg = np.linspace(self.d, n_y - self.d, n_y * self.ratio)
-
-        return (x, z, x_sg, z_sg)
-
-    def interpolate_to_sub_grid(self, field, coords):
-        x, z, x_sg, z_sg = coords
-        interp_f = interpolate.RectBivariateSpline(x,
-                                                   z,
-                                                   field,
-                                                   kx=self.interpolation,
-                                                   ky=self.interpolation)
-        f_i = interp_f(x_sg, z_sg)
-        return f_i
-
-    def update_magnetic(self):
-
-        # Copy previous time step magnetic field values to the previous
-        # time step variables
-        self.update_previous_timestep_fields(self.fn_m)
-
-        for obj in self.magnetic_slices:
-
-            # Grab the main grid fields used to interpolate across the IS
-            # f = self.Hi[slice]
-            f_1, f_2 = self.f1_f2_magnetic(obj)
-            w = obj[-2]
-            c1, c2 = calculate_weighting_coefficients(w, self.ratio)
-            # transverse interpolated h field
-            f_t = c1 * f_1 + c2 * f_2
-
-            # interpolate over a fine grid
-            f_i = self.interpolate_to_sub_grid(f_t, obj[1])
-
-            if f_i == f_t:
-                raise ValueError
-            # discard the outer nodes only required for interpolation
-            f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
-            setattr(self, obj[0], f)
-
-    def update_electric(self):
-
-        self.update_previous_timestep_fields(self.fn_e)
-
-        for obj in self.electric_slices:
-            f_m = self.f_m_electric(obj)
-            f_i = self.interpolate_to_sub_grid(f_m, obj[1])
-            # discard the outer nodes only required for interpolation
-            f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
-            setattr(self, obj[0], f)
-
-
-class PrecursorNodesFiltered(PrecusorNodesBase):
-
-    def __init__(self, fdtd_grid, sub_grid):
-
-        super().__init__(fdtd_grid, sub_grid)
-        self.initialize_magnetic_slices_array()
-        self.initialize_electric_slices_array()
-
     def initialize_magnetic_slices_array(self):
 
         # Array contains the indices at which the main grid should be sliced
@@ -283,99 +447,59 @@ class PrecursorNodesFiltered(PrecusorNodesBase):
 
         # Extend the surface so that the outer fields can be interpolated
         # more accurately
-        i0 = self.i0 - 1
+        #i0 = self.i0 - 1
-        j0 = self.j0 - 1
+        #j0 = self.j0 - 1
-        k0 = self.k0 - 1
+        #k0 = self.k0 - 1
-        i1 = self.i1 + 1
+        #i1 = self.i1 + 1
-        j1 = self.j1 + 1
+        #j1 = self.j1 + 1
-        k1 = self.k1 + 1
+        #k1 = self.k1 + 1
+
+        # not extended
+        i0 = self.i0
+        j0 = self.j0
+        k0 = self.k0
+        i1 = self.i1
+        j1 = self.j1
+        k1 = self.k1
 
         slices = [
             ['hy_left_1', False,
-                (self.i0 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
                 (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
-                (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                self.l_weight,
-                self.Hy],
             ['hy_right_1', False,
-                (self.i1 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
                 (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
-                (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                self.r_weight,
-                self.Hy],
             ['hz_left_1', True,
-                (self.i0 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
                 (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
-                (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                self.l_weight,
-                self.Hz],
             ['hz_right_1', True,
-                (self.i1 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
                 (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
-                (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                self.r_weight,
-                self.Hz],
             ['hx_front_1', False,
-                (slice(i0, i1 + 1, 1), self.j0 - 2, slice(k0, k1, 1)),
                 (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)), self.Hx],
-                (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)),
-                self.l_weight,
-                self.Hx],
             ['hx_back_1', False,
-                (slice(i0, i1 + 1, 1), self.j1 - 2, slice(k0, k1, 1)),
                 (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)), self.Hx],
-                (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)),
-                self.r_weight,
-                self.Hx],
             ['hz_front_1', True,
-                (slice(i0, i1, 1), self.j0 - 2, slice(k0, k1 + 1, 1)),
                 (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
-                (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)), self.Hz],
-                (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)),
-                self.l_weight,
-                self.Hz],
             ['hz_back_1', True,
-                (slice(i0, i1, 1), self.j1 - 2, slice(k0, k1 + 1, 1)),
                 (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
-                (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)), self.Hz],
-                (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)),
-                self.r_weight,
-                self.Hz],
             ['hx_bottom_1', False,
                 # check these indexes
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 2),
                 (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0), self.Hx],
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1),
-                self.l_weight,
-                self.Hx],
             ['hx_top_1', False,
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 2),
                 (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1), self.Hx],
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1),
-                self.r_weight,
-                self.Hx],
             ['hy_bottom_1', True,
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 2),
                 (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0), self.Hy],
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1),
-                self.l_weight,
-                self.Hy],
             ['hy_top_1', True,
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 2),
                 (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1), self.Hy]
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1),
-                self.r_weight,
-                self.Hy]
         ]
 
         for obj in slices:
@@ -387,208 +511,20 @@ class PrecursorNodesFiltered(PrecusorNodesBase):
     def initialize_electric_slices_array(self):
         # Extend the region sliced from the main grid by 1 cell.
         # this allows more accurate interpolation for the outernodes
-        i0 = self.i0 - 1
+        #i0 = self.i0 - 1
-        j0 = self.j0 - 1
+        #j0 = self.j0 - 1
-        k0 = self.k0 - 1
+        #k0 = self.k0 - 1
-        i1 = self.i1 + 1
+        #i1 = self.i1 + 1
-        j1 = self.j1 + 1
+        #j1 = self.j1 + 1
-        k1 = self.k1 + 1
+        #k1 = self.k1 + 1
 
-        # Spatially interpolate nodes
+        # not extended
-        slices = [
+        i0 = self.i0
-            ['ex_front_1', True,
+        j0 = self.j0
-                (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
+        k0 = self.k0
-                (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
+        i1 = self.i1
-                (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)),
+        j1 = self.j1
-                self.Ex],
+        k1 = self.k1
-            ['ex_back_1', True,
-                (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
-                (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
-                (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)),
-                self.Ex],
-            ['ez_front_1', False,
-                (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)),
-                self.Ez],
-            ['ez_back_1', False,
-                (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)),
-                self.Ez],
-            ['ey_left_1', True,
-                (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                self.Ey],
-            ['ey_right_1', True,
-                (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                self.Ey],
-            ['ez_left_1', False,
-                (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                self.Ez],
-            ['ez_right_1', False,
-                (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                self.Ez],
-
-            ['ex_bottom_1', True,
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1),
-                self.Ex],
-            ['ex_top_1', True,
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1),
-                self.Ex],
-            ['ey_bottom_1', False,
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1),
-                self.Ey],
-            ['ey_top_1', False,
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1),
-                self.Ey]
-        ]
-
-        for obj in slices:
-            sliced_field = obj[-1][obj[2]]
-            obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
-
-        self.electric_slices = slices
-
-    def f1_f2_magnetic(self, obj):
-        f_u_1 = obj[-1][obj[2]]
-        f_u_2 = obj[-1][obj[3]]
-        f_u_3 = obj[-1][obj[4]]
-        f_u_4 = obj[-1][obj[5]]
-
-        f_1 = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
-        f_2 = 0.25 * f_u_2 + 0.5 * f_u_3 + 0.25 * f_u_4
-
-        return f_1, f_2
-
-    def f_m_electric(self, obj):
-
-        f_u_1 = obj[-1][obj[2]]
-        f_u_2 = obj[-1][obj[3]]
-        f_u_3 = obj[-1][obj[4]]
-
-        f_m = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
-
-        return f_m
-
-
-class PrecursorNodes(PrecusorNodesBase):
-
-    def __init__(self, fdtd_grid, sub_grid):
-
-        super().__init__(fdtd_grid, sub_grid)
-        self.initialize_magnetic_slices_array()
-        self.initialize_electric_slices_array()
-
-    def initialize_magnetic_slices_array(self):
-
-        # Array contains the indices at which the main grid should be sliced
-        # to obtain the 2 2d array of H nodes required for interpolation
-        # across the IS boundary for each h field on each face of the subgrid
-
-        # Extend the surface so that the outer fields can be interpolated
-        # more accurately
-        i0 = self.i0 - 1
-        j0 = self.j0 - 1
-        k0 = self.k0 - 1
-        i1 = self.i1 + 1
-        j1 = self.j1 + 1
-        k1 = self.k1 + 1
-
-        slices = [
-            ['hy_left_1', False,
-                (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                self.l_weight,
-                self.Hy],
-            ['hy_right_1', False,
-                (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
-                self.r_weight,
-                self.Hy],
-            ['hz_left_1', True,
-                (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                self.l_weight,
-                self.Hz],
-            ['hz_right_1', True,
-                (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
-                self.r_weight,
-                self.Hz],
-            ['hx_front_1', False,
-                (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
-                self.l_weight,
-                self.Hx],
-            ['hx_back_1', False,
-                (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
-                (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
-                self.r_weight,
-                self.Hx],
-            ['hz_front_1', True,
-                (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
-                (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
-                self.l_weight,
-                self.Hz],
-            ['hz_back_1', True,
-                (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
-                (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
-                self.r_weight,
-                self.Hz],
-            ['hx_bottom_1', False,
-                # check these indexes
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
-                self.l_weight,
-                self.Hx],
-            ['hx_top_1', False,
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
-                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
-                self.r_weight,
-                self.Hx],
-            ['hy_bottom_1', True,
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
-                self.l_weight,
-                self.Hy],
-            ['hy_top_1', True,
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
-                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
-                self.r_weight,
-                self.Hy]
-        ]
-
-        for obj in slices:
-            sliced_field = obj[-1][obj[2]]
-            obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
-
-        self.magnetic_slices = slices
-
-    def initialize_electric_slices_array(self):
-        # Extend the region sliced from the main grid by 1 cell.
-        # this allows more accurate interpolation for the outernodes
-        i0 = self.i0 - 1
-        j0 = self.j0 - 1
-        k0 = self.k0 - 1
-        i1 = self.i1 + 1
-        j1 = self.j1 + 1
-        k1 = self.k1 + 1
 
         # Spatially interpolate nodes
         slices = [
@@ -614,10 +550,84 @@ class PrecursorNodes(PrecusorNodesBase):
 
         self.electric_slices = slices
 
-    def f1_f2_magnetic(self, obj):
+    def create_interpolated_coords(self, mid, field):
-        f_1 = obj[-1][obj[2]]
-        f_2 = obj[-1][obj[3]]
-        return f_1, f_2
 
-    def f_m_electric(self, obj):
+        n_x = field.shape[0]
-        return obj[-1][obj[2]]
+        n_y = field.shape[1]
+
+        if mid:
+            x = np.arange(0.5, n_x, 1.0)
+            z = np.arange(0, n_y, 1.0)
+
+            # Coordinates that require interpolated values
+            x_sg = np.linspace(self.d, n_x - self.d, n_x * self.ratio)
+            z_sg = np.linspace(0, n_y - 1, (n_y - 1) * self.ratio + 1)
+
+        else:
+            x = np.arange(0, n_x, 1.0)
+            z = np.arange(0.5, n_y, 1.0)
+
+            # Coordinates that require interpolated values
+            x_sg = np.linspace(0, n_x - 1, (n_x - 1) * self.ratio + 1)
+            z_sg = np.linspace(self.d, n_y - self.d, n_y * self.ratio)
+
+        return (x, z, x_sg, z_sg)
+
+    def update_magnetic(self):
+
+        # Copy previous time step magnetic field values to the previous
+        # time step variables
+        self.update_previous_timestep_fields(self.fn_m)
+
+        for obj in self.magnetic_slices:
+
+            # Grab the main grid fields used to interpolate across the IS
+            # f = self.Hi[slice]
+            f_1 = obj[-1][obj[2]]
+            f_2 = obj[-1][obj[3]]
+            if ('left' in obj[0] or
+                'bottom' in obj[0] or
+                    'front' in obj[0]):
+                w = self.l_weight
+            else:
+                w = self.r_weight
+            c1, c2 = calculate_weighting_coefficients(w, self.ratio)
+            # transverse interpolated h field
+            f_t = c1 * f_1 + c2 * f_2
+
+            # interpolate over a fine grid
+            f_i = self.interpolate_to_sub_grid(f_t, obj[1])
+
+            if f_i == f_t:
+                raise ValueError
+
+            # discard the outer nodes only required for interpolation
+            #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+            f = f_i
+            setattr(self, obj[0], f)
+
+    def update_electric(self):
+
+        self.update_previous_timestep_fields(self.fn_e)
+
+        for obj in self.electric_slices:
+            f_m = obj[-1][obj[2]]
+            f_i = self.interpolate_to_sub_grid(f_m, obj[1])
+            f = f_i
+            #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+            setattr(self, obj[0], f)
+
+    def interpolate_to_sub_grid(self, field, coords):
+        x, z, x_sg, z_sg = coords
+        ex_t = np.transpose(field)
+        f = interpolate.interp2d(x, z, ex_t, kind=self.interpolation)
+        #f = interpolate.RectBivariateSpline(x, z, field)
+        ex_inter_t = f(x_sg, z_sg)
+        ex_inter = np.transpose(ex_inter_t)
+        #ex_inter = ex_inter_t
+
+        return ex_inter
+
+
+class PlaneError(Exception):
+    pass

gprMax/subgrids/precursor_nodes_filtered.py

@@ -0,0 +1,738 @@
+import numpy as np
+from scipy import interpolate
+import sys
+
+
+def calculate_weighting_coefficients(x1, x):
+    c1 = (x - x1) / x
+    c2 = x1 / x
+    return (c1, c2)
+
+
+class PrecusorNodes2dBase(object):
+
+    def __init__(self, fdtd_grid, sub_grid):
+        self.G = fdtd_grid
+        self.ratio = sub_grid.ratio
+        self.nwx = sub_grid.nwx
+        self.nwy = sub_grid.nwy
+        self.sub_grid = sub_grid
+        self.interpolation = sub_grid.interpolation
+
+        self.Hx = fdtd_grid.Hx
+        self.Hy = fdtd_grid.Hy
+        self.Hz = fdtd_grid.Hz
+        self.Ex = fdtd_grid.Ex
+        self.Ey = fdtd_grid.Ey
+        self.Ez = fdtd_grid.Ez
+
+        # Main grid indices of subgrids
+        self.i0 = sub_grid.i0
+        self.j0 = sub_grid.j0
+        self.k0 = sub_grid.k0
+        self.i1 = sub_grid.i1
+        self.j1 = sub_grid.j1
+        self.k1 = sub_grid.k1
+
+        # dl / 2 sub cell
+        self.d = 1 / (2 * self.ratio)
+
+        self._initialize_fields()
+        self._initialize_field_names()
+
+        self.l_weight = self.ratio // 2
+        self.r_weight = self.ratio - self.l_weight
+        #self.l_weight = 1
+        #self.r_weight = 2
+
+    def _initialize_fields(self):
+        print('dont call me')
+        sys.exit()
+        pass
+
+    def _initialize_field_names(self):
+        print('dont call me')
+        sys.exit()
+        pass
+
+    def interpolate_magnetic_in_time(self, m):
+        self.weight_pre_and_current_fields(m, self.fn_m)
+
+    def interpolate_electric_in_time(self, m):
+        self.weight_pre_and_current_fields(m, self.fn_e)
+
+    def weight_pre_and_current_fields(self, m, field_names):
+        c1, c2 = calculate_weighting_coefficients(m, self.ratio)
+
+        for f in field_names:
+            try:
+                val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
+            except ValueError:
+                print(self.ex_front_0.shape)
+                print(self.ex_front_1.shape)
+                raise Exception(f)
+            setattr(self, f, val)
+
+    def calc_exact_field(self, field_names):
+        """
+            Function to set the fields used in update calculations to the
+            values at the current main time step.
+            i.e. ey_left = copy.ey_left_1
+        """
+        for f in field_names:
+            val = np.copy(getattr(self, f + '_1'))
+            setattr(self, f, val)
+
+    def calc_exact_magnetic_in_time(self):
+        self.calc_exact_field(self.fn_m)
+
+    def calc_exact_electric_in_time(self):
+        self.calc_exact_field(self.fn_e)
+
+
+class PrecursorNodes2dTM(PrecusorNodes2dBase):
+
+    def __init__(self, fdtd_grid, sub_grid):
+        super().__init__(fdtd_grid, sub_grid)
+
+    def _initialize_fields(self):
+
+        # H precursors
+        self.hx_front_0 = np.zeros(self.nwx + 1)
+        self.hx_front_1 = np.zeros(self.nwx + 1)
+
+        self.hx_back_0 = np.zeros(self.nwx + 1)
+        self.hx_back_1 = np.zeros(self.nwx + 1)
+
+        self.hy_left_0 = np.zeros(self.nwy + 1)
+        self.hy_left_1 = np.zeros(self.nwy + 1)
+
+        self.hy_right_0 = np.zeros(self.nwy + 1)
+        self.hy_right_1 = np.zeros(self.nwy + 1)
+
+        # E precursors
+        self.ez_front_0 = np.zeros(self.nwx + 1)
+        self.ez_front_1 = np.zeros(self.nwx + 1)
+
+        self.ez_back_0 = np.zeros(self.nwx + 1)
+        self.ez_back_1 = np.zeros(self.nwx + 1)
+
+        self.ez_left_0 = np.zeros(self.nwy + 1)
+        self.ez_left_1 = np.zeros(self.nwy + 1)
+
+        self.ez_right_0 = np.zeros(self.nwy + 1)
+        self.ez_right_1 = np.zeros(self.nwy + 1)
+
+    def _initialize_field_names(self):
+        self.fn_m = ['hy_left', 'hy_right', 'hx_back', 'hx_front']
+        self.fn_e = ['ez_left', 'ez_right', 'ez_back', 'ez_front']
+
+    def interpolate_across_sub_cells(self, y, n0, n1):
+        n = n1 - n0
+        x = np.arange(0, y.shape[0], 1.0)
+        x_new = np.linspace(1, y.shape[0] - 2, n * self.ratio + 1)
+        f = interpolate.interp1d(x, y, kind='linear')
+        y_new = f(x_new)
+        """plt.plot(x, y, 'x', x_new, y_new, '.')
+        plt.show()
+        sys.exit()"""
+        return y_new
+
+    def update_electric(self):
+        # f = np.sin(np.arange(0, 13))
+
+        # line of Ez nodes along the left edge of the IS
+        self.ez_left_0 = np.copy(self.ez_left_1)
+        # interpolate nodes between two Ez nodes 1 cell beyond/infront the IS
+        f = self.Ez[self.i0, self.j0 - 1:self.j1 + 2, 0]
+        self.ez_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        self.ez_right_0 = np.copy(self.ez_right_1)
+        f = self.Ez[self.i1, self.j0 - 1:self.j1 + 2, 0]
+        self.ez_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        self.ez_front_0 = np.copy(self.ez_front_1)
+        f = self.Ez[self.i0 - 1:self.i1 + 2, self.j0, 0]
+        self.ez_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+        self.ez_back_0 = np.copy(self.ez_back_1)
+        f = self.Ez[self.i0 - 1:self.i1 + 2, self.j1, 0]
+        self.ez_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+    def update_hy_left_1(self):
+
+        self.hy_left_0 = np.copy(self.hy_left_1)
+        c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+        l = self.Hy[self.i0 - 1, self.j0 - 1:self.j1 + 2, 0]
+        r = self.Hy[self.i0, self.j0 - 1:self.j1 + 2, 0]
+
+        # Tranverse interpolated hz
+        hz_t_i = c1 * l + c2 * r
+        # Tangential interpolated hz
+        self.hy_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hy_right_1(self):
+
+        self.hy_right_0 = np.copy(self.hy_right_1)
+        c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+        l = self.Hy[self.i1 - 1, self.j0 - 1:self.j1 + 2, 0]
+        r = self.Hy[self.i1, self.j0 - 1:self.j1 + 2, 0]
+        hz_t_i = c1 * l + c2 * r
+        self.hy_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hx_back_1(self):
+        self.hx_back_0 = np.copy(self.hx_back_1)
+        c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+        l = self.Hx[self.i0 - 1:self.i1 + 2, self.j1 - 1, 0]
+        r = self.Hx[self.i0 - 1:self.i1 + 2, self.j1, 0]
+        hz_t_i = c1 * l + c2 * r
+        self.hx_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_hx_front_1(self):
+
+        self.hx_front_0 = np.copy(self.hx_front_1)
+        c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+        l = self.Hx[self.i0 - 1:self.i1 + 2, self.j0 - 1, 0]
+        r = self.Hx[self.i0 - 1:self.i1 + 2, self.j0, 0]
+        hz_t_i = c1 * l + c2 * r
+        self.hx_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_magnetic(self):
+        self.update_hy_left_1()
+        self.update_hy_right_1()
+        self.update_hx_back_1()
+        self.update_hx_front_1()
+
+
+class PrecursorNodes2d(PrecusorNodes2dBase):
+
+    def __init__(self, fdtd_grid, sub_grid):
+        super().__init__(fdtd_grid, sub_grid)
+
+    def _initialize_fields(self):
+        # E precursors
+        self.ex_front_0 = np.zeros(self.nwx)
+        self.ex_front_1 = np.zeros(self.nwx)
+
+        self.ex_back_0 = np.zeros(self.nwx)
+        self.ex_back_1 = np.zeros(self.nwx)
+
+        self.ey_left_0 = np.zeros(self.nwy)
+        self.ey_left_1 = np.zeros(self.nwy)
+
+        self.ey_right_0 = np.zeros(self.nwy)
+        self.ey_right_1 = np.zeros(self.nwy)
+
+        # H precursors
+        self.hz_front_0 = np.zeros(self.nwx)
+        self.hz_front_1 = np.zeros(self.nwx)
+
+        self.hz_back_0 = np.zeros(self.nwx)
+        self.hz_back_1 = np.zeros(self.nwx)
+
+        self.hz_left_0 = np.zeros(self.nwy)
+        self.hz_left_1 = np.zeros(self.nwy)
+
+        self.hz_right_0 = np.zeros(self.nwy)
+        self.hz_right_1 = np.zeros(self.nwy)
+
+    def _initialize_field_names(self):
+        # Field names
+        self.fn_m = ['hz_left', 'hz_right', 'hz_back', 'hz_front']
+        self.fn_e = ['ey_left', 'ey_right', 'ex_back', 'ex_front']
+
+    def update_hz_left_1(self):
+        self.hz_left_0 = np.copy(self.hz_left_1)
+        c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+        l = self.Hz[self.i0 - 1, self.j0 - 1:self.j1 + 1, 1]
+        r = self.Hz[self.i0, self.j0 - 1:self.j1 + 1, 1]
+
+        # Tranverse interpolated hz
+        hz_t_i = c1 * l + c2 * r
+        # Tangential interpolated hz
+        self.hz_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hz_right_1(self):
+
+        self.hz_right_0 = np.copy(self.hz_right_1)
+        c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+        l = self.Hz[self.i1 - 1, self.j0 - 1:self.j1 + 1, 1]
+        r = self.Hz[self.i1, self.j0 - 1:self.j1 + 1, 1]
+        hz_t_i = c1 * l + c2 * r
+        self.hz_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+    def update_hz_back_1(self):
+        self.hz_back_0 = np.copy(self.hz_back_1)
+        c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+        l = self.Hz[self.i0 - 1:self.i1 + 1, self.j1 - 1, 1]
+        r = self.Hz[self.i0 - 1:self.i1 + 1, self.j1, 1]
+        hz_t_i = c1 * l + c2 * r
+        self.hz_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_hz_front_1(self):
+        self.hz_front_0 = np.copy(self.hz_front_1)
+        c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+        l = self.Hz[self.i0 - 1:self.i1 + 1, self.j0 - 1, 1]
+        r = self.Hz[self.i0 - 1:self.i1 + 1, self.j0, 1]
+        hz_t_i = c1 * l + c2 * r
+        self.hz_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+    def update_magnetic(self):
+        self.update_hz_left_1()
+        self.update_hz_right_1()
+        self.update_hz_back_1()
+        self.update_hz_front_1()
+
+    def update_electric(self):
+
+        """Function to calculate the precursor nodes at the next main
+            grid time-step
+        """
+
+        # LEFT BOUNDARY
+        # Save the last precursor node values - these will be used interpolate
+        # field values at sub grid time step values.
+        self.ey_left_0 = np.copy(self.ey_left_1)
+        # line of Ex nodes along the left edge of the IS
+        f = self.Ey[self.i0, self.j0 - 1:self.j1 + 1, 1]
+        self.ey_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        # RIGHT BOUNDARY
+        self.ey_right_0 = np.copy(self.ey_right_1)
+        f = self.Ey[self.i1, self.j0 - 1:self.j1 + 1, 1]
+        self.ey_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+        # BACK BOUNDARY
+        self.ex_back_0 = np.copy(self.ex_back_1)
+        f = self.Ex[self.i0 - 1:self.i1 + 1, self.j1, 1]
+        self.ex_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+        # FRONT BOUNDARY
+        self.ex_front_0 = np.copy(self.ex_front_1)
+        f = self.Ex[self.i0 - 1:self.i1 + 1, self.j0, 1]
+        self.ex_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+    def interpolate_across_sub_cells(self, y, n0, n1):
+        n = n1 - n0
+        x = np.arange(0.5, y.shape[0], 1.0)
+        x_new = np.linspace(1 + self.d, (1 + n) - self.d, n * self.ratio)
+        f = interpolate.interp1d(x, y, kind='linear')
+        y_new = f(x_new)
+
+        return y_new
+
+
+class PrecursorNodes(PrecusorNodes2dBase):
+
+    def __init__(self, fdtd_grid, sub_grid):
+
+        self.nwz = sub_grid.nwz
+        super().__init__(fdtd_grid, sub_grid)
+        self.initialize_magnetic_slices_array()
+        self.initialize_electric_slices_array()
+
+    def build_coefficient_matrix(self, update_coefficients, lookup_id, inc_field, name):
+        """Function builds a 2d matrix of update coefficients for each face and field.
+            This allows the is nodes to be updated via slicing rather than iterating which is much faster
+        """
+        pass
+
+    def _initialize_fields(self):
+
+        # Initialise the precursor arrays
+
+        # The precursors are divided up into the 6. Each represent 1
+        # face of a huygens cube surface. We represent each face as a 2d array
+        # containing a field in a particular direction.
+
+        # _1 are the fields at the current main grid timestep
+        # _0 are the fields at the previous main grid timestep
+        # We store both fields so we can do different interpolations between
+        # them on the fly.
+
+        # Front face
+        self.ex_front_1 = np.zeros((self.nwx, self.nwz + 1))
+        self.ex_front_0 = np.copy(self.ex_front_1)
+        self.ez_front_1 = np.zeros((self.nwx + 1, self.nwz))
+        self.ez_front_0 = np.copy(self.ez_front_1)
+
+        # The same as the opposite face
+        self.ex_back_1 = np.copy(self.ex_front_1)
+        self.ex_back_0 = np.copy(self.ex_front_1)
+        self.ez_back_1 = np.copy(self.ez_front_1)
+        self.ez_back_0 = np.copy(self.ez_front_1)
+
+        self.ey_left_1 = np.zeros((self.nwy, self.nwz + 1))
+        self.ey_left_0 = np.copy(self.ey_left_1)
+        self.ez_left_1 = np.zeros((self.nwy + 1, self.nwz))
+        self.ez_left_0 = np.copy(self.ez_left_1)
+
+        self.ey_right_1 = np.copy(self.ey_left_1)
+        self.ey_right_0 = np.copy(self.ey_left_1)
+        self.ez_right_1 = np.copy(self.ez_left_1)
+        self.ez_right_0 = np.copy(self.ez_left_1)
+
+        self.ex_bottom_1 = np.zeros((self.nwx, self.nwy + 1))
+        self.ex_bottom_0 = np.copy(self.ex_bottom_1)
+        self.ey_bottom_1 = np.zeros((self.nwx + 1, self.nwy))
+        self.ey_bottom_0 = np.copy(self.ey_bottom_1)
+
+        self.ex_top_1 = np.copy(self.ex_bottom_1)
+        self.ex_top_0 = np.copy(self.ex_bottom_1)
+        self.ey_top_1 = np.copy(self.ey_bottom_1)
+        self.ey_top_0 = np.copy(self.ey_bottom_1)
+
+        # Initialize the H precursor fields
+        self.hx_front_1 = np.copy(self.ez_front_1)
+        self.hx_front_0 = np.copy(self.ez_front_1)
+        self.hz_front_1 = np.copy(self.ex_front_1)
+        self.hz_front_0 = np.copy(self.ex_front_1)
+
+        self.hx_back_1 = np.copy(self.hx_front_1)
+        self.hx_back_0 = np.copy(self.hx_front_1)
+        self.hz_back_1 = np.copy(self.hz_front_1)
+        self.hz_back_0 = np.copy(self.hz_front_1)
+
+        self.hy_left_1 = np.copy(self.ez_left_1)
+        self.hy_left_0 = np.copy(self.ez_left_1)
+        self.hz_left_1 = np.copy(self.ey_left_1)
+        self.hz_left_0 = np.copy(self.ey_left_1)
+
+        self.hy_right_1 = np.copy(self.hy_left_1)
+        self.hy_right_0 = np.copy(self.hy_left_1)
+        self.hz_right_1 = np.copy(self.hz_left_1)
+        self.hz_right_0 = np.copy(self.hz_left_1)
+
+        self.hx_top_1 = np.copy(self.ey_top_1)
+        self.hx_top_0 = np.copy(self.ey_top_1)
+        self.hy_top_1 = np.copy(self.ex_top_1)
+        self.hy_top_0 = np.copy(self.ex_top_1)
+
+        self.hx_bottom_1 = np.copy(self.hx_top_1)
+        self.hx_bottom_0 = np.copy(self.hx_top_1)
+        self.hy_bottom_1 = np.copy(self.hy_top_1)
+        self.hy_bottom_0 = np.copy(self.hy_top_1)
+
+    def _initialize_field_names(self):
+
+        self.fn_m = [
+            'hx_front', 'hz_front',
+            'hx_back', 'hz_back',
+            'hy_left', 'hz_left',
+            'hy_right', 'hz_right',
+            'hx_top', 'hy_top',
+            'hx_bottom', 'hy_bottom'
+        ]
+
+        self.fn_e = [
+            'ex_front', 'ez_front',
+            'ex_back', 'ez_back',
+            'ey_left', 'ez_left',
+            'ey_right', 'ez_right',
+            'ex_top', 'ey_top',
+            'ex_bottom', 'ey_bottom'
+        ]
+
+    def update_previous_timestep_fields(self, field_names):
+        for fn in field_names:
+            val = getattr(self, fn + '_1')
+            val_c = np.copy(val)
+            setattr(self, fn + '_0', val_c)
+
+    def initialize_magnetic_slices_array(self):
+
+        # Array contains the indices at which the main grid should be sliced
+        # to obtain the 2 2d array of H nodes required for interpolation
+        # across the IS boundary for each h field on each face of the subgrid
+
+        # Extend the surface so that the outer fields can be interpolated
+        # more accurately
+        #i0 = self.i0 - 1
+        #j0 = self.j0 - 1
+        #k0 = self.k0 - 1
+        #i1 = self.i1 + 1
+        #j1 = self.j1 + 1
+        #k1 = self.k1 + 1
+
+        # not extended
+        i0 = self.i0
+        j0 = self.j0
+        k0 = self.k0
+        i1 = self.i1
+        j1 = self.j1
+        k1 = self.k1
+
+        slices = [
+            ['hy_left_1', False,
+                (self.i0 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
+            ['hy_right_1', False,
+                (self.i1 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
+            ['hz_left_1', True,
+                (self.i0 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
+            ['hz_right_1', True,
+                (self.i1 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
+            ['hx_front_1', False,
+                (slice(i0, i1 + 1, 1), self.j0 - 2, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)), self.Hx],
+            ['hx_back_1', False,
+                (slice(i0, i1 + 1, 1), self.j1 - 2, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)), self.Hx],
+            ['hz_front_1', True,
+                (slice(i0, i1, 1), self.j0 - 2, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)), self.Hz],
+            ['hz_back_1', True,
+                (slice(i0, i1, 1), self.j1 - 2, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)), self.Hz],
+            ['hx_bottom_1', False,
+                # check these indexes
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 2),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1), self.Hx],
+            ['hx_top_1', False,
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 2),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1), self.Hx],
+            ['hy_bottom_1', True,
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 2),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1), self.Hy],
+            ['hy_top_1', True,
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 2),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1), self.Hy]
+        ]
+
+        for obj in slices:
+            sliced_field = obj[-1][obj[2]]
+            obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
+
+        self.magnetic_slices = slices
+
+    def initialize_electric_slices_array(self):
+        # Extend the region sliced from the main grid by 1 cell.
+        # this allows more accurate interpolation for the outernodes
+        #i0 = self.i0 - 1
+        #j0 = self.j0 - 1
+        #k0 = self.k0 - 1
+        #i1 = self.i1 + 1
+        #j1 = self.j1 + 1
+        #k1 = self.k1 + 1
+
+        # not extended
+        i0 = self.i0
+        j0 = self.j0
+        k0 = self.k0
+        i1 = self.i1
+        j1 = self.j1
+        k1 = self.k1
+
+        # Spatially interpolate nodes
+        slices = [
+            ['ex_front_1', True,
+                (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)),
+                self.Ex],
+            ['ex_back_1', True,
+                (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
+                (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)),
+                self.Ex],
+            ['ez_front_1', False,
+                (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)),
+                self.Ez],
+            ['ez_back_1', False,
+                (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
+                (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)),
+                self.Ez],
+            ['ey_left_1', True,
+                (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                self.Ey],
+            ['ey_right_1', True,
+                (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+                self.Ey],
+            ['ez_left_1', False,
+                (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                self.Ez],
+            ['ez_right_1', False,
+                (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+                self.Ez],
+
+            ['ex_bottom_1', True,
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1),
+                self.Ex],
+            ['ex_top_1', True,
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
+                (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1),
+                self.Ex],
+            ['ey_bottom_1', False,
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1),
+                self.Ey],
+            ['ey_top_1', False,
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
+                (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1),
+                self.Ey]
+        ]
+
+        for obj in slices:
+            sliced_field = obj[-1][obj[2]]
+            obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
+
+        self.electric_slices = slices
+
+    def create_interpolated_coords(self, mid, field):
+
+        n_x = field.shape[0]
+        n_y = field.shape[1]
+
+        if mid:
+            x = np.arange(0.5, n_x, 1.0)
+            z = np.arange(0, n_y, 1.0)
+
+            # Coordinates that require interpolated values
+            x_sg = np.linspace(self.d, n_x - self.d, n_x * self.ratio)
+            z_sg = np.linspace(0, n_y - 1, (n_y - 1) * self.ratio + 1)
+
+        else:
+            x = np.arange(0, n_x, 1.0)
+            z = np.arange(0.5, n_y, 1.0)
+
+            # Coordinates that require interpolated values
+            x_sg = np.linspace(0, n_x - 1, (n_x - 1) * self.ratio + 1)
+            z_sg = np.linspace(self.d, n_y - self.d, n_y * self.ratio)
+
+        return (x, z, x_sg, z_sg)
+
+    def update_magnetic(self):
+
+        # Copy previous time step magnetic field values to the previous
+        # time step variables
+        self.update_previous_timestep_fields(self.fn_m)
+
+        for obj in self.magnetic_slices:
+
+            # Grab the main grid fields used to interpolate across the IS
+            # f = self.Hi[slice]
+
+            f_u_1 = obj[-1][obj[2]]
+            f_u_2 = obj[-1][obj[3]]
+            f_u_3 = obj[-1][obj[4]]
+            f_u_4 = obj[-1][obj[5]]
+
+            f_1 = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
+            f_2 = 0.25 * f_u_2 + 0.5 * f_u_3 + 0.25 * f_u_4
+
+            #f_1 = obj[-1][obj[2]]
+            #f_2 = obj[-1][obj[3]]
+            if ('left' in obj[0] or
+                'bottom' in obj[0] or
+                    'front' in obj[0]):
+                w = self.l_weight
+            else:
+                w = self.r_weight
+            c1, c2 = calculate_weighting_coefficients(w, self.ratio)
+            # transverse interpolated h field
+            f_t = c1 * f_1 + c2 * f_2
+
+            # interpolate over a fine grid
+            f_i = self.interpolate_to_sub_grid(f_t, obj[1])
+
+            if f_i == f_t:
+                raise ValueError
+
+            if obj[0] == 'hz_front_1' or obj[0] == 'hz_back_1':
+
+                pre_0 = int(f_i.shape[0] / 2 - 1 - self.ratio // 2)
+                pre_1 = int(f_i.shape[0] / 2 + self.ratio // 2)
+
+                f_i[pre_0 + 1:pre_0 + 1 + self.ratio // 2, :] = np.stack((f_i[pre_0, :], f_i[pre_0, :]))
+                f_i[pre_1 - self.ratio // 2:pre_1, :] = np.stack((f_i[pre_1, :], f_i[pre_1, :]))
+
+
+            # discard the outer nodes only required for interpolation
+            #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+            f = f_i
+            setattr(self, obj[0], f)
+
+    def update_electric(self):
+
+        self.update_previous_timestep_fields(self.fn_e)
+
+        for obj in self.electric_slices:
+
+            f_u_1 = obj[-1][obj[2]]
+            f_u_2 = obj[-1][obj[3]]
+            f_u_3 = obj[-1][obj[4]]
+
+            f_m = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
+
+            f_i = self.interpolate_to_sub_grid(f_m, obj[1])
+
+            # correction interpolated points at the pec traversal point
+            if obj[0] == 'ex_front_1' or obj[0] == 'ex_back_1':
+
+                pre_0 = int(f_i.shape[0] / 2 - 1 - self.ratio // 2)
+                pre_1 = int(f_i.shape[0] / 2 + self.ratio // 2)
+
+                f_i[pre_0 + 1:pre_0 + 1 + self.ratio // 2, :] = np.stack((f_i[pre_0, :], f_i[pre_0, :]))
+                f_i[pre_1 - self.ratio // 2:pre_1, :] = np.stack((f_i[pre_1, :], f_i[pre_1, :]))
+
+            f = f_i
+            #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+            setattr(self, obj[0], f)
+
+    def interpolate_to_sub_grid(self, field, coords):
+        x, z, x_sg, z_sg = coords
+        ex_t = np.transpose(field)
+        f = interpolate.interp2d(x, z, ex_t, kind=self.interpolation)
+        #f = interpolate.RectBivariateSpline(x, z, field)
+        ex_inter_t = f(x_sg, z_sg)
+        ex_inter = np.transpose(ex_inter_t)
+        #ex_inter = ex_inter_t
+
+        return ex_inter
+
+
+class PlaneError(Exception):
+    pass

gprMax/subgrids/solver.py

@@ -0,0 +1,355 @@
+from ..cython.fields_updates_normal import update_electric
+from ..cython.fields_updates_normal import update_magnetic
+from ..fields_outputs import store_outputs
+from ..utilities import get_terminal_width
+from ..exceptions import GeneralError
+
+from .subgrid_hsg import SubGridHSG
+from .precursor_nodes import PrecursorNodes as PrecursorNodesHSG
+from .precursor_nodes_filtered import PrecursorNodes as PrecursorNodesFilteredHSG
+
+from tqdm import tqdm
+from time import perf_counter
+
+import os
+import sys
+
+
+def create_solver(G):
+    """Return the solver for the given subgrids."""
+    updaters = []
+
+    for sg in G.subgrids:
+        print(sg)
+        sg_type = type(sg)
+        if sg_type == SubGridHSG and sg.filter:
+            precursors = PrecursorNodesFilteredHSG(G, sg)
+        elif sg_type == SubGridHSG and not sg.filter:
+            precursors = PrecursorNodesHSG(G, sg)
+        else:
+            raise GeneralError(str(sg) + ' is not a subgrid type')
+
+        sgu = SubgridUpdater(sg, precursors, G)
+        updaters.append(sgu)
+
+    solver = SubGridSolver(G, updaters)
+    return solver
+
+
+class SubGridSolver:
+    """Solver for subgridding simulations."""
+
+    """Class to call the various update methods required for an HSG-Subgrid simulation.
+    Multiple subgrids can be updated by adding more subgrid_updater objects to the subgrid_updater
+    array.
+    """
+
+    def __init__(self, G, subgrid_updaters, hsg=True):
+        """
+        Args:
+            G (G): Grid class instance - holds essential parameters
+            describing the model.
+            subgrid_updaters: (list): list of subgrid_updaters used for updating
+            the subgrids
+            iterations (int): number of iterations for the simulation.
+            hsg (bool): HSG methods for subgrids will not be called if False.
+
+        """
+
+        self.G = G
+        self.grids = [self.G] + G.subgrids
+        self.subgrid_updaters = subgrid_updaters
+        self.materials = G.materials
+        self.abs_time = 0
+        self.hsg = hsg
+
+    def store_outputs(self):
+        """Method to store field outputs for all grid for each main grid iteration."""
+        for grid in self.grids:
+            store_outputs(self.G.iteration, grid.Ex, grid.Ey, grid.Ez,
+                          grid.Hx, grid.Hy, grid.Hz, grid)
+
+    def store_snapshots(self):
+        """Store any snapshots."""
+        for snap in self.G.snapshots:
+            if snap.time == self.G.iteration + 1:
+                snap.store(self.G)
+
+    def hsg_1(self):
+        """Method to update the subgrids over the first phase."""
+        for sg_updater in self.subgrid_updaters:
+            sg_updater.hsg_1()
+
+    def hsg_2(self):
+        """Method to update the subgrids over the second phase."""
+        for sg_updater in self.subgrid_updaters:
+            sg_updater.hsg_2()
+
+    def solve(self, iterations):
+        """Run timestepping."""
+        tsolvestart = perf_counter()
+        self.iterations = iterations
+        # for time step in range(self.G.iterations):
+
+        # The main grid FDTD loop
+        for iteration in self.iterations:
+
+            # Keep track of the index. Required for saving output correctly
+            self.G.iteration = iteration
+
+            # Write any snapshots of the E, H, I fields/currents
+            self.write_snapshots(iteration)
+
+            self.store_outputs()
+            # Update main grid electric field components including sources
+            self.update_magnetic()
+            # Update the fields in the subgrids / main grid
+            if self.hsg:
+                self.hsg_2()
+
+            # Update main grid electric field components including sources
+            self.update_electric()
+            # Update the fields in the subgrids / main grid
+            if self.hsg:
+                self.hsg_1()
+
+        # Return the elapsed time
+        tsolve = perf_counter() - tsolvestart
+
+        return tsolve
+
+    def update_electric(self):
+        """Method to update E fields, PML and electric sources."""
+        # All materials are non-dispersive so do standard update
+        G = self.G
+
+        update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID,
+                        G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+        # Update electric field components with the PML correction
+        for pml in G.pmls:
+            pml.update_electric(G)
+        # Update electric field components from sources (update any Hertzian dipole sources last)
+        for source in G.voltagesources + G.transmissionlines + G.hertziandipoles:
+            source.update_electric(G.iteration, G.updatecoeffsE, G.ID, G.Ex,
+                                   G.Ey, G.Ez, G)
+
+    def update_magnetic(self):
+        """Method to update H fields, PML and magnetic sources."""
+        # Update magnetic field components
+        G = self.G
+
+        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)
+
+        # Update magnetic field components with the PML correction
+        for pml in G.pmls:
+            pml.update_magnetic(G)
+
+        # Update magnetic field components from sources
+        for source in G.transmissionlines + G.magneticdipoles:
+            source.update_magnetic(G.iteration, G.updatecoeffsH, G.ID,
+                                   G.Hx, G.Hy, G.Hz, G)
+
+    def write_snapshots(self, iteration):
+        # Write any snapshots to file
+        for i, snap in enumerate(self.G.snapshots):
+            if snap.time == iteration + 1:
+                snapiters = 36 * (((snap.xf - snap.xs) / snap.dx) * ((snap.yf - snap.ys) / snap.dy) * ((snap.zf - snap.zs) / snap.dz))
+                pbar = tqdm(total=snapiters, leave=False, unit='byte', unit_scale=True, desc='  Writing snapshot file {} of {}, {}'.format(i + 1, len(self.G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=self.G.tqdmdisable)
+
+                # Use this call to print out main grid and subgrids
+                snap.write_vtk_imagedata(self.G.Ex, self.G.Ey, self.G.Ez, self.G.Hx, self.G.Hy, self.G.Hz, self.G, pbar, sub_grids=self.G.subgrids)
+
+                # Use this call to print out the standard grid without subgrid
+                # snap.write_vtk_imagedata(self.G.Ex, self.G.Ey, self.G.Ez, self.G.Hx, self.G.Hy, self.G.Hz, self.G, pbar)
+
+                # Use this call to print out only the subgrid - use in combination with commented code in .multi_cmds/snapshots.py
+                # snap.write_vtk_imagedata_fast(self.grid)
+                pbar.close()
+
+
+class SubgridUpdater:
+    """Class to handle updating the electric and magnetic fields of an HSG
+    subgrid. The IS, OS, subgrid region and the electric/magnetic sources are updated
+    using the precursor regions.
+    """
+
+    def __init__(self, subgrid, precursors, G):
+        """
+        Args:
+            subgrid (SubGrid3d): Subgrid to be updated
+            precursors (PrecursorNodes): Precursor nodes associated with
+            the subgrid
+            G (class): Grid class instance - holds essential parameters
+            describing the model.
+        """
+        self.subgrid = subgrid
+        self.precursors = precursors
+        self.G = G
+        self.source_iteration = 0
+
+    def hsg_1(self):
+        """This is the first half of the subgrid update. Takes the time step
+        up to the main grid magnetic update"""
+        G = self.G
+        sub_grid = self.subgrid
+        precursors = self.precursors
+
+        precursors.update_electric()
+
+        upper_m = int(sub_grid.ratio / 2 - 0.5)
+
+        for m in range(1, upper_m + 1):
+
+            # store_outputs(self.grid)
+            # STD update, interpolate inc. field in time, apply correction
+            update_electric(sub_grid.nx,
+                            sub_grid.ny,
+                            sub_grid.nz,
+                            G.nthreads,
+                            sub_grid.updatecoeffsE,
+                            sub_grid.ID,
+                            sub_grid.Ex,
+                            sub_grid.Ey,
+                            sub_grid.Ez,
+                            sub_grid.Hx,
+                            sub_grid.Hy,
+                            sub_grid.Hz)
+            for pml in sub_grid.pmls:
+                pml.update_electric(sub_grid)
+            precursors.interpolate_magnetic_in_time(int(m + sub_grid.ratio / 2 - 0.5))
+            sub_grid.update_electric_is(precursors)
+
+            self.update_sub_grid_electric_sources()
+
+            # STD update, interpolate inc. field in time, apply correction
+            update_magnetic(sub_grid.nx,
+                            sub_grid.ny,
+                            sub_grid.nz,
+                            G.nthreads,
+                            sub_grid.updatecoeffsH,
+                            sub_grid.ID,
+                            sub_grid.Ex,
+                            sub_grid.Ey,
+                            sub_grid.Ez,
+                            sub_grid.Hx,
+                            sub_grid.Hy,
+                            sub_grid.Hz)
+            for pml in sub_grid.pmls:
+                pml.update_magnetic(sub_grid)
+
+            precursors.interpolate_electric_in_time(m)
+            sub_grid.update_magnetic_is(precursors)
+            self.update_sub_grid_magnetic_sources()
+
+        # store_outputs(self.grid)
+        update_electric(sub_grid.nx,
+                        sub_grid.ny,
+                        sub_grid.nz,
+                        G.nthreads,
+                        sub_grid.updatecoeffsE,
+                        sub_grid.ID,
+                        sub_grid.Ex,
+                        sub_grid.Ey,
+                        sub_grid.Ez,
+                        sub_grid.Hx,
+                        sub_grid.Hy,
+                        sub_grid.Hz)
+        for pml in sub_grid.pmls:
+            pml.update_electric(sub_grid)
+        precursors.calc_exact_magnetic_in_time()
+        sub_grid.update_electric_is(precursors)
+        self.update_sub_grid_electric_sources()
+
+        sub_grid.update_electric_os(G)
+
+    def hsg_2(self):
+        """This is the first half of the subgrid update. Takes the time step
+        up to the main grid electric update"""
+        G = self.G
+        sub_grid = self.subgrid
+        precursors = self.precursors
+
+        precursors.update_magnetic()
+
+        upper_m = int(sub_grid.ratio / 2 - 0.5)
+
+        for m in range(1, upper_m + 1):
+
+            update_magnetic(sub_grid.nx,
+                            sub_grid.ny,
+                            sub_grid.nz,
+                            G.nthreads,
+                            sub_grid.updatecoeffsH,
+                            sub_grid.ID,
+                            sub_grid.Ex,
+                            sub_grid.Ey,
+                            sub_grid.Ez,
+                            sub_grid.Hx,
+                            sub_grid.Hy,
+                            sub_grid.Hz)
+
+            for pml in sub_grid.pmls:
+                pml.update_magnetic(sub_grid)
+
+            precursors.interpolate_electric_in_time(int(m + sub_grid.ratio / 2 - 0.5))
+            sub_grid.update_magnetic_is(precursors)
+            self.update_sub_grid_magnetic_sources()
+
+            # store_outputs(self.grid)
+            update_electric(sub_grid.nx,
+                            sub_grid.ny,
+                            sub_grid.nz,
+                            G.nthreads,
+                            sub_grid.updatecoeffsE,
+                            sub_grid.ID,
+                            sub_grid.Ex,
+                            sub_grid.Ey,
+                            sub_grid.Ez,
+                            sub_grid.Hx,
+                            sub_grid.Hy,
+                            sub_grid.Hz)
+
+            for pml in sub_grid.pmls:
+                pml.update_electric(sub_grid)
+
+            precursors.interpolate_magnetic_in_time(m)
+            sub_grid.update_electric_is(precursors)
+            self.update_sub_grid_electric_sources()
+    
+        update_magnetic(sub_grid.nx,
+                        sub_grid.ny,
+                        sub_grid.nz,
+                        G.nthreads,
+                        sub_grid.updatecoeffsH,
+                        sub_grid.ID,
+                        sub_grid.Ex,
+                        sub_grid.Ey,
+                        sub_grid.Ez,
+                        sub_grid.Hx,
+                        sub_grid.Hy,
+                        sub_grid.Hz)
+        for pml in sub_grid.pmls:
+            pml.update_magnetic(sub_grid)
+
+        precursors.calc_exact_electric_in_time()
+        sub_grid.update_magnetic_is(precursors)
+        self.update_sub_grid_magnetic_sources()
+
+        sub_grid.update_magnetic_os(G)
+
+    def update_sub_grid_electric_sources(self):
+        """Update any electric sources in the subgrid"""
+        sg = self.subgrid
+        for source in sg.voltagesources + sg.transmissionlines + sg.hertziandipoles:
+            source.update_electric(self.source_iteration, sg.updatecoeffsE, sg.ID,
+                                   sg.Ex, sg.Ey, sg.Ez, sg)
+        self.source_iteration += 1
+
+    def update_sub_grid_magnetic_sources(self):
+        """Update any magnetic sources in the subgrid"""
+        sg = self.subgrid
+        for source in sg.transmissionlines + sg.magneticdipoles:
+            source.update_magnetic(self.source_iteration, sg.updatecoeffsH, sg.ID,
+                                   sg.Hx, sg.Hy, sg.Hz, sg)

gprMax/subgrids/subgrid_hsg.py

@@ -1,25 +1,7 @@
-# Copyright (C) 2015-2019: The University of Edinburgh
-#                 Authors: Craig Warren and Antonis Giannopoulos
-#
-# This file is part of gprMax.
-#
-# gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# gprMax is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
-
 from .base import SubGridBase
 from ..cython.fields_updates_hsg import cython_update_is
-from ..cython.fields_updates_hsg import cython_update_electric_os
 from ..cython.fields_updates_hsg import cython_update_magnetic_os
+from ..cython.fields_updates_hsg import cython_update_electric_os
 from ..utilities import human_size
 
 from colorama import init, Fore, Style
@@ -34,89 +16,98 @@ class SubGridHSG(SubGridBase):
         super().__init__(**kwargs)
         self.gridtype = SubGridHSG.gridtype
 
-    def memory_usage(self):
-        pass
-
     def update_magnetic_is(self, precursors):
-        """Apply the incident field corrections to the subgrid fields at the IS.
+        """Update the subgrid nodes at the IS with the currents derived
-
+        from the main grid.
-        Args:
+        Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
-              precursors (PrecursorNodes): PrecursorNodes are the incident
-              fields calculated at the correct time and position.
         """
-        # bottom and top faces
+        # Hz = c0Hz - c1Ey + c2Ex
+        # Hy = c0Hy - c3Ex + c1Ez
+        # Hx = c0Hx - c2Ez + c3Ey
+        # bottom and top
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwy + 1, self.nwz, 1, self.Hy, precursors.ex_bottom, precursors.ex_top, self.IDlookup['Hy'], 1, -1, 3, self.nthreads)
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwy, self.nwz, 1, self.Hx, precursors.ey_bottom, precursors.ey_top, self.IDlookup['Hx'], -1, 1, 3, self.nthreads)
 
-        # left and right faces
+        # left and right
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy, self.nwz + 1, self.nwx, 2, self.Hz, precursors.ey_left, precursors.ey_right, self.IDlookup['Hz'], 1, -1, 1, self.nthreads)
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy + 1, self.nwz, self.nwx, 2, self.Hy, precursors.ez_left, precursors.ez_right, self.IDlookup['Hy'], -1, 1, 1, self.nthreads)
 
-        # front and back faces
+        # front and back
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwz + 1, self.nwy, 3, self.Hz, precursors.ex_front, precursors.ex_back, self.IDlookup['Hz'], -1, 1, 2, self.nthreads)
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, self.nthreads)
 
     def update_electric_is(self, precursors):
-        """Apply the incident field corrections to the subgrid fields at the IS.
+        # Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
 
-        Args:
+        # Ex = c0(Ex) + c2(dHz) - c3(dHy)
-              precursors (PrecursorNodes): PrecursorNodes are the incident
+        # Ey = c0(Ey) + c3(dHx) - c1(dHz)
-              fields calculated at the correct time and position.
+        # Ez = c0(Ez) + c1(dHy) - c2(dHx)
-        """
+
-        # bottom and top faces
+        # bottom and top
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwy + 1, self.nwz, 1, self.Ex, precursors.hy_bottom, precursors.hy_top, self.IDlookup['Ex'], 1, -1, 3, self.nthreads)
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwy, self.nwz, 1, self.Ey, precursors.hx_bottom, precursors.hx_top, self.IDlookup['Ey'], -1, 1, 3, self.nthreads)
 
-        # left and right faces
+        # left and right
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy, self.nwz + 1, self.nwx, 2, self.Ey, precursors.hz_left, precursors.hz_right, self.IDlookup['Ey'], 1, -1, 1, self.nthreads)
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy + 1, self.nwz, self.nwx, 2, self.Ez, precursors.hy_left, precursors.hy_right, self.IDlookup['Ez'], -1, 1, 1, self.nthreads)
 
-        # front and back faces
+        # front and back
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwz + 1, self.nwy, 3, self.Ex, precursors.hz_front, precursors.hz_back, self.IDlookup['Ex'], -1, 1, 2, self.nthreads)
         cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwz, self.nwy, 3, self.Ez, precursors.hx_front, precursors.hx_back, self.IDlookup['Ez'], 1, -1, 2, self.nthreads)
 
     def update_electric_os(self, main_grid):
-        """Apply the incident field corrections to the main grid fields at the OS.
+        i_l = self.i0 - self.is_os_sep
+        i_u = self.i1 + self.is_os_sep
+        j_l = self.j0 - self.is_os_sep
+        j_u = self.j1 + self.is_os_sep
+        k_l = self.k0 - self.is_os_sep
+        k_u = self.k1 + self.is_os_sep
 
-        Main grid fields are collocated with subgrid fields therefore
+        # Args: sub_grid, normal, l_l, l_u, m_l, m_u, n_l, n_u, nwn, lookup_id, field, inc_field, co, sign_n, sign_f
-        precursor fields are not required.
 
-        Args:
+        # Form of FDTD update equations for E
-              main_grid (FDTDGrid): main fdtd solver grid
+        # Ex = c0(Ex) + c2(dHz) - c3(dHy)
-        """
+        # Ey = c0(Ey) + c3(dHx) - c1(dHz)
-        # front and back
+        # Ez = c0(Ez) + c1(dHy) - c2(dHx)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, self.i_l, self.i_u, self.k_l, self.k_u + 1, self.j_l, self.j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, self.i_l, self.i_u + 1, self.k_l, self.k_u, self.j_l, self.j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l, j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l, j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
 
         # Left and Right
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, self.j_l, self.j_u, self.k_l, self.k_u + 1, self.i_l, self.i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l, i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, self.j_l, self.j_u + 1, self.k_l, self.k_u, self.i_l, self.i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l, i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
 
         # Bottom and Top
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, self.i_l, self.i_u, self.j_l, self.j_u + 1, self.k_l, self.k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l, k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, self.i_l, self.i_u + 1, self.j_l, self.j_u, self.k_l, self.k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l, k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
 
     def update_magnetic_os(self, main_grid):
-        """Apply the incident field corrections to the main grid fields at the OS.
 
-        Main grid fields are collocated with subgrid fields therefore
+        i_l = self.i0 - self.is_os_sep
-        precursor fields are not required.
+        i_u = self.i1 + self.is_os_sep
+        j_l = self.j0 - self.is_os_sep
+        j_u = self.j1 + self.is_os_sep
+        k_l = self.k0 - self.is_os_sep
+        k_u = self.k1 + self.is_os_sep
 
-        Args:
+        # Form of FDTD update equations for H
-              main_grid (FDTDGrid): main fdtd solver grid
+        # Hz = c0Hz - c1Ey + c2Ex
-        """
+        # Hy = c0Hy - c3Ex + c1Ez
+        # Hx = c0Hx - c2Ez + c3Ey
+
+        # Args: sub_grid, normal, l_l, l_u, m_l, m_u, n_l, n_u, nwn, lookup_id, field, inc_field, co, sign_n, sign_f):
 
         # Front and back
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, self.i_l, self.i_u, self.k_l, self.k_u + 1, self.j_l - 1, self.j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, self.i_l, self.i_u + 1, self.k_l, self.k_u, self.j_l - 1, self.j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
 
         # Left and Right
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, self.j_l, self.j_u, self.k_l, self.k_u + 1, self.i_l - 1, self.i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, self.j_l, self.j_u + 1, self.k_l, self.k_u, self.i_l - 1, self.i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
 
         # bottom and top
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, self.i_l, self.i_u, self.j_l, self.j_u + 1, self.k_l - 1, self.k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, self.i_l, self.i_u + 1, self.j_l, self.j_u, self.k_l - 1, self.k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
 
     def __str__(self):
 

gprMax/subgrids/updates.py

@@ -16,6 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 from ..updates import CPUUpdates
+from gprMax.utilities import timer
+from ..cython.fields_updates_normal import update_electric
+from ..cython.fields_updates_normal import update_magnetic
 
 
 def create_subgrid_updates(G):
@@ -41,6 +44,78 @@ def create_subgrid_updates(G):
     return updates
 
 
+class SubgridSolver:
+
+    """Generic solver for Update objects"""
+
+    def __init__(self, updates):
+        """Context for the model to run in. Sub-class this with contexts
+        i.e. an MPI context.
+
+        Args:
+            updates (Updates): updates contains methods to run FDTD algorithm
+            iterator (iterator): can be range() or tqdm()
+        """
+        self.updates = updates
+
+    def get_G(self):
+        return self.updates.grid
+
+    def update_electric_1(self):
+        """Method to update E fields, PML and electric sources"""
+        G = self.get_G()
+        update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+        # Update electric field components with the PML correction
+        for pml in G.pmls:
+            pml.update_electric(G)
+        # Update electric field components from sources (update any Hertzian dipole sources last)
+        for source in G.voltagesources + G.transmissionlines + G.hertziandipoles:
+            source.update_electric(G.iteration, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G)
+
+    def update_electric_2(self):
+        return
+
+    def update_magnetic(self):
+        """Method to update H fields, PML and magnetic sources"""
+        # Update magnetic field components
+        G = self.get_G()
+
+        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)
+
+        # Update magnetic field components with the PML correction
+        for pml in G.pmls:
+            pml.update_magnetic(G)
+
+        # Update magnetic field components from sources
+        for source in G.transmissionlines + G.magneticdipoles:
+            source.update_magnetic(G.iteration, G.updatecoeffsH, G.ID, G.Hx, G.Hy, G.Hz, G)
+
+    def solve(self, iterator):
+        """Time step the FDTD model."""
+        tsolvestart = timer()
+
+        for iteration in iterator:
+            # Update main grid electric field components including sources
+            self.updates.store_outputs(iteration)
+            self.update_electric_1()
+            #self.updates.update_electric_pml()
+            #self.updates.update_electric_sources(iteration)
+            # Update the fields in the subgrids / main grid
+            self.updates.hsg_1()
+            # apply 2nd dispersive update after OS updates
+            #self.updates.update_electric_b()
+            # Update main grid electric field components including sources
+            self.update_magnetic()
+            #self.updates.update_magnetic_pml()
+            #self.updates.update_magnetic_sources(iteration)
+            # Update the fields in the subgrids / main grid
+            self.updates.hsg_2()
+
+        tsolve = timer() - tsolvestart
+        return tsolve
+
+
 class SubgridUpdates(CPUUpdates):
     """Top level subgrid updater. Manages the collection of subgrids."""
 
@@ -48,6 +123,11 @@ class SubgridUpdates(CPUUpdates):
         super().__init__(G)
         self.subgrid_updaters = subgrid_updaters
 
+    def store_outputs(self, iteration):
+        super().store_outputs(iteration)
+        for updater in self.subgrid_updaters:
+            updater.store_outputs(iteration)
+
     def hsg_1(self):
         """Method to update the subgrids over the first phase"""
         for sg_updater in self.subgrid_updaters:
@@ -58,14 +138,6 @@ class SubgridUpdates(CPUUpdates):
         for sg_updater in self.subgrid_updaters:
             sg_updater.hsg_2()
 
-    def update_electric_a(self):
-        super().update_electric_a()
-        self.hsg_1()
-
-    def update_magnetic(self):
-        super().update_magnetic()
-        self.hsg_2()
-
     def set_dispersive_updates(self, props):
         # set the dispersive update functions for the main grid updates
         super().set_dispersive_updates(props)
@@ -94,33 +166,51 @@ class SubgridUpdater(CPUUpdates):
         self.source_iteration = 0
         self.G = G
 
+    def update_electric_2(self):
+        return
+
+    def update_electric_1(self):
+        subgrid = self.grid
+        # Update electric field components
+        # All materials are non-dispersive so do standard update
+        update_electric(subgrid.nx, subgrid.ny, subgrid.nz, subgrid.nthreads, subgrid.updatecoeffsE, subgrid.ID, subgrid.Ex, subgrid.Ey, subgrid.Ez, subgrid.Hx, subgrid.Hy, subgrid.Hz)
+
     def hsg_1(self):
-        """This is the first half of the subgrid update. Takes the time step
+        G = self.G
-        up to the main grid magnetic update"""
+        sub_grid = self.grid
+
         precursors = self.precursors
 
         precursors.update_electric()
 
-        # the subgrid
-        sub_grid = self.grid
-
         upper_m = int(sub_grid.ratio / 2 - 0.5)
 
         for m in range(1, upper_m + 1):
 
             # store_outputs(self.grid)
             # STD update, interpolate inc. field in time, apply correction
-            self.update_electric_a()
+            self.update_electric_1()
 
             for pml in sub_grid.pmls:
                 pml.update_electric(sub_grid)
             self.update_sub_grid_electric_sources()
             precursors.interpolate_magnetic_in_time(int(m + sub_grid.ratio / 2 - 0.5))
             sub_grid.update_electric_is(precursors)
-            self.update_electric_b()
+            self.update_electric_2()
 
             # STD update, interpolate inc. field in time, apply correction
-            self.update_magnetic()
+            update_magnetic(sub_grid.nx,
+                            sub_grid.ny,
+                            sub_grid.nz,
+                            G.nthreads,
+                            sub_grid.updatecoeffsH,
+                            sub_grid.ID,
+                            sub_grid.Ex,
+                            sub_grid.Ey,
+                            sub_grid.Ez,
+                            sub_grid.Hx,
+                            sub_grid.Hy,
+                            sub_grid.Hz)
 
             for pml in sub_grid.pmls:
                 pml.update_magnetic(sub_grid)
@@ -129,7 +219,7 @@ class SubgridUpdater(CPUUpdates):
             self.update_sub_grid_magnetic_sources()
 
         # store_outputs(self.grid)
-        self.update_electric_a()
+        self.update_electric_1()
         for pml in sub_grid.pmls:
             pml.update_electric(sub_grid)
 
@@ -137,13 +227,12 @@ class SubgridUpdater(CPUUpdates):
 
         precursors.calc_exact_magnetic_in_time()
         sub_grid.update_electric_is(precursors)
-        self.update_electric_b()
+        self.update_electric_2()
 
-        sub_grid.update_electric_os(self.G)
+        sub_grid.update_electric_os(G)
 
     def hsg_2(self):
-        """This is the first half of the subgrid update. Takes the time step
+        G = self.G
-        up to the main grid electric update"""
         sub_grid = self.grid
         precursors = self.precursors
 
@@ -153,7 +242,18 @@ class SubgridUpdater(CPUUpdates):
 
         for m in range(1, upper_m + 1):
 
-            self.update_magnetic()
+            update_magnetic(sub_grid.nx,
+                            sub_grid.ny,
+                            sub_grid.nz,
+                            G.nthreads,
+                            sub_grid.updatecoeffsH,
+                            sub_grid.ID,
+                            sub_grid.Ex,
+                            sub_grid.Ey,
+                            sub_grid.Ez,
+                            sub_grid.Hx,
+                            sub_grid.Hy,
+                            sub_grid.Hz)
 
             for pml in sub_grid.pmls:
                 pml.update_magnetic(sub_grid)
@@ -163,7 +263,7 @@ class SubgridUpdater(CPUUpdates):
             self.update_sub_grid_magnetic_sources()
 
             # store_outputs(self.grid)
-            self.update_electric_a()
+            self.update_electric_1()
 
             for pml in sub_grid.pmls:
                 pml.update_electric(sub_grid)
@@ -171,9 +271,22 @@ class SubgridUpdater(CPUUpdates):
 
             precursors.interpolate_magnetic_in_time(m)
             sub_grid.update_electric_is(precursors)
-            self.update_electric_b()
+            self.update_electric_2()
 
-        self.update_magnetic()
+        #sub_grid.update_magnetic_os(G)
+
+        update_magnetic(sub_grid.nx,
+                        sub_grid.ny,
+                        sub_grid.nz,
+                        G.nthreads,
+                        sub_grid.updatecoeffsH,
+                        sub_grid.ID,
+                        sub_grid.Ex,
+                        sub_grid.Ey,
+                        sub_grid.Ez,
+                        sub_grid.Hx,
+                        sub_grid.Hy,
+                        sub_grid.Hz)
 
         for pml in sub_grid.pmls:
             pml.update_magnetic(sub_grid)
@@ -182,7 +295,7 @@ class SubgridUpdater(CPUUpdates):
         sub_grid.update_magnetic_is(precursors)
         self.update_sub_grid_magnetic_sources()
 
-        sub_grid.update_magnetic_os(self.G)
+        sub_grid.update_magnetic_os(G)
 
     def update_sub_grid_electric_sources(self):
         """Update any electric sources in the subgrid"""

gprMax/subgrids/user_objects.py

@@ -1,11 +1,10 @@
-from .hsg import SubGridHSG as SubGridHSGUser
+from .subgrid_hsg import SubGridHSG as SubGridHSGUser
 from .multi import ReferenceRx as ReferenceRxUser
 from ..exceptions import CmdInputError
 from ..cmds_multiple import UserObjectMulti
 from ..cmds_geometry.cmds_geometry import UserObjectGeometry
 from ..cmds_multiple import Rx
-import gprMax.config as config
+from gprMax import config
-
 
 from copy import copy
 
@@ -39,23 +38,33 @@ class SubGridBase(UserObjectMulti):
 
     def set_discretisation(self, sg, grid):
         """Set the spatial discretisation."""
-        sg.dl = grid.dl / sg.ratio
+        sg.dx = grid.dx / sg.ratio
-        sg.dx, sg.dy, sg.dz = grid.dl
+        sg.dy = grid.dy / sg.ratio
+        sg.dz = grid.dz / sg.ratio
+        sg.dl = np.array([sg.dx, sg.dy, sg.dz])
 
     def set_main_grid_indices(self, sg, grid, uip, p1, p2):
         """Set subgrid indices related to main grid placement."""
+        # Main grid indices of the sub grid. These are dummy indices. They are
+        # not user internal except for printing to the user
+        sg.i0_u, sg.j0_u, sg.k0_u = p1
+        sg.i1_u, sg.j1_u, sg.k1_u = p2
 
-        # IS indices
+        # The actual sub gridded area (IS index) is 4 cells in
-        sg.i0, sg.j0, sg.k0 = p1
+        sg.i0, sg.j0, sg.k0 = np.add([sg.i0_u, sg.j0_u, sg.k0_u], sg.is_os_sep)
-        sg.i1, sg.j1, sg.k1 = p2
+        sg.i1, sg.j1, sg.k1 = np.subtract([sg.i1_u, sg.j1_u, sg.k1_u], sg.is_os_sep)
 
-        # OS indices
+        # Main grid indices of the sub grid. These are dummy indices. They are
-        sg.i_l, sg.j_l, sg.k_l = p1 - sg.is_os_sep
+        # not user internal except for printing to the user
-        sg.i_u, sg.j_u, sg.k_u = p2 + sg.is_os_sep
+        sg.x1_u, sg.y1_u, sg.z1_u = uip.round_to_grid(p1)
+        sg.x2_u, sg.y2_u, sg.z2_u = uip.round_to_grid(p2)
+
+        sg.x1, sg.y1, sg.z1 = np.add([sg.x1_u, sg.y1_u, sg.z1_u], sg.is_os_sep * sg.dx)
+        sg.x2, sg.y2, sg.z2 = np.subtract([sg.x2_u, sg.y2_u, sg.z2_u], sg.is_os_sep * sg.dx)
+
+    def set_name(self):
+        self.name = self.kwargs['id']
 
-        # discretisted coordinates of the IS
-        sg.x1, sg.y1, sg.z1 = uip.round_to_grid(p1)
-        sg.x2, sg.y2, sg.z2 = uip.round_to_grid(p2)
 
     def set_working_region_cells(self, sg):
         """Number of cells in each dimension for the working region."""
@@ -73,9 +82,6 @@ class SubGridBase(UserObjectMulti):
         """Set number of iterations that will take place in the subgrid."""
         sg.iterations = main.iterations * sg.ratio
 
-    def set_name(self, sg):
-        sg.name = self.kwargs['ID']
-
     def setup(self, sg, grid, uip):
         """"Common setup to both all subgrid types."""
         p1 = self.kwargs['p1']
@@ -83,8 +89,6 @@ class SubGridBase(UserObjectMulti):
 
         p1, p2 = uip.check_box_points(p1, p2, self.__str__())
 
-        self.set_name(sg)
-
         self.check_filters(grid)
 
         self.set_discretisation(sg, grid)
@@ -98,9 +102,20 @@ class SubGridBase(UserObjectMulti):
         # set the indices related to the subgrids main grid placement
         self.set_main_grid_indices(sg, grid, uip, p1, p2)
 
+        """
+        try:
+            uip.check_box_points([sg.i0, sg.j0, sg.k0],
+                                 [sg.i1, sg.j1, sg.k1], cmd_str)
+        except CmdInputError:
+            es_f = 'The subgrid should extend at least {} cells'
+            es = es_f.format(sg.is_os_sep * 2)
+            raise CmdInputError(cmd_str, es)
+        """
+
         self.set_working_region_cells(sg)
         self.set_total_cells(sg)
         self.set_iterations(sg, grid)
+        self.set_name()
 
         # Copy a reference for the main grid to the sub grid
         sg.parent_grid = grid
@@ -133,9 +148,11 @@ class SubGridHSG(SubGridBase):
                  ratio=3,
                  ID='',
                  is_os_sep=3,
+                 pml_separation=4,
                  subgrid_pml_thickness=6,
-                 interpolation=3,
+                 interpolation='linear',
                  loss_mechanism=False,
+                 loss_factor=False,
                  filter=True,
                  **kwargs):
         """Constructor."""
@@ -146,10 +163,12 @@ class SubGridHSG(SubGridBase):
         kwargs['ratio'] = ratio
         kwargs['ID'] = ID
         kwargs['is_os_sep'] = is_os_sep
-        kwargs['pml_separation'] = ratio // 2 + 2
+        kwargs['pml_separation'] = pml_separation
         kwargs['subgrid_pml_thickness'] = subgrid_pml_thickness
         kwargs['interpolation'] = interpolation
         kwargs['filter'] = filter
+        kwargs['loss_mechanism'] = loss_mechanism
+        kwargs['loss_factor'] = loss_factor
 
         super().__init__(**kwargs)
         self.order = 18
@@ -158,10 +177,12 @@ class SubGridHSG(SubGridBase):
     def create(self, grid, uip):
         sg = SubGridHSGUser(**self.kwargs)
         self.setup(sg, grid, uip)
-        if config.general['messages']:
+        if config.is_messages():
             print(sg)
         return sg
 
+
+
 class ReferenceRx(Rx):
     """ReferenceRx User Object."""
 

gprMax/utilities.py

@@ -38,6 +38,7 @@ import numpy as np
 from .exceptions import GeneralError
 
 
+
 def get_terminal_width():
     """Get/set width of terminal being used.
 
@@ -421,8 +422,8 @@ def timer():
 
 class Printer():
 
-    def __init__(self, sim_config):
+    def __init__(self, config):
-        self.printing = sim_config.general['messages']
+        self.printing = config.is_messages()
 
     def print(self, str):
         if self.printing: