From 4dfe017d609cb52650dcd8bce20255f7b7fcb8f7 Mon Sep 17 00:00:00 2001 From: John Hartley Date: Mon, 29 Jul 2019 15:08:07 +0100 Subject: [PATCH] transfer existing subgridding modules --- gprMax/__init__.py | 45 + gprMax/cmds_geometry/__init__.py | 0 gprMax/cmds_geometry/add_grass.py | 164 +++ gprMax/cmds_geometry/add_surface_roughness.py | 126 ++ gprMax/cmds_geometry/add_surface_water.py | 118 ++ gprMax/cmds_geometry/box.py | 92 ++ gprMax/cmds_geometry/cmds_geometry.py | 28 + gprMax/cmds_geometry/cylinder.py | 97 ++ gprMax/cmds_geometry/cylindrical_sector.py | 137 +++ gprMax/cmds_geometry/edge.py | 64 ++ gprMax/cmds_geometry/fractal_box.py | 92 ++ gprMax/cmds_geometry/fractal_box_builder.py | 286 +++++ gprMax/cmds_geometry/geometry_objects_read.py | 89 ++ gprMax/cmds_geometry/plate.py | 112 ++ gprMax/cmds_geometry/sphere.py | 91 ++ gprMax/cmds_geometry/triangle.py | 132 +++ gprMax/cmds_multiple.py | 1022 +++++++++++++++++ gprMax/cmds_single_use.py | 507 ++++++++ gprMax/cython/fields_updates_hsg.pyx | 269 +++++ gprMax/scene.py | 122 ++ gprMax/subgrids/base.py | 61 + gprMax/subgrids/factory.py | 28 + gprMax/subgrids/hsg.py | 142 +++ gprMax/subgrids/precursor_nodes.py | 608 ++++++++++ gprMax/user_inputs.py | 147 +++ temp/gprMax | 1 + todo | 1 + 27 files changed, 4581 insertions(+) create mode 100644 gprMax/cmds_geometry/__init__.py create mode 100644 gprMax/cmds_geometry/add_grass.py create mode 100644 gprMax/cmds_geometry/add_surface_roughness.py create mode 100644 gprMax/cmds_geometry/add_surface_water.py create mode 100644 gprMax/cmds_geometry/box.py create mode 100644 gprMax/cmds_geometry/cmds_geometry.py create mode 100644 gprMax/cmds_geometry/cylinder.py create mode 100644 gprMax/cmds_geometry/cylindrical_sector.py create mode 100644 gprMax/cmds_geometry/edge.py create mode 100644 gprMax/cmds_geometry/fractal_box.py create mode 100644 gprMax/cmds_geometry/fractal_box_builder.py create mode 100644 gprMax/cmds_geometry/geometry_objects_read.py create mode 100644 gprMax/cmds_geometry/plate.py create mode 100644 gprMax/cmds_geometry/sphere.py create mode 100644 gprMax/cmds_geometry/triangle.py create mode 100644 gprMax/cmds_multiple.py create mode 100644 gprMax/cmds_single_use.py create mode 100644 gprMax/cython/fields_updates_hsg.pyx create mode 100644 gprMax/scene.py create mode 100644 gprMax/subgrids/base.py create mode 100644 gprMax/subgrids/factory.py create mode 100644 gprMax/subgrids/hsg.py create mode 100644 gprMax/subgrids/precursor_nodes.py create mode 100644 gprMax/user_inputs.py create mode 160000 temp/gprMax create mode 100644 todo diff --git a/gprMax/__init__.py b/gprMax/__init__.py index ef348ef1..ee721632 100644 --- a/gprMax/__init__.py +++ b/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' diff --git a/gprMax/cmds_geometry/__init__.py b/gprMax/cmds_geometry/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py new file mode 100644 index 00000000..efc052a8 --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py new file mode 100644 index 00000000..65289fea --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/add_surface_water.py b/gprMax/cmds_geometry/add_surface_water.py new file mode 100644 index 00000000..75e80629 --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/box.py b/gprMax/cmds_geometry/box.py new file mode 100644 index 00000000..bb8636f1 --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/cmds_geometry.py b/gprMax/cmds_geometry/cmds_geometry.py new file mode 100644 index 00000000..d9cc5f6a --- /dev/null +++ b/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.""" diff --git a/gprMax/cmds_geometry/cylinder.py b/gprMax/cmds_geometry/cylinder.py new file mode 100644 index 00000000..d6f21d55 --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/cylindrical_sector.py b/gprMax/cmds_geometry/cylindrical_sector.py new file mode 100644 index 00000000..9fa100f9 --- /dev/null +++ b/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))) diff --git a/gprMax/cmds_geometry/edge.py b/gprMax/cmds_geometry/edge.py new file mode 100644 index 00000000..e5158f2b --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py new file mode 100644 index 00000000..1d96bd07 --- /dev/null +++ b/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) diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py new file mode 100644 index 00000000..ebce587c --- /dev/null +++ b/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) diff --git a/gprMax/cmds_geometry/geometry_objects_read.py b/gprMax/cmds_geometry/geometry_objects_read.py new file mode 100644 index 00000000..bd5c1ee2 --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/plate.py b/gprMax/cmds_geometry/plate.py new file mode 100644 index 00000000..eea23c47 --- /dev/null +++ b/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))) diff --git a/gprMax/cmds_geometry/sphere.py b/gprMax/cmds_geometry/sphere.py new file mode 100644 index 00000000..983c6e06 --- /dev/null +++ b/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)) diff --git a/gprMax/cmds_geometry/triangle.py b/gprMax/cmds_geometry/triangle.py new file mode 100644 index 00000000..c742f7cb --- /dev/null +++ b/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))) diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py new file mode 100644 index 00000000..558bdf7b --- /dev/null +++ b/gprMax/cmds_multiple.py @@ -0,0 +1,1022 @@ +"""Object which can be created multiple times.""" +from gprMax.constants import z0 +from .constants import floattype +from .utilities import round_value +from .cmds_geometry.cmds_geometry import UserObjectGeometry +from .waveforms import Waveform as WaveformUser +from .sources import VoltageSource as VoltageSourceUser +from .sources import HertzianDipole as HertzianDipoleUser +from .sources import MagneticDipole as MagneticDipoleUser +from .sources import TransmissionLine as TransmissionLineUser +from .snapshots import Snapshot as SnapshotUser +from .receivers import Rx as RxUser +from .materials import Material as MaterialUser +from .materials import PeplinskiSoil +from .geometry_outputs import GeometryObjects as GeometryObjectsUser +from .pml import CFSParameter +from .pml import CFS +from .subgrids.base import SubGridBase + +from .exceptions import CmdInputError + + +import numpy as np +from tqdm import tqdm + + +class UserObjectMulti: + """Specific multiobject object.""" + + def __init__(self, **kwargs): + """Constructor.""" + self.kwargs = kwargs + self.order = None + 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 create(self, grid, uip): + """Create the object and add it to the grid.""" + pass + + def params_str(self): + """Readble string of parameters given to object.""" + return self.hash + ': ' + str(self.kwargs) + + +class Waveform(UserObjectMulti): + """Create Waveform different excitation types.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 0 + self.hash = '#waveform' + + def create(self, grid, uip): + """Create the waveform and add it to the grid.""" + try: + wavetype = self.kwargs['wave_type'].lower() + amp = self.kwargs['amp'] + freq = self.kwargs['freq'] + ID = self.kwargs['id'] + + except KeyError: + raise CmdInputError(self.params_str() + ' requires exactly four parameters') + + if wavetype not in WaveformUser.types: + raise CmdInputError(self.__str__() + ' must have one of the following types {}'.format(','.join(WaveformUser.types))) + if freq <= 0: + raise CmdInputError(self.__str__() + ' requires an excitation frequency value of greater than zero') + if any(x.ID == ID for x in grid.waveforms): + raise CmdInputError(self.__str__() + ' with ID {} already exists'.format(ID)) + + w = WaveformUser() + w.ID = ID + w.type = wavetype + w.amp = amp + w.freq = freq + + if grid.messages: + print('Waveform {} of type {} with maximum amplitude scaling {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq)) + + grid.waveforms.append(w) + + +class VoltageSource(UserObjectMulti): + """User Object for a voltage source.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 1 + self.hash = '#voltage_source' + + def create(self, grid, uip): + """Create voltage source and add it to the grid.""" + try: + polarisation = self.kwargs['polarisation'].lower() + p1 = self.kwargs['p1'] + resistance = self.kwargs['resistance'] + waveform_id = self.kwargs['waveform_id'] + + except KeyError: + raise CmdInputError(self.__str__() + " requires at least six parameters") + + # Check polarity & position parameters + if polarisation not in ('x', 'y', 'z'): + raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__())) + if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__())) + elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__())) + elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): + raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__())) + + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + + if resistance < 0: + raise CmdInputError("'{}' requires a source resistance of zero or greater".format(self.__str__())) + + # Check if there is a waveformID in the waveforms list + if not any(x.ID == waveform_id for x in grid.waveforms): + raise CmdInputError("'{}' there is no waveform with the identifier {}'.format(tmp[5]".format(self.__str__())) + + v = VoltageSourceUser() + v.polarisation = polarisation + v.xcoord = xcoord + v.ycoord = ycoord + v.zcoord = zcoord + v.ID = v.__class__.__name__ + '(' + str(v.xcoord) + ',' + str(v.ycoord) + ',' + str(v.zcoord) + ')' + v.resistance = resistance + v.waveformID = waveform_id + + try: + start = self.kwargs['start'] + stop = self.kwargs['stop'] + # Check source start & source remove time parameters + if start < 0: + raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__())) + if stop < 0: + raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__())) + if stop - start <= 0: + raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__())) + v.start = start + if stop > grid.timewindow: + v.stop = grid.timewindow + else: + v.stop = stop + startstop = ' start time {:g} secs, finish time {:g} secs '.format(v.start, v.stop) + except KeyError: + v.start = 0 + v.stop = grid.timewindow + startstop = ' ' + + v.calculate_waveform_values(grid) + + if grid.messages: + print('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * grid.dx, v.ycoord * grid.dy, v.zcoord * grid.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID)) + + grid.voltagesources.append(v) + + +class HertzianDipole(UserObjectMulti): + """User Object for HertzianDipole.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 2 + self.hash = '#hertzian_dipole' + + def create(self, grid, uip): + """Create HertzianDipole and add it to the grid.""" + try: + polarisation = self.kwargs['polarisation'].lower() + p1 = self.kwargs['p1'] + waveform_id = self.kwargs['waveform_id'] + + except KeyError: + raise CmdInputError("'{}' requires at least 3 parameters".format(self.params_str())) + + # Check polarity & position parameters + if polarisation not in ('x', 'y', 'z'): + raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__())) + if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__())) + elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__())) + elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): + raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__())) + + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + + # Check if there is a waveformID in the waveforms list + if not any(x.ID == waveform_id for x in grid.waveforms): + raise CmdInputError("'{}' there is no waveform with the identifier {}'.format(tmp[4]".format(self.__str__())) + + h = HertzianDipoleUser() + h.polarisation = polarisation + + # Set length of dipole to grid size in polarisation direction + if h.polarisation == 'x': + h.dl = grid.dx + elif h.polarisation == 'y': + h.dl = grid.dy + elif h.polarisation == 'z': + h.dl = grid.dz + + h.xcoord = xcoord + h.ycoord = ycoord + h.zcoord = zcoord + h.xcoordorigin = xcoord + h.ycoordorigin = ycoord + h.zcoordorigin = zcoord + h.ID = h.__class__.__name__ + '(' + str(h.xcoord) + ',' + str(h.ycoord) + ',' + str(h.zcoord) + ')' + h.waveformID = waveform_id + + try: + # Check source start & source remove time parameters + start = self.kwargs['start'] + stop = self.kwargs['stop'] + if start < 0: + raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__())) + if stop < 0: + raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__())) + if stop - start <= 0: + raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__())) + h.start = start + if stop > grid.timewindow: + h.stop = grid.timewindow + else: + h.stop = stop + startstop = ' start time {:g} secs, finish time {:g} secs '.format(h.start, h.stop) + except KeyError: + h.start = 0 + h.stop = grid.timewindow + print(grid.timewindow) + startstop = ' ' + + h.calculate_waveform_values(grid) + + if grid.messages: + if grid.mode == '2D': + print('Hertzian dipole is a line source in 2D with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * grid.dx, h.ycoord * grid.dy, h.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(h.waveformID)) + else: + print('Hertzian dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * grid.dx, h.ycoord * grid.dy, h.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(h.waveformID)) + + grid.hertziandipoles.append(h) + + +class MagneticDipole(UserObjectMulti): + """Magnetic Dipole User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 3 + self.hash = '#magnetic_dipole' + + def create(self, grid, uip): + """Create Magnetic Dipole and add it the grid.""" + try: + polarisation = self.kwargs['polarisation'].lower() + p1 = self.kwargs['p1'] + waveform_id = self.kwargs['waveform_id'] + except KeyError: + raise CmdInputError("'{}' requires at least five parameters".format(self.__str__())) + + # Check polarity & position parameters + if polarisation not in ('x', 'y', 'z'): + raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__())) + if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__())) + elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__())) + elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): + raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__())) + + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + + # Check if there is a waveformID in the waveforms list + if not any(x.ID == waveform_id for x in grid.waveforms): + raise CmdInputError("'{}' there is no waveform with the identifier {}".format(self.__str__(), waveform_id)) + + m = MagneticDipoleUser() + m.polarisation = polarisation + m.xcoord = xcoord + m.ycoord = ycoord + m.zcoord = zcoord + m.xcoordorigin = xcoord + m.ycoordorigin = ycoord + m.zcoordorigin = zcoord + m.ID = m.__class__.__name__ + '(' + str(m.xcoord) + ',' + str(m.ycoord) + ',' + str(m.zcoord) + ')' + m.waveformID = waveform_id + + try: + # Check source start & source remove time parameters + start = self.kwargs['start'] + stop = self.kwargs['stop'] + if start < 0: + raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__())) + if stop < 0: + raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__())) + if stop - start <= 0: + raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__())) + m.start = start + if stop > grid.timewindow: + m.stop = grid.timewindow + else: + m.stop = stop + startstop = ' start time {:g} secs, finish time {:g} secs '.format(m.start, m.stop) + except KeyError: + m.start = 0 + m.stop = grid.timewindow + startstop = ' ' + + m.calculate_waveform_values(grid) + + if grid.messages: + print('Magnetic dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(m.polarisation, m.xcoord * grid.dx, m.ycoord * grid.dy, m.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(m.waveformID)) + + grid.magneticdipoles.append(m) + + +class TransmissionLine(UserObjectMulti): + """Magnetic Dipole User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 4 + self.hash = '#transmission_line' + + def create(self, grid, uip): + + try: + polarisation = self.kwargs['polarisation'].lower() + p1 = self.kwargs['p1'] + waveform_id = self.kwargs['waveform_id'] + resistance = self.kwargs['resistance'] + + except KeyError: + raise CmdInputError("'{}' requires at least six parameters".format(self.params_str())) + + # Warn about using a transmission line on GPU + if grid.gpu is not None: + raise CmdInputError("'{}' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.".format(self.__str__())) + + # Check polarity & position parameters + if polarisation not in ('x', 'y', 'z'): + raise CmdInputError("'{}' polarisation must be x, y, or z".format(self.__str__())) + if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be x in 2D TMx mode".format(self.__str__())) + elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): + raise CmdInputError("'{}' polarisation must be y in 2D TMy mode".format(self.__str__())) + elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): + raise CmdInputError("'{}' polarisation must be z in 2D TMz mode".format(self.__str__())) + + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + + if resistance <= 0 or resistance >= z0: + raise CmdInputError("'{}' requires a resistance greater than zero and less than the impedance of free space, i.e. 376.73 Ohms".format(self.__str__())) + + # Check if there is a waveformID in the waveforms list + if not any(x.ID == waveform_id for x in grid.waveforms): + raise CmdInputError("'{}' there is no waveform with the identifier {}'.format(tmp[5]".format(self.__str__())) + + t = TransmissionLineUser(grid) + t.polarisation = polarisation + t.xcoord = xcoord + t.ycoord = ycoord + t.zcoord = zcoord + t.ID = t.__class__.__name__ + '(' + str(t.xcoord) + ',' + str(t.ycoord) + ',' + str(t.zcoord) + ')' + t.resistance = resistance + t.waveformID = waveform_id + + try: + # Check source start & source remove time parameters + start = self.kwargs['start'] + stop = self.kwargs['stop'] + if start < 0: + raise CmdInputError("'{}' delay of the initiation of the source should not be less than zero".format(self.__str__())) + if stop < 0: + raise CmdInputError("'{}' time to remove the source should not be less than zero".format(self.__str__())) + if stop - start <= 0: + raise CmdInputError("'{}' duration of the source should not be zero or less".format(self.__str__())) + t.start = start + if stop > grid.timewindow: + t.stop = grid.timewindow + else: + t.stop = stop + startstop = ' start time {:g} secs, finish time {:g} secs '.format(t.start, t.stop) + except KeyError: + t.start = 0 + t.stop = grid.timewindow + startstop = ' ' + + t.calculate_waveform_values(grid) + t.calculate_incident_V_I(grid) + + if grid.messages: + print('Transmission line with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(t.polarisation, t.xcoord * grid.dx, t.ycoord * grid.dy, t.zcoord * grid.dz, t.resistance) + startstop + 'using waveform {} created.'.format(t.waveformID)) + + grid.transmissionlines.append(t) + + +class Rx(UserObjectMulti): + """Magnetic Dipole User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 5 + self.hash = '#rx' + self.constructor = RxUser + + def create(self, grid, uip): + try: + p1 = self.kwargs['p1'] + except KeyError: + raise CmdInputError("'{}' has an incorrect number of parameters".format(self.__str__())) + + p = uip.check_src_rx_point(p1, self.__str__()) + + r = self.constructor() + r.xcoord, r.ycoord, r.zcoord = p + r.xcoordorigin, r.ycoordorigin, r.zcoordorigin = p + + try: + r.ID = self.kwargs['ID'] + outputs = self.kwargs['outputs'] + # Get allowable outputs + if grid.gpu is not None: + allowableoutputs = RxUser.gpu_allowableoutputs + else: + allowableoutputs = RxUser.allowableoutputs + # Check and add field output names + for field in outputs: + if field in allowableoutputs: + r.outputs[field] = np.zeros(grid.iterations, dtype=floattype) + else: + raise CmdInputError("{} contains an output type that is not allowable. Allowable outputs in current context are {}".format(self.__str__(), allowableoutputs)) + + # If no ID or outputs are specified, use default + except KeyError: + r.ID = r.__class__.__name__ + '(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')' + for key in RxUser.defaultoutputs: + r.outputs[key] = np.zeros(grid.iterations, dtype=floattype) + + if grid.messages: + print('Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs))) + + grid.rxs.append(r) + + return r + +class RxArray(UserObjectMulti): + """Receiver Array User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 6 + self.hash = '#rx_array' + + def create(self, grid, uip): + + try: + p1 = self.kwargs['p1'] + p2 = self.kwargs['p2'] + dl = self.kwargs['dl'] + + except KeyError: + raise CmdInputError("'{}' requires exactly 9 parameters".format(self.__str__())) + + xs, ys, zs = uip.check_src_rx_point(p1, self.__str__(), 'lower') + xf, yf, zf = uip.check_src_rx_point(p2, self.__str__(), 'upper') + dx, dy, dz = uip.discretise_point(dl) + + if xs > xf or ys > yf or zs > zf: + raise CmdInputError("'{}' the lower coordinates should be less than the upper coordinates".format(self.__str__())) + if dx < 0 or dy < 0 or dz < 0: + raise CmdInputError("'{}' the step size should not be less than zero".format(self.__str__())) + if dx < 1: + if dx == 0: + dx = 1 + else: + raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__())) + if dy < 1: + if dy == 0: + dy = 1 + else: + raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__())) + if dz < 1: + if dz == 0: + dz = 1 + else: + raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__())) + + if grid.messages: + print('Receiver array {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g}m'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz)) + + for x in range(xs, xf + 1, dx): + for y in range(ys, yf + 1, dy): + for z in range(zs, zf + 1, dz): + r = RxUser() + r.xcoord = x + r.ycoord = y + r.zcoord = z + r.xcoordorigin = x + r.ycoordorigin = y + r.zcoordorigin = z + r.ID = r.__class__.__name__ + '(' + str(x) + ',' + str(y) + ',' + str(z) + ')' + for key in RxUser.defaultoutputs: + r.outputs[key] = np.zeros(grid.iterations, dtype=floattype) + if grid.messages: + print(' Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs))) + grid.rxs.append(r) + + +class Snapshot(UserObjectMulti): + """Snapshot User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 19 + self.hash = '#snapshot' + + def create(self, grid, uip): + + if isinstance(grid, SubGridBase): + raise CmdInputError("'{}' Do not add Snapshots to Subgrids.".format(self.__str__())) + try: + p1 = self.kwargs['p1'] + p2 = self.kwargs['p2'] + dl = self.kwargs['dl'] + filename = self.kwargs['filename'] + except KeyError: + raise CmdInputError("'{}' requires exactly 11 parameters".format(self.__str__())) + + p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + xs, ys, zs = p1 + xf, yf, zf = p2 + dx, dy, dz = uip.discretise_point(dl) + + # If number of iterations given + try: + iterations = self.kwargs['iterations'] + # If real floating point value given + except KeyError: + try: + time = self.kwargs['time'] + except KeyError: + raise CmdInputError("'{}' requires exactly 5 parameters".format(self.__str__())) + if time > 0: + iterations = round_value((time / grid.dt)) + 1 + else: + raise CmdInputError("'{}' time value must be greater than zero".format(self.__str__())) + + if dx < 0 or dy < 0 or dz < 0: + raise CmdInputError("'{}' the step size should not be less than zero".format(self.__str__())) + if dx < 1 or dy < 1 or dz < 1: + raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__())) + if iterations <= 0 or iterations > grid.iterations: + raise CmdInputError("'{}' time value is not valid".format(self.__str__())) + + # Replace with old style snapshots if there are subgrids + if grid.subgrids: + from .snapshot_subgrid import Snapshot as SnapshotSub + s = SnapshotSub(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, filename) + else: + s = SnapshotUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, filename) + + if grid.messages: + print('Snapshot from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, at {:g} secs with filename {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, s.time * grid.dt, s.basefilename)) + + grid.snapshots.append(s) + +class Material(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 8 + self.hash = '#material' + + def create(self, grid, uip): + try: + er = self.kwargs['er'] + se = self.kwargs['se'] + mr = self.kwargs['mr'] + sm = self.kwargs['sm'] + material_id = self.kwargs['id'] + except KeyError: + raise CmdInputError('{} requires exactly five parameters'.format(self.params_str())) + + if er < 1: + raise CmdInputError('{} requires a positive value of one or greater for static (DC) permittivity'.format(self.__str__())) + if se != 'inf': + se = float(se) + if se < 0: + raise CmdInputError('{} requires a positive value for conductivity'.format(self.__str__())) + else: + se = float('inf') + if mr < 1: + raise CmdInputError('{} requires a positive value of one or greater for permeability'.format(self.__str__())) + if sm < 0: + raise CmdInputError('{} requires a positive value for magnetic conductivity'.format(self.__str__())) + if any(x.ID == material_id for x in grid.materials): + raise CmdInputError('{} with ID {} already exists'.format(material_id).format(self.__str__())) + + # Create a new instance of the Material class material (start index after pec & free_space) + m = MaterialUser(len(grid.materials), material_id) + m.er = er + m.se = se + m.mr = mr + m.sm = sm + + # Set material averaging to False if infinite conductivity, i.e. pec + if m.se == float('inf'): + m.averagable = False + + if grid.messages: + tqdm.write('Material {} with eps_r={:g}, sigma={:g} S/m; mu_r={:g}, sigma*={:g} Ohm/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm)) + + # Append the new material object to the materials list + grid.materials.append(m) + + +class AddDebyeDispersion(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 9 + self.hash = '#add_dispersion_debye' + + def create(self, grid, uip): + + try: + poles = self.kwargs['n_poles'] + er_delta = self.kwargs['er_delta'] + tau = self.kwargs['tau'] + material_ids = self.kwargs['material_ids'] + + except KeyError: + raise CmdInputError(self.__str__() + ' requires at least four parameters') + + if poles < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for number of poles') + + # Look up requested materials in existing list of material instances + materials = [y for x in material_ids for y in grid.materials if y.ID == x] + + if len(materials) != len(material_ids): + notfound = [x for x in material_ids if x not in materials] + raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound)) + + for material in materials: + material.type = 'debye' + material.poles = poles + material.averagable = False + for i in range(0, poles): + # N.B Not checking if relaxation times are greater than time-step + if tau[i] > 0: + material.deltaer.append(er_delta[i]) + material.tau.append(tau[i]) + else: + raise CmdInputError(self.__str__() + ' requires positive values for the permittivity difference.') + if material.poles > MaterialUser.maxpoles: + MaterialUser.maxpoles = material.poles + + if grid.messages: + tqdm.write('Debye disperion added to {} with delta_eps_r={}, and tau={} secs created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau))) + + +class AddLorentzDispersion(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 10 + self.hash = '#add_dispersion_lorentz' + + def create(self, grid, uip): + + try: + poles = self.kwargs['n_poles'] + er_delta = self.kwargs['er_delta'] + tau = self.kwargs['omega'] + alpha = self.kwargs['delta'] + material_ids = self.kwargs['material_ids'] + except KeyError: + raise CmdInputError(self.__str__() + ' requires at least five parameters') + + if poles < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for number of poles') + + # Look up requested materials in existing list of material instances + materials = [y for x in material_ids for y in grid.materials if y.ID == x] + + if len(materials) != len(material_ids): + notfound = [x for x in material_ids if x not in materials] + raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound)) + + for material in materials: + material.type = 'lorentz' + material.poles = poles + material.averagable = False + for i in range(0, poles): + if er_delta[i] > 0 and tau[i] > grid.dt and alpha[i] > grid.dt: + material.deltaer.append(er_delta[i]) + material.tau.append(tau[i]) + material.alpha.append(alpha[i]) + else: + raise CmdInputError(self.__str__() + ' requires positive values for the permittivity difference and frequencies, and associated times that are greater than the time step for the model.') + if material.poles > MaterialUser.maxpoles: + MaterialUser.maxpoles = material.poles + + if grid.messages: + tqdm.write('Lorentz disperion added to {} with delta_eps_r={}, omega={} secs, and gamma={} created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha))) + + +class AddDrudeDispersion(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 11 + self.hash = '#add_dispersion_Drude' + + def create(self, grid, uip): + + try: + poles = self.kwargs['n_poles'] + tau = self.kwargs['tau'] + alpha = self.kwargs['alpha'] + material_ids = self.kwargs['material_ids'] + except KeyError: + raise CmdInputError(self.__str__() + ' requires at least four parameters') + if poles < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for number of poles') + + # Look up requested materials in existing list of material instances + materials = [y for x in material_ids for y in grid.materials if y.ID == x] + + if len(materials) != len(material_ids): + notfound = [x for x in material_ids if x not in materials] + raise CmdInputError(self.__str__() + ' material(s) {} do not exist'.format(notfound)) + + for material in materials: + material.type = 'drude' + material.poles = poles + material.averagable = False + for i in range(0, poles): + if tau[i] > 0 and alpha[i] > grid.dt: + material.tau.append(tau[i]) + material.alpha.append(alpha[i]) + else: + raise CmdInputError(self.__str__() + ' requires positive values for the frequencies, and associated times that are greater than the time step for the model.') + if material.poles > MaterialUser.maxpoles: + MaterialUser.maxpoles = material.poles + + if grid.messages: + tqdm.write('Drude disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha))) + + +class SoilPeplinski(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 12 + self.hash = '#soil_peplinski' + + def create(self, grid, uip): + + try: + sand_fraction = self.kwargs['sand_fraction'] + clay_fraction = self.kwargs['clay_fraction'] + bulk_density = self.kwargs['bulk_density'] + sand_density = self.kwargs['sand_density'] + water_fraction_lower = self.kwargs['water_fraction_lower'] + water_fraction_upper = self.kwargs['water_fraction_upper'] + ID = self.kwargs['ID'] + + except KeyError: + raise CmdInputError(self.__str__() + ' requires at exactly seven parameters') + + if sand_fraction < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for the sand fraction') + if clay_fraction < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for the clay fraction') + if bulk_density < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for the bulk density') + if sand_density < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for the sand particle density') + if water_fraction_lower < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for the lower limit of the water volumetric fraction') + if water_fraction_upper < 0: + raise CmdInputError(self.__str__() + ' requires a positive value for the upper limit of the water volumetric fraction') + if any(x.ID == ID for x in grid.mixingmodels): + raise CmdInputError(self.__str__() + ' with ID {} already exists'.format(ID)) + + # Create a new instance of the Material class material (start index after pec & free_space) + s = PeplinskiSoil(ID, sand_fraction, clay_fraction, bulk_density, sand_density, (water_fraction_lower, water_fraction_upper)) + + if grid.messages: + print('Mixing model (Peplinski) used to create {} with sand fraction {:g}, clay fraction {:g}, bulk density {:g}g/cm3, sand particle density {:g}g/cm3, and water volumetric fraction {:g} to {:g} created.'.format(s.ID, s.S, s.C, s.rb, s.rs, s.mu[0], s.mu[1])) + + # Append the new material object to the materials list + grid.mixingmodels.append(s) + + +class GeometryView(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 18 + self.hash = '#geometry_view' + self.multi_grid = False + + def geometry_view_constructor(self, grid, output_type): + # Check if user want geometry output for all grids + try: + self.kwargs['multi_grid'] + # there is no voxel output for multi grid output + if isinstance(grid, SubGridBase): + raise CmdInputError("'{}' Do not add multi_grid output to subgrid user object. Please add to Scene".format(self.__str__())) + if output_type == 'n': + raise CmdInputError("'{}' Voxel output type (n) is not supported for multigrid output :(".format(self.__str__())) + # Change constructor to the multi grid output + from .geometry_outputs import GeometryViewFineMultiGrid as GeometryViewUser + self.multi_grid = True + except KeyError: + self.multi_grid = False + from .geometry_outputs import GeometryView as GeometryViewUser + + return GeometryViewUser + + def create(self, grid, uip): + try: + p1 = self.kwargs['p1'] + p2 = self.kwargs['p2'] + dl = self.kwargs['dl'] + output_type = self.kwargs['output_type'].lower() + filename = self.kwargs['filename'] + except KeyError: + raise CmdInputError("'{}' requires exactly eleven parameters".format(self.__str__())) + + GeometryViewUser = self.geometry_view_constructor(grid, output_type) + + p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + xs, ys, zs = p1 + xf, yf, zf = p2 + + dx, dy, dz = uip.discretise_point(dl) + + if dx < 0 or dy < 0 or dz < 0: + raise CmdInputError("'{}' the step size should not be less than zero".format(self.__str__())) + if dx > grid.nx or dy > grid.ny or dz > grid.nz: + raise CmdInputError("'{}' the step size should be less than the domain size".format(self.__str__())) + if dx < 1 or dy < 1 or dz < 1: + raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__())) + if output_type != 'n' and output_type != 'f': + raise CmdInputError("'{}' requires type to be either n (normal) or f (fine)".format(self.__str__())) + if output_type == 'f' and (dx * grid.dx != grid.dx or dy * grid.dy != grid.dy or dz * grid.dz != grid.dz): + raise CmdInputError("'{}' requires the spatial discretisation for the geometry view to be the same as the model for geometry view of type f (fine)".format(self.__str__())) + + # Set type of geometry file + if output_type == 'n': + fileext = '.vti' + else: + fileext = '.vtp' + + g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, grid) + + if grid.messages: + print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, multi_grid {}, grid={}, with filename base {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, self.multi_grid, grid.name, g.basefilename)) + + # Append the new GeometryView object to the geometry views list + grid.geometryviews.append(g) + + +class GeometryObjectsWrite(UserObjectMulti): + """Material User Object.""" + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 14 + self.hash = '#geometry_objects_write' + + def create(self, grid, uip): + try: + p1 = self.kwargs['p1'] + p2 = self.kwargs['p2'] + filename = self.kwargs['filename'] + except KeyError: + raise CmdInputError("'{}' requires exactly seven parameters".format(self.__str__())) + + p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + x0, y0, z0 = p1 + x1, y1, z1 = p2 + + g = GeometryObjectsUser(x0, y0, z0, x1, y1, z1, filename) + + if grid.messages: + print('Geometry objects in the volume from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, will be written to {}, with materials written to {}'.format(p1[0] * grid.dx, p1[1] * grid.dy, p1[2] * grid.dz, p2[0] * grid.dx, p2[1] * grid.dy, p2[2] * grid.dz, g.filename, g.materialsfilename)) + + # Append the new GeometryView object to the geometry objects to write list + grid.geometryobjectswrite.append(g) + + +class PMLCFS(UserObjectMulti): + """Material User Object.""" + + count = 0 + + def __init__(self, **kwargs): + """Constructor.""" + super().__init__(**kwargs) + self.order = 15 + self.hash = '#pml_cfs' + PMLCFS.count += 1 + if PMLCFS.count == 2: + raise CmdInputError(self.__str__() + ' can only be used up to two times, for up to a 2nd order PML') + + def create(self, grid, uip): + + try: + alphascalingprofile = self.kwargs['alphascalingprofile'] + alphascalingdirection = self.kwargs['alphascalingdirection'] + alphamin = self.kwargs['alphamin'] + alphamax = self.kwargs['alphamax'] + kappascalingprofile = self.kwargs['kappascalingprofile'] + kappascalingdirection = self.kwargs['kappascalingdirection'] + kappamin = self.kwargs['kappamin'] + kappamax = self.kwargs['kappamax'] + sigmascalingprofile = self.kwargs['sigmascalingprofile'] + sigmascalingdirection = self.kwargs['sigmascalingdirection'] + sigmamin = self.kwargs['sigmamin'] + sigmamax = self.kwargs['sigmamax'] + + except KeyError: + raise CmdInputError(self.__str__() + ' requires exactly twelve parameters') + + if alphascalingprofile not in CFSParameter.scalingprofiles.keys() or kappascalingprofile not in CFSParameter.scalingprofiles.keys() or sigmascalingprofile not in CFSParameter.scalingprofiles.keys(): + raise CmdInputError(self.__str__() + ' must have scaling type {}'.format(','.join(CFSParameter.scalingprofiles.keys()))) + if alphascalingdirection not in CFSParameter.scalingdirections or kappascalingdirection not in CFSParameter.scalingdirections or sigmascalingdirection not in CFSParameter.scalingdirections: + raise CmdInputError(self.__str__() + ' must have scaling type {}'.format(','.join(CFSParameter.scalingdirections))) + if float(alphamin) < 0 or float(alphamax) < 0 or float(kappamin) < 0 or float(kappamax) < 0 or float(sigmamin) < 0: + raise CmdInputError(self.__str__() + ' minimum and maximum scaling values must be greater than zero') + if float(kappamin) < 1: + raise CmdInputError(self.__str__() + ' minimum scaling value for kappa must be greater than or equal to one') + + cfsalpha = CFSParameter() + cfsalpha.ID = 'alpha' + cfsalpha.scalingprofile = alphascalingprofile + cfsalpha.scalingdirection = alphascalingdirection + cfsalpha.min = float(alphamin) + cfsalpha.max = float(alphamax) + cfskappa = CFSParameter() + cfskappa.ID = 'kappa' + cfskappa.scalingprofile = kappascalingprofile + cfskappa.scalingdirection = kappascalingdirection + cfskappa.min = float(kappamin) + cfskappa.max = float(kappamax) + cfssigma = CFSParameter() + cfssigma.ID = 'sigma' + cfssigma.scalingprofile = sigmascalingprofile + cfssigma.scalingdirection = sigmascalingdirection + cfssigma.min = float(sigmamin) + if sigmamax == 'None': + cfssigma.max = None + else: + cfssigma.max = float(sigmamax) + cfs = CFS() + cfs.alpha = cfsalpha + cfs.kappa = cfskappa + cfs.sigma = cfssigma + + if grid.messages: + print('PML CFS parameters: alpha (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), kappa (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), sigma (scaling: {}, scaling direction: {}, min: {:g}, max: {}) created.'.format(cfsalpha.scalingprofile, cfsalpha.scalingdirection, cfsalpha.min, cfsalpha.max, cfskappa.scalingprofile, cfskappa.scalingdirection, cfskappa.min, cfskappa.max, cfssigma.scalingprofile, cfssigma.scalingdirection, cfssigma.min, cfssigma.max)) + + grid.cfs.append(cfs) + +class Subgrid(UserObjectMulti): + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.children_multiple = [] + self.children_geometry = [] + + def add(self, node): + if isinstance(node, UserObjectMulti): + self.children_multiple.append(node) + elif isinstance(node, UserObjectGeometry): + self.children_geometry.append(node) + else: + raise Exception('This Object is Unknown to gprMax') + +class SubgridHSG(UserObjectMulti): + + def __init__(self, **kwargs): + super().__init__(**kwargs) diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py new file mode 100644 index 00000000..dd9b6d4e --- /dev/null +++ b/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') diff --git a/gprMax/cython/fields_updates_hsg.pyx b/gprMax/cython/fields_updates_hsg.pyx new file mode 100644 index 00000000..33caadbd --- /dev/null +++ b/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 diff --git a/gprMax/scene.py b/gprMax/scene.py new file mode 100644 index 00000000..cc2a4983 --- /dev/null +++ b/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 diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py new file mode 100644 index 00000000..a73ab474 --- /dev/null +++ b/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 diff --git a/gprMax/subgrids/factory.py b/gprMax/subgrids/factory.py new file mode 100644 index 00000000..cd086464 --- /dev/null +++ b/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 diff --git a/gprMax/subgrids/hsg.py b/gprMax/subgrids/hsg.py new file mode 100644 index 00000000..5b37ee00 --- /dev/null +++ b/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 diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py new file mode 100644 index 00000000..3fc7c6d5 --- /dev/null +++ b/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]] diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py new file mode 100644 index 00000000..62ad67c0 --- /dev/null +++ b/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 diff --git a/temp/gprMax b/temp/gprMax new file mode 160000 index 00000000..411a3f9e --- /dev/null +++ b/temp/gprMax @@ -0,0 +1 @@ +Subproject commit 411a3f9ec8da465271cf8e75057a597fe52369d4 diff --git a/todo b/todo new file mode 100644 index 00000000..f9902507 --- /dev/null +++ b/todo @@ -0,0 +1 @@ +memory update for all grid (model_build_run.py)