transfer existing subgridding modules

gprMax/__init__.py

@@ -8,6 +8,51 @@ Electromagnetic wave propagation simulation software.
 """
 
 from ._version import __version__
+from .cmds_single_use import Discretisation
+from .cmds_single_use import Domain
+from .cmds_single_use import TimeWindow
+from .cmds_single_use import Messages
+from .cmds_single_use import Title
+from .cmds_single_use import NumThreads
+from .cmds_single_use import TimeStepStabilityFactor
+from .cmds_single_use import PMLCells
+from .cmds_single_use import SrcSteps
+from .cmds_single_use import RxSteps
+from .cmds_single_use import ExcitationFile
+from .cmds_multiple import Waveform
+from .cmds_multiple import VoltageSource
+from .cmds_multiple import HertzianDipole
+from .cmds_multiple import MagneticDipole
+from .cmds_multiple import TransmissionLine
+from .cmds_multiple import Rx
+from .cmds_multiple import RxArray
+from .cmds_multiple import Snapshot
+from .cmds_multiple import Material
+from .cmds_multiple import AddDebyeDispersion
+from .cmds_multiple import AddLorentzDispersion
+from .cmds_multiple import AddDrudeDispersion
+from .cmds_multiple import SoilPeplinski
+from .cmds_multiple import GeometryView
+from .cmds_multiple import GeometryObjectsWrite
+from .cmds_multiple import PMLCFS
+from .subgrids.user_objects import SubGridHSG
+from .subgrids.user_objects import SubGridSHSG
+from .subgrids.user_objects import ReferenceRx
+
+# import geometry commands
+from .cmds_geometry.edge import Edge
+from .cmds_geometry.plate import Plate
+from .cmds_geometry.triangle import Triangle
+from .cmds_geometry.box import Box
+from .cmds_geometry.cylinder import Cylinder
+from .cmds_geometry.cylindrical_sector import CylindricalSector
+from .cmds_geometry.sphere import Sphere
+from .cmds_geometry.fractal_box import FractalBox
+from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness
+from .cmds_geometry.add_surface_water import AddSurfaceWater
+from .cmds_geometry.add_grass import AddGrass
+
+from .scene import Scene
 from .gprMax import api as run
 
 __name__ = 'gprMax'

gprMax/cmds_geometry/add_grass.py

@@ -0,0 +1,164 @@
+"""Class for add_grass command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..utilities import round_value
+from ..materials import Material
+from ..fractals import FractalSurface
+from ..fractals import Grass
+
+from tqdm import tqdm
+import numpy as np
+
+
+class AddGrass(UserObjectGeometry):
+    """User class for Grass command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 12
+        self.hash = '#add_grass'
+
+    def create(self, grid, uip):
+        """Add Grass to fractal box."""
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+            fractal_box_id = self.kwargs['fractal_box_id']
+            frac_dim = self.kwargs['frac_dim']
+            limits = self.kwargs['limits']
+            n_blades = self.kwargs['n_blades']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' requires at least eleven parameters')
+
+        try:
+            seed = self.kwargs['seed']
+        except KeyError:
+            seed = None
+
+        # grab the correct fractal volume
+        volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
+        if volumes:
+            volume = volumes[0]
+        else:
+            raise CmdInputError(self.__str__() + ' Cant find FractalBox {}'.format(fractal_box_id))
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        if frac_dim < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal dimension')
+        if limits[0] < 0 or limits[1] < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the minimum and maximum heights for grass blades')
+
+        # Check for valid orientations
+        if xs == xf:
+            if ys == yf or zs == zf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if xs != volume.xs and xs != volume.xf:
+                raise CmdInputError(self.__str__() + ' must specify external surfaces on a fractal box')
+            fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
+            # xminus surface
+            if xs == volume.xs:
+                raise CmdInputError(self.__str__() + ' grass can only be specified on surfaces in the positive axis direction')
+            # xplus surface
+            elif xf == volume.xf:
+                if fractalrange[1] > grid.nx:
+                    raise CmdInputError(self.__str__() + ' cannot apply grass to fractal box as it would exceed the domain size in the x direction')
+                requestedsurface = 'xplus'
+
+        elif ys == yf:
+            if xs == xf or zs == zf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if ys != volume.ys and ys != volume.yf:
+                raise CmdInputError(self.__str__() + ' must specify external surfaces on a fractal box')
+            fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
+            # yminus surface
+            if ys == volume.ys:
+                raise CmdInputError(self.__str__() + ' grass can only be specified on surfaces in the positive axis direction')
+            # yplus surface
+            elif yf == volume.yf:
+                if fractalrange[1] > grid.ny:
+                    raise CmdInputError(self.__str__() + ' cannot apply grass to fractal box as it would exceed the domain size in the y direction')
+                requestedsurface = 'yplus'
+
+        elif zs == zf:
+            if xs == xf or ys == yf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if zs != volume.zs and zs != volume.zf:
+                raise CmdInputError(self.__str__() + ' must specify external surfaces on a fractal box')
+            fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
+            # zminus surface
+            if zs == volume.zs:
+                raise CmdInputError(self.__str__() + ' grass can only be specified on surfaces in the positive axis direction')
+            # zplus surface
+            elif zf == volume.zf:
+                if fractalrange[1] > grid.nz:
+                    raise CmdInputError(self.__str__() + ' cannot apply grass to fractal box as it would exceed the domain size in the z direction')
+                requestedsurface = 'zplus'
+
+        else:
+            raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+
+        surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim)
+        surface.ID = 'grass'
+        surface.surfaceID = requestedsurface
+        surface.seed = seed
+
+        # Set the fractal range to scale the fractal distribution between zero and one
+        surface.fractalrange = (0, 1)
+        surface.operatingonID = volume.ID
+        surface.generate_fractal_surface(grid)
+        if n_blades > surface.fractalsurface.shape[0] * surface.fractalsurface.shape[1]:
+            raise CmdInputError(self.__str__() + ' the specified surface is not large enough for the number of grass blades/roots specified')
+
+        # Scale the distribution so that the summation is equal to one, i.e. a probability distribution
+        surface.fractalsurface = surface.fractalsurface / np.sum(surface.fractalsurface)
+
+        # Set location of grass blades using probability distribution
+        # Create 1D vector of probability values from the 2D surface
+        probability1D = np.cumsum(np.ravel(surface.fractalsurface))
+
+        # Create random numbers between zero and one for the number of blades of grass
+        R = np.random.RandomState(surface.seed)
+        A = R.random_sample(n_blades)
+
+        # Locate the random numbers in the bins created by the 1D vector of probability values, and convert the 1D index back into a x, y index for the original surface.
+        bladesindex = np.unravel_index(np.digitize(A, probability1D), (surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
+
+        # Set the fractal range to minimum and maximum heights of the grass blades
+        surface.fractalrange = fractalrange
+
+        # Set the fractal surface using the pre-calculated spatial distribution and a random height
+        surface.fractalsurface = np.zeros((surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
+        for i in range(len(bladesindex[0])):
+                surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(surface.fractalrange[0], surface.fractalrange[1], size=1)
+
+        # Create grass geometry parameters
+        g = Grass(n_blades)
+        g.seed = surface.seed
+        surface.grass.append(g)
+
+        # Check to see if grass has been already defined as a material
+        if not any(x.ID == 'grass' for x in grid.materials):
+            m = Material(len(grid.materials), 'grass')
+            m.averagable = False
+            m.type = 'builtin, debye'
+            m.er = Material.grasseri
+            m.deltaer.append(Material.grassdeltaer)
+            m.tau.append(Material.grasstau)
+            grid.materials.append(m)
+            if Material.maxpoles == 0:
+                Material.maxpoles = 1
+
+        # Check if time step for model is suitable for using grass
+        grass = next((x for x in grid.materials if x.ID == 'grass'))
+        testgrass = next((x for x in grass.tau if x < grid.dt), None)
+        if testgrass:
+            raise CmdInputError(self.__str__() + ' requires the time step for the model to be less than the relaxation time required to model grass.')
+
+        volume.fractalsurfaces.append(surface)
+
+        if grid.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

@@ -0,0 +1,126 @@
+"""Class for surface roughness command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..fractals import FractalSurface
+from ..utilities import round_value
+
+from tqdm import tqdm
+import numpy as np
+
+
+class AddSurfaceRoughness(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 10
+        self.hash = '#add_surface_roughness'
+
+    def create(self, grid, uip):
+
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+            frac_dim = self.kwargs['frac_dim']
+            weighting = np.array(self.kwargs['weighting'])
+            limits = self.kwargs['limits']
+            fractal_box_id = self.kwargs['fractal_box_id']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' Incorrect parameters')
+
+        try:
+            seed = self.kwargs['seed']
+        except KeyError:
+            seed = None
+
+        # grab the correct fractal volume
+        volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
+        if volumes:
+            volume = volumes[0]
+        else:
+            raise CmdInputError(self.__str__() + ' Cant find FractalBox {}'.format(fractal_box_id))
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        if frac_dim < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal dimension')
+        if weighting[0] < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal weighting in the first direction of the surface')
+        if weighting[1] < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal weighting in the second direction of the surface')
+
+        # Check for valid orientations
+        if xs == xf:
+            if ys == yf or zs == zf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if xs != volume.xs and xs != volume.xf:
+                raise CmdInputError(self.__str__() + ' can only be used on the external surfaces of a fractal box')
+            fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
+            # xminus surface
+            if xs == volume.xs:
+                if fractalrange[0] < 0 or fractalrange[1] > volume.xf:
+                    raise CmdInputError(self.__str__() + ' cannot apply fractal surface to fractal box as it would exceed either the upper coordinates of the fractal box or the domain in the x direction')
+                requestedsurface = 'xminus'
+            # xplus surface
+            elif xf == volume.xf:
+                if fractalrange[0] < volume.xs or fractalrange[1] > grid.nx:
+                    raise CmdInputError(self.__str__() + ' cannot apply fractal surface to fractal box as it would exceed either the lower coordinates of the fractal box or the domain in the x direction')
+                requestedsurface = 'xplus'
+
+        elif ys == yf:
+            if xs == xf or zs == zf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if ys != volume.ys and ys != volume.yf:
+                raise CmdInputError(self.__str__() + ' can only be used on the external surfaces of a fractal box')
+            fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
+            # yminus surface
+            if ys == volume.ys:
+                if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
+                    raise CmdInputError(self.__str__() + ' cannot apply fractal surface to fractal box as it would exceed either the upper coordinates of the fractal box or the domain in the y direction')
+                requestedsurface = 'yminus'
+            # yplus surface
+            elif yf == volume.yf:
+                if fractalrange[0] < volume.ys or fractalrange[1] > grid.ny:
+                    raise CmdInputError(self.__str__() + ' cannot apply fractal surface to fractal box as it would exceed either the lower coordinates of the fractal box or the domain in the y direction')
+                requestedsurface = 'yplus'
+
+        elif zs == zf:
+            if xs == xf or ys == yf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if zs != volume.zs and zs != volume.zf:
+                raise CmdInputError(self.__str__() + ' can only be used on the external surfaces of a fractal box')
+            fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
+            # zminus surface
+            if zs == volume.zs:
+                if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
+                    raise CmdInputError(self.__str__() + ' cannot apply fractal surface to fractal box as it would exceed either the upper coordinates of the fractal box or the domain in the x direction')
+                requestedsurface = 'zminus'
+            # zplus surface
+            elif zf == volume.zf:
+                if fractalrange[0] < volume.zs or fractalrange[1] > grid.nz:
+                    raise CmdInputError(self.__str__() + ' cannot apply fractal surface to fractal box as it would exceed either the lower coordinates of the fractal box or the domain in the z direction')
+                requestedsurface = 'zplus'
+
+        else:
+            raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+
+        surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim)
+        surface.surfaceID = requestedsurface
+        surface.fractalrange = fractalrange
+        surface.operatingonID = volume.ID
+        surface.seed = seed
+        surface.weighting = weighting
+
+        # List of existing surfaces IDs
+        existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
+        if surface.surfaceID in existingsurfaceIDs:
+            raise CmdInputError(self.__str__() + ' has already been used on the {} surface'.format(surface.surfaceID))
+
+        surface.generate_fractal_surface(grid)
+        volume.fractalsurfaces.append(surface)
+
+        if grid.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

@@ -0,0 +1,118 @@
+"""Class for surface water command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..utilities import round_value
+from ..materials import Material
+
+from tqdm import tqdm
+
+
+class AddSurfaceWater(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 11
+        self.hash = '#add_surface_water'
+
+    def create(self, grid, uip):
+        """"Create surface water on fractal box."""
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+            fractal_box_id = self.kwargs['fractal_box_id']
+            depth = self.kwargs['depth']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' requires exactly eight parameters')
+
+        # grab the correct fractal volume
+        volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
+        if volumes:
+            volume = volumes[0]
+        else:
+            raise CmdInputError(self.__str__() + ' Cant find FractalBox {}'.format(fractal_box_id))
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        if depth <= 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the depth of water')
+
+        # Check for valid orientations
+        if xs == xf:
+            if ys == yf or zs == zf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if xs != volume.xs and xs != volume.xf:
+                raise CmdInputError(self.__str__() + ' can only be used on the external surfaces of a fractal box')
+            # xminus surface
+            if xs == volume.xs:
+                requestedsurface = 'xminus'
+            # xplus surface
+            elif xf == volume.xf:
+                requestedsurface = 'xplus'
+            filldepthcells = round_value(depth / grid.dx)
+            filldepth = filldepthcells * grid.dx
+
+        elif ys == yf:
+            if xs == xf or zs == zf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if ys != volume.ys and ys != volume.yf:
+                raise CmdInputError(self.__str__() + ' can only be used on the external surfaces of a fractal box')
+            # yminus surface
+            if ys == volume.ys:
+                requestedsurface = 'yminus'
+            # yplus surface
+            elif yf == volume.yf:
+                requestedsurface = 'yplus'
+            filldepthcells = round_value(depth / grid.dy)
+            filldepth = filldepthcells * grid.dy
+
+        elif zs == zf:
+            if xs == xf or ys == yf:
+                raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+            if zs != volume.zs and zs != volume.zf:
+                raise CmdInputError(self.__str__() + ' can only be used on the external surfaces of a fractal box')
+            # zminus surface
+            if zs == volume.zs:
+                requestedsurface = 'zminus'
+            # zplus surface
+            elif zf == volume.zf:
+                requestedsurface = 'zplus'
+            filldepthcells = round_value(depth / grid.dz)
+            filldepth = filldepthcells * grid.dz
+
+        else:
+            raise CmdInputError(self.__str__() + ' dimensions are not specified correctly')
+
+        surface = next((x for x in volume.fractalsurfaces if x.surfaceID == requestedsurface), None)
+        if not surface:
+            raise CmdInputError(self.__str__() + ' specified surface {} does not have a rough surface applied'.format(requestedsurface))
+
+        surface.filldepth = filldepthcells
+
+        # Check that requested fill depth falls within range of surface roughness
+        if surface.filldepth < surface.fractalrange[0] or surface.filldepth > surface.fractalrange[1]:
+            raise CmdInputError(self.__str__() + ' requires a value for the depth of water that lies with the range of the requested surface roughness')
+
+        # Check to see if water has been already defined as a material
+        if not any(x.ID == 'water' for x in grid.materials):
+            m = Material(len(grid.materials), 'water')
+            m.averagable = False
+            m.type = 'builtin, debye'
+            m.er = Material.watereri
+            m.deltaer.append(Material.waterdeltaer)
+            m.tau.append(Material.watertau)
+            grid.materials.append(m)
+            if Material.maxpoles == 0:
+                Material.maxpoles = 1
+
+        # Check if time step for model is suitable for using water
+        water = next((x for x in grid.materials if x.ID == 'water'))
+        testwater = next((x for x in water.tau if x < grid.dt), None)
+        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:
+            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

@@ -0,0 +1,92 @@
+"""Class for triangle command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..materials import Material
+from ..geometry_primitives_ext import build_box
+
+from tqdm import tqdm
+import numpy as np
+
+
+class Box(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 5
+        self.hash = '#box'
+
+    def create(self, grid, uip):
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' Please specify two points.')
+        # check materials have been specified
+        # isotropic case
+        try:
+            materialsrequested = [self.kwargs['material_id']]
+        except KeyError:
+            # Anisotropic case
+            try:
+                materialsrequested = self.kwargs['material_ids']
+            except KeyError:
+                raise CmdInputError(self.__str__() + ' No materials have been specified')
+
+        # check averaging
+        try:
+            # go with user specified averaging
+            averagebox = self.kwargs['averaging']
+        except KeyError:
+            # if they havent specfied - go with the grid default
+            averagebox = grid.averagevolumeobjects
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        # Look up requested materials in existing list of material instances
+        materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+        if len(materials) != len(materialsrequested):
+            notfound = [x for x in materialsrequested if x not in materials]
+            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+
+        # Isotropic case
+        if len(materials) == 1:
+            averaging = materials[0].averagable and averagebox
+            numID = numIDx = numIDy = numIDz = materials[0].numID
+
+        # Uniaxial anisotropic case
+        elif len(materials) == 3:
+            averaging = False
+            numIDx = materials[0].numID
+            numIDy = materials[1].numID
+            numIDz = materials[2].numID
+            requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+            averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+            if averagedmaterial:
+                numID = averagedmaterial.numID
+            else:
+                numID = len(grid.materials)
+                m = Material(numID, requiredID)
+                m.type = 'dielectric-smoothed'
+                # Create dielectric-smoothed constituents for material
+                m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
+                m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
+                m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+
+                # Append the new material object to the materials list
+                grid.materials.append(m)
+
+        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 averaging:
+                dielectricsmoothing = 'on'
+            else:
+                dielectricsmoothing = 'off'
+            tqdm.write('Box from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created, dielectric smoothing is {}.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, ', '.join(materialsrequested), dielectricsmoothing))

gprMax/cmds_geometry/cmds_geometry.py

@@ -0,0 +1,28 @@
+
+class UserObjectGeometry:
+    """Specific Geometry object."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        self.kwargs = kwargs
+        # define the order of priority for calling create()
+        self.order = None
+        # hash command
+        self.hash = '#example'
+
+    def __str__(self):
+        """Readble user string as per hash commands."""
+        s = ''
+        for k, v in self.kwargs.items():
+            if isinstance(v, tuple) or isinstance(v, list):
+                v = ' '.join([str(el) for el in v])
+            s += str(v) + ' '
+
+        return '{}: {}'.format(self.hash, s[:-1])
+
+    def params_str(self):
+        """Readble string of parameters given to object."""
+        return self.hash + ': ' + str(self.kwargs)
+
+    def create(self, grid, uip):
+        """Create the object and add it to the grid."""

gprMax/cmds_geometry/cylinder.py

@@ -0,0 +1,97 @@
+"""Class for cylinder command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..materials import Material
+from ..geometry_primitives_ext import build_cylinder
+
+from tqdm import tqdm
+import numpy as np
+
+
+class Cylinder(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 6
+        self.hash = '#cylinder'
+
+    def create(self, grid, uip):
+
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+            r = self.kwargs['r']
+
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' Please specify 2 points and a radius')
+
+        # check averaging
+        try:
+            # go with user specified averaging
+            averagecylinder = self.kwargs['averaging']
+        except KeyError:
+            # if they havent specfied - go with the grid default
+            averagecylinder = grid.averagevolumeobjects
+
+        # check materials have been specified
+        # isotropic case
+        try:
+            materialsrequested = [self.kwargs['material_id']]
+        except KeyError:
+            # Anisotropic case
+            try:
+                materialsrequested = self.kwargs['material_ids']
+            except KeyError:
+                raise CmdInputError(self.__str__() + ' No materials have been specified')
+
+        x1, y1, z1 = uip.round_to_grid(p1)
+        x2, y2, z2 = uip.round_to_grid(p2)
+
+        if r <= 0:
+            raise CmdInputError(self.__str__() + ' the radius {:g} should be a positive value.'.format(r))
+
+        # Look up requested materials in existing list of material instances
+        materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+        if len(materials) != len(materialsrequested):
+            notfound = [x for x in materialsrequested if x not in materials]
+            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+
+        # Isotropic case
+        if len(materials) == 1:
+            averaging = materials[0].averagable and averagecylinder
+            numID = numIDx = numIDy = numIDz = materials[0].numID
+
+        # Uniaxial anisotropic case
+        elif len(materials) == 3:
+            averaging = False
+            numIDx = materials[0].numID
+            numIDy = materials[1].numID
+            numIDz = materials[2].numID
+            requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+            averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+            if averagedmaterial:
+                numID = averagedmaterial.numID
+            else:
+                numID = len(grid.materials)
+                m = Material(numID, requiredID)
+                m.type = 'dielectric-smoothed'
+                # Create dielectric-smoothed constituents for material
+                m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
+                m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
+                m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+
+                # Append the new material object to the materials list
+                grid.materials.append(m)
+
+        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 averaging:
+                dielectricsmoothing = 'on'
+            else:
+                dielectricsmoothing = 'off'
+            tqdm.write('Cylinder with face centres {:g}m, {:g}m, {:g}m and {:g}m, {:g}m, {:g}m, with radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(x1, y1, z1, x2, y2, z2, r, ', '.join(materialsrequested), dielectricsmoothing))

gprMax/cmds_geometry/cylindrical_sector.py

@@ -0,0 +1,137 @@
+"""Class for cylinder command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..materials import Material
+from ..geometry_primitives_ext import build_cylindrical_sector
+
+from tqdm import tqdm
+import numpy as np
+
+
+class CylindricalSector(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 7
+        self.hash = '#cylindrical_sector'
+
+    def create(self, grid, uip):
+
+        try:
+            normal = self.kwargs['normal'].lower()
+            ctr1 = self.kwargs['ctr1']
+            ctr2 = self.kwargs['ctr2']
+            extent1 = self.kwargs['extent1']
+            extent2 = self.kwargs['extent2']
+            start = self.kwargs['start']
+            end = self.kwargs['end']
+            r = self.kwargs['r']
+            thickness = extent2 - extent1
+        except KeyError:
+            raise CmdInputError(self.__str__())
+
+        # check averaging
+        try:
+            # go with user specified averaging
+            averagecylindricalsector = self.kwargs['averaging']
+        except KeyError:
+            # if they havent specfied - go with the grid default
+            averagecylindricalsector = grid.averagevolumeobjects
+
+        # check materials have been specified
+        # isotropic case
+        try:
+            materialsrequested = [self.kwargs['material_id']]
+        except KeyError:
+            # Anisotropic case
+            try:
+                materialsrequested = self.kwargs['material_ids']
+            except KeyError:
+                raise CmdInputError(self.__str__() + ' No materials have been specified')
+
+        sectorstartangle = 2 * np.pi * (start / 360)
+        sectorangle = 2 * np.pi * (end / 360)
+
+        if normal != 'x' and normal != 'y' and normal != 'z':
+            raise CmdInputError(self.__str__() + ' the normal direction must be either x, y or z.')
+        if r <= 0:
+            raise CmdInputError(self.__str__() + ' the radius {:g} should be a positive value.'.format(r))
+        if sectorstartangle < 0 or sectorangle <= 0:
+            raise CmdInputError(self.__str__() + ' the starting angle and sector angle should be a positive values.')
+        if sectorstartangle >= 2 * np.pi or sectorangle >= 2 * np.pi:
+            raise CmdInputError(self.__str__() + ' the starting angle and sector angle must be less than 360 degrees.')
+
+        # Look up requested materials in existing list of material instances
+        materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+        if len(materials) != len(materialsrequested):
+            notfound = [x for x in materialsrequested if x not in materials]
+            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+
+        if thickness > 0:
+            # Isotropic case
+            if len(materials) == 1:
+                averaging = materials[0].averagable and averagecylindricalsector
+                numID = numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 3:
+                averaging = False
+                numIDx = materials[0].numID
+                numIDy = materials[1].numID
+                numIDz = materials[2].numID
+                requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+                averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+                if averagedmaterial:
+                    numID = averagedmaterial.numID
+                else:
+                    numID = len(grid.materials)
+                    m = Material(numID, requiredID)
+                    m.type = 'dielectric-smoothed'
+                    # Create dielectric-smoothed constituents for material
+                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
+                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
+                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+
+                    # Append the new material object to the materials list
+                    grid.materials.append(m)
+        else:
+            averaging = False
+            # Isotropic case
+            if len(materials) == 1:
+                numID = numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 3:
+                # numID requires a value but it will not be used
+                numID = None
+                numIDx = materials[0].numID
+                numIDy = materials[1].numID
+                numIDz = materials[2].numID
+
+        # yz-plane cylindrical sector
+        if normal == 'x':
+            level, ctr1, ctr2 = uip.round_to_grid((extent1, ctr1, ctr2))
+
+        # xz-plane cylindrical sector
+        elif normal == 'y':
+            ctr1, level, ctr2 = uip.round_to_grid((ctr1, extent1, ctr2))
+
+        # xy-plane cylindrical sector
+        elif normal == 'z':
+            ctr1, ctr2, level = uip.round_to_grid((ctr1, ctr2, extent1))
+
+        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 thickness > 0:
+                if averaging:
+                    dielectricsmoothing = 'on'
+                else:
+                    dielectricsmoothing = 'off'
+                tqdm.write('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, thickness {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, thickness, ', '.join(materialsrequested), dielectricsmoothing))
+            else:
+                tqdm.write('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, of material(s) {} created.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, ', '.join(materialsrequested)))

gprMax/cmds_geometry/edge.py

@@ -0,0 +1,64 @@
+"""Class for edge command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..geometry_primitives_ext import build_edge_x
+from ..geometry_primitives_ext import build_edge_y
+from ..geometry_primitives_ext import build_edge_z
+
+from tqdm import tqdm
+
+
+class Edge(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 2
+        self.hash = '#edge'
+
+    def create(self, grid, uip):
+        """Create edge and add it to the grid."""
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+            material_id = self.kwargs['material_id']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' requires exactly 3 parameters')
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        material = next((x for x in grid.materials if x.ID == material_id), None)
+
+        if not material:
+            raise CmdInputError('Material with ID {} does not exist'.format(material_id))
+
+        # Check for valid orientations
+        # x-orientated wire
+        if xs != xf:
+            if ys != yf or zs != zf:
+                raise CmdInputError(self.__str__() + ' the edge is not specified correctly')
+            else:
+                for i in range(xs, xf):
+                    build_edge_x(i, ys, zs, material.numID, grid.rigidE, grid.rigidH, grid.ID)
+
+        # y-orientated wire
+        elif ys != yf:
+            if xs != xf or zs != zf:
+                raise CmdInputError(self.__str__() + ' the edge is not specified correctly')
+            else:
+                for j in range(ys, yf):
+                    build_edge_y(xs, j, zs, material.numID, grid.rigidE, grid.rigidH, grid.ID)
+
+        # z-orientated wire
+        elif zs != zf:
+            if xs != xf or ys != yf:
+                raise CmdInputError(self.__str__() + ' the edge is not specified correctly')
+            else:
+                for k in range(zs, zf):
+                    build_edge_z(xs, ys, k, material.numID, grid.rigidE, grid.rigidH, grid.ID)
+
+        if grid.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

@@ -0,0 +1,92 @@
+"""Class for surface roughness command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..fractals import FractalVolume
+
+from tqdm import tqdm
+import numpy as np
+
+
+class FractalBox(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 9
+        self.hash = '#fractal_box'
+
+    def create(self, grid, uip):
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+            frac_dim = self.kwargs['frac_dim']
+            weighting = np.array(self.kwargs['weighting'])
+            n_materials = self.kwargs['n_materials']
+            mixing_model_id = self.kwargs['mixing_model_id']
+            ID = self.kwargs['ID']
+
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' Incorrect parameters')
+
+        try:
+            seed = self.kwargs['seed']
+        except KeyError:
+            seed = None
+
+        # Default is no dielectric smoothing for a fractal box
+        averagefractalbox = False
+
+        # check averaging
+        try:
+            # go with user specified averaging
+            averagefractalbox = self.kwargs['averaging']
+        except KeyError:
+            # if they havent specfied - go with the grid default
+            averagefractalbox = grid.averagevolumeobjects
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        if frac_dim < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal dimension')
+        if weighting[0] < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal weighting in the x direction')
+        if weighting[1] < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal weighting in the y direction')
+        if weighting[2] < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the fractal weighting in the z direction')
+        if n_materials < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for the number of bins')
+
+        # Find materials to use to build fractal volume, either from mixing models or normal materials
+        mixingmodel = next((x for x in grid.mixingmodels if x.ID == mixing_model_id), None)
+        material = next((x for x in grid.materials if x.ID == mixing_model_id), None)
+        nbins = n_materials
+
+        if mixingmodel:
+            if nbins == 1:
+                raise CmdInputError(self.__str__() + ' must be used with more than one material from the mixing model.')
+            # Create materials from mixing model as number of bins now known from fractal_box command
+            mixingmodel.calculate_debye_properties(nbins, grid)
+        elif not material:
+            raise CmdInputError(self.__str__() + ' mixing model or material with ID {} does not exist'.format(mixing_model_id))
+
+        volume = FractalVolume(xs, xf, ys, yf, zs, zf, frac_dim)
+        volume.ID = ID
+        volume.operatingonID = mixing_model_id
+        volume.nbins = nbins
+        volume.seed = seed
+        volume.weighting = weighting
+        volume.averaging = averagefractalbox
+        volume.mixingmodel = mixingmodel
+
+        if grid.messages:
+            if volume.averaging:
+                dielectricsmoothing = 'on'
+            else:
+                dielectricsmoothing = 'off'
+            tqdm.write('Fractal box {} from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with {}, fractal dimension {:g}, fractal weightings {:g}, {:g}, {:g}, fractal seeding {}, with {} material(s) created, dielectric smoothing is {}.'.format(volume.ID, xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, volume.operatingonID, volume.dimension, volume.weighting[0], volume.weighting[1], volume.weighting[2], volume.seed, volume.nbins, dielectricsmoothing))
+
+        grid.fractalvolumes.append(volume)

gprMax/cmds_geometry/fractal_box_builder.py

@@ -0,0 +1,286 @@
+"""Class for fractal box modifications command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..constants import floattype
+from gprMax.geometry_primitives_ext import build_voxels_from_array
+from gprMax.geometry_primitives_ext import build_voxels_from_array_mask
+
+import numpy as np
+
+
+class FractalBoxBuilder(UserObjectGeometry):
+    """Internal class for fractal box modifications. This class should be used
+    internally only when surface modification have been made to a fractal box"""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 13
+        self.hash = '#fractal_box_modifications'
+
+    def create(self, grid, uip):
+
+        for volume in grid.fractalvolumes:
+
+            if volume.fractalsurfaces:
+                volume.originalxs = volume.xs
+                volume.originalxf = volume.xf
+                volume.originalys = volume.ys
+                volume.originalyf = volume.yf
+                volume.originalzs = volume.zs
+                volume.originalzf = volume.zf
+
+                # Extend the volume to accomodate any rough surfaces, grass, or roots
+                for surface in volume.fractalsurfaces:
+                    if surface.surfaceID == 'xminus':
+                        if surface.fractalrange[0] < volume.xs:
+                            volume.nx += volume.xs - surface.fractalrange[0]
+                            volume.xs = surface.fractalrange[0]
+                    elif surface.surfaceID == 'xplus':
+                        if surface.fractalrange[1] > volume.xf:
+                            volume.nx += surface.fractalrange[1] - volume.xf
+                            volume.xf = surface.fractalrange[1]
+                    elif surface.surfaceID == 'yminus':
+                        if surface.fractalrange[0] < volume.ys:
+                            volume.ny += volume.ys - surface.fractalrange[0]
+                            volume.ys = surface.fractalrange[0]
+                    elif surface.surfaceID == 'yplus':
+                        if surface.fractalrange[1] > volume.yf:
+                            volume.ny += surface.fractalrange[1] - volume.yf
+                            volume.yf = surface.fractalrange[1]
+                    elif surface.surfaceID == 'zminus':
+                        if surface.fractalrange[0] < volume.zs:
+                            volume.nz += volume.zs - surface.fractalrange[0]
+                            volume.zs = surface.fractalrange[0]
+                    elif surface.surfaceID == 'zplus':
+                        if surface.fractalrange[1] > volume.zf:
+                            volume.nz += surface.fractalrange[1] - volume.zf
+                            volume.zf = surface.fractalrange[1]
+
+                # If there is only 1 bin then a normal material is being used, otherwise a mixing model
+                if volume.nbins == 1:
+                    volume.fractalvolume = np.ones((volume.nx, volume.ny, volume.nz), dtype=floattype)
+                    materialnumID = next(x.numID for x in grid.materials if x.ID == volume.operatingonID)
+                    volume.fractalvolume *= materialnumID
+                else:
+                    volume.generate_fractal_volume(grid)
+                    volume.fractalvolume += volume.mixingmodel.startmaterialnum
+
+                volume.generate_volume_mask()
+
+                # Apply any rough surfaces and add any surface water to the 3D mask array
+                for surface in volume.fractalsurfaces:
+                    if surface.surfaceID == 'xminus':
+                        for i in range(surface.fractalrange[0], surface.fractalrange[1]):
+                            for j in range(surface.ys, surface.yf):
+                                for k in range(surface.zs, surface.zf):
+                                    if i > surface.fractalsurface[j - surface.ys, k - surface.zs]:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
+                                    elif surface.filldepth > 0 and i > surface.filldepth:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
+                                    else:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
+
+                    elif surface.surfaceID == 'xplus':
+                        if not surface.ID:
+                            for i in range(surface.fractalrange[0], surface.fractalrange[1]):
+                                for j in range(surface.ys, surface.yf):
+                                    for k in range(surface.zs, surface.zf):
+                                        if i < surface.fractalsurface[j - surface.ys, k - surface.zs]:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
+                                        elif surface.filldepth > 0 and i < surface.filldepth:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
+                                        else:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
+                        elif surface.ID == 'grass':
+                            g = surface.grass[0]
+                            # Build the blades of the grass
+                            blade = 0
+                            for j in range(surface.ys, surface.yf):
+                                for k in range(surface.zs, surface.zf):
+                                    if surface.fractalsurface[j - surface.ys, k - surface.zs] > 0:
+                                        height = 0
+                                        for i in range(volume.xs, surface.fractalrange[1]):
+                                            if i < surface.fractalsurface[j - surface.ys, k - surface.zs] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
+                                                y, z = g.calculate_blade_geometry(blade, height)
+                                                # Add y, z coordinates to existing location
+                                                yy = int(j - volume.ys + y)
+                                                zz = int(k - volume.zs + z)
+                                                # If these coordinates are outwith fractal volume stop building the blade, otherwise set the mask for grass
+                                                if yy < 0 or yy >= volume.mask.shape[1] or zz < 0 or zz >= volume.mask.shape[2]:
+                                                    break
+                                                else:
+                                                    volume.mask[i - volume.xs, yy, zz] = 3
+                                                    height += 1
+                                        blade += 1
+
+                            # Build the roots of the grass
+                            root = 0
+                            for j in range(surface.ys, surface.yf):
+                                for k in range(surface.zs, surface.zf):
+                                    if surface.fractalsurface[j - surface.ys, k - surface.zs] > 0:
+                                        depth = 0
+                                        i = volume.xf - 1
+                                        while i > volume.xs:
+                                            if i > volume.originalxf - (surface.fractalsurface[j - surface.ys, k - surface.zs] - volume.originalxf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
+                                                y, z = g.calculate_root_geometry(root, depth)
+                                                # Add y, z coordinates to existing location
+                                                yy = int(j - volume.ys + y)
+                                                zz = int(k - volume.zs + z)
+                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
+                                                if yy < 0 or yy >= volume.mask.shape[1] or zz < 0 or zz >= volume.mask.shape[2]:
+                                                    break
+                                                else:
+                                                    volume.mask[i - volume.xs, yy, zz] = 3
+                                                    depth += 1
+                                            i -= 1
+                                        root += 1
+
+                    elif surface.surfaceID == 'yminus':
+                        for i in range(surface.xs, surface.xf):
+                            for j in range(surface.fractalrange[0], surface.fractalrange[1]):
+                                for k in range(surface.zs, surface.zf):
+                                    if j > surface.fractalsurface[i - surface.xs, k - surface.zs]:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
+                                    elif surface.filldepth > 0 and j > surface.filldepth:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
+                                    else:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
+
+                    elif surface.surfaceID == 'yplus':
+                        if not surface.ID:
+                            for i in range(surface.xs, surface.xf):
+                                for j in range(surface.fractalrange[0], surface.fractalrange[1]):
+                                    for k in range(surface.zs, surface.zf):
+                                        if j < surface.fractalsurface[i - surface.xs, k - surface.zs]:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
+                                        elif surface.filldepth > 0 and j < surface.filldepth:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
+                                        else:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
+                        elif surface.ID == 'grass':
+                            g = surface.grass[0]
+                            # Build the blades of the grass
+                            blade = 0
+                            for i in range(surface.xs, surface.xf):
+                                for k in range(surface.zs, surface.zf):
+                                    if surface.fractalsurface[i - surface.xs, k - surface.zs] > 0:
+                                        height = 0
+                                        for j in range(volume.ys, surface.fractalrange[1]):
+                                            if j < surface.fractalsurface[i - surface.xs, k - surface.zs] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
+                                                x, z = g.calculate_blade_geometry(blade, height)
+                                                # Add x, z coordinates to existing location
+                                                xx = int(i - volume.xs + x)
+                                                zz = int(k - volume.zs + z)
+                                                # If these coordinates are outwith fractal volume stop building the blade, otherwise set the mask for grass
+                                                if xx < 0 or xx >= volume.mask.shape[0] or zz < 0 or zz >= volume.mask.shape[2]:
+                                                    break
+                                                else:
+                                                    volume.mask[xx, j - volume.ys, zz] = 3
+                                                    height += 1
+                                        blade += 1
+
+                            # Build the roots of the grass
+                            root = 0
+                            for i in range(surface.xs, surface.xf):
+                                for k in range(surface.zs, surface.zf):
+                                    if surface.fractalsurface[i - surface.xs, k - surface.zs] > 0:
+                                        depth = 0
+                                        j = volume.yf - 1
+                                        while j > volume.ys:
+                                            if j > volume.originalyf - (surface.fractalsurface[i - surface.xs, k - surface.zs] - volume.originalyf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
+                                                x, z = g.calculate_root_geometry(root, depth)
+                                                # Add x, z coordinates to existing location
+                                                xx = int(i - volume.xs + x)
+                                                zz = int(k - volume.zs + z)
+                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
+                                                if xx < 0 or xx >= volume.mask.shape[0] or zz < 0 or zz >= volume.mask.shape[2]:
+                                                    break
+                                                else:
+                                                    volume.mask[xx, j - volume.ys, zz] = 3
+                                                    depth += 1
+                                            j -= 1
+                                        root += 1
+
+                    elif surface.surfaceID == 'zminus':
+                        for i in range(surface.xs, surface.xf):
+                            for j in range(surface.ys, surface.yf):
+                                for k in range(surface.fractalrange[0], surface.fractalrange[1]):
+                                    if k > surface.fractalsurface[i - surface.xs, j - surface.ys]:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
+                                    elif surface.filldepth > 0 and k > surface.filldepth:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
+                                    else:
+                                        volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
+
+                    elif surface.surfaceID == 'zplus':
+                        if not surface.ID:
+                            for i in range(surface.xs, surface.xf):
+                                for j in range(surface.ys, surface.yf):
+                                    for k in range(surface.fractalrange[0], surface.fractalrange[1]):
+                                        if k < surface.fractalsurface[i - surface.xs, j - surface.ys]:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 1
+                                        elif surface.filldepth > 0 and k < surface.filldepth:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 2
+                                        else:
+                                            volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] = 0
+                        elif surface.ID == 'grass':
+                            g = surface.grass[0]
+                            # Build the blades of the grass
+                            blade = 0
+                            for i in range(surface.xs, surface.xf):
+                                for j in range(surface.ys, surface.yf):
+                                    if surface.fractalsurface[i - surface.xs, j - surface.ys] > 0:
+                                        height = 0
+                                        for k in range(volume.zs, surface.fractalrange[1]):
+                                            if k < surface.fractalsurface[i - surface.xs, j - surface.ys] and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] != 1:
+                                                x, y = g.calculate_blade_geometry(blade, height)
+                                                # Add x, y coordinates to existing location
+                                                xx = int(i - volume.xs + x)
+                                                yy = int(j - volume.ys + y)
+                                                # If these coordinates are outwith the fractal volume stop building the blade, otherwise set the mask for grass
+                                                if xx < 0 or xx >= volume.mask.shape[0] or yy < 0 or yy >= volume.mask.shape[1]:
+                                                    break
+                                                else:
+                                                    volume.mask[xx, yy, k - volume.zs] = 3
+                                                    height += 1
+                                        blade += 1
+
+                            # Build the roots of the grass
+                            root = 0
+                            for i in range(surface.xs, surface.xf):
+                                for j in range(surface.ys, surface.yf):
+                                    if surface.fractalsurface[i - surface.xs, j - surface.ys] > 0:
+                                        depth = 0
+                                        k = volume.zf - 1
+                                        while k > volume.zs:
+                                            if k > volume.originalzf - (surface.fractalsurface[i - surface.xs, j - surface.ys] - volume.originalzf) and volume.mask[i - volume.xs, j - volume.ys, k - volume.zs] == 1:
+                                                x, y = g.calculate_root_geometry(root, depth)
+                                                # Add x, y coordinates to existing location
+                                                xx = int(i - volume.xs + x)
+                                                yy = int(j - volume.ys + y)
+                                                # If these coordinates are outwith the fractal volume stop building the root, otherwise set the mask for grass
+                                                if xx < 0 or xx >= volume.mask.shape[0] or yy < 0 or yy >= volume.mask.shape[1]:
+                                                    break
+                                                else:
+                                                    volume.mask[xx, yy, k - volume.zs] = 3
+                                                    depth += 1
+                                            k -= 1
+                                        root += 1
+
+                # Build voxels from any true values of the 3D mask array
+                waternumID = next((x.numID for x in grid.materials if x.ID == 'water'), 0)
+                grassnumID = next((x.numID for x in grid.materials if x.ID == 'grass'), 0)
+                data = volume.fractalvolume.astype('int16', order='C')
+                mask = volume.mask.copy(order='C')
+                build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs, waternumID, grassnumID, volume.averaging, mask, data, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+
+            else:
+                if volume.nbins == 1:
+                    raise CmdInputError(self.__str__() + ' is being used with a single material and no modifications, therefore please use a #box command instead.')
+                else:
+                    volume.generate_fractal_volume(grid)
+                    volume.fractalvolume += volume.mixingmodel.startmaterialnum
+
+                data = volume.fractalvolume.astype('int16', order='C')
+                build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0, volume.averaging, data, grid.solid, grid.rigidE, grid.rigidH, grid.ID)

gprMax/cmds_geometry/geometry_objects_read.py

@@ -0,0 +1,89 @@
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+import os
+from tqdm import tqdm
+
+
+class GeometryObjectsRead(UserObjectGeometry):
+    """MORE WORK REQUIRED HERE."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        self.order = 1
+        self.hash = '#geometry_objects_read'
+
+    def create(self, G, uip):
+        """Create the object and add it to the grid."""
+
+        try:
+            p1 = self.kwargs['p1']
+            geofile = self.kwargs['geofile']
+            matfile = self.kwargs['matfile']
+        except:
+            raise CmdInputError(self.__str__() + 'requires exactly five parameters')
+
+        # discretise the point using uip object. This has different behaviour
+        # depending on the type of uip object. So we can use it for
+        # the main grid or the subgrid.
+        xs, ys, zs = uip.discretise_point(p1)
+
+        # See if material file exists at specified path and if not try input file directory
+        if not os.path.isfile(matfile):
+            matfile = os.path.abspath(os.path.join(G.inputdirectory, matfile))
+
+        matstr = os.path.splitext(os.path.split(matfile)[1])[0]
+        numexistmaterials = len(G.materials)
+
+        # Read materials from file
+        with open(matfile, 'r') as f:
+            # Strip out any newline characters and comments that must begin with double hashes
+            materials = [line.rstrip() + '{' + matstr + '}\n' for line in f if(not line.startswith('##') and line.rstrip('\n'))]
+
+        # Check validity of command names
+        singlecmdsimport, multicmdsimport, geometryimport = check_cmd_names(materials, checkessential=False)
+
+        # Process parameters for commands that can occur multiple times in the model
+        process_multicmds(multicmdsimport, G)
+
+        # Update material type
+        for material in G.materials:
+            if material.numID >= numexistmaterials:
+                if material.type:
+                    material.type += ',\nimported'
+                else:
+                    material.type = 'imported'
+
+        # See if geometry object file exists at specified path and if not try input file directory
+        if not os.path.isfile(geofile):
+            geofile = os.path.abspath(os.path.join(G.inputdirectory, geofile))
+
+        # Open geometry object file and read/check spatial resolution attribute
+        f = h5py.File(geofile, 'r')
+        dx_dy_dz = f.attrs['dx_dy_dz']
+        if round_value(dx_dy_dz[0] / G.dx) != 1 or round_value(dx_dy_dz[1] / G.dy) != 1 or round_value(dx_dy_dz[2] / G.dz) != 1:
+            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the spatial resolution of the geometry objects file to match the spatial resolution of the model')
+
+        data = f['/data'][:]
+
+        # Should be int16 to allow for -1 which indicates background, i.e.
+        # don't build anything, but AustinMan/Woman maybe uint16
+        if data.dtype != 'int16':
+            data = data.astype('int16')
+
+        # Look to see if rigid and ID arrays are present (these should be
+        # present if the original geometry objects were written from gprMax)
+        try:
+            rigidE = f['/rigidE'][:]
+            rigidH = f['/rigidH'][:]
+            ID = f['/ID'][:]
+            G.solid[xs:xs + data.shape[0], ys:ys + data.shape[1], zs:zs + data.shape[2]] = data + numexistmaterials
+            G.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE
+            G.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH
+            G.ID[:, xs:xs + ID.shape[1], ys:ys + ID.shape[2], zs:zs + ID.shape[3]] = ID + numexistmaterials
+            if G.messages:
+                tqdm.write('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))
+        except KeyError:
+            averaging = False
+            build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, data, G.solid, G.rigidE, G.rigidH, G.ID)
+            if G.messages:
+                tqdm.write('Geometry objects from file (voxels only) {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))

gprMax/cmds_geometry/plate.py

@@ -0,0 +1,112 @@
+"""Class for edge command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..geometry_primitives_ext import build_face_yz
+from ..geometry_primitives_ext import build_face_xz
+from ..geometry_primitives_ext import build_face_xy
+
+from tqdm import tqdm
+
+
+class Plate(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 3
+        self.hash = '#plate'
+
+    def create(self, grid, uip):
+
+        try:
+            p1 = self.kwargs['p1']
+            p2 = self.kwargs['p2']
+
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' 2 points must be specified')
+
+        # isotropic
+        try:
+            materialsrequested = [self.kwargs['material_id']]
+        except KeyError:
+            # Anisotropic case
+            try:
+                materialsrequested = self.kwargs['material_ids']
+            except KeyError:
+                raise CmdInputError(self.__str__() + ' No materials have been specified')
+
+        p1, p2 = uip.check_box_points(p1, p2, self.__str__())
+        xs, ys, zs = p1
+        xf, yf, zf = p2
+
+        # Check for valid orientations
+        if xs == xf:
+            if ys == yf or zs == zf:
+                raise CmdInputError(self.__str__() + ' the plate is not specified correctly')
+
+        elif ys == yf:
+            if xs == xf or zs == zf:
+                raise CmdInputError(self.__str__() + ' the plate is not specified correctly')
+
+        elif zs == zf:
+            if xs == xf or ys == yf:
+                raise CmdInputError(self.__str__() + ' the plate is not specified correctly')
+
+        else:
+            raise CmdInputError(self.__str__() + ' the plate is not specified correctly')
+
+        # Look up requested materials in existing list of material instances
+        materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+        if len(materials) != len(materialsrequested):
+            notfound = [x for x in materialsrequested if x not in materials]
+            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+
+        # yz-plane plate
+        if xs == xf:
+            # Isotropic case
+            if len(materials) == 1:
+                numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 2:
+                numIDy = materials[0].numID
+                numIDz = materials[1].numID
+
+            for j in range(ys, yf):
+                for k in range(zs, zf):
+                    build_face_yz(xs, j, k, numIDy, numIDz, grid.rigidE, grid.rigidH, grid.ID)
+
+        # xz-plane plate
+        elif ys == yf:
+            # Isotropic case
+            if len(materials) == 1:
+                numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 2:
+                numIDx = materials[0].numID
+                numIDz = materials[1].numID
+
+            for i in range(xs, xf):
+                for k in range(zs, zf):
+                    build_face_xz(i, ys, k, numIDx, numIDz, grid.rigidE, grid.rigidH, grid.ID)
+
+        # xy-plane plate
+        elif zs == zf:
+            # Isotropic case
+            if len(materials) == 1:
+                numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 2:
+                numIDx = materials[0].numID
+                numIDy = materials[1].numID
+
+            for i in range(xs, xf):
+                for j in range(ys, yf):
+                    build_face_xy(i, j, zs, numIDx, numIDy, grid.rigidE, grid.rigidH, grid.ID)
+
+        if grid.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

@@ -0,0 +1,91 @@
+"""Class for sphere command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..materials import Material
+from ..geometry_primitives_ext import build_sphere
+
+from tqdm import tqdm
+import numpy as np
+
+
+class Sphere(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 8
+        self.hash = '#sphere'
+
+    def create(self, grid, uip):
+        try:
+            p1 = self.kwargs['p1']
+            r = self.kwargs['r']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' Please specify a point and a radius.')
+
+        # check averaging
+        try:
+            # go with user specified averaging
+            averagesphere = self.kwargs['averaging']
+        except KeyError:
+            # if they havent specfied - go with the grid default
+            averagesphere = grid.averagevolumeobjects
+
+        # check materials have been specified
+        # isotropic case
+        try:
+            materialsrequested = [self.kwargs['material_id']]
+        except KeyError:
+            # Anisotropic case
+            try:
+                materialsrequested = self.kwargs['material_ids']
+            except KeyError:
+                raise CmdInputError(self.__str__() + ' No materials have been specified')
+
+        # Centre of sphere
+        xc, yc, zc = uip.round_to_grid(p1)
+
+        # Look up requested materials in existing list of material instances
+        materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+        if len(materials) != len(materialsrequested):
+            notfound = [x for x in materialsrequested if x not in materials]
+            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+
+        # Isotropic case
+        if len(materials) == 1:
+            averaging = materials[0].averagable and averagesphere
+            numID = numIDx = numIDy = numIDz = materials[0].numID
+
+        # Uniaxial anisotropic case
+        elif len(materials) == 3:
+            averaging = False
+            numIDx = materials[0].numID
+            numIDy = materials[1].numID
+            numIDz = materials[2].numID
+            requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+            averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+            if averagedmaterial:
+                numID = averagedmaterial.numID
+            else:
+                numID = len(grid.materials)
+                m = Material(numID, requiredID)
+                m.type = 'dielectric-smoothed'
+                # Create dielectric-smoothed constituents for material
+                m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
+                m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
+                m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+
+                # Append the new material object to the materials list
+                grid.materials.append(m)
+
+        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 averaging:
+                dielectricsmoothing = 'on'
+            else:
+                dielectricsmoothing = 'off'
+            tqdm.write('Sphere with centre {:g}m, {:g}m, {:g}m, radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(xc * grid.dx, yc * grid.dy, zc * grid.dz, r, ', '.join(materialsrequested), dielectricsmoothing))

gprMax/cmds_geometry/triangle.py

@@ -0,0 +1,132 @@
+"""Class for triangle command."""
+from .cmds_geometry import UserObjectGeometry
+from ..exceptions import CmdInputError
+from ..materials import Material
+from ..geometry_primitives_ext import build_triangle
+
+from tqdm import tqdm
+import numpy as np
+
+
+class Triangle(UserObjectGeometry):
+    """User class for edge command."""
+
+    def __init__(self, **kwargs):
+        """Constructor."""
+        super().__init__(**kwargs)
+        self.order = 4
+        self.hash = '#triangle'
+
+    def create(self, grid, uip):
+
+        try:
+            up1 = self.kwargs['p1']
+            up2 = self.kwargs['p2']
+            up3 = self.kwargs['p3']
+            thickness = self.kwargs['thickness']
+
+        except KeyError:
+            raise CmdInputError(self.params_str() + ' Specify 3 points and a thickness')
+
+        # check averaging
+        try:
+            # go with user specified averaging
+            averagetriangularprism = self.kwargs['averaging']
+        except KeyError:
+            # if they havent specfied - go with the grid default
+            averagetriangularprism = grid.averagevolumeobjects
+
+        # check materials have been specified
+        # isotropic case
+        try:
+            materialsrequested = [self.kwargs['material_id']]
+        except KeyError:
+            # Anisotropic case
+            try:
+                materialsrequested = self.kwargs['material_ids']
+            except KeyError:
+                raise CmdInputError(self.__str__() + ' No materials have been specified')
+
+        # Check whether points are valid against grid
+        uip.check_tri_points(up1, up2, up3, object)
+        # Convert points to metres
+        x1, y1, z1 = uip.round_to_grid(up1)
+        x2, y2, z2 = uip.round_to_grid(up2)
+        x3, y3, z3 = uip.round_to_grid(up3)
+
+        if thickness < 0:
+            raise CmdInputError(self.__str__() + ' requires a positive value for thickness')
+
+        # Check for valid orientations
+        # yz-plane triangle
+        if x1 == x2 and x2 == x3:
+            normal = 'x'
+        # xz-plane triangle
+        elif y1 == y2 and y2 == y3:
+            normal = 'y'
+        # xy-plane triangle
+        elif z1 == z2 and z2 == z3:
+            normal = 'z'
+        else:
+            raise CmdInputError(self.__str__() + ' the triangle is not specified correctly')
+
+        # Look up requested materials in existing list of material instances
+        materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+        if len(materials) != len(materialsrequested):
+            notfound = [x for x in materialsrequested if x not in materials]
+            raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound))
+
+        if thickness > 0:
+            # Isotropic case
+            if len(materials) == 1:
+                averaging = materials[0].averagable and averagetriangularprism
+                numID = numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 3:
+                averaging = False
+                numIDx = materials[0].numID
+                numIDy = materials[1].numID
+                numIDz = materials[2].numID
+                requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+                averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+                if averagedmaterial:
+                    numID = averagedmaterial.numID
+                else:
+                    numID = len(grid.materials)
+                    m = Material(numID, requiredID)
+                    m.type = 'dielectric-smoothed'
+                    # Create dielectric-smoothed constituents for material
+                    m.er = np.mean((materials[0].er, materials[1].er, materials[2].er), axis=0)
+                    m.se = np.mean((materials[0].se, materials[1].se, materials[2].se), axis=0)
+                    m.mr = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+                    m.sm = np.mean((materials[0].mr, materials[1].mr, materials[2].mr), axis=0)
+
+                    # Append the new material object to the materials list
+                    grid.materials.append(m)
+        else:
+            averaging = False
+            # Isotropic case
+            if len(materials) == 1:
+                numID = numIDx = numIDy = numIDz = materials[0].numID
+
+            # Uniaxial anisotropic case
+            elif len(materials) == 3:
+                # numID requires a value but it will not be used
+                numID = None
+                numIDx = materials[0].numID
+                numIDy = materials[1].numID
+                numIDz = materials[2].numID
+
+        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 thickness > 0:
+                if averaging:
+                    dielectricsmoothing = 'on'
+                else:
+                    dielectricsmoothing = 'off'
+                tqdm.write('Triangle with coordinates {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m and thickness {:g}m of material(s) {} created, dielectric smoothing is {}.'.format(x1, y1, z1, x2, y2, z2, x3, y3, z3, thickness, ', '.join(materialsrequested), dielectricsmoothing))
+            else:
+                tqdm.write('Triangle with coordinates {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m of material(s) {} created.'.format(x1, y1, z1, x2, y2, z2, x3, y3, z3, ', '.join(materialsrequested)))

gprMax/cmds_single_use.py

@@ -0,0 +1,507 @@
+import os
+import sys
+import inspect
+
+from .constants import c
+from .constants import floattype
+from .exceptions import CmdInputError
+from .waveforms import Waveform
+
+import numpy as np
+from scipy import interpolate
+
+from colorama import init
+from colorama import Fore
+from colorama import Style
+init()
+
+
+class UserObjectSingle:
+
+    """
+        Specific geometry object
+    """
+
+    def __init__(self, **kwargs):
+        # each single command has an order to specify the order in which
+        # the commands are constructed. IE. discretisation must be
+        # created before the domain
+        self.order = None
+        self.kwargs = kwargs
+
+    def __str__(self):
+        pass
+
+    def create(self, grid, uip):
+        pass
+
+
+class DomainSingle(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 2
+
+    def __str__(self):
+        try:
+            s = '#domain: {} {} {}'.format(self.kwargs['p1'][0],
+                                           self.kwargs['p1'][1],
+                                           self.kwargs['p1'][2])
+        except KeyError:
+            print('error message')
+
+        return s
+
+    def create(self, G, uip):
+
+        # code to create the gprMax domain as per input_cmds_singleuse.py
+        try:
+            G.nx, G.ny, G.nz = uip.discretise_point(self.kwargs['p1'])
+        except KeyError:
+            raise CmdInputError(self.__str__ + ' Please specify a point')
+
+        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 G.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
+        if G.nx == 1:
+            G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
+            G.mode = '2D TMx'
+            G.pmlthickness['x0'] = 0
+            G.pmlthickness['xmax'] = 0
+        elif G.ny == 1:
+            G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dz) * (1 / G.dz)))
+            G.mode = '2D TMy'
+            G.pmlthickness['y0'] = 0
+            G.pmlthickness['ymax'] = 0
+        elif G.nz == 1:
+            G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy)))
+            G.mode = '2D TMz'
+            G.pmlthickness['z0'] = 0
+            G.pmlthickness['zmax'] = 0
+        else:
+            G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
+            G.mode = '3D'
+
+        # Round down time step to nearest float with precision one less than hardware maximum. Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
+        G.round_time_step()
+
+        if G.messages:
+            print('Mode: {}'.format(G.mode))
+            print('Time step (at CFL limit): {:g} secs'.format(G.dt))
+
+        # do threads here
+        from .utilities import get_host_info
+        hostinfo = get_host_info()
+
+        # Number of threads (OpenMP) to use
+        if sys.platform == 'darwin':
+            os.environ['OMP_WAIT_POLICY'] = 'ACTIVE'  # Should waiting threads consume CPU power (can drastically effect performance)
+        os.environ['OMP_DYNAMIC'] = 'FALSE'  # Number of threads may be adjusted by the run time environment to best utilize system resources
+        os.environ['OMP_PLACES'] = 'cores'  # Each place corresponds to a single core (having one or more hardware threads)
+        os.environ['OMP_PROC_BIND'] = 'TRUE'  # Bind threads to physical cores
+        # os.environ['OMP_DISPLAY_ENV'] = 'TRUE' # Prints OMP version and environment variables (useful for debug)
+
+        # Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785)
+        if 'Microsoft' in G.hostinfo['osversion']:
+            os.environ['KMP_AFFINITY'] = 'disabled'
+            del os.environ['OMP_PLACES']
+            del os.environ['OMP_PROC_BIND']
+
+        if os.environ.get('OMP_NUM_THREADS'):
+            G.nthreads = int(os.environ.get('OMP_NUM_THREADS'))
+        else:
+            # Set number of threads to number of physical CPU cores
+            G.nthreads = hostinfo['physicalcores']
+            os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
+
+        if G.messages:
+            print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
+        if G.nthreads > G.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)
+
+
+
+class Domain:
+
+    """Restrict user object so there can only be one instance
+    https://python-3-patterns-idioms-test.readthedocs.io/en/latest/index.html
+    """
+
+    instance = None
+
+    def __new__(cls, **kwargs):  # __new__ always a classmethod
+        if not Domain.instance:
+            Domain.instance = DomainSingle(**kwargs)
+        return Domain.instance
+
+
+class Discretisation(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 1
+
+    def __str__(self):
+        try:
+            s = '#dx_dy_dz: {} {} {}'.format(self.kwargs['p1'][0],
+                                             self.kwargs['p1'][1],
+                                             self.kwargs['p1'][2])
+        except KeyError:
+            print('error message')
+
+        return s
+
+    def create(self, G, uip):
+
+        try:
+            G.dx, G.dy, G.dz = self.kwargs['p1']
+
+        except KeyError:
+            raise CmdInputError('Discretisation requires a point')
+
+        if G.dx <= 0:
+            raise CmdInputError('Discretisation requires the x-direction spatial step to be greater than zero')
+        if G.dy <= 0:
+            raise CmdInputError(' Discretisation requires the y-direction spatial step to be greater than zero')
+        if G.dz <= 0:
+            raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero')
+
+        if G.messages:
+            print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz))
+
+
+class TimeWindow(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 4
+
+    def __str__(self):
+        try:
+            s = '#time_window: {}'.format(self.kwargs['time'])
+        except KeyError:
+            try:
+                s = '#time_window: {}'.format(self.kwargs['iterations'])
+            except KeyError:
+                print('time window error')
+
+        return s
+
+    def create(self, G, uip):
+        # If number of iterations given
+        # The +/- 1 used in calculating the number of iterations is to account for
+        # the fact that the solver (iterations) loop runs from 0 to < G.iterations
+        try:
+            iterations = int(self.kwargs['iterations'])
+            G.timewindow = (iterations - 1) * G.dt
+            G.iterations = iterations
+
+        except KeyError:
+            pass
+
+        try:
+            tmp = float(self.kwargs['time'])
+            if tmp > 0:
+                G.timewindow = tmp
+                G.iterations = int(np.ceil(tmp / G.dt)) + 1
+            else:
+                raise CmdInputError(self.__str__() + ' must have a value greater than zero')
+
+        except KeyError:
+            pass
+
+        if not G.timewindow:
+            raise CmdInputError('TimeWindow: Specify a time or number of iterations')
+
+        if G.messages:
+            print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
+
+
+class Messages(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 0
+
+    def __str__(self):
+        try:
+            s = '#messages: {}'.format(self.kwargs['yn'])
+        except KeyError:
+            print('messages problem')
+
+    def create(self, G, uip):
+        try:
+            yn = self.kwargs['yn']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' requires exactly one parameter')
+
+        if yn.lower() == 'y':
+            G.messages = True
+        elif yn.lower() == 'n':
+            G.messages = False
+        else:
+            raise CmdInputError(self.__str__() + ' requires input values of either y or n')
+
+
+class Title(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 5
+
+    def create(self, G, uip):
+        # Title
+        try:
+            title = self.kwargs['name']
+            G.title = title
+        except KeyError:
+            pass
+
+        if G.messages:
+            print('Model title: {}'.format(G.title))
+
+class NumThreads(UserObjectSingle):
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 6
+
+    def __str__(self):
+        try:
+            return '#n_threads: {}'.format(self.kwargs['n'])
+        except KeyError:
+            return '#n_threads:'
+
+    def create(self, G, uip):
+        # Get information about host machine
+        from .utilities import get_host_info
+        hostinfo = get_host_info()
+
+        try:
+            n = self.kwargs['n']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' requires exactly one parameter to specify the number of threads to use')
+
+        if n < 1:
+            raise CmdInputError(self.__str__() + ' requires the value to be an integer not less than one')
+
+        G.nthreads = n
+        os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
+
+        if G.messages:
+            print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
+        if G.nthreads > G.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 G.messages:
+            if G.gpu is not None:
+                print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))
+
+
+# Time step stability factor
+class TimeStepStabilityFactor(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 7
+
+    def __str__(self):
+
+        try:
+            return '#time_step_stability_factor: {}'.format(self.kwargs['f'])
+        except KeyError:
+            return '#time_step_stability_factor:'
+
+    def create(self, G, uip):
+
+        try:
+            f = self.kwargs['f']
+        except KeyError:
+            raise CmdInputError(self.__str__() + ' requires exactly one parameter')
+
+        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 G.messages:
+            print('Time step (modified): {:g} secs'.format(G.dt))
+
+
+class PMLCells(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 8
+
+    def create(self, G, uip):
+
+        try:
+            thickness = self.kwargs['thickness']
+
+            for key in G.pmlthickness.keys():
+                G.pmlthickness[key] = int(thickness)
+
+        except KeyError:
+            try:
+                G.pmlthickness['x0'] = int(self.kwargs['x0'])
+                G.pmlthickness['y0'] = int(self.kwargs['y0'])
+                G.pmlthickness['z0'] = int(self.kwargs['z0'])
+                G.pmlthickness['xmax'] = int(self.kwargs['xmax'])
+                G.pmlthickness['ymax'] = int(self.kwargs['ymax'])
+                G.pmlthickness['zmax'] = int(self.kwargs['zmax'])
+            except KeyError:
+                raise CmdInputError('#pml_cells: requires either one or six parameter(s)')
+
+        if (2 * G.pmlthickness['x0'] >= G.nx or
+            2 * G.pmlthickness['y0'] >= G.ny or
+            2 * G.pmlthickness['z0'] >= G.nz or
+            2 * G.pmlthickness['xmax'] >= G.nx or
+            2 * G.pmlthickness['ymax'] >= G.ny or
+            2 * G.pmlthickness['zmax'] >= G.nz):
+                raise CmdInputError('#pml_thickness: has too many cells for the domain size')
+
+
+class SrcSteps(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 9
+
+    def create(self, G, uip):
+        try:
+            G.srcsteps = uip.discretise_point(self.kwargs['p1'])
+        except KeyError:
+            raise CmdInputError('#src_steps: requires exactly three parameters')
+        # src_steps
+        if G.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))
+
+
+class RxSteps(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 10
+
+    def create(self, G, uip):
+        try:
+            G.rxsteps = uip.discretise_point(self.kwargs['p1'])
+        except KeyError:
+            raise CmdInputError('#rx_steps: requires exactly three parameters')
+        if G.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))
+
+
+class ExcitationFile(UserObjectSingle):
+
+    def create(self, G, uip):
+    # Excitation file for user-defined source waveforms
+        try:
+            kwargs = dict()
+            excitationfile = self.kwargs['filepath']
+            kwargs['kind'] = self.kwargs['kind']
+            kwargs['fill_value'] = self.kwargs['fill_value']
+
+        except KeyError:
+            try:
+                excitationfile = self.kwargs['filepath']
+                args, varargs, keywords, defaults = inspect.getargspec(interpolate.interp1d)
+                kwargs = dict(zip(reversed(args), reversed(defaults)))
+            except KeyError:
+                raise CmdInputError('#excitation_file: requires either one or three parameter(s)')
+
+            # See if file exists at specified path and if not try input file directory
+            if not os.path.isfile(excitationfile):
+                excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
+
+            if G.messages:
+                print('\nExcitation file: {}'.format(excitationfile))
+
+            # Get waveform names
+            with open(excitationfile, 'r') as f:
+                waveformIDs = f.readline().split()
+
+            # Read all waveform values into an array
+            waveformvalues = np.loadtxt(excitationfile, skiprows=1, dtype=floattype)
+
+            # Time array (if specified) for interpolation, otherwise use simulation time
+            if waveformIDs[0].lower() == 'time':
+                waveformIDs = waveformIDs[1:]
+                waveformtime = waveformvalues[:, 0]
+                waveformvalues = waveformvalues[:, 1:]
+                timestr = 'user-defined time array'
+            else:
+                waveformtime = np.arange(0, G.timewindow + G.dt, G.dt)
+                timestr = 'simulation time array'
+
+            for waveform in range(len(waveformIDs)):
+                if any(x.ID == waveformIDs[waveform] for x in G.waveforms):
+                    raise CmdInputError('Waveform with ID {} already exists'.format(waveformIDs[waveform]))
+                w = Waveform()
+                w.ID = waveformIDs[waveform]
+                w.type = 'user'
+
+                # Select correct column of waveform values depending on array shape
+                singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, waveform]
+
+                # Truncate waveform array if it is longer than time array
+                if len(singlewaveformvalues) > len(waveformtime):
+                    singlewaveformvalues = singlewaveformvalues[:len(waveformtime)]
+                # Zero-pad end of waveform array if it is shorter than time array
+                elif len(singlewaveformvalues) < len(waveformtime):
+                    singlewaveformvalues = np.lib.pad(singlewaveformvalues, (0, len(singlewaveformvalues) - len(waveformvalues)), 'constant', constant_values=0)
+
+                # Interpolate waveform values
+                w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)
+
+                if G.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)
+
+
+class OutputDir(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        # dont need to define parameters in advance. Just catch errors
+        # when they occur
+        super().__init__(**kwargs)
+        self.order = 11
+
+    def create(self, grid, uip):
+        grid.outputdirectory = self.kwargs['dir']
+
+
+class NumberOfModelRuns(UserObjectSingle):
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.order = 12
+
+    def create(self, grid, uip):
+        try:
+            grid.numberofmodelruns = self.kwargs['n']
+        except KeyError:
+            raise CmdInputError('#numberofmodelruns: requires exactly one parameter')

gprMax/cython/fields_updates_hsg.pyx

@@ -0,0 +1,269 @@
+from __future__ import print_function
+from gprMax.constants cimport floattype_t
+from gprMax.constants cimport complextype_t
+
+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,
+        int face,
+        int l_l,
+        int l_u,
+        int m_l,
+        int m_u,
+        size_t n_l,
+        size_t n_u,
+        int nwn,
+        size_t lookup_id,
+        np.float64_t[:, :, :] field,
+        np.float64_t[:, :, :] inc_field,
+        size_t co,
+        int sign_n,
+        int sign_f,
+        int mid,
+        int r,
+        int s,
+        int nb,
+        int nthreads
+):
+        """
+            Args:
+
+            subgrid: (Subgrid)
+            n: (String) the normal to the face to update
+            nwn: (Int) number of working cell in the normal direction
+            to the face
+            lookup_id: (Int) id of the H component we wish to update at
+            each node
+            field: (Numpy array) main grid field to be updated
+            inc_field: (Numpy array) incident sub_grid field
+            co: (Int) Coefficient used by gprMax update equations which
+            is specific to the field component being updated.
+            sign_n: (Int) 1 or -1 sign of the incident field on the near face.
+            sign_f: (Int) 1 or -1 sign of the incident field on the far face.
+            mid: (Bool) is the H node midway along the lower edge?
+            r = self.ratio
+            s = self.is_os_sep
+            nb = self.n_boundary_cells
+        """
+        # Comments here as as per left and right face
+
+        cdef Py_ssize_t l, m, l_s, m_s, n_s_l, n_s_r, material_e_l, material_e_r, i0, j0, k0, i1, j1, k1, i2, j2, k2, i3, j3, k3
+        cdef int os
+        cdef double inc_n, inc_f
+
+        # surface normal index for the subgrid near face h nodes (left i index)
+        n_s_l = nb - s * r - r + r // 2
+        # surface normal index for the subgrid far face h nodes (right i index)
+        n_s_r = nb + nwn + s * r + r // 2
+        # OS at the left face
+        os = nb - r * s
+
+        # Iterate over a slice of the main grid using dummy indices
+        for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
+
+            # Calculate the subgrid j component of the H nodes
+            # i.e. Hz node of the left or right face
+            if mid == 1:
+                l_s = os + (l - l_l) * r + r // 2
+            # i.e. the Hy node of the left or right face
+            else:
+                l_s = os + (l - l_l) * r
+
+            for m in range(m_l, m_u):
+
+                # Calculate the subgrid k component of the H nodes
+                if mid == 1:
+                    m_s = os + (m - m_l) * r
+                else:
+                    m_s = os + (m - m_l) * r + r // 2
+
+                # left and right
+                if face == 2:
+                    # main grid index
+                    i0, j0, k0 = n_l, l, m
+                    # equivalent subgrid index
+                    i1, j1, k1 = n_s_l, l_s, m_s
+                    i2, j2, k2 = n_u, l, m
+                    i3, j3, k3 = n_s_r, l_s, m_s
+                # front and back
+                if face == 3:
+                    i0, j0, k0 = l, n_l, m
+                    i1, j1, k1 = l_s, n_s_l, m_s
+                    i2, j2, k2 = l, n_u, m
+                    i3, j3, k3 = l_s, n_s_r, m_s
+                # top bottom
+                if face == 1:
+                    i0, j0, k0 = l, m, n_l
+                    i1, j1, k1 = l_s, m_s, n_s_l
+                    i2, j2, k2 = l, m, n_u
+                    i3, j3, k3 = l_s, m_s, n_s_r
+                # Update the left face
+
+                # Get the material at main grid index
+                material_e_l = ID[lookup_id, i0, j0, k0]
+                # Get the associated indident field from the subgrid
+                inc_n = inc_field[i1, j1, k1] * sign_n
+                # Update the main grid E field with the corrected H field
+                field[i0, j0, k0] += updatecoeffsE[material_e_l, co] * inc_n
+
+                # Update the right face
+                material_e_r = ID[lookup_id, i2, j2, k2]
+                inc_f = inc_field[i3, j3, k3] * sign_f
+                field[i2, j2, k2] += updatecoeffsE[material_e_r, co] * inc_f
+
+cpdef void cython_update_magnetic_os(
+        np.float64_t[:, :] updatecoeffsH,
+        np.uint32_t[:, :, :, :] ID,
+        int face,
+        int l_l,
+        int l_u,
+        int m_l,
+        int m_u,
+        size_t n_l,
+        size_t n_u,
+        int nwn,
+        size_t lookup_id,
+        np.float64_t[:, :, :] field,
+        np.float64_t[:, :, :] inc_field,
+        size_t co,
+        int sign_n,
+        int sign_f,
+        int mid,
+        int r,
+        int s,
+        int nb,
+        int nthreads
+):
+        """
+        int r ratio,
+        int s is_os_sep,
+        int nb n_boundary_cells
+        """
+
+        cdef Py_ssize_t l, m, l_s, m_s, n_s_l, n_s_r, material_e_l, material_e_r, i0, j0, k0, i1, j1, k1, i2, j2, k2, i3, j3, k3
+        cdef int os
+        cdef double inc_n, inc_f
+
+        # i index (normal to os) for the subgrid near face e node
+        n_s_l = nb - r * s
+        # Normal index for the subgrid far face e node
+        n_s_r = nb + nwn + s * r
+
+        # os inner index for the sub grid
+        os = nb - r * s
+
+        for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
+
+            # y coord of the Ex field component
+            if mid == 1:
+                l_s = os + (l - l_l) * r + r // 2
+            # y coord of the Ez field component
+            else:
+                l_s = os + (l - l_l) * r
+
+            for m in range(m_l, m_u):
+
+                # z coordinate of the Ex node in the subgrid
+                if mid == 1:
+                    m_s = os + (m - m_l) * r
+                else:
+                    m_s = os + (m - m_l) * r + r // 2
+
+                # associate the given indices with their i, j, k values
+
+                # left and right
+                if face == 2:
+                    # main grid index
+                    i0, j0, k0 = n_l, l, m
+                    # equivalent subgrid index
+                    i1, j1, k1 = n_s_l, l_s, m_s
+                    i2, j2, k2 = n_u, l, m
+                    i3, j3, k3 = n_s_r, l_s, m_s
+                # front and back
+                if face == 3:
+                    i0, j0, k0 = l, n_l, m
+                    i1, j1, k1 = l_s, n_s_l, m_s
+                    i2, j2, k2 = l, n_u, m
+                    i3, j3, k3 = l_s, n_s_r, m_s
+                # top bottom
+                if face == 1:
+                    i0, j0, k0 = l, m, n_l
+                    i1, j1, k1 = l_s, m_s, n_s_l
+                    i2, j2, k2 = l, m, n_u
+                    i3, j3, k3 = l_s, m_s, n_s_r
+
+                material_e_l = ID[lookup_id, i0, j0, k0]
+                inc_n = inc_field[i1, j1, k1] * sign_n
+
+                # make sure these are the correct grid
+                field[i0, j0, k0] += updatecoeffsH[material_e_l, co] * inc_n
+
+                # Far face
+                material_e_r = ID[lookup_id, i2, j2, k2]
+                inc_f = inc_field[i3, j3, k3] * sign_f
+                field[i2, j2, k2] += updatecoeffsH[material_e_r, co] * inc_f
+
+cpdef void cython_update_is(
+        int nwx,
+        int nwy,
+        int nwz,
+        floattype_t[:, :] updatecoeffsE,
+        np.uint32_t[:, :, :, :] ID,
+        int n,
+        int offset,
+        int nwl,
+        int nwm,
+        int nwn,
+        int face,
+        floattype_t[:, :, :] field,
+        floattype_t[:, :] inc_field_l,
+        floattype_t[:, :] inc_field_u,
+        Py_ssize_t lookup_id,
+        int sign_l,
+        int sign_u,
+        Py_ssize_t co,
+        int nthreads
+):
+
+        cdef Py_ssize_t l, m, i1, j1, k1, i2, j2, k2, field_material_l, field_material_u, inc_i, inc_j
+        cdef double inc_l, inc_u, f_l, f_u
+        # for inner faces H nodes are 1 cell before n boundary cells
+        cdef int n_o = n + offset
+
+        for l in prange(n, nwl + n, nogil=True, schedule='static', num_threads=nthreads):
+            for m in range(n, nwm + n):
+
+                # bottom and top
+                if face == 1:
+                    i1, j1, k1 = l, m, n_o
+                    i2, j2, k2 = l, m, n + nwz
+                # left and right
+                if face == 2:
+                    i1, j1, k1 = n_o, l, m
+                    i2, j2, k2 = n + nwx, l, m
+                # front and back
+                if face == 3:
+                    i1, j1, k1 = l, n_o, m
+                    i2, j2, k2 = l, n + nwy, m
+
+                inc_i = l - n
+                inc_j = m - n
+
+                field_material_l = ID[lookup_id, i1, j1, k1]
+                inc_l = inc_field_l[inc_i, inc_j]
+                # Additional field at i, j, k
+                f_l = updatecoeffsE[field_material_l, co] * inc_l * sign_l
+                # Set the new value
+                field[i1, j1, k1] += f_l
+
+                field_material_u = ID[lookup_id, i2, j2, k2]
+                inc_u = inc_field_u[inc_i, inc_j]
+                # Additional field at i, j, k
+                f_u = updatecoeffsE[field_material_u, co] * inc_u * sign_u
+                # Set the new value
+                field[i2, j2, k2] += f_u

gprMax/scene.py

@@ -0,0 +1,122 @@
+from .user_inputs import create_user_input_points
+from .materials import create_built_in_materials
+from .cmds_single_use import UserObjectSingle
+from .cmds_single_use import DomainSingle
+from .cmds_single_use import Discretisation
+from .cmds_single_use import TimeWindow
+from .cmds_multiple import UserObjectMulti
+from .subgrids.user_objects import SubGridUserBase
+from .cmds_geometry.cmds_geometry import UserObjectGeometry
+from .exceptions import CmdInputError
+from .cmds_geometry.fractal_box_builder import FractalBoxBuilder
+from .utilities import human_size
+from .grid import memory_check
+from .snapshots import memory_check as snapshot_memory_check
+
+
+class Scene:
+
+    """Scene stores all of the user created objects
+    """
+
+    def __init__(self):
+        self.multiple_cmds = []
+        self.single_cmds = []
+        self.geometry_cmds = []
+        self.essential_cmds = [DomainSingle, TimeWindow, Discretisation]
+
+        # fractal box commands have an additional nonuser object which
+        # process modifications
+        fbb = FractalBoxBuilder()
+        self.add(fbb)
+
+    def add(self, node):
+        if isinstance(node, UserObjectMulti):
+            self.children_multiple.append(node)
+        elif isinstance(node, UserObjectGeometry):
+            self.children_geometry.append(node)
+        elif isinstance(node, UserObjectSingle):
+            self.children_single.append(node)
+        else:
+            raise Exception('This Object is Unknown to gprMax')
+
+    def process_subgrid_commands(self, subgrids):
+
+        # check for subgrid user objects
+        def func(obj):
+            if isinstance(obj, SubGridUserBase):
+                return True
+            else:
+                return False
+
+        subgrid_cmds = list(filter(func, self.children_multiple))
+
+        # iterate through the user command objects under the subgrid user object
+        for sg_cmd in subgrid_cmds:
+            # when the subgrid is created its reference is attached to its user
+            # object. this reference allows the multi and geo user objects
+            # to build in the correct subgrid.
+            sg = sg_cmd.subgrid
+            self.process_cmds(sg_cmd.children_multiple, sg)
+            self.process_cmds(sg_cmd.children_geometry, sg)
+
+    def process_cmds(self, commands, grid):
+        cmds_sorted = sorted(commands, key=lambda cmd: cmd.order)
+        for obj in cmds_sorted:
+            # in the first level all objects belong to the main grid
+            uip = create_user_input_points(grid)
+            # Create an instance to check the geometry points provided by the
+            # user. The way the point are checked depends on which grid the
+            # points belong to.
+            obj.create(grid, uip)
+
+        return self
+
+    def process_singlecmds(self, G):
+
+        # check for duplicate commands and warn user if they exist
+        cmds_unique = list(set(self.children_single))
+        if len(cmds_unique) != len(self.children_single):
+            raise CmdInputError('Duplicate Single Commands exist in the input.')
+
+        # check essential cmds and warn user if missing
+        for cmd_type in self.essential_cmds:
+            d = any([isinstance(cmd, cmd_type) for cmd in cmds_unique])
+            if not d:
+                raise CmdInputError('Your input file is missing essential commands required to run a model. Essential commands are: Domain, Discretisation, Time Window')
+
+        self.process_cmds(cmds_unique, G)
+
+
+    def create_internal_objects(self, G):
+
+        # gprMax API presents the user with UserObjects in order to build
+        # the internal Rx(), Cylinder() etc... objects. This function
+        # essentially calls the UserObject.create() function in the correct
+        # way
+
+        # Traverse all the user objects in the correct order and create them.
+        create_built_in_materials(G)
+
+        # process commands that can onlyhave a single instance
+        self.process_singlecmds(G)
+
+        # Process main grid multiple commands
+        self.process_cmds(self.children_multiple, G)
+
+        # Estimate and check memory (RAM) usage
+        memory_check(G)
+        snapshot_memory_check(G)
+
+        # Initialise an array for volumetric material IDs (solid), boolean
+        # arrays for specifying materials not to be averaged (rigid),
+        # an array for cell edge IDs (ID)
+        G.initialise_grids()
+
+        # Process the main grid geometry commands
+        self.process_cmds(self.children_geometry, G)
+
+        # Process all the commands for the subgrid
+        self.process_subgrid_commands(G.subgrids)
+
+        return self

gprMax/subgrids/base.py

@@ -0,0 +1,61 @@
+from ..grid import FDTDGrid
+
+from scipy.constants import c
+import numpy as np
+
+
+class SubGridBase(FDTDGrid):
+
+    def __init__(self, **kwargs):
+        super().__init__()
+
+        self.mode = '3D'
+
+        # subgridding ratio
+        self.ratio = kwargs['ratio']
+
+        if self.ratio % 2 == 0:
+            raise ValueError('Subgrid Error: Only odd ratios are supported')
+
+        # Name of the grid
+        self.name = kwargs['ID']
+
+        self.filter = kwargs['filter']
+
+        # Number of main grid cells between the IS and OS
+        self.is_os_sep = kwargs['is_os_sep']
+        # Number of subgrid grid cells between the IS and OS
+        self.s_is_os_sep = self.is_os_sep * self.ratio
+
+        # Distance from OS to pml or the edge of the grid when pml is off
+        self.pml_separation = kwargs['pml_separation']
+
+        self.pmlthickness['x0'] = kwargs['subgrid_pml_thickness']
+        self.pmlthickness['y0'] = kwargs['subgrid_pml_thickness']
+        self.pmlthickness['z0'] = kwargs['subgrid_pml_thickness']
+        self.pmlthickness['xmax'] = kwargs['subgrid_pml_thickness']
+        self.pmlthickness['ymax'] = kwargs['subgrid_pml_thickness']
+        self.pmlthickness['zmax'] = kwargs['subgrid_pml_thickness']
+
+        # Number of sub cells to extend the sub grid beyond the IS boundary
+        d_to_pml = self.s_is_os_sep + self.pml_separation
+        self.n_boundary_cells = d_to_pml + self.pmlthickness['x0']
+        self.n_boundary_cells_x = d_to_pml + self.pmlthickness['x0']
+        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):
+        """
+        Return the equivalent spatial index of the global (i, j, k) in the subgrid.
+        Args:
+              p (numpy.array): i, j, k indices of a point in the main grid.
+        """
+        return self.n_boundary_cells_p + (p - self.p0) * self.ratio

gprMax/subgrids/factory.py

@@ -0,0 +1,28 @@
+from .subgrid_solver import SubGridSolver
+from .cpu_solvers import CPUSolver
+from .subgrid_updaters import SubgridUpdater
+from ..config import get_iterations
+from ..config import filter
+from precursor_nodes import PrecursorNodes
+from precursor_nodes import PrecursorNodesFiltered
+
+
+def create_solver(G, config):
+    """Factory method to return a solver instance for a CPU or subgrid simulation"""
+
+    iterations = get_iterations(config, G)
+
+    if filter:
+        from precursor_nodes import PrecursorNodesFiltered as PrecursorNodes
+    else:
+        from precursor_nodes import PrecursorNodes
+
+    updaters = []
+    for sg in G.subgrids:
+        precursors = PrecursorNodes(G, sg)
+        sgu = SubgridUpdater(sg, precursors, G)
+        updaters.append(sgu)
+
+    solver = SubGridSolver(G, updaters, iterations)
+
+    return solver

gprMax/subgrids/hsg.py

@@ -0,0 +1,142 @@
+from .base import SubGridBase
+from .hsg_fields_updates_ext import cython_update_is
+from .hsg_fields_updates_ext import cython_update_electric_os
+from .hsg_fields_updates_ext import cython_update_magnetic_os
+from ..utilities import human_size
+from ..utilities import memory_usage
+
+from colorama import init, Fore, Style
+init()
+
+
+class SubGridHSG(SubGridBase):
+
+    gridtype = '3DSUBGRID'
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.gridtype = SubGridHSG.gridtype
+
+        # upper and lower indices for the OS.
+        self.i_l = self.i0 - self.is_os_sep
+        self.i_u = self.i1 + self.is_os_sep
+        self.j_l = self.j0 - self.is_os_sep
+        self.j_u = self.j1 + self.is_os_sep
+        self.k_l = self.k0 - self.is_os_sep
+        self.k_u = self.k1 + self.is_os_sep
+
+    def update_magnetic_is(self, precursors):
+        """Apply the incident field corrections to the subgrid fields at the IS.
+
+        Args:
+              precursors (PrecursorNodes): PrecursorNodes are the incident
+              fields calculated at the correct time and position.
+        """
+        # bottom and top faces
+        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
+        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
+        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:
+              precursors (PrecursorNodes): PrecursorNodes are the incident
+              fields calculated at the correct time and position.
+        """
+        # bottom and top faces
+        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
+        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
+        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.
+
+        Main grid fields are collocated with subgrid fields therefore
+        precursor fields are not required.
+
+        Args:
+              main_grid (FDTDGrid): main fdtd solver grid
+        """
+        # front and back
+        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)
+
+        # 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, 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)
+
+        # 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, 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)
+
+    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
+        precursor fields are not required.
+
+        Args:
+              main_grid (FDTDGrid): main fdtd solver grid
+        """
+
+        # 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, 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)
+
+        # 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, 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)
+
+        # 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, 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)
+
+    def __str__(self):
+        s = '\n'
+        s += Fore.CYAN
+        s += 'Sub Grid HSG\n'
+        s += 'Name: {}\n'.format(self.name)
+        s += 'dx, dy, dz: {}m {}m {}m\n'.format(self.dx, self.dy, self.dz)
+        s += 'dt: {}s\n'.format(self.dt)
+        s += 'Memory Estimate: {}\n'.format(human_size(memory_usage(self)))
+        s += 'Position: ({}m, {}m, {}m), ({}m, {}m, {}m)\n'.format(self.x1,
+                                                                   self.y1,
+                                                                   self.z1,
+                                                                   self.x2,
+                                                                   self.y2,
+                                                                   self.z2)
+        s += 'Main Grid Indices: lower left({}, {}, {}), upper right({}, {}, {})\n'.format(self.i0, self.j0, self.k0, self.i1, self.j1, self.k1)
+        s += 'Total Cells: {} {} {}\n'.format(self.nx, self.ny, self.nz)
+        s += 'Working Region Cells: {} {} {}\n'.format(self.nwx,
+                                                       self.nwy,
+                                                       self.nwz)
+        for h in self.hertziandipoles:
+            s += 'Hertizian dipole: {} {} {}\n'.format(h.xcoord,
+                                                       h.ycoord,
+                                                       h.zcoord)
+            s += str([x for x in self.waveforms
+                      if x.ID == h.waveformID][0]) + '\n'
+        for r in self.rxs:
+            s += 'Receiver: {} {} {}\n'.format(r.xcoord, r.ycoord, r.zcoord)
+
+        for tl in self.transmissionlines:
+            s += 'Transmission Line: {} {} {}\n'.format(tl.xcoord, tl.ycoord, tl.zcoord)
+            s += str([x for x in self.waveforms
+                      if x.ID == tl.waveformID][0]) + '\n'
+        s += Style.RESET_ALL
+        return s

gprMax/subgrids/precursor_nodes.py

@@ -0,0 +1,608 @@
+from .config import sg_interp_d
+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 PrecusorNodesBase:
+
+    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 = sg_interp_d
+
+        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
+
+    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 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]
+            f1, f2 = 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(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
+        # 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 - 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.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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 1),
+                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 = [
+            ['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 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 = [
+            ['ex_front_1', True,  (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)), self.Ex],
+            ['ex_back_1', True,   (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)), self.Ex],
+            ['ez_front_1', False, (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)), self.Ez],
+            ['ez_back_1', False,   (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)), self.Ez],
+
+            ['ey_left_1', True,   (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Ey],
+            ['ey_right_1', True,  (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Ey],
+            ['ez_left_1', False,   (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Ez],
+            ['ez_right_1', False,  (self.i1, 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), self.Ex],
+            ['ex_top_1', True,    (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1), self.Ex],
+            ['ey_bottom_1', False, (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0), self.Ey],
+            ['ey_top_1', False,    (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1), 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_1 = obj[-1][obj[2]]
+        f_2 = obj[-1][obj[3]]
+        return f1, f2
+
+    def f_m_electric(self, obj):
+        return obj[-1][obj[2]]

gprMax/user_inputs.py

@@ -0,0 +1,147 @@
+from .exceptions import CmdInputError
+from .subgrids.base import SubGridBase
+from .utilities import round_value
+
+import numpy as np
+from colorama import init
+from colorama import Fore
+from colorama import Style
+
+init()
+
+"""Module contains classes to handle points supplied by a user. The
+classes implement a common interface such that geometry building objects
+such as Box or triangle do not need to have any knowledge which grid to which they
+are rounding continuous points or checking the point is within the grid.
+Additionally all logic related to rounding points etc is encapulsated here.
+"""
+
+
+def create_user_input_points(grid):
+    """Return a point checker class based on the grid supplied."""
+    if isinstance(grid, SubGridBase):
+        return SubgridUserInput(grid)
+    else:
+        return MainGridUserInput(grid)
+
+
+class UserInput:
+
+    """Base class to handle (x, y, z) points supplied by the user."""
+
+    def __init__(self, grid):
+        self.grid = grid
+
+    def point_within_bounds(self, p, cmd_str, name):
+        try:
+            self.grid.within_bounds(p)
+        except ValueError as err:
+            v = ['x', 'y', 'z']
+            # Discretisation
+            dl = getattr(self.grid, 'd' + err.args[0])
+            # Incorrect index
+            i = p[v.index(err.args[0])]
+            if name:
+                fs = "'{}' the {} {}-coordinate {:g} is not within the model domain"
+                s = fs.format(cmd_str, err.args[0], name, i * dl)
+            else:
+                fs = "'{}' {}-coordinate {:g} is not within the model domain"
+                s = fs.format(cmd_str, err.args[0], i * dl)
+            raise CmdInputError(s)
+
+    def discretise_point(self, p):
+        """Function to get the index of a continuous point with the grid."""
+        rv = np.vectorize(round_value)
+        return rv(p / self.grid.dl)
+
+    def round_to_grid(self, p):
+        """Function to get the nearest continuous point on the grid from a continuous point in space."""
+        return self.discretise_point(p) * self.grid.dl
+
+
+class MainGridUserInput(UserInput):
+    """Class to handle (x, y, z) points supplied by the user in the main grid."""
+
+    def __init__(self, grid):
+        super().__init__(grid)
+
+    def check_point(self, p, cmd_str, name=''):
+        p = self.discretise_point(p)
+        self.point_within_bounds(p, cmd_str, name)
+        return p
+
+    def check_src_rx_point(self, p, cmd_str, name=''):
+        p = self.check_point(p, cmd_str, name)
+
+        if self.grid.within_pml(p):
+            s = """WARNING: '{}' sources and receivers should not normally be \
+            positioned within the PML."""
+            fs = s.format(cmd_str)
+            print(Fore.RED + fs + Style.RESET_ALL)
+
+        return p
+
+    def check_box_points(self, p1, p2, cmd_str):
+        p1 = self.check_point(p1, cmd_str, name='lower')
+        p2 = self.check_point(p2, cmd_str, name='upper')
+
+        if np.greater(p1, p2).any():
+            s = """'{}' the lower coordinates should be less than the upper \
+            coordinates"""
+            raise CmdInputError(s.format(cmd_str))
+
+        return p1, p2
+
+    def check_tri_points(self, p1, p2, p3, cmd_str):
+        p1 = self.check_point(p1, cmd_str, name='vertex_1')
+        p2 = self.check_point(p2, cmd_str, name='vertex_2')
+        p3 = self.check_point(p3, cmd_str, name='vertex_3')
+
+        return p1, p2, p3
+
+
+class SubgridUserInput(MainGridUserInput):
+    """Class to handle (x, y, z) points supplied by the user in the sub grid."""
+
+    def __init__(self, grid):
+        super().__init__(grid)
+
+        # defines the region exposed to the user
+        self.inner_bound = np.array([
+                                    grid.pmlthickness['x0'] + grid.pml_separation,
+                                    grid.pmlthickness['y0'] + grid.pml_separation,
+                                    grid.pmlthickness['z0'] + grid.pml_separation])
+
+        self.outer_bound = np.subtract([grid.nx, grid.ny, grid.nz],
+                                       self.inner_bound)
+
+    def translate_to_gap(self, p):
+        """Function to translate the user input point to the real point in the
+        subgrid"""
+        return np.add(p, self.inner_bound)
+
+    def discretise_point(self, p):
+        """"Function to discretise a point. Does not provide any checks. The
+        user enters coordinates relative to self.inner_bound. This function
+        translate the user point to the correct index for building objects"""
+        p = super().discretise_point(p)
+        p_t = self.translate_to_gap(p)
+        return p_t
+
+    def round_to_grid(self, p):
+        p_t = self.discretise_point(p)
+        p_m = p_t * self.grid.dl
+        return p_m
+
+    def check_point(self, p, cmd_str, name=''):
+        p_t = super().check_point(p, cmd_str, name)
+
+        # Provide user within a warning if they have placed objects within
+        # the OS non-working region.
+        if (np.less(p_t, self.inner_bound).any() or
+            np.greater(p_t, self.outer_bound).any()):
+                s = """WARNING: '{}' This object traverses the Outer Surface. This \
+                is an advanced feature."""
+                print(Fore.RED + s.format(cmd_str) + Style.RESET_ALL)
+
+        return p_t

todo

@@ -0,0 +1 @@
+memory update for all grid (model_build_run.py)