publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-07 23:14:03 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Updates to docstrings and formatting.

浏览代码
这个提交包含在：
Craig Warren
2022-11-04 16:17:20 +00:00
父节点 c56ab30876
当前提交 69e3fe55c5
共有 29 个文件被更改，包括 1301 次插入 和 1011 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

34
gprMax/cmds_geometry/add_grass.py
查看文件

@@ -29,20 +29,20 @@ logger = logging.getLogger(__name__)
class AddGrass(UserObjectGeometry):
"""Allows you to add grass with roots to a FractalBox class in the model.
"""Adds grass with roots to a FractalBox class in the model.
Attributes:
p1: a list of the lower left (x,y,z) coordinates of a surface on a
p1: list of the lower left (x,y,z) coordinates of a surface on a
FractalBox class.
p2: a list of the upper right (x,y,z) coordinates of a surface on a
p2: list of the upper right (x,y,z) coordinates of a surface on a
FractalBox class.
frac_dim: a float for the fractal dimension which, for an orthogonal
frac_dim:float for the fractal dimension which, for an orthogonal
parallelepiped, should take values between zero and three.
limits: a list to define lower and upper limits for a range over which
limits: list to define lower and upper limits for a range over which
the height of the blades of grass can vary.
n_blades: an int for the number of blades of grass that should be
n_blades:int for the number of blades of grass that should be
applied to the surface area.
fractal_box_id: a string identifier for the FractalBox class that the
fractal_box_id: string identifier for the FractalBox class that the
grass should be applied to.
"""
@@ -191,7 +191,8 @@ class AddGrass(UserObjectGeometry):
'enough for the number of grass blades/roots specified')
raise ValueError
# Scale the distribution so that the summation is equal to one, i.e. a probability distribution
# 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
@@ -202,16 +203,23 @@ class AddGrass(UserObjectGeometry):
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]))
# 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]))
# 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)
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)

26
gprMax/cmds_geometry/add_surface_roughness.py
查看文件

@@ -28,20 +28,20 @@ logger = logging.getLogger(__name__)
class AddSurfaceRoughness(UserObjectGeometry):
"""Allows you to add grass with roots to a FractalBox class in the model.
"""Adds surface roughness to a FractalBox class in the model.
Attributes:
p1: a list of the lower left (x,y,z) coordinates of a surface on a
p1: list of the lower left (x,y,z) coordinates of a surface on a
FractalBox class.
p2: a list of the upper right (x,y,z) coordinates of a surface on a
p2: list of the upper right (x,y,z) coordinates of a surface on a
FractalBox class.
frac_dim: a float for the fractal dimension which, for an orthogonal
frac_dim: float for the fractal dimension which, for an orthogonal
parallelepiped, should take values between zero and three.
weighting: a list with weightings in the first and second direction of
weighting: list with weightings in the first and second direction of
the surface.
limits: a list to define lower and upper limits for a range over which
limits: ist to define lower and upper limits for a range over which
the surface roughness can vary.
fractal_box_id: a string identifier for the FractalBox class
fractal_box_id: string identifier for the FractalBox class
that the surface roughness should be applied to.
seed: (optional) float parameter which controls the seeding of the random
number generator used to create the fractals.
@@ -119,7 +119,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
logger.exception(self.__str__() + ' can only be used on the external ' +
'surfaces of a fractal box')
raise ValueError
fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
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:
@@ -147,7 +148,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
logger.exception(self.__str__() + ' can only be used on the external ' +
'surfaces of a fractal box')
raise ValueError
fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
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:
@@ -175,7 +177,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
logger.exception(self.__str__() + ' can only be used on the external ' +
'surfaces of a fractal box')
raise ValueError
fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
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:
@@ -212,7 +215,8 @@ class AddSurfaceRoughness(UserObjectGeometry):
# List of existing surfaces IDs
existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
if surface.surfaceID in existingsurfaceIDs:
logger.exception(self.__str__() + f' has already been used on the {surface.surfaceID} surface')
logger.exception(self.__str__() + f' has already been used on the ' +
f'{surface.surfaceID} surface')
raise ValueError
surface.generate_fractal_surface()

16
gprMax/cmds_geometry/add_surface_water.py
查看文件

@@ -28,17 +28,17 @@ logger = logging.getLogger(__name__)
class AddSurfaceWater(UserObjectGeometry):
"""Allows you to add surface water to a FractalBox class in the model.
"""Adds surface water to a FractalBox class in the model.
Attributes:
p1: a list of the lower left (x,y,z) coordinates of a surface on a
p1: list of the lower left (x,y,z) coordinates of a surface on a
FractalBox class.
p2: a list of the upper right (x,y,z) coordinates of a surface on a
p2: list of the upper right (x,y,z) coordinates of a surface on a
FractalBox class.
depth: a float that defines the depth of the water, which should be
depth: float that defines the depth of the water, which should be
specified relative to the dimensions of the #fractal_box that
the surface water is being applied.
fractal_box_id: a string identifier for the FractalBox class that the
fractal_box_id: string identifier for the FractalBox class that the
surface water should be applied to.
"""
@@ -87,7 +87,8 @@ class AddSurfaceWater(UserObjectGeometry):
xf, yf, zf = p2
if depth <= 0:
logger.exception(self.__str__() + ' requires a positive value for the depth of water')
logger.exception(self.__str__() + ' requires a positive value for ' +
'the depth of water')
raise ValueError
# Check for valid orientations
@@ -148,7 +149,8 @@ class AddSurfaceWater(UserObjectGeometry):
surface = next((x for x in volume.fractalsurfaces if x.surfaceID == requestedsurface), None)
if not surface:
logger.exception(self.__str__() + f' specified surface {requestedsurface} does not have a rough surface applied')
logger.exception(self.__str__() + f' specified surface {requestedsurface} ' +
'does not have a rough surface applied')
raise ValueError
surface.filldepth = filldepthcells

29
gprMax/cmds_geometry/box.py
查看文件

@@ -28,16 +28,16 @@ logger = logging.getLogger(__name__)
class Box(UserObjectGeometry):
"""Allows you to introduce an orthogonal parallelepiped with specific
properties into the model.
"""Introduces an orthogonal parallelepiped with specific properties into
the model.
Attributes:
p1: a list of the lower left (x,y,z) coordinates of the parallelepiped.
p2: a list of the upper right (x,y,z) coordinates of the parallelepiped.
material_id: a string for the material identifier that must correspond
p1: list of the lower left (x,y,z) coordinates of the parallelepiped.
p2: list of the upper right (x,y,z) coordinates of the parallelepiped.
material_id: string for the material identifier that must correspond
to material that has already been defined.
material_ids: a list of material identifiers in the x, y, z directions.
averaging: a string (y or n) used to switch on and off dielectric smoothing.
material_ids: list of material identifiers in the x, y, z directions.
averaging: string (y or n) used to switch on and off dielectric smoothing.
"""
def __init__(self, **kwargs):
@@ -124,15 +124,20 @@ class Box(UserObjectGeometry):
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)
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)
build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz,
averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
dielectricsmoothing = 'on' if averaging else 'off'

42
gprMax/cmds_geometry/cmds_geometry.py
查看文件

@@ -45,15 +45,15 @@ class UserObjectGeometry:
return f'{self.hash}: {s[:-1]}'
def create(self, grid, uip):
"""Create object and add it to the grid."""
"""Creates object and adds it to the grid."""
pass
def rotate(self, axis, angle, origin=None):
"""Rotate object - specialised for each object."""
"""Rotates object - specialised for each object."""
pass
def grid_name(self, grid):
"""Return subgrid name for use with logging info. Return an empty
"""Returns subgrid name for use with logging info. Returns an empty
string if the grid is the main grid.
"""
@@ -64,16 +64,16 @@ class UserObjectGeometry:
def rotate_point(p, axis, angle, origin=(0, 0, 0)):
"""Rotate a point.
"""Rotates a point.
Args:
p: an array of coordinates of point (x, y, z).
axis: a string which defines the axis about which to perform rotation (x, y, or z).
angle: an int specifying the angle of rotation (degrees).
origin: a tuple defining the point about which to perform rotation (x, y, z).
p: array of coordinates of point (x, y, z).
axis: string which defines the axis about which to perform rotation (x, y, or z).
angle: int specifying the angle of rotation (degrees).
origin: tuple defining the point about which to perform rotation (x, y, z).
Returns:
p: an array of coordinates of rotated point (x, y, z)
p: array of coordinates of rotated point (x, y, z)
"""
origin = np.array(origin)
@@ -97,13 +97,13 @@ def rotate_2point_object(pts, axis, angle, origin=None):
"""Rotate a geometry object that is defined by 2 points.
Args:
pts: an array ofcoordinates of points of object to be rotated.
axis: a string which defines the axis about which to perform rotation (x, y, or z).
angle: an int specifying the angle of rotation (degrees).
origin: a tuple defining the point about which to perform rotation (x, y, z).
pts: array ofcoordinates of points of object to be rotated.
axis: string which defines the axis about which to perform rotation (x, y, or z).
angle: int specifying the angle of rotation (degrees).
origin: tuple defining the point about which to perform rotation (x, y, z).
Returns:
new_pts: an array of coordinates of points of rotated object.
new_pts: array of coordinates of points of rotated object.
"""
# Use origin at centre of object if not given
@@ -152,18 +152,18 @@ def rotate_2point_object(pts, axis, angle, origin=None):
def rotate_polarisation(p, polarisation, axis, angle, G):
"""Rotate a geometry object that is defined by a point and polarisation.
"""Rotates a geometry object that is defined by a point and polarisation.
Args:
p: an array of coordinates of point (x, y, z).
polarisation: a string defining the current polarisation (x, y, or z).
axis: a string which defines the axis about which to perform rotation (x, y, or z).
angle: an int specifying the angle of rotation (degrees).
p: array of coordinates of point (x, y, z).
polarisation: string defining the current polarisation (x, y, or z).
axis: string which defines the axis about which to perform rotation (x, y, or z).
angle: int specifying the angle of rotation (degrees).
G: FDTDGrid class describing a grid in a model.
Returns:
pts: an array of coordinates of points of rotated object.
new_polarisation: a string defining the new polarisation (x, y, or z).
pts: array of coordinates of points of rotated object.
new_polarisation: string defining the new polarisation (x, y, or z).
"""
if polarisation.lower() == 'x':

24
gprMax/cmds_geometry/cylinder.py
查看文件

@@ -28,18 +28,18 @@ logger = logging.getLogger(__name__)
class Cylinder(UserObjectGeometry):
"""Allows you to introduce a circular cylinder into the model.
"""Introduces a circular cylinder into the model.
Attributes:
p1: a list of the coordinates (x,y,z) of the centre of the first face
p1: list of the coordinates (x,y,z) of the centre of the first face
of the cylinder.
p2: a list of the coordinates (x,y,z) of the centre of the second face
p2: list of the coordinates (x,y,z) of the centre of the second face
of the cylinder.
r: float of the radius of the cylinder.
material_id: a string for the material identifier that must correspond
material_id: string for the material identifier that must correspond
to material that has already been defined.
material_ids: a list of material identifiers in the x, y, z directions.
averaging: a string (y or n) used to switch on and off dielectric smoothing.
material_ids: list of material identifiers in the x, y, z directions.
averaging: string (y or n) used to switch on and off dielectric smoothing.
"""
def __init__(self, **kwargs):
@@ -113,10 +113,14 @@ class Cylinder(UserObjectGeometry):
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)
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)

37
gprMax/cmds_geometry/cylindrical_sector.py
查看文件

@@ -28,28 +28,27 @@ logger = logging.getLogger(__name__)
class CylindricalSector(UserObjectGeometry):
"""Allows you to introduce a cylindrical sector (shaped like a slice of pie)
into the model.
"""Introduces a cylindrical sector (shaped like a slice of pie) into the model.
Attributes:
normal: a string for the direction of the axis of the cylinder from which
normal: string for the direction of the axis of the cylinder from which
the sector is defined and can be x, y, or z.
ctr1: a float for the first coordinate of the centre of the cylindrical
ctr1: float for the first coordinate of the centre of the cylindrical
sector.
ctr2: a float for the second coordinate of the centre of the cylindrical
ctr2: float for the second coordinate of the centre of the cylindrical
sector.
extent1: a float for the first thickness from the centre of the
extent1: float for the first thickness from the centre of the
cylindrical sector.
extent2: a float for the second thickness from the centre of the
extent2: float for the second thickness from the centre of the
cylindrical sector.
r: a float for the radius of the cylindrical sector.
start: a float for the starting angle (in degrees) for the cylindrical
r: float for the radius of the cylindrical sector.
start: float for the starting angle (in degrees) for the cylindrical
sector.
end: a float for the angle (in degrees) swept by the cylindrical sector.
material_id: a string for the material identifier that must correspond
end: float for the angle (in degrees) swept by the cylindrical sector.
material_id: string for the material identifier that must correspond
to material that has already been defined.
material_ids: a list of material identifiers in the x, y, z directions.
averaging: a string (y or n) used to switch on and off dielectric smoothing.
material_ids: list of material identifiers in the x, y, z directions.
averaging: string (y or n) used to switch on and off dielectric smoothing.
"""
def __init__(self, **kwargs):
@@ -138,10 +137,14 @@ class CylindricalSector(UserObjectGeometry):
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)
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)

12
gprMax/cmds_geometry/edge.py
查看文件

@@ -28,12 +28,12 @@ logger = logging.getLogger(__name__)
class Edge(UserObjectGeometry):
"""Allows you to introduce a wire with specific properties into the model.
"""Introduces a wire with specific properties into the model.
Attributes:
p1: a list of the coordinates (x,y,z) of the starting point of the edge.
p2: a list of the coordinates (x,y,z) of the ending point of the edge.
material_id: a string for the material identifier that must correspond
p1: list of the coordinates (x,y,z) of the starting point of the edge.
p2: list of the coordinates (x,y,z) of the ending point of the edge.
material_id: string for the material identifier that must correspond
to material that has already been defined.
"""
@@ -49,14 +49,14 @@ class Edge(UserObjectGeometry):
self.dorotate = True
def __dorotate(self):
"""Perform rotation."""
"""Performs rotation."""
pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
self.kwargs['p1'] = tuple(rot_pts[0, :])
self.kwargs['p2'] = tuple(rot_pts[1, :])
def create(self, grid, uip):
"""Create edge and add it to the grid."""
"""Creates edge and adds it to the grid."""
try:
p1 = self.kwargs['p1']
p2 = self.kwargs['p2']

37
gprMax/cmds_geometry/fractal_box.py
查看文件

@@ -27,27 +27,27 @@ logger = logging.getLogger(__name__)
class FractalBox(UserObjectGeometry):
"""Allows you to introduce an orthogonal parallelepiped with fractal
distributed properties which are related to a mixing model or
normal material into the model.
"""Introduces an orthogonal parallelepiped with fractal distributed
properties which are related to a mixing model or normal material into
the model.
Attributes:
p1: a list of the lower left (x,y,z) coordinates of the parallelepiped.
p2: a list of the upper right (x,y,z) coordinates of the parallelepiped.
frac_dim: a float for the fractal dimension which, for an orthogonal
p1: list of the lower left (x,y,z) coordinates of the parallelepiped.
p2: list of the upper right (x,y,z) coordinates of the parallelepiped.
frac_dim: float for the fractal dimension which, for an orthogonal
parallelepiped, should take values between zero and three.
weighting: a list of the weightings in the x, y, z direction of the
weighting: list of the weightings in the x, y, z direction of the
parallelepiped.
n_materials: an int of the number of materials to use for the fractal
n_materials: int of the number of materials to use for the fractal
distribution (defined according to the associated
mixing model). This should be set to one if using a
normal material instead of a mixing model.
mixing_model_id: a string identifier for the associated mixing model or
mixing_model_id: string identifier for the associated mixing model or
material.
id: a string identifier for the fractal box itself.
seed: seed: (optional) float parameter which controls the seeding of the random
number generator used to create the fractals.
averaging: a string (y or n) used to switch on and off dielectric smoothing.
id: string identifier for the fractal box itself.
seed: (optional) float parameter which controls the seeding of the
random number generator used to create the fractals.
averaging: string (y or n) used to switch on and off dielectric smoothing.
"""
def __init__(self, **kwargs):
@@ -62,7 +62,7 @@ class FractalBox(UserObjectGeometry):
self.dorotate = True
def __dorotate(self):
"""Perform rotation."""
"""Performs rotation."""
pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
self.kwargs['p1'] = tuple(rot_pts[0, :])
@@ -94,7 +94,8 @@ class FractalBox(UserObjectGeometry):
# Go with user specified averaging
averagefractalbox = self.kwargs['averaging']
except KeyError:
# If they havent specfied - default is no dielectric smoothing for a fractal box
# If they havent specified - default is no dielectric smoothing for
# a fractal box.
averagefractalbox = False
p3 = uip.round_to_grid_static_point(p1)
@@ -124,7 +125,8 @@ class FractalBox(UserObjectGeometry):
'number of bins')
raise ValueError
# Find materials to use to build fractal volume, either from mixing models or normal materials
# 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
@@ -134,7 +136,8 @@ class FractalBox(UserObjectGeometry):
logger.exception(self.__str__() + ' must be used with more than ' +
'one material from the mixing model.')
raise ValueError
# Create materials from mixing model as number of bins now known from fractal_box command
# Create materials from mixing model as number of bins now known
# from fractal_box command.
mixingmodel.calculate_debye_properties(nbins, grid)
elif not material:
logger.exception(self.__str__() + f' mixing model or material with ' +

18
gprMax/cmds_geometry/fractal_box_builder.py
查看文件

@@ -124,7 +124,8 @@ class FractalBoxBuilder(UserObjectGeometry):
# 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 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:
@@ -146,7 +147,8 @@ class FractalBoxBuilder(UserObjectGeometry):
# 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 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:
@@ -192,7 +194,8 @@ class FractalBoxBuilder(UserObjectGeometry):
# 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 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:
@@ -214,7 +217,8 @@ class FractalBoxBuilder(UserObjectGeometry):
# 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 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:
@@ -260,7 +264,8 @@ class FractalBoxBuilder(UserObjectGeometry):
# 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 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:
@@ -282,7 +287,8 @@ class FractalBoxBuilder(UserObjectGeometry):
# 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 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:

5
gprMax/cmds_geometry/geometry_objects_read.py
查看文件

@@ -40,7 +40,7 @@ class GeometryObjectsRead(UserObjectGeometry):
pass
def create(self, grid, uip):
"""Create the object and add it to the grid."""
"""Creates the object and adds it to the grid."""
try:
p1 = self.kwargs['p1']
geofile = self.kwargs['geofile']
@@ -88,7 +88,8 @@ class GeometryObjectsRead(UserObjectGeometry):
else:
material.type = 'imported'
# See if geometry object file exists at specified path and if not try input file directory
# See if geometry object file exists at specified path and if not try
# input file directory.
geofile = Path(geofile)
if not geofile.exists():
geofile = Path(config.sim_config.input_file_path.parent, geofile)

12
gprMax/cmds_geometry/plate.py
查看文件

@@ -28,14 +28,14 @@ logger = logging.getLogger(__name__)
class Plate(UserObjectGeometry):
"""Allows you to introduce a plate with specific properties into the model.
"""Introduces a plate with specific properties into the model.
Attributes:
p1: a list of the lower left (x,y,z) coordinates of the plate.
p2: a list of the upper right (x,y,z) coordinates of the plate.
material_id: a string for the material identifier that must correspond
p1: list of the lower left (x,y,z) coordinates of the plate.
p2: list of the upper right (x,y,z) coordinates of the plate.
material_id: string for the material identifier that must correspond
to material that has already been defined.
material_ids: a list of material identifiers in the x, y, z directions.
material_ids: list of material identifiers in the x, y, z directions.
"""
def __init__(self, **kwargs):
@@ -50,7 +50,7 @@ class Plate(UserObjectGeometry):
self.dorotate = True
def __dorotate(self):
"""Perform rotation."""
"""Performs rotation."""
pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
self.kwargs['p1'] = tuple(rot_pts[0, :])

29
gprMax/cmds_geometry/sphere.py
查看文件

@@ -28,16 +28,15 @@ logger = logging.getLogger(__name__)
class Sphere(UserObjectGeometry):
"""Allows you to introduce a spherical object with specific parameters
into the model.
"""Introduces a spherical object with specific parameters into the model.
Attributes:
p1: a list of the coordinates (x,y,z) of the centre of the sphere.
r: a float of radius of the sphere.
material_id: a string for the material identifier that must correspond
p1: list of the coordinates (x,y,z) of the centre of the sphere.
r: float of radius of the sphere.
material_id: string for the material identifier that must correspond
to material that has already been defined.
material_ids: a list of material identifiers in the x, y, z directions.
averaging: a string (y or n) used to switch on and off dielectric smoothing.
material_ids: list of material identifiers in the x, y, z directions.
averaging: string (y or n) used to switch on and off dielectric smoothing.
"""
def __init__(self, **kwargs):
@@ -104,15 +103,21 @@ class Sphere(UserObjectGeometry):
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)
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)
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)
dielectricsmoothing = 'on' if averaging else 'off'
logger.info(self.grid_name(grid) + f"Sphere with centre {p2[0]:g}m, " +

34
gprMax/cmds_geometry/triangle.py
查看文件

@@ -28,19 +28,19 @@ logger = logging.getLogger(__name__)
class Triangle(UserObjectGeometry):
"""Allows you to introduce a triangular patch or a triangular prism with
specific properties into the model.
"""Introduces a triangular patch or a triangular prism with specific
properties into the model.
Attributes:
p1: a list of the coordinates (x,y,z) of the first apex of the triangle.
p2: a list of the coordinates (x,y,z) of the second apex of the triangle.
p3: a list of the coordinates (x,y,z) of the third apex of the triangle.
thickness: a float for the thickness of the triangular prism. If the
p1: list of the coordinates (x,y,z) of the first apex of the triangle.
p2: list of the coordinates (x,y,z) of the second apex of the triangle.
p3: list of the coordinates (x,y,z) of the third apex of the triangle.
thickness: float for the thickness of the triangular prism. If the
thickness is zero then a triangular patch is created.
material_id: a string for the material identifier that must correspond
material_id: string for the material identifier that must correspond
to material that has already been defined.
material_ids: a list of material identifiers in the x, y, z directions.
averaging: a string (y or n) used to switch on and off dielectric smoothing.
material_ids: list of material identifiers in the x, y, z directions.
averaging: string (y or n) used to switch on and off dielectric smoothing.
"""
def __init__(self, **kwargs):
@@ -48,14 +48,14 @@ class Triangle(UserObjectGeometry):
self.hash = '#triangle'
def rotate(self, axis, angle, origin=None):
"""Set parameters for rotation."""
"""Sets parameters for rotation."""
self.axis = axis
self.angle = angle
self.origin = origin
self.dorotate = True
def __dorotate(self):
"""Perform rotation."""
"""Performs rotation."""
p1 = rotate_point(self.kwargs['p1'], self.axis, self.angle, self.origin)
p2 = rotate_point(self.kwargs['p2'], self.axis, self.angle, self.origin)
p3 = rotate_point(self.kwargs['p3'], self.axis, self.angle, self.origin)
@@ -154,10 +154,14 @@ class Triangle(UserObjectGeometry):
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)
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)

144
gprMax/cython/fields_updates_dispersive_template.jinja
查看文件

@@ -39,26 +39,26 @@ cdef extern from "complex.h" nogil:
{% for item in functions %}
cpdef void {{ item.name_a }}(
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.field_type }}[:, ::1] updatecoeffsE,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez,
{{ item.field_type }}[:, :, ::1] Hx,
{{ item.field_type }}[:, :, ::1] Hy,
{{ item.field_type }}[:, :, ::1] Hz
):
"""This function updates the electric field components when dispersive
materials (with multiple poles) are present.
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.field_type }}[:, ::1] updatecoeffsE,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez,
{{ item.field_type }}[:, :, ::1] Hx,
{{ item.field_type }}[:, :, ::1] Hy,
{{ item.field_type }}[:, :, ::1] Hz
):
"""Updates the electric field components when dispersive materials
(with multiple poles) are present.
Args:
nx, ny, nz: int for grid size in cells.
@@ -158,22 +158,22 @@ cpdef void {{ item.name_a }}(
{% for item in functions %}
cpdef void {{ item.name_b }}(
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez
):
"""This function updates a temporary dispersive material array when
disperisive materials (with multiple poles) are present.
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez
):
"""Updates a temporary dispersive material array when disperisive materials
(with multiple poles) are present.
Args:
nx, ny, nz: int for grid size in cells.
@@ -228,26 +228,26 @@ cpdef void {{ item.name_b }}(
{% for item in functions %}
cpdef void {{ item.name_a_1 }}(
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.field_type }}[:, ::1] updatecoeffsE,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez,
{{ item.field_type }}[:, :, ::1] Hx,
{{ item.field_type }}[:, :, ::1] Hy,
{{ item.field_type }}[:, :, ::1] Hz
):
"""This function updates the electric field components when dispersive
materials (with 1 pole) are present.
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.field_type }}[:, ::1] updatecoeffsE,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez,
{{ item.field_type }}[:, :, ::1] Hx,
{{ item.field_type }}[:, :, ::1] Hy,
{{ item.field_type }}[:, :, ::1] Hz
):
"""Updates the electric field components when dispersive materials
(with 1 pole) are present.
Args:
nx, ny, nz: int for grid size in cells.
@@ -339,22 +339,22 @@ cpdef void {{ item.name_a_1 }}(
{% for item in functions %}
cpdef void {{ item.name_b_1 }}(
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez
):
"""This function updates a temporary dispersive material array when
disperisive materials (with 1 pole) are present.
int nx,
int ny,
int nz,
int nthreads,
int maxpoles,
{{ item.dispersive_type }}[:, ::1] updatecoeffsdispersive,
np.uint32_t[:, :, :, ::1] ID,
{{ item.dispersive_type }}[:, :, :, ::1] Tx,
{{ item.dispersive_type }}[:, :, :, ::1] Ty,
{{ item.dispersive_type }}[:, :, :, ::1] Tz,
{{ item.field_type }}[:, :, ::1] Ex,
{{ item.field_type }}[:, :, ::1] Ey,
{{ item.field_type }}[:, :, ::1] Ez
):
"""Updates a temporary dispersive material array when disperisive materials
(with 1 pole) are present.
Args:
nx, ny, nz: int for grid size in cells.

449
gprMax/cython/fields_updates_hsg.pyx
查看文件

@@ -22,262 +22,267 @@ 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
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:
"""
Args:
subgrid: Subgrid class instance.
n: string for the normal to the face to update.
nwn: int for the number of working cell in the normal direction
to the face.
lookup_id: int id of the H component to update at each node.
field: nparray of main grid field to be updated.
inc_field: nparray of 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: boolean for the H node midway along the lower edge?
r: int for subgrid ratio.
s: int for separation of inner and outer surfaces.
nb: int for number of boundary cells.
"""
# Comments here as as per left and right face
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
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
# 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):
# 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
# Calculate the subgrid j component of the H nodes
# i.e. Hz node of the left or right face
for m in range(m_l, m_u):
# Calculate the subgrid k component of the H nodes
if mid == 1:
l_s = os + (l - l_l) * r + r // 2
# i.e. the Hy node of the left or right face
m_s = os + (m - m_l) * r
else:
l_s = os + (l - l_l) * r
m_s = os + (m - m_l) * r + r // 2
for m in range(m_l, m_u):
# 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
# 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
# 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
# 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
# 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
# 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
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
"""
"""
Args:
r: int for subgrid ratio.
s: int for separation of inner and outer surfaces.
nb: int for number of 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
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
# 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
# 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):
for l in prange(l_l, l_u, nogil=True, schedule='static', num_threads=nthreads):
# y coord of the Ex field component
# 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:
l_s = os + (l - l_l) * r + r // 2
# y coord of the Ez field component
m_s = os + (m - m_l) * r
else:
l_s = os + (l - l_l) * r
m_s = os + (m - m_l) * r + r // 2
for m in range(m_l, m_u):
# associate the given indices with their i, j, k values
# 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
# 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
# associate the given indices with their i, j, k values
material_e_l = ID[lookup_id, i0, j0, k0]
inc_n = inc_field[i1, j1, k1] * sign_n
# 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
# make sure these are the correct grid
field[i0, j0, k0] += updatecoeffsH[material_e_l, co] * inc_n
material_e_l = ID[lookup_id, i0, j0, k0]
inc_n = inc_field[i1, j1, k1] * sign_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
# 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,
np.float64_t[:, :] updatecoeffsE,
np.uint32_t[:, :, :, :] ID,
int n,
int offset,
int nwl,
int nwm,
int nwn,
int face,
np.float64_t[:, :, :] field,
np.float64_t[:, :] inc_field_l,
np.float64_t[:, :] inc_field_u,
Py_ssize_t lookup_id,
int sign_l,
int sign_u,
Py_ssize_t co,
int nthreads
int nwx,
int nwy,
int nwz,
np.float64_t[:, :] updatecoeffsE,
np.uint32_t[:, :, :, :] ID,
int n,
int offset,
int nwl,
int nwm,
int nwn,
int face,
np.float64_t[:, :, :] field,
np.float64_t[:, :] inc_field_l,
np.float64_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
"""
Args:
"""
for l in prange(n, nwl + n, nogil=True, schedule='static', num_threads=nthreads):
for m in range(n, nwm + n):
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
# 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
for l in prange(n, nwl + n, nogil=True, schedule='static', num_threads=nthreads):
for m in range(n, nwm + n):
inc_i = l - n
inc_j = m - 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
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
inc_i = l - n
inc_j = m - n
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
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

136
gprMax/cython/fields_updates_normal.pyx
查看文件

@@ -23,29 +23,27 @@ from cython.parallel import prange
from gprMax.config cimport float_or_double
###############################################
# Electric field updates - standard materials #
###############################################
cpdef void update_electric(
int nx,
int ny,
int nz,
int nthreads,
float_or_double[:, ::1] updatecoeffsE,
np.uint32_t[:, :, :, ::1] ID,
float_or_double[:, :, ::1] Ex,
float_or_double[:, :, ::1] Ey,
float_or_double[:, :, ::1] Ez,
float_or_double[:, :, ::1] Hx,
float_or_double[:, :, ::1] Hy,
float_or_double[:, :, ::1] Hz
):
"""This function updates the electric field components.
int nx,
int ny,
int nz,
int nthreads,
float_or_double[:, ::1] updatecoeffsE,
np.uint32_t[:, :, :, ::1] ID,
float_or_double[:, :, ::1] Ex,
float_or_double[:, :, ::1] Ey,
float_or_double[:, :, ::1] Ez,
float_or_double[:, :, ::1] Hx,
float_or_double[:, :, ::1] Hy,
float_or_double[:, :, ::1] Hz
):
"""Updates the electric field components (for standard materials).
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
nx, ny, nz: ints for grid size in cells.
nthreads: int for number of threads to use.
updatecoeffs, ID, E, H: memoryviews to access update coefficients,
ID and field component arrays
"""
cdef Py_ssize_t i, j, k
@@ -57,7 +55,9 @@ cpdef void update_electric(
for j in range(1, ny):
for k in range(1, nz):
materialEx = ID[0, i, j, k]
Ex[i, j, k] = updatecoeffsE[materialEx, 0] * Ex[i, j, k] + updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1])
Ex[i, j, k] = (updatecoeffsE[materialEx, 0] * Ex[i, j, k] +
updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) -
updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]))
# 2D - Ey component
elif ny == 1:
@@ -65,7 +65,9 @@ cpdef void update_electric(
for j in range(0, ny):
for k in range(1, nz):
materialEy = ID[1, i, j, k]
Ey[i, j, k] = updatecoeffsE[materialEy, 0] * Ey[i, j, k] + updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
Ey[i, j, k] = (updatecoeffsE[materialEy, 0] * Ey[i, j, k] +
updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) -
updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]))
# 2D - Ez component
elif nz == 1:
@@ -73,7 +75,9 @@ cpdef void update_electric(
for j in range(1, ny):
for k in range(0, nz):
materialEz = ID[2, i, j, k]
Ez[i, j, k] = updatecoeffsE[materialEz, 0] * Ez[i, j, k] + updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
Ez[i, j, k] = (updatecoeffsE[materialEz, 0] * Ez[i, j, k] +
updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) -
updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]))
# 3D
else:
@@ -83,52 +87,62 @@ cpdef void update_electric(
materialEx = ID[0, i, j, k]
materialEy = ID[1, i, j, k]
materialEz = ID[2, i, j, k]
Ex[i, j, k] = updatecoeffsE[materialEx, 0] * Ex[i, j, k] + updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1])
Ey[i, j, k] = updatecoeffsE[materialEy, 0] * Ey[i, j, k] + updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
Ez[i, j, k] = updatecoeffsE[materialEz, 0] * Ez[i, j, k] + updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
Ex[i, j, k] = (updatecoeffsE[materialEx, 0] * Ex[i, j, k] +
updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) -
updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]))
Ey[i, j, k] = (updatecoeffsE[materialEy, 0] * Ey[i, j, k] +
updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) -
updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]))
Ez[i, j, k] = (updatecoeffsE[materialEz, 0] * Ez[i, j, k] +
updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) -
updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k]))
# Ex components at i = 0
for j in prange(1, ny, nogil=True, schedule='static', num_threads=nthreads):
for k in range(1, nz):
materialEx = ID[0, 0, j, k]
Ex[0, j, k] = updatecoeffsE[materialEx, 0] * Ex[0, j, k] + updatecoeffsE[materialEx, 2] * (Hz[0, j, k] - Hz[0, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[0, j, k] - Hy[0, j, k - 1])
Ex[0, j, k] = (updatecoeffsE[materialEx, 0] * Ex[0, j, k] +
updatecoeffsE[materialEx, 2] * (Hz[0, j, k] - Hz[0, j - 1, k]) -
updatecoeffsE[materialEx, 3] * (Hy[0, j, k] - Hy[0, j, k - 1]))
# Ey components at j = 0
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for k in range(1, nz):
materialEy = ID[1, i, 0, k]
Ey[i, 0, k] = updatecoeffsE[materialEy, 0] * Ey[i, 0, k] + updatecoeffsE[materialEy, 3] * (Hx[i, 0, k] - Hx[i, 0, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, 0, k] - Hz[i - 1, 0, k])
Ey[i, 0, k] = (updatecoeffsE[materialEy, 0] * Ey[i, 0, k] +
updatecoeffsE[materialEy, 3] * (Hx[i, 0, k] - Hx[i, 0, k - 1]) -
updatecoeffsE[materialEy, 1] * (Hz[i, 0, k] - Hz[i - 1, 0, k]))
# Ez components at k = 0
for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(1, ny):
materialEz = ID[2, i, j, 0]
Ez[i, j, 0] = updatecoeffsE[materialEz, 0] * Ez[i, j, 0] + updatecoeffsE[materialEz, 1] * (Hy[i, j, 0] - Hy[i - 1, j, 0]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, 0] - Hx[i, j - 1, 0])
Ez[i, j, 0] = (updatecoeffsE[materialEz, 0] * Ez[i, j, 0] +
updatecoeffsE[materialEz, 1] * (Hy[i, j, 0] - Hy[i - 1, j, 0]) -
updatecoeffsE[materialEz, 2] * (Hx[i, j, 0] - Hx[i, j - 1, 0]))
##########################
# Magnetic field updates #
##########################
cpdef void update_magnetic(
int nx,
int ny,
int nz,
int nthreads,
float_or_double[:, ::1] updatecoeffsH,
np.uint32_t[:, :, :, ::1] ID,
float_or_double[:, :, ::1] Ex,
float_or_double[:, :, ::1] Ey,
float_or_double[:, :, ::1] Ez,
float_or_double[:, :, ::1] Hx,
float_or_double[:, :, ::1] Hy,
float_or_double[:, :, ::1] Hz
):
"""This function updates the magnetic field components.
int nx,
int ny,
int nz,
int nthreads,
float_or_double[:, ::1] updatecoeffsH,
np.uint32_t[:, :, :, ::1] ID,
float_or_double[:, :, ::1] Ex,
float_or_double[:, :, ::1] Ey,
float_or_double[:, :, ::1] Ez,
float_or_double[:, :, ::1] Hx,
float_or_double[:, :, ::1] Hy,
float_or_double[:, :, ::1] Hz
):
"""Updates the magnetic field components.
Args:
nx, ny, nz (int): Grid size in cells
nthreads (int): Number of threads to use
updatecoeffs, ID, E, H (memoryviews): Access to update coeffients, ID and field component arrays
nx, ny, nz: ints for grid size in cells.
nthreads: int for number of threads to use.
updatecoeffs, ID, E, H: memoryviews to access update coefficients,
ID and field component arrays
"""
cdef Py_ssize_t i, j, k
@@ -142,7 +156,9 @@ cpdef void update_magnetic(
for j in range(0, ny):
for k in range(0, nz):
materialHx = ID[3, i, j, k]
Hx[i, j, k] = updatecoeffsH[materialHx, 0] * Hx[i, j, k] - updatecoeffsH[materialHx, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) + updatecoeffsH[materialHx, 3] * (Ey[i, j, k + 1] - Ey[i, j, k])
Hx[i, j, k] = (updatecoeffsH[materialHx, 0] * Hx[i, j, k] -
updatecoeffsH[materialHx, 2] * (Ez[i, j + 1, k] - Ez[i, j, k]) +
updatecoeffsH[materialHx, 3] * (Ey[i, j, k + 1] - Ey[i, j, k]))
# Hy component
if nx == 1 or nz == 1:
@@ -150,7 +166,9 @@ cpdef void update_magnetic(
for j in range(1, ny):
for k in range(0, nz):
materialHy = ID[4, i, j, k]
Hy[i, j, k] = updatecoeffsH[materialHy, 0] * Hy[i, j, k] - updatecoeffsH[materialHy, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) + updatecoeffsH[materialHy, 1] * (Ez[i + 1, j, k] - Ez[i, j, k])
Hy[i, j, k] = (updatecoeffsH[materialHy, 0] * Hy[i, j, k] -
updatecoeffsH[materialHy, 3] * (Ex[i, j, k + 1] - Ex[i, j, k]) +
updatecoeffsH[materialHy, 1] * (Ez[i + 1, j, k] - Ez[i, j, k]))
# Hz component
if nx == 1 or ny == 1:
@@ -158,7 +176,9 @@ cpdef void update_magnetic(
for j in range(0, ny):
for k in range(1, nz):
materialHz = ID[5, i, j, k]
Hz[i, j, k] = updatecoeffsH[materialHz, 0] * Hz[i, j, k] - updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k] - Ey[i, j, k]) + updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k] - Ex[i, j, k])
Hz[i, j, k] = (updatecoeffsH[materialHz, 0] * Hz[i, j, k] -
updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k] - Ey[i, j, k]) +
updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k] - Ex[i, j, k]))
# 3D
else:
for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
@@ -167,6 +187,12 @@ cpdef void update_magnetic(
materialHx = ID[3, i + 1, j, k]
materialHy = ID[4, i, j + 1, k]
materialHz = ID[5, i, j, k + 1]
Hx[i + 1, j, k] = updatecoeffsH[materialHx, 0] * Hx[i + 1, j, k] - updatecoeffsH[materialHx, 2] * (Ez[i + 1, j + 1, k] - Ez[i + 1, j, k]) + updatecoeffsH[materialHx, 3] * (Ey[i + 1, j, k + 1] - Ey[i + 1, j, k])
Hy[i, j + 1, k] = updatecoeffsH[materialHy, 0] * Hy[i, j + 1, k] - updatecoeffsH[materialHy, 3] * (Ex[i, j + 1, k + 1] - Ex[i, j + 1, k]) + updatecoeffsH[materialHy, 1] * (Ez[i + 1, j + 1, k] - Ez[i, j + 1, k])
Hz[i, j, k + 1] = updatecoeffsH[materialHz, 0] * Hz[i, j, k + 1] - updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k + 1] - Ey[i, j, k + 1]) + updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k + 1] - Ex[i, j, k + 1])
Hx[i + 1, j, k] = (updatecoeffsH[materialHx, 0] * Hx[i + 1, j, k] -
updatecoeffsH[materialHx, 2] * (Ez[i + 1, j + 1, k] - Ez[i + 1, j, k]) +
updatecoeffsH[materialHx, 3] * (Ey[i + 1, j, k + 1] - Ey[i + 1, j, k]))
Hy[i, j + 1, k] = (updatecoeffsH[materialHy, 0] * Hy[i, j + 1, k] -
updatecoeffsH[materialHy, 3] * (Ex[i, j + 1, k + 1] - Ex[i, j + 1, k]) +
updatecoeffsH[materialHy, 1] * (Ez[i + 1, j + 1, k] - Ez[i, j + 1, k]))
Hz[i, j, k + 1] = (updatecoeffsH[materialHz, 0] * Hz[i, j, k + 1] -
updatecoeffsH[materialHz, 1] * (Ey[i + 1, j, k + 1] - Ey[i, j, k + 1]) +
updatecoeffsH[materialHz, 2] * (Ex[i, j + 1, k + 1] - Ex[i, j, k + 1]))

61
gprMax/cython/fractals_generate.pyx
查看文件

@@ -23,17 +23,29 @@ from cython.parallel import prange
from gprMax.config cimport float_or_double_complex
cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, ::1] A, float_or_double_complex[:, ::1] fractalsurface):
"""This function generates a fractal surface for a 2D array.
cpdef void generate_fractal2D(
int nx,
int ny,
int nthreads,
int b,
np.float64_t[:] weighting,
np.float64_t[:] v1,
np.complex128_t[:, ::1] A,
float_or_double_complex[:, ::1] fractalsurface
):
"""Generates a fractal surface for a 2D array.
Args:
nx, ny (int): Fractal surface size in cells
nthreads (int): Number of threads to use
b (int): Constant related to fractal dimension
weighting (memoryview): Access to weighting vector
v1 (memoryview): Access to positional vector at centre of array, scaled by weighting
A (memoryview): Access to array containing random numbers (to be convolved with fractal function)
fractalsurface (memoryview): Access to array containing fractal surface data
nx, ny: int for fractal surface size in cells.
nthreads: int for number of threads to use
b: int for constant related to fractal dimension.
weighting: memoryview for access to weighting vector.
v1: memoryview for access to positional vector at centre of array,
scaled by weighting.
A: memoryview for access to array containing random numbers
(to be convolved with fractal function).
fractalsurface: memoryview for access to array containing fractal
surface data.
"""
cdef Py_ssize_t i, j
@@ -54,17 +66,30 @@ cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[
fractalsurface[i, j] = A[i, j] / B
cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, :, ::1] A, float_or_double_complex[:, :, ::1] fractalvolume):
"""This function generates a fractal volume for a 3D array.
cpdef void generate_fractal3D(
int nx,
int ny,
int nz,
int nthreads,
int b,
np.float64_t[:] weighting,
np.float64_t[:] v1,
np.complex128_t[:, :, ::1] A,
float_or_double_complex[:, :, ::1] fractalvolume
):
"""Generates a fractal volume for a 3D array.
Args:
nx, ny, nz (int): Fractal volume size in cells
nthreads (int): Number of threads to use
b (int): Constant related to fractal dimension
weighting (memoryview): Access to weighting vector
v1 (memoryview): Access to positional vector at centre of array, scaled by weighting
A (memoryview): Access to array containing random numbers (to be convolved with fractal function)
fractalvolume (memoryview): Access to array containing fractal volume data
nx, ny, nz: int for fractal volume size in cells.
nthreads: int for number of threads to use
b: int for constant related to fractal dimension.
weighting: memoryview for access to weighting vector.
v1: memoryview for access to positional vector at centre of array,
scaled by weighting.
A: memoryview for access to array containing random numbers
(to be convolved with fractal function).
fractalsurface: memoryview for access to array containing fractal
volume data.
"""
cdef Py_ssize_t i, j, k

31
gprMax/cython/geometry_outputs.pyx
查看文件

@@ -20,21 +20,30 @@ import numpy as np
cimport numpy as np
cpdef write_lines(float xs, float ys, float zs, int nx, int ny, int nz,
float dx, float dy, float dz, np.uint32_t[:, :, :, :] ID):
"""This function generates arrays with to be written as lines (cell edges)
to a VTK file.
cpdef write_lines(
float xs,
float ys,
float zs,
int nx,
int ny,
int nz,
float dx,
float dy,
float dz,
np.uint32_t[:, :, :, :] ID
):
"""Generates arrays with to be written as lines (cell edges) to a VTK file.
Args:
xs, ys, zs (float): Starting coordinates of geometry view in metres
nx, ny, nz (int): Size of the volume in cells
dx, dy, dz (float): Spatial discretisation of geometry view in metres
ID (nparray): Sampled ID array according to geometry view spatial
discretisation
xs, ys, zs: float for starting coordinates of geometry view in metres.
nx, ny, nz: int for size of the volume in cells.
dx, dy, dz: float for spatial discretisation of geometry view in metres.
ID: memoryview of sampled ID array according to geometry view spatial
discretisation.
Returns:
x, y, z (nparray): 1D arrays with coordinates of the vertex of the lines
lines (nparray): array of material IDs for each line (cell edge) required
x, y, z: 1D nparrays with coordinates of the vertex of the lines.
lines: nparray of material IDs for each line (cell edge) required.
"""
cdef Py_ssize_t i, j, k

712
gprMax/cython/geometry_primitives.pyx
查看文件

@@ -34,15 +34,15 @@ from gprMax.utilities.utilities import round_value
cpdef bint are_clockwise(
float v1x,
float v1y,
float v2x,
float v2y
):
float v1x,
float v1y,
float v2x,
float v2y
):
"""Find if vector 2 is clockwise relative to vector 1.
Args:
v1x, v1y, v2x, v2y (float): Coordinates of vectors.
v1x, v1y, v2x, v2y: floats of coordinates of vectors.
Returns:
(boolean)
@@ -52,15 +52,15 @@ cpdef bint are_clockwise(
cpdef bint is_within_radius(
float vx,
float vy,
float radius
):
float vx,
float vy,
float radius
):
"""Check if the point is within a given radius of the centre of the circle.
Args:
vx, vy (float): Coordinates of vector.
radius (float): Radius.
vx, vy: floats of coordinates of vector.
radius: float for radius.
Returns:
(boolean)
@@ -70,14 +70,14 @@ cpdef bint is_within_radius(
cpdef bint is_inside_sector(
float px,
float py,
float ctrx,
float ctry,
float sectorstartangle,
float sectorangle,
float radius
):
float px,
float py,
float ctrx,
float ctry,
float sectorstartangle,
float sectorangle,
float radius
):
"""For a point to be inside a circular sector, it has to meet the following tests:
It has to be positioned anti-clockwise from the start "arm" of the sector
It has to be positioned clockwise from the end arm of the sector
@@ -86,11 +86,11 @@ cpdef bint is_inside_sector(
i.e. sector defined in an anti-clockwise direction
Args:
px, py (float): Coordinates of point.
ctrx, ctry (float): Coordinates of centre of circle.
sectorstartangle (float): Angle (in radians) of start of sector.
sectorangle (float): Angle (in radians) that sector makes.
radius (float): Radius.
px, py: floats for coordinates of point.
ctrx, ctry: floats for coordinates of centre of circle.
sectorstartangle: float for angle (in radians) of start of sector.
sectorangle: float for angle (in radians) that sector makes.
radius: float for radius.
Returns:
(boolean)
@@ -111,18 +111,18 @@ cpdef bint is_inside_sector(
cpdef bint point_in_polygon(
float px,
float py,
list polycoords
):
float px,
float py,
list polycoords
):
"""Calculates, using a ray casting algorithm, whether a point lies within a polygon.
Args:
px, py (float): Coordinates of point to test.
polycoords (list): x, y tuples of coordinates that define the polygon.
px, py: float for coordinates of point to test.
polycoords: list of x, y tuples of coordinates that define the polygon.
Returns:
inside (boolean)
inside: boolean
"""
cdef int i
@@ -164,20 +164,20 @@ cpdef bint point_in_polygon(
cpdef void build_edge_x(
int i,
int j,
int k,
int numIDx,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numIDx,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set x-orientated edges in the rigid and ID arrays for a Yee voxel.
Args:
i, j, k (int): Cell coordinates of edge.
numIDz (int): Numeric ID of material.
rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
i, j, k: ints for cell coordinates of edge.
numIDx: int for numeric ID of material.
rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
"""
set_rigid_Ex(i, j, k, rigidE)
@@ -185,20 +185,20 @@ cpdef void build_edge_x(
cpdef void build_edge_y(
int i,
int j,
int k,
int numIDy,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numIDy,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set y-orientated edges in the rigid and ID arrays for a Yee voxel.
Args:
i, j, k (int): Cell coordinates of edge.
numIDz (int): Numeric ID of material.
rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
i, j, k: ints for cell coordinates of edge.
numIDy: int for numeric ID of material.
rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
"""
set_rigid_Ey(i, j, k, rigidE)
@@ -206,20 +206,20 @@ cpdef void build_edge_y(
cpdef void build_edge_z(
int i,
int j,
int k,
int numIDz,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numIDz,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set z-orientated edges in the rigid and ID arrays for a Yee voxel.
Args:
i, j, k (int): Cell coordinates of edge.
numIDz (int): Numeric ID of material.
rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
i, j, k: ints for cell coordinates of edge.
numIDz: int for numeric ID of material.
rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
"""
set_rigid_Ez(i, j, k, rigidE)
@@ -227,21 +227,21 @@ cpdef void build_edge_z(
cpdef void build_face_yz(
int i,
int j,
int k,
int numIDy,
int numIDz,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numIDy,
int numIDz,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set the edges of the yz-plane face of a Yell cell in the rigid and ID arrays.
Args:
i, j, k (int): Cell coordinates of the face.
numIDx, numIDy (int): Numeric ID of material.
rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
i, j, k: ints for cell coordinates of the face.
numIDy, numIDz: int for numeric ID of material.
rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
"""
set_rigid_Ey(i, j, k, rigidE)
@@ -255,21 +255,21 @@ cpdef void build_face_yz(
cpdef void build_face_xz(
int i,
int j,
int k,
int numIDx,
int numIDz,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numIDx,
int numIDz,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set the edges of the xz-plane face of a Yell cell in the rigid and ID arrays.
Args:
i, j, k (int): Cell coordinates of the face.
numIDx, numIDy (int): Numeric ID of material.
rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
i, j, k: ints for cell coordinates of the face.
numIDx, numIDz: int for numeric ID of material.
rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
"""
set_rigid_Ex(i, j, k, rigidE)
@@ -283,21 +283,21 @@ cpdef void build_face_xz(
cpdef void build_face_xy(
int i,
int j,
int k,
int numIDx,
int numIDy,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numIDx,
int numIDy,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set the edges of the xy-plane face of a Yell cell in the rigid and ID arrays.
Args:
i, j, k (int): Cell coordinates of the face.
numIDx, numIDy (int): Numeric ID of material.
rigidE, rigidH, ID (memoryviews): Access to rigid and ID arrays.
i, j, k: ints for cell coordinates of the face.
numIDx, numIDy: int for numeric ID of material.
rigidE, rigidH, ID: memoryviews to access rigid and ID arrays.
"""
set_rigid_Ex(i, j, k, rigidE)
@@ -311,26 +311,27 @@ cpdef void build_face_xy(
cpdef void build_voxel(
int i,
int j,
int k,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int i,
int j,
int k,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Set values in the solid, rigid and ID arrays for a Yee voxel.
Args:
i, j, k (int): Cell coordinates of voxel.
numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
averaging (bint): Whether material property averaging will occur for the object.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
i, j, k: ints for cell coordinates of voxel.
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
averaging: bint for whether material property averaging will occur for
the object.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
if averaging:
@@ -375,43 +376,44 @@ cpdef void build_voxel(
cpdef void build_triangle(
float x1,
float y1,
float z1,
float x2,
float y2,
float z2,
float x3,
float y3,
float z3,
str normal,
float thickness,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
float x1,
float y1,
float z1,
float x2,
float y2,
float z2,
float x3,
float y3,
float z3,
str normal,
float thickness,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""
Builds #triangle and #triangular_prism commands which sets values in the
solid, rigid and ID arrays for a Yee voxel.
Builds triangles and triangular prisms which sets values in the solid,
rigid and ID arrays for a Yee voxel.
Args:
x1, y1, z1, x2, y2, z2, x3, y3, z3 (float): Coordinates of the vertices
of the triangular prism.
normal (char): Normal direction to the plane of the triangular prism.
thickness (float): Thickness of the triangular prism.
dx, dy, dz (float): Spatial discretisation.
numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
averaging (bint): Whether material property averaging will occur for the object.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
x1, y1, z1, x2, y2, z2, x3, y3, z3: floats of coordinates of the vertices
of the triangular prism.
normal: string for normal direction to the plane of the triangular prism.
thickness: float for thickness of the triangular prism.
dx, dy, dz: floats for spatial discretisation.
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
averaging: bint for whether material property averaging will occur for
the object.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t i, j, k
@@ -449,7 +451,8 @@ cpdef void build_triangle(
for i in range(i1, i2):
for j in range(j1, j2):
# Calculate the areas of the 3 triangles defined by the 3 vertices of the main triangle and the point under test
# Calculate the areas of the 3 triangles defined by the 3 vertices
# of the main triangle and the point under test.
if normal == 'x':
ir = (i + 0.5) * dy
jr = (j + 0.5) * dz
@@ -470,70 +473,78 @@ cpdef void build_triangle(
if s > 0 and t > 0 and (s + t) < 2 * area * sign:
if thicknesscells == 0:
if normal == 'x':
build_face_yz(level, i, j, numIDy, numIDz, rigidE, rigidH, ID)
build_face_yz(level, i, j, numIDy, numIDz,
rigidE, rigidH, ID)
elif normal == 'y':
build_face_xz(i, level, j, numIDx, numIDz, rigidE, rigidH, ID)
build_face_xz(i, level, j, numIDx, numIDz,
rigidE, rigidH, ID)
elif normal == 'z':
build_face_xy(i, j, level, numIDx, numIDy, rigidE, rigidH, ID)
build_face_xy(i, j, level, numIDx, numIDy,
rigidE, rigidH, ID)
else:
for k in range(level, level + thicknesscells):
if normal == 'x':
build_voxel(k, i, j, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(k, i, j, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
elif normal == 'y':
build_voxel(i, k, j, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, k, j, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
elif normal == 'z':
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
cpdef void build_cylindrical_sector(
float ctr1,
float ctr2,
int level,
float sectorstartangle,
float sectorangle,
float radius,
str normal,
float thickness,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
float ctr1,
float ctr2,
int level,
float sectorstartangle,
float sectorangle,
float radius,
str normal,
float thickness,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""
Builds #cylindrical_sector commands which sets values in the solid, rigid
and ID arrays for a Yee voxel. It defines a sector of cylinder given by the
direction of the axis of the coordinates of the cylinder face centre, depth
coordinates, sector start point, sector angle, and sector radius. N.B
Assumes sector start is always clockwise from sector end, i.e. sector
defined in an anti-clockwise direction.
Builds cylindrical sectors which sets values in the solid, rigid and ID
arrays for a Yee voxel. It defines a sector of cylinder given by the
direction of the axis of the coordinates of the cylinder face centre,
depth coordinates, sector start point, sector angle, and sector radius.
N.B Assumes sector start is always clockwise from sector end,
i.e. sector defined in an anti-clockwise direction.
Args:
ctr1, ctr2 (float): Coordinates of centre of circle.
level (int): Third dimensional coordinate.
sectorstartangle (float): Angle (in radians) of start of sector.
sectorangle (float): Angle (in radians) that sector makes.
radius (float): Radius of the cylindrical sector.
normal (char): Normal direction to the plane of the cylindrical sector.
thickness (float): Thickness of the cylindrical sector.
dx, dy, dz (float): Spatial discretisation.
numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
averaging (bint): Whether material property averaging will occur for the object.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
ctr1, ctr2: floats for coordinates of centre of circle.
level: int for the third dimensional coordinate.
sectorstartangle: float for angle (in radians) of start of sector.
sectorangle: float for angle (in radians) that sector makes.
radius: float for radius of the cylindrical sector.
normal: string for normal direction to the plane of the triangular prism.
thickness: float for thickness of the triangular prism.
dx, dy, dz: floats for spatial discretisation.
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
averaging: bint for whether material property averaging will occur for
the object.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t x, y, z
cdef int x1, x2, y1, y2, z1, z2, thicknesscells
if normal == 'x':
# Angles are defined from zero degrees on the positive y-axis going towards positive z-axis
# Angles are defined from zero degrees on the positive y-axis going
# towards positive z-axis.
y1 = round_value((ctr1 - radius)/dy)
y2 = round_value((ctr1 + radius)/dy)
z1 = round_value((ctr2 - radius)/dz)
@@ -552,15 +563,19 @@ cpdef void build_cylindrical_sector(
for y in range(y1, y2):
for z in range(z1, z2):
if is_inside_sector(y * dy + 0.5 * dy, z * dz + 0.5 * dz, ctr1, ctr2, sectorstartangle, sectorangle, radius):
if is_inside_sector(y * dy + 0.5 * dy, z * dz + 0.5 * dz, ctr1,
ctr2, sectorstartangle, sectorangle, radius):
if thicknesscells == 0:
build_face_yz(level, y, z, numIDy, numIDz, rigidE, rigidH, ID)
build_face_yz(level, y, z, numIDy, numIDz,
rigidE, rigidH, ID)
else:
for x in range(level, level + thicknesscells):
build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(x, y, z, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
elif normal == 'y':
# Angles are defined from zero degrees on the positive x-axis going towards positive z-axis
# Angles are defined from zero degrees on the positive x-axis going
# towards positive z-axis.
x1 = round_value((ctr1 - radius)/dx)
x2 = round_value((ctr1 + radius)/dx)
z1 = round_value((ctr2 - radius)/dz)
@@ -579,15 +594,19 @@ cpdef void build_cylindrical_sector(
for x in range(x1, x2):
for z in range(z1, z2):
if is_inside_sector(x * dx + 0.5 * dx, z * dz + 0.5 * dz, ctr1, ctr2, sectorstartangle, sectorangle, radius):
if is_inside_sector(x * dx + 0.5 * dx, z * dz + 0.5 * dz, ctr1,
ctr2, sectorstartangle, sectorangle, radius):
if thicknesscells == 0:
build_face_xz(x, level, z, numIDx, numIDz, rigidE, rigidH, ID)
build_face_xz(x, level, z, numIDx, numIDz,
rigidE, rigidH, ID)
else:
for y in range(level, level + thicknesscells):
build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(x, y, z, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
elif normal == 'z':
# Angles are defined from zero degrees on the positive x-axis going towards positive y-axis
# Angles are defined from zero degrees on the positive x-axis going
# towards positive y-axis.
x1 = round_value((ctr1 - radius)/dx)
x2 = round_value((ctr1 + radius)/dx)
y1 = round_value((ctr2 - radius)/dy)
@@ -606,38 +625,42 @@ cpdef void build_cylindrical_sector(
for x in range(x1, x2):
for y in range(y1, y2):
if is_inside_sector(x * dx + 0.5 * dx, y * dy + 0.5 * dy, ctr1, ctr2, sectorstartangle, sectorangle, radius):
if is_inside_sector(x * dx + 0.5 * dx, y * dy + 0.5 * dy, ctr1,
ctr2, sectorstartangle, sectorangle, radius):
if thicknesscells == 0:
build_face_xy(x, y, level, numIDx, numIDy, rigidE, rigidH, ID)
build_face_xy(x, y, level, numIDx, numIDy,
rigidE, rigidH, ID)
else:
for z in range(level, level + thicknesscells):
build_voxel(x, y, z, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(x, y, z, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
cpdef void build_box(
int xs,
int xf,
int ys,
int yf,
int zs,
int zf,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds #box commands which sets values in the solid, rigid and ID arrays.
int xs,
int xf,
int ys,
int yf,
int zs,
int zf,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds boxes which sets values in the solid, rigid and ID arrays.
Args:
xs, xf, ys, yf, zs, zf (int): Cell coordinates of entire box.
numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
averaging (bint): Whether material property averaging will occur for the object.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
xs, xf, ys, yf, zs, zf: ints for cell coordinates of entire box.
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
averaging: bint for whether material property averaging will occur for
the object.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t i, j, k
@@ -689,35 +712,38 @@ cpdef void build_box(
cpdef void build_cylinder(
float x1,
float y1,
float z1,
float x2,
float y2,
float z2,
float r,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds #cylinder commands which sets values in the solid, rigid and ID arrays for a Yee voxel.
float x1,
float y1,
float z1,
float x2,
float y2,
float z2,
float r,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds cylinders which sets values in the solid, rigid and ID arrays for
a Yee voxel.
Args:
x1, y1, z1, x2, y2, z2 (float): Coordinates of the centres of cylinder faces.
r (float): Radius of the cylinder.
dx, dy, dz (float): Spatial discretisation.
numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
averaging (bint): Whether material property averaging will occur for the object.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
x1, y1, z1, x2, y2, z2: floats for coordinates of the centres of cylinder
faces.
r: float for radius of the cylinder.
dx, dy, dz: floats for spatial discretisation.
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
averaging: bint for whether material property averaging will occur for
the object.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t i, j, k
@@ -729,15 +755,18 @@ cpdef void build_cylinder(
# Check if cylinder is aligned with an axis
x_align = y_align = z_align = 0
# x-aligned
if round_value(y1 / dy) == round_value(y2 / dy) and round_value(z1 / dz) == round_value(z2 / dz):
if (round_value(y1 / dy) == round_value(y2 / dy) and
round_value(z1 / dz) == round_value(z2 / dz)):
x_align = 1
# y-aligned
elif round_value(x1 / dx) == round_value(x2 / dx) and round_value(z1 / dz) == round_value(z2 / dz):
elif (round_value(x1 / dx) == round_value(x2 / dx) and
round_value(z1 / dz) == round_value(z2 / dz)):
y_align = 1
# z-aligned
elif round_value(x1 / dx) == round_value(x2 / dx) and round_value(y1 / dy) == round_value(y2 / dy):
elif (round_value(x1 / dx) == round_value(x2 / dx) and
round_value(y1 / dy) == round_value(y2 / dy):)
z_align = 1
# Calculate a bounding box for the cylinder
@@ -804,21 +833,24 @@ cpdef void build_cylinder(
for k in range(zs, zf):
if np.sqrt((j * dy + 0.5 * dy - y1)**2 + (k * dz + 0.5 * dz - z1)**2) <= r:
for i in range(xs, xf):
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
# y-aligned cylinder
elif y_align:
for i in range(xs, xf):
for k in range(zs, zf):
if np.sqrt((i * dx + 0.5 * dx - x1)**2 + (k * dz + 0.5 * dz - z1)**2) <= r:
for j in range(ys, yf):
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
# z-aligned cylinder
elif z_align:
for i in range(xs, xf):
for j in range(ys, yf):
if np.sqrt((i * dx + 0.5 * dx - x1)**2 + (j * dy + 0.5 * dy - y1)**2) <= r:
for k in range(zs, zf):
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
# Not aligned with any axis
else:
@@ -836,9 +868,13 @@ cpdef void build_cylinder(
# Build flag - default false, set to True if point is in cylinder
build = 0
# Vector from centre of first cylinder face to test point
f1pt = np.array([i * dx + 0.5 * dx - x1, j * dy + 0.5 * dy - y1, k * dz + 0.5 * dz - z1], dtype=np.float32)
f1pt = np.array([i * dx + 0.5 * dx - x1,
j * dy + 0.5 * dy - y1,
k * dz + 0.5 * dz - z1], dtype=np.float32)
# Vector from centre of second cylinder face to test point
f2pt = np.array([i * dx + 0.5 * dx - x2, j * dy + 0.5 * dy - y2, k * dz + 0.5 * dz - z2], dtype=np.float32)
f2pt = np.array([i * dx + 0.5 * dx - x2,
j * dy + 0.5 * dy - y2,
k * dz + 0.5 * dz - z2], dtype=np.float32)
# Magnitudes
f1ptmag = np.sqrt((f1pt*f1pt).sum(axis=0))
f2ptmag = np.sqrt((f2pt*f2pt).sum(axis=0))
@@ -862,40 +898,44 @@ cpdef void build_cylinder(
theta2 = 0
distance1 = f1ptmag * np.sin(theta1)
distance2 = f2ptmag * np.sin(theta2)
if (distance1 <= r or distance2 <= r) and theta1 <= np.pi/2 and theta2 <= np.pi/2:
if ((distance1 <= r or distance2 <= r) and
theta1 <= np.pi/2 and theta2 <= np.pi/2):
build = 1
if build:
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
cpdef void build_sphere(
int xc,
int yc,
int zc,
float r,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds #sphere commands which sets values in the solid, rigid and ID arrays for a Yee voxel.
int xc,
int yc,
int zc,
float r,
float dx,
float dy,
float dz,
int numID,
int numIDx,
int numIDy,
int numIDz,
bint averaging,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds spheres which sets values in the solid, rigid and ID arrays for
a Yee voxel.
Args:
xc, yc, zc (int): Cell coordinates of the centre of the sphere.
r (float): Radius of the sphere.
dx, dy, dz (float): Spatial discretisation.
numID, numIDx, numIDy, numIDz (int): Numeric ID of material.
averaging (bint): Whether material property averaging will occur for the object.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
xc, yc, zc: ints for cell coordinates of the centre of the sphere.
r: float for radius of the sphere.
dx, dy, dz: floats for spatial discretisation.
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
averaging: bint for whether material property averaging will occur for
the object.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t i, j, k
@@ -926,30 +966,36 @@ cpdef void build_sphere(
for i in range(xs, xf):
for j in range(ys, yf):
for k in range(zs, zf):
if np.sqrt((i + 0.5 - xc)**2 * dx**2 + (j + 0.5 - yc)**2 * dy**2 + (k + 0.5 - zc)**2 * dz**2) <= r:
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
if (np.sqrt((i + 0.5 - xc)**2 * dx**2 +
(j + 0.5 - yc)**2 * dy**2 +
(k + 0.5 - zc)**2 * dz**2) <= r):
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
cpdef void build_voxels_from_array(
int xs,
int ys,
int zs,
int numexistmaterials,
bint averaging,
np.int16_t[:, :, ::1] data,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int xs,
int ys,
int zs,
int numexistmaterials,
bint averaging,
np.int16_t[:, :, ::1] data,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds Yee voxels by reading integers from an array.
Args:
xs, ys, zs (int): Cell coordinates of position of start of array in domain.
numexistmaterials (int): Number of existing materials in model prior to building voxels.
averaging (bint): Whether material property averaging will occur for the object.
data (memoryview): Access to array containing numeric IDs of voxels to create.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
xs, ys, zs: ints for cell coordinates of position of start of array in
domain.
numexistmaterials: int for number of existing materials in model prior
to building voxels.
averaging: bint for whether material property averaging will occur for
the object.
data: memoryview to access array containing numeric IDs of voxels to create.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t i, j, k
@@ -987,28 +1033,29 @@ cpdef void build_voxels_from_array(
cpdef void build_voxels_from_array_mask(
int xs,
int ys,
int zs,
int waternumID,
int grassnumID,
bint averaging,
np.int8_t[:, :, ::1] mask,
np.int16_t[:, :, ::1] data,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
int xs,
int ys,
int zs,
int waternumID,
int grassnumID,
bint averaging,
np.int8_t[:, :, ::1] mask,
np.int16_t[:, :, ::1] data,
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID
):
"""Builds Yee voxels by reading integers from an array.
Args:
xs, ys, zs (int): Cell coordinates of position of start of array in domain.
waternumID, grassnumID (int): Numeric ID of water and grass materials.
averaging (bint): Whether material property averaging will occur for the object.
data (memoryview): Access to array containing numeric IDs of voxels to create.
mask (memoryview): Access to array containing a mask of voxels to create.
solid, rigidE, rigidH, ID (memoryviews): Access to solid, rigid and ID arrays.
xs, ys, zs: ints for cell coordinates of position of start of array in domain.
waternumID, grassnumID: ints for numeric ID of water and grass materials.
averaging: bint for whether material property averaging will occur for
the object.
data: memoryview to access array containing numeric IDs of voxels to create.
mask: memoryview to access to array containing a mask of voxels to create.
solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
"""
cdef Py_ssize_t i, j, k
@@ -1024,10 +1071,13 @@ cpdef void build_voxels_from_array_mask(
for k in range(zs, zf):
if mask[i - xs, j - ys, k - zs] == 1:
numID = numIDx = numIDy = numIDz = data[i - xs, j - ys, k - zs]
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
elif mask[i - xs, j - ys, k - zs] == 2:
numID = numIDx = numIDy = numIDz = waternumID
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)
elif mask[i - xs, j - ys, k - zs] == 3:
numID = numIDx = numIDy = numIDz = grassnumID
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz, averaging, solid, rigidE, rigidH, ID)
build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
averaging, solid, rigidE, rigidH, ID)

44
gprMax/cython/snapshots.pyx
查看文件

@@ -20,31 +20,29 @@ from gprMax.config cimport float_or_double
cpdef void calculate_snapshot_fields(
int nx,
int ny,
int nz,
float_or_double[:, :, ::1] sliceEx,
float_or_double[:, :, ::1] sliceEy,
float_or_double[:, :, ::1] sliceEz,
float_or_double[:, :, ::1] sliceHx,
float_or_double[:, :, ::1] sliceHy,
float_or_double[:, :, ::1] sliceHz,
float_or_double[:, :, ::1] snapEx,
float_or_double[:, :, ::1] snapEy,
float_or_double[:, :, ::1] snapEz,
float_or_double[:, :, ::1] snapHx,
float_or_double[:, :, ::1] snapHy,
float_or_double[:, :, ::1] snapHz
):
"""This function calculates electric and magnetic values at points from
averaging values in cells.
int nx,
int ny,
int nz,
float_or_double[:, :, ::1] sliceEx,
float_or_double[:, :, ::1] sliceEy,
float_or_double[:, :, ::1] sliceEz,
float_or_double[:, :, ::1] sliceHx,
float_or_double[:, :, ::1] sliceHy,
float_or_double[:, :, ::1] sliceHz,
float_or_double[:, :, ::1] snapEx,
float_or_double[:, :, ::1] snapEy,
float_or_double[:, :, ::1] snapEz,
float_or_double[:, :, ::1] snapHx,
float_or_double[:, :, ::1] snapHy,
float_or_double[:, :, ::1] snapHz
):
"""Calculates electric and magnetic values at points from averaging values
in cells.
Args:
nx, ny, nz (int): Size of snapshot array
sliceEx, sliceEy, sliceEz,
sliceHx, sliceHy, sliceHz (memoryview): Access to slices of field arrays
snapEx, snapEy, snapEz,
snapHx, snapHy, snapHz (memoryview): Access to snapshot arrays
nx, ny, nz: ints for size of snapshot array.
slice: memoryviews to access slices of field arrays.
snap: memoryviews to access snapshot arrays.
"""
cdef Py_ssize_t i, j, k

103
gprMax/cython/yee_cell_build.pyx
查看文件

@@ -28,18 +28,30 @@ from gprMax.cython.yee_cell_setget_rigid cimport get_rigid_Hz
from gprMax.materials import Material
cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2, int numID3, int numID4, int componentID, G):
"""This function creates a new material by averaging the dielectric properties of the surrounding cells.
cpdef void create_electric_average(
int i,
int j,
int k,
int numID1,
int numID2,
int numID3,
int numID4,
int componentID,
G
):
"""Creates a new material by averaging the dielectric properties of the
surrounding cells.
Args:
i, j, k (int): Cell coordinates.
numID1, numID2, numID3, numID4 (int): Numeric IDs for materials in surrounding cells.
componentID (int): Numeric ID for electric field component.
G (class): Grid class instance - holds essential parameters describing the model.
i, j, k: ints for cell coordinates.
numID: ints for numeric IDs for materials in surrounding cells.
componentID: int for numeric ID for electric field component.
G: FDTDGrid class describing a grid in a model.
"""
# Make an ID composed of the names of the four materials that will be averaged
requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' + G.materials[numID3].ID + '+' + G.materials[numID4].ID
requiredID = (G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' +
G.materials[numID3].ID + '+' + G.materials[numID4].ID)
# Check if this material already exists
tmp = requiredID.split('+')
@@ -57,10 +69,14 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
m = Material(newNumID, requiredID)
m.type = 'dielectric-smoothed'
# Create averaged constituents for material
m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, G.materials[numID3].er, G.materials[numID4].er), axis=0)
m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, G.materials[numID3].se, G.materials[numID4].se), axis=0)
m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er,
G.materials[numID3].er, G.materials[numID4].er), axis=0)
m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se,
G.materials[numID3].se, G.materials[numID4].se), axis=0)
m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr,
G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm,
G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
# Append the new material object to the materials list
G.materials.append(m)
@@ -68,14 +84,23 @@ cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2,
G.ID[componentID, i, j, k] = newNumID
cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2, int componentID, G):
"""This function creates a new material by averaging the dielectric properties of the surrounding cells.
cpdef void create_magnetic_average(
int i,
int j,
int k,
int numID1,
int numID2,
int componentID,
G
):
"""Creates a new material by averaging the dielectric properties of the
surrounding cells.
Args:
i, j, k (int): Cell coordinates.
numID1, numID2 (int): Numeric IDs for materials in surrounding cells.
componentID (int): Numeric ID for electric field component.
G (class): Grid class instance - holds essential parameters describing the model.
i, j, k: ints for cell coordinates.
numID: ints for numeric IDs for materials in surrounding cells.
componentID: int for numeric ID for magnetic field component.
G: FDTDGrid class describing a grid in a model.
"""
# Make an ID composed of the names of the two materials that will be averaged
@@ -107,12 +132,17 @@ cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2,
G.ID[componentID, i, j, k] = newNumID
cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.uint32_t[:, :, :, ::1] ID, G):
"""This function builds the electric field components in the ID array.
cpdef void build_electric_components(
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidE,
np.uint32_t[:, :, :, ::1] ID,
G
):
"""Builds the electric field components in the ID array.
Args:
solid, rigid, ID (memoryviews): Access to solid, rigid and ID arrays
G (class): Grid class instance - holds essential parameters describing the model.
solid, rigid, ID: memoryviews to access solid, rigid and ID arrays.
G: FDTDGrid class describing a grid in a model.
"""
cdef Py_ssize_t i, j, k
@@ -138,7 +168,8 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
ID[componentID, i, j, k] = numID1
else:
# Averaging is required
create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
create_electric_average(i, j, k, numID1, numID2,
numID3, numID4, componentID, G)
# Ey component
componentID = G.IDlookup['Ey']
@@ -160,7 +191,8 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
ID[componentID, i, j, k] = numID1
else:
# Averaging is required
create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
create_electric_average(i, j, k, numID1, numID2,
numID3, numID4, componentID, G)
# Ez component
componentID = G.IDlookup['Ez']
@@ -182,15 +214,21 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
ID[componentID, i, j, k] = numID1
else:
# Averaging is required
create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
create_electric_average(i, j, k, numID1, numID2,
numID3, numID4, componentID, G)
cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID, G):
"""This function builds the magnetic field components in the ID array.
cpdef void build_magnetic_components(
np.uint32_t[:, :, ::1] solid,
np.int8_t[:, :, :, ::1] rigidH,
np.uint32_t[:, :, :, ::1] ID,
G
):
"""Builds the magnetic field components in the ID array.
Args:
solid, rigid, ID (memoryviews): Access to solid, rigid and ID arrays
G (class): Grid class instance - holds essential parameters describing the model.
solid, rigid, ID: memoryviews to access solid, rigid and ID arrays.
G: FDTDGrid class describing a grid in a model.
"""
cdef Py_ssize_t i, j, k
@@ -214,7 +252,8 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
ID[componentID, i, j, k] = numID1
else:
# Averaging is required
create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
create_magnetic_average(i, j, k, numID1, numID2,
componentID, G)
# Hy component
componentID = G.IDlookup['Hy']
@@ -234,7 +273,8 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
ID[4, i, j, k] = numID1
else:
# Averaging is required
create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
create_magnetic_average(i, j, k, numID1, numID2,
componentID, G)
# Hz component
componentID = G.IDlookup['Hz']
@@ -254,4 +294,5 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
ID[5, i, j, k] = numID1
else:
# Averaging is required
create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
create_magnetic_average(i, j, k, numID1, numID2,
componentID, G)

127
gprMax/cython/yee_cell_setget_rigid.pyx
查看文件

@@ -23,7 +23,13 @@ cimport numpy as np
# Get and set functions for the rigid electric component array. The rigid array
# is 4D with the 1st dimension holding the 12 electric edge components of a
# cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef bint get_rigid_Ex(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
cdef bint result
result = False
if rigidE[0, i, j, k]:
@@ -39,7 +45,13 @@ cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
result = True
return result
cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef bint get_rigid_Ey(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
cdef bint result
result = False
if rigidE[4, i, j, k]:
@@ -55,7 +67,13 @@ cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
result = True
return result
cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef bint get_rigid_Ez(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
cdef bint result
result = False
if rigidE[8, i, j, k]:
@@ -71,7 +89,13 @@ cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
result = True
return result
cdef void set_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef void set_rigid_Ex(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
rigidE[0, i, j, k] = True
if j != 0:
rigidE[1, i, j - 1, k] = True
@@ -80,7 +104,13 @@ cdef void set_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
if j != 0 and k != 0:
rigidE[2, i, j - 1, k - 1] = True
cdef void set_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef void set_rigid_Ey(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
rigidE[4, i, j, k] = True
if i != 0:
rigidE[7, i - 1, j, k] = True
@@ -89,7 +119,13 @@ cdef void set_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
if i != 0 and k != 0:
rigidE[6, i - 1, j, k - 1] = True
cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef void set_rigid_Ez(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
rigidE[8, i, j, k] = True
if i != 0:
rigidE[9, i - 1, j, k] = True
@@ -98,16 +134,33 @@ cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
if i != 0 and j != 0:
rigidE[10, i - 1, j - 1, k] = True
cdef void set_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef void set_rigid_E(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
rigidE[:, i, j, k] = True
cdef void unset_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE):
cdef void unset_rigid_E(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidE
):
rigidE[:, i, j, k] = False
# Get and set functions for the rigid magnetic component array. The rigid array
# is 4D with the 1st dimension holding the 6 magnetic edge components - Hx1,
# Hx2, Hy1, Hy2, Hz1, Hz2
cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef bint get_rigid_Hx(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
cdef bint result
result = False
if rigidH[0, i, j, k]:
@@ -117,7 +170,13 @@ cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
result = True
return result
cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef bint get_rigid_Hy(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
cdef bint result
result = False
if rigidH[2, i, j, k]:
@@ -127,7 +186,13 @@ cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
result = True
return result
cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef bint get_rigid_Hz(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
cdef bint result
result = False
if rigidH[4, i, j, k]:
@@ -137,23 +202,53 @@ cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
result = True
return result
cdef void set_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef void set_rigid_Hx(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
rigidH[0, i, j, k] = True
if i != 0:
rigidH[1, i - 1, j, k] = True
cdef void set_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef void set_rigid_Hy(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
rigidH[2, i, j, k] = True
if j != 0:
rigidH[3, i, j - 1, k] = True
cdef void set_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef void set_rigid_Hz(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
rigidH[4, i, j, k] = True
if k != 0:
rigidH[5, i, j, k - 1] = True
cdef void set_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef void set_rigid_H(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
rigidH[:, i, j, k] = True
cdef void unset_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH):
cdef void unset_rigid_H(
int i,
int j,
int k,
np.int8_t[:, :, :, ::1] rigidH
):
rigidH[:, i, j, k] = False

3
gprMax/subgrids/base.py
查看文件

@@ -17,7 +17,6 @@
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import logging
import numpy as np
from ..grid import FDTDGrid
@@ -57,7 +56,7 @@ class SubGridBase(FDTDGrid):
# Number of sub cells to extend the sub grid beyond the IS boundary
d_to_pml = self.s_is_os_sep + self.pml_separation
# index of the is
# Index of the is
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']

44
gprMax/subgrids/multi.py
查看文件

@@ -34,14 +34,14 @@ class ReferenceRx(Rx):
logger.debug('ReferenceRx has no offset member.')
def get_field(self, str_id, field):
"""Return the field value at the equivalent coarse yee cell.
"""Field value at the equivalent coarse yee cell.
Args:
str_id (str): 'Ex' etc...
field (nparray): Numpy array of grid.Ez
str_id: string of 'Ex' etc...
field: nparray of grid.Ez
Returns:
e (float): Field value.
e: float of field value.
"""
d = {
@@ -58,13 +58,13 @@ class ReferenceRx(Rx):
return e
def get_Ex_from_field(self, Ex):
"""Return the Ex field value from the equivalent coarse yee cell.
"""Ex field value from the equivalent coarse yee cell.
Args:
Ex (nparray): Numpy array of Ex field.
Ex: nparray of Ex field.
Returns:
e (float): Ex field value.
e: float of Ex field value.
"""
# offset = ratio // 2
@@ -72,62 +72,62 @@ class ReferenceRx(Rx):
return e
def get_Ey_from_field(self, Ey):
"""Return the Ey field value from the equivalent coarse yee cell.
"""Ey field value from the equivalent coarse yee cell.
Args:
Ey (nparray): Numpy array of Ey field.
Ey: nparray of Ey field.
Returns:
e (float): Ey field value.
e: float of Ey field value.
"""
e = Ey[self.xcoord, self.ycoord + self.offset, self.zcoord]
return e
def get_Ez_from_field(self, Ez):
"""Return the Ez field value from the equivalent coarse yee cell.
"""Ez field value from the equivalent coarse yee cell.
Args:
Ez (nparray): Numpy array of Ez field.
Ez: nparray of Ez field.
Returns:
e (float): Ez field value.
e: float of Ez field value.
"""
e = Ez[self.xcoord, self.ycoord, self.zcoord + self.offset]
return e
def get_Hx_from_field(self, Hx):
"""Return the Hx field value from the equivalent coarse yee cell.
"""Hx field value from the equivalent coarse yee cell.
Args:
Hx (nparray): Numpy array of Hx field.
Hx: nparray of Hx field.
Returns:
e (float): Hx field value.
e: float of Hx field value.
"""
e = Hx[self.xcoord, self.ycoord + self.offset, self.zcoord + self.offset]
return e
def get_Hy_from_field(self, Hy):
"""Return the Hy field value from the equivalent coarse yee cell.
"""Hy field value from the equivalent coarse yee cell.
Args:
Hy (nparray): Numpy array of Hy field.
Hy: nparray of Hy field.
Returns:
e (float): Hy field value.
e: float of Hy field value.
"""
e = Hy[self.xcoord + self.offset, self.ycoord, self.zcoord + self.offset]
return e
def get_Hz_from_field(self, Hz):
"""Return the Hz field value from the equivalent coarse yee cell.
"""Hz field value from the equivalent coarse yee cell.
Args:
Hz (nparray): Numpy array of Hz field.
Hz: nparray of Hz field.
Returns:
e (float): Hz field value.
e: float of Hz field value.
"""
e = Hz[self.xcoord + self.offset, self.ycoord + self.offset, self.zcoord]
return e

5
gprMax/subgrids/precursor_nodes.py
查看文件

@@ -180,9 +180,8 @@ class PrecursorNodesBase:
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
"""Sets 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'))

6
gprMax/subgrids/subgrid_hsg.py
查看文件

@@ -36,7 +36,7 @@ class SubGridHSG(SubGridBase):
self.gridtype = SubGridHSG.gridtype
def update_magnetic_is(self, precursors):
"""Update the subgrid nodes at the IS with the currents derived
"""Updates the subgrid nodes at the IS with the currents derived
from the main grid.
Args:
@@ -59,7 +59,7 @@ class SubGridHSG(SubGridBase):
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, config.get_model_config().ompthreads)
def update_electric_is(self, precursors):
"""Update the subgrid nodes at the IS with the currents derived
"""Updates the subgrid nodes at the IS with the currents derived
from the main grid.
Args:
@@ -142,7 +142,7 @@ class SubGridHSG(SubGridBase):
cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
def print_info(self):
"""Print information about the subgrid."""
"""Prints information about the subgrid."""
logger.info('')
logger.info(f'[{self.name}] Type: {self.gridtype}')

30
gprMax/subgrids/updates.py
查看文件

@@ -47,37 +47,37 @@ def create_updates(G):
class SubgridUpdates(CPUUpdates):
"""Provides update functions for the Sub gridding simulation."""
"""Update functions for the Sub gridding simulation."""
def __init__(self, G, updaters):
super().__init__(G)
self.updaters = updaters
def hsg_1(self):
"""Method to update the subgrids over the first phase."""
# updaters update each subgrid
"""Updates the subgrids over the first phase."""
for sg_updater in self.updaters:
sg_updater.hsg_1()
def hsg_2(self):
"""Method to update the subgrids over the second phase."""
"""Updates the subgrids over the second phase."""
for sg_updater in self.updaters:
sg_updater.hsg_2()
class SubgridUpdater(CPUUpdates):
"""Class to handle updating the electric and magnetic fields of an HSG
subgrid. The IS, OS, subgrid region and the electric/magnetic sources
are updated using the precursor regions.
"""Handles updating the electric and magnetic fields of an HSG subgrid.
The IS, OS, subgrid region and the electric/magnetic sources are updated
using the precursor regions.
"""
def __init__(self, subgrid, precursors, G):
"""
Args:
subgrid (SubGrid3d): Subgrid to be updated.
precursors (PrecursorNodes): Precursor nodes associated with
the subgrid - contain interpolated fields.
G (FDTDGrid): Parameters describing a grid in a model.
subgrid: SubGrid3d instance to be updated.
precursors (PrecursorNodes): PrecursorNodes instance nodes associated
with the subgrid - contain interpolated
fields.
G: FDTDGrid class describing a grid in a model.
"""
super().__init__(subgrid)
self.precursors = precursors
@@ -85,8 +85,8 @@ class SubgridUpdater(CPUUpdates):
self.source_iteration = 0
def hsg_1(self):
"""This is the first half of the subgrid update. Takes the time step
up to the main grid magnetic update.
"""First half of the subgrid update. Takes the time step up to the main
grid magnetic update.
"""
G = self.G
@@ -126,8 +126,8 @@ class SubgridUpdater(CPUUpdates):
sub_grid.update_electric_os(G)
def hsg_2(self):
"""This is the first half of the subgrid update. Takes the time step
up to the main grid electric update.
"""Second half of the subgrid update. Takes the time step up to the main
grid electric update.
"""
G = self.G

58
gprMax/subgrids/user_objects.py
查看文件

@@ -30,18 +30,17 @@ logger = logging.getLogger(__name__)
class SubGridBase(UserObjectMulti):
"""Class to allow UserObjectMulti and UserObjectGeometry to be nested
in SubGrid type user objects.
"""Allows UserObjectMulti and UserObjectGeometry to be nested in SubGrid
type user objects.
"""
def __init__(self, **kwargs):
"""Constructor."""
super().__init__(**kwargs)
self.children_multiple = []
self.children_geometry = []
def add(self, node):
"""Function to add other user objects. Geometry and multi only."""
"""Adds other user objects. Geometry and multi only."""
if isinstance(node, UserObjectMulti):
self.children_multiple.append(node)
elif isinstance(node, UserObjectGeometry):
@@ -51,14 +50,13 @@ class SubGridBase(UserObjectMulti):
raise ValueError
def set_discretisation(self, sg, grid):
"""Set the spatial discretisation."""
sg.dx = grid.dx / sg.ratio
sg.dy = grid.dy / sg.ratio
sg.dz = grid.dz / sg.ratio
sg.dl = np.array([sg.dx, sg.dy, sg.dz])
def set_main_grid_indices(self, sg, grid, uip, p1, p2):
"""Set subgrid indices related to main grid placement."""
"""Sets subgrid indices related to main grid placement."""
# location of the IS
sg.i0, sg.j0, sg.k0 = p1
sg.i1, sg.j1, sg.k1 = p2
@@ -82,7 +80,7 @@ class SubGridBase(UserObjectMulti):
sg.nz = 2 * sg.n_boundary_cells_z + sg.nwz
def set_iterations(self, sg, main):
"""Set number of iterations that will take place in the subgrid."""
"""Sets number of iterations that will take place in the subgrid."""
sg.iterations = main.iterations * sg.ratio
def setup(self, sg, grid, uip):
@@ -97,7 +95,7 @@ class SubGridBase(UserObjectMulti):
# Set the temporal discretisation
sg.calculate_dt()
# set the indices related to the subgrids main grid placement
# Set the indices related to the subgrids main grid placement
self.set_main_grid_indices(sg, grid, uip, p1, p2)
"""
@@ -121,8 +119,8 @@ class SubGridBase(UserObjectMulti):
sg.timewindow = grid.timewindow
# Copy a subgrid reference to self so that children.create(grid, uip) can access
# the correct grid
# Copy a subgrid reference to self so that children.create(grid, uip)
# can access the correct grid.
self.subgrid = sg
# Copy over built in materials
@@ -141,25 +139,25 @@ class SubGridBase(UserObjectMulti):
class SubGridHSG(SubGridBase):
"""Huygens Surface subgridding (HSG) user object.
:param p1: Position of the lower left corner of the Inner Surface (x, y, z) in the main grid.
:type p1: list, non-optional
:param p2: Position of the upper right corner of the Inner Surface (x, y, z) in the main grid.
:type p2: list, non-optional
:param ratio: Ratio of the main grid spatial step to the sub-grid spatial step. Must be an odd integer.
:type ratio: int, non-optional
:param id: Identifier for the sub-grid.
:type id: str, non-optional
:param is_os_sep: Number of main grid cells between the Inner Surface and the Outer Surface. Defaults to 3.
:type is_os_sep: str, optional
:param pml_separation: Number of sub-grid cells between the Outer Surface and the PML. Defaults to ratio // 2 + 2
:type pml_separation: int, optional
:param subgrid_pml_thickness: Thickness of the PML on each of the 6 sides of the sub-grid. Defaults to 6.
:type subgrid_pml_thickness: int, optional
:param interpolation: Degree of the interpolation scheme used for spatial interpolation of the fields at the Inner Surface. Defaults to Linear
:type interpolation: str, optional
:param filter: Turn on the 3-pole filter. Increases numerical stability. Defaults to True
:type filter: bool, optional
Attributes:
p1: list of the position of the lower left corner of the Inner Surface
(x, y, z) in the main grid.
p2: list of the position of the upper right corner of the Inner Surface
(x, y, z) in the main grid.
ratio: int of the ratio of the main grid spatial step to the sub-grid
spatial step. Must be an odd integer.
id: string identifier for the sub-grid.
is_os_sep: int for the number of main grid cells between the Inner
Surface and the Outer Surface. Defaults to 3.
pml_separation: int for the number of sub-grid cells between the Outer
Surface and the PML. Defaults to ratio // 2 + 2
subgrid_pml_thickness: int for the thickness of the PML on each of the
6 sides of the sub-grid. Defaults to 6.
interpolation: string for the degree of the interpolation scheme used
for spatial interpolation of the fields at the Inner
Surface. Defaults to Linear.
filter: boolean to turn on the 3-pole filter. Increases numerical
stability. Defaults to True.
"""
def __init__(self,
p1=None,
@@ -175,7 +173,7 @@ class SubGridHSG(SubGridBase):
pml_separation = ratio // 2 + 2
# copy over the optional parameters
# Copy over the optional parameters
kwargs['p1'] = p1
kwargs['p2'] = p2
kwargs['ratio'] = ratio