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

Changed name of rvalue function to roundvalue.

浏览代码
这个提交包含在：
Craig Warren
2016-01-11 13:33:41 +00:00
父节点 e193dfd5da
当前提交 e123e66817
共有 5 个文件被更改，包括 192 次插入 和 192 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

104
gprMax/geometry_primitives.pyx
查看文件

@@ -19,7 +19,7 @@
import numpy as np
cimport numpy as np
np.seterr(divide='raise')
from gprMax.utilities import rvalue
from gprMax.utilities import roundvalue
from gprMax.yee_cell_setget_rigid cimport set_rigid_Ex, set_rigid_Ey, set_rigid_Ez, set_rigid_Hx, set_rigid_Hy, set_rigid_Hz, set_rigid_E, unset_rigid_E, set_rigid_H, unset_rigid_H
@@ -271,28 +271,28 @@ cpdef build_triangle(float x1, float y1, float z1, float x2, float y2, float z2,
# Calculate a bounding box for the triangle
if normal == 'x':
area = 0.5 * (-z2 * y3 + z1 * (-y2 + y3) + y1 * (z2 - z3) + y2 * z3)
i1 = rvalue(np.amin([y1, y2, y3]) / dy) - 1
i2 = rvalue(np.amax([y1, y2, y3]) / dy) + 1
j1 = rvalue(np.amin([z1, z2, z3]) / dz) - 1
j2 = rvalue(np.amax([z1, z2, z3]) / dz) + 1
level = rvalue(x1 / dx)
thicknesscells = rvalue(thickness / dx)
i1 = roundvalue(np.amin([y1, y2, y3]) / dy) - 1
i2 = roundvalue(np.amax([y1, y2, y3]) / dy) + 1
j1 = roundvalue(np.amin([z1, z2, z3]) / dz) - 1
j2 = roundvalue(np.amax([z1, z2, z3]) / dz) + 1
level = roundvalue(x1 / dx)
thicknesscells = roundvalue(thickness / dx)
elif normal == 'y':
area = 0.5 * (-z2 * x3 + z1 * (-x2 + x3) + x1 * (z2 - z3) + x2 * z3)
i1 = rvalue(np.amin([x1, x2, x3]) / dx) - 1
i2 = rvalue(np.amax([x1, x2, x3]) / dx) + 1
j1 = rvalue(np.amin([z1, z2, z3]) / dz) - 1
j2 = rvalue(np.amax([z1, z2, z3]) / dz) + 1
level = rvalue(y1 /dy)
thicknesscells = rvalue(thickness / dy)
i1 = roundvalue(np.amin([x1, x2, x3]) / dx) - 1
i2 = roundvalue(np.amax([x1, x2, x3]) / dx) + 1
j1 = roundvalue(np.amin([z1, z2, z3]) / dz) - 1
j2 = roundvalue(np.amax([z1, z2, z3]) / dz) + 1
level = roundvalue(y1 /dy)
thicknesscells = roundvalue(thickness / dy)
elif normal == 'z':
area = 0.5 * (-y2 * x3 + y1 * (-x2 + x3) + x1 * (y2 - y3) + x2 * y3)
i1 = rvalue(np.amin([x1, x2, x3]) / dx) - 1
i2 = rvalue(np.amax([x1, x2, x3]) / dx) + 1
j1 = rvalue(np.amin([y1, y2, y3]) / dy) - 1
j2 = rvalue(np.amax([y1, y2, y3]) / dy) + 1
level = rvalue(z1 / dz)
thicknesscells = rvalue(thickness / dz)
i1 = roundvalue(np.amin([x1, x2, x3]) / dx) - 1
i2 = roundvalue(np.amax([x1, x2, x3]) / dx) + 1
j1 = roundvalue(np.amin([y1, y2, y3]) / dy) - 1
j2 = roundvalue(np.amax([y1, y2, y3]) / dy) + 1
level = roundvalue(z1 / dz)
thicknesscells = roundvalue(thickness / dz)
sign = np.sign(area)
@@ -356,11 +356,11 @@ cpdef build_cylindrical_sector(float ctr1, float ctr2, int level, float sectorst
if normal == 'x':
# Angles are defined from zero degrees on the positive y-axis going towards positive z-axis
y1 = rvalue((ctr1 - radius)/dy)
y2 = rvalue((ctr1 + radius)/dy)
z1 = rvalue((ctr2 - radius)/dz)
z2 = rvalue((ctr2 + radius)/dz)
thicknesscells = rvalue(thickness/dx)
y1 = roundvalue((ctr1 - radius)/dy)
y2 = roundvalue((ctr1 + radius)/dy)
z1 = roundvalue((ctr2 - radius)/dz)
z2 = roundvalue((ctr2 + radius)/dz)
thicknesscells = roundvalue(thickness/dx)
for y in range(y1, y2):
for z in range(z1, z2):
@@ -373,11 +373,11 @@ cpdef build_cylindrical_sector(float ctr1, float ctr2, int level, float sectorst
elif normal == 'y':
# Angles are defined from zero degrees on the positive x-axis going towards positive z-axis
x1 = rvalue((ctr1 - radius)/dx)
x2 = rvalue((ctr1 + radius)/dx)
z1 = rvalue((ctr2 - radius)/dz)
z2 = rvalue((ctr2 + radius)/dz)
thicknesscells = rvalue(thickness/dy)
x1 = roundvalue((ctr1 - radius)/dx)
x2 = roundvalue((ctr1 + radius)/dx)
z1 = roundvalue((ctr2 - radius)/dz)
z2 = roundvalue((ctr2 + radius)/dz)
thicknesscells = roundvalue(thickness/dy)
for x in range(x1, x2):
for z in range(z2, z2):
@@ -390,11 +390,11 @@ cpdef build_cylindrical_sector(float ctr1, float ctr2, int level, float sectorst
elif normal == 'z':
# Angles are defined from zero degrees on the positive x-axis going towards positive y-axis
x1 = rvalue((ctr1 - radius)/dx)
x2 = rvalue((ctr1 + radius)/dx)
y1 = rvalue((ctr2 - radius)/dy)
y2 = rvalue((ctr2 + radius)/dy)
thicknesscells = rvalue(thickness/dz)
x1 = roundvalue((ctr1 - radius)/dx)
x2 = roundvalue((ctr1 + radius)/dx)
y1 = roundvalue((ctr2 - radius)/dy)
y2 = roundvalue((ctr2 + radius)/dy)
thicknesscells = roundvalue(thickness/dz)
for x in range(x1, x2):
for y in range(y1, y2):
@@ -483,23 +483,23 @@ cpdef build_cylinder(float x1, float y1, float z1, float x2, float y2, float z2,
# Calculate a bounding box for the cylinder
if x1 < x2:
xs = rvalue((x1 - r) / dx) - 1
xf = rvalue((x2 + r) / dx) + 1
xs = roundvalue((x1 - r) / dx) - 1
xf = roundvalue((x2 + r) / dx) + 1
else:
xs = rvalue((x2 - r) / dx) - 1
xf = rvalue((x1 + r) / dx) + 1
xs = roundvalue((x2 - r) / dx) - 1
xf = roundvalue((x1 + r) / dx) + 1
if y1 < y2:
ys = rvalue((y1 - r) / dy) - 1
yf = rvalue((y2 + r) / dy) + 1
ys = roundvalue((y1 - r) / dy) - 1
yf = roundvalue((y2 + r) / dy) + 1
else:
ys = rvalue((y2 - r) / dy) - 1
yf = rvalue((y1 + r) / dy) + 1
ys = roundvalue((y2 - r) / dy) - 1
yf = roundvalue((y1 + r) / dy) + 1
if z1 < z2:
zs = rvalue((z1 - r) / dz) - 1
zf = rvalue((z2 + r) / dz) + 1
zs = roundvalue((z1 - r) / dz) - 1
zf = roundvalue((z2 + r) / dz) + 1
else:
zs = rvalue((z2 - r) / dz) - 1
zf = rvalue((z1 + r) / dz) + 1
zs = roundvalue((z2 - r) / dz) - 1
zf = roundvalue((z1 + r) / dz) + 1
# Set bounds to domain if they outside
if xs < 0:
@@ -577,12 +577,12 @@ cpdef build_sphere(int xc, int yc, int zc, float r, float dx, float dy, float dz
cdef int i, j, k, xs, xf, ys, yf, zs, zf
# Calculate a bounding box for sphere
xs = rvalue(((xc * dx) - r) / dx) - 1
xf = rvalue(((xc * dx) + r) / dx) + 1
ys = rvalue(((yc * dy) - r) / dy) - 1
yf = rvalue(((yc * dy) + r) / dy) + 1
zs = rvalue(((zc * dz) - r) / dz) - 1
zf = rvalue(((zc * dz) + r) / dz) + 1
xs = roundvalue(((xc * dx) - r) / dx) - 1
xf = roundvalue(((xc * dx) + r) / dx) + 1
ys = roundvalue(((yc * dy) - r) / dy) - 1
yf = roundvalue(((yc * dy) + r) / dy) + 1
zs = roundvalue(((zc * dz) - r) / dz) - 1
zf = roundvalue(((zc * dz) + r) / dz) + 1
# Set bounds to domain if they outside
if xs < 0:

16
gprMax/geometry_views.py
查看文件

@@ -20,7 +20,7 @@ import sys
import numpy as np
from struct import pack
from gprMax.utilities import rvalue
from gprMax.utilities import roundvalue
class GeometryView:
@@ -75,12 +75,12 @@ class GeometryView:
self.filename += '.vti'
# Calculate number of cells according to requested sampling
self.vtk_xscells = rvalue(self.xs / self.dx)
self.vtk_xfcells = rvalue(self.xf / self.dx)
self.vtk_yscells = rvalue(self.ys / self.dy)
self.vtk_yfcells = rvalue(self.yf / self.dy)
self.vtk_zscells = rvalue(self.zs / self.dz)
self.vtk_zfcells = rvalue(self.zf / self.dz)
self.vtk_xscells = roundvalue(self.xs / self.dx)
self.vtk_xfcells = roundvalue(self.xf / self.dx)
self.vtk_yscells = roundvalue(self.ys / self.dy)
self.vtk_yfcells = roundvalue(self.yf / self.dy)
self.vtk_zscells = roundvalue(self.zs / self.dz)
self.vtk_zfcells = roundvalue(self.zf / self.dz)
with open(self.filename, 'wb') as f:
f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
f.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
@@ -91,7 +91,7 @@ class GeometryView:
f.write('</CellData>\n</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
# Calculate number of bytes of appended data section
datasize = rvalue(np.dtype(np.uint32).itemsize * (self.vtk_nx / self.dx) * (self.vtk_ny / self.dy) * (self.vtk_nz / self.dz))
datasize = roundvalue(np.dtype(np.uint32).itemsize * (self.vtk_nx / self.dx) * (self.vtk_ny / self.dy) * (self.vtk_nz / self.dz))
# Write number of bytes of appended data as UInt32
f.write(pack('I', datasize))
for k in range(self.zs, self.zf, self.dz):

162
gprMax/input_cmds_geometry.py
查看文件

@@ -24,7 +24,7 @@ from gprMax.exceptions import CmdInputError
from gprMax.fractals import FractalSurface, FractalVolume, Grass
from gprMax.geometry_primitives import build_edge_x, build_edge_y, build_edge_z, build_face_yz, build_face_xz, build_face_xy, build_triangle, build_voxel, build_box, build_cylinder, build_cylindrical_sector, build_sphere
from gprMax.materials import Material
from gprMax.utilities import rvalue
from gprMax.utilities import roundvalue
def process_geometrycmds(geometry, G):
@@ -42,12 +42,12 @@ def process_geometrycmds(geometry, G):
if len(tmp) != 8:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters')
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
if xs < 0 or xs > G.nx:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs * G.dx))
@@ -113,12 +113,12 @@ def process_geometrycmds(geometry, G):
else:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
if xs < 0 or xs > G.nx:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs * G.dx))
@@ -233,15 +233,15 @@ def process_geometrycmds(geometry, G):
else:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
x1 = rvalue(float(tmp[1])/G.dx) * G.dx
y1 = rvalue(float(tmp[2])/G.dy) * G.dy
z1 = rvalue(float(tmp[3])/G.dz) * G.dz
x2 = rvalue(float(tmp[4])/G.dx) * G.dx
y2 = rvalue(float(tmp[5])/G.dy) * G.dy
z2 = rvalue(float(tmp[6])/G.dz) * G.dz
x3 = rvalue(float(tmp[7])/G.dx) * G.dx
y3 = rvalue(float(tmp[8])/G.dy) * G.dy
z3 = rvalue(float(tmp[9])/G.dz) * G.dz
x1 = roundvalue(float(tmp[1])/G.dx) * G.dx
y1 = roundvalue(float(tmp[2])/G.dy) * G.dy
z1 = roundvalue(float(tmp[3])/G.dz) * G.dz
x2 = roundvalue(float(tmp[4])/G.dx) * G.dx
y2 = roundvalue(float(tmp[5])/G.dy) * G.dy
z2 = roundvalue(float(tmp[6])/G.dz) * G.dz
x3 = roundvalue(float(tmp[7])/G.dx) * G.dx
y3 = roundvalue(float(tmp[8])/G.dy) * G.dy
z3 = roundvalue(float(tmp[9])/G.dz) * G.dz
thickness = float(tmp[10])
if x1 < 0 or x2 < 0 or x3 < 0 or x1 > G.nx or x2 > G.nx or x3 > G.nx:
@@ -353,12 +353,12 @@ def process_geometrycmds(geometry, G):
else:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
if xs < 0 or xs > G.nx:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs * G.dx))
@@ -445,12 +445,12 @@ def process_geometrycmds(geometry, G):
else:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
x1 = rvalue(float(tmp[1])/G.dx) * G.dx
y1 = rvalue(float(tmp[2])/G.dy) * G.dy
z1 = rvalue(float(tmp[3])/G.dz) * G.dz
x2 = rvalue(float(tmp[4])/G.dx) * G.dx
y2 = rvalue(float(tmp[5])/G.dy) * G.dy
z2 = rvalue(float(tmp[6])/G.dz) * G.dz
x1 = roundvalue(float(tmp[1])/G.dx) * G.dx
y1 = roundvalue(float(tmp[2])/G.dy) * G.dy
z1 = roundvalue(float(tmp[3])/G.dz) * G.dz
x2 = roundvalue(float(tmp[4])/G.dx) * G.dx
y2 = roundvalue(float(tmp[5])/G.dy) * G.dy
z2 = roundvalue(float(tmp[6])/G.dz) * G.dz
r = float(tmp[7])
if r <= 0:
@@ -595,21 +595,21 @@ def process_geometrycmds(geometry, G):
# yz-plane cylindrical sector
if normal == 'x':
ctr1 = rvalue(ctr1/G.dy) * G.dy
ctr2 = rvalue(ctr2/G.dz) * G.dz
level = rvalue(extent1/G.dx)
ctr1 = roundvalue(ctr1/G.dy) * G.dy
ctr2 = roundvalue(ctr2/G.dz) * G.dz
level = roundvalue(extent1/G.dx)
# xz-plane cylindrical sector
elif normal == 'y':
ctr1 = rvalue(ctr1/G.dx) * G.dx
ctr2 = rvalue(ctr2/G.dz) * G.dz
level = rvalue(extent1/G.dy)
ctr1 = roundvalue(ctr1/G.dx) * G.dx
ctr2 = roundvalue(ctr2/G.dz) * G.dz
level = roundvalue(extent1/G.dy)
# xy-plane cylindrical sector
elif normal == 'z':
ctr1 = rvalue(ctr1/G.dx) * G.dx
ctr2 = rvalue(ctr2/G.dy) * G.dy
level = rvalue(extent1/G.dz)
ctr1 = roundvalue(ctr1/G.dx) * G.dx
ctr2 = roundvalue(ctr2/G.dy) * G.dy
level = roundvalue(extent1/G.dz)
build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
@@ -651,9 +651,9 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
# Centre of sphere
xc = rvalue(float(tmp[1])/G.dx)
yc = rvalue(float(tmp[2])/G.dy)
zc = rvalue(float(tmp[3])/G.dz)
xc = roundvalue(float(tmp[1])/G.dx)
yc = roundvalue(float(tmp[2])/G.dy)
zc = roundvalue(float(tmp[3])/G.dz)
r = float(tmp[4])
# Look up requested materials in existing list of material instances
@@ -710,12 +710,12 @@ def process_geometrycmds(geometry, G):
else:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
if xs < 0 or xs > G.nx:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs * G.dx))
@@ -739,13 +739,13 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the y direction')
if float(tmp[10]) < 0:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the fractal weighting in the z direction')
if rvalue(tmp[11]) < 0:
if roundvalue(tmp[11]) < 0:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires a positive value for the number of bins')
# Find materials to use to build fractal volume, either from mixing models or normal materials
mixingmodel = next((x for x in G.mixingmodels if x.ID == tmp[12]), None)
material = next((x for x in G.materials if x.ID == tmp[12]), None)
nbins = rvalue(tmp[11])
nbins = roundvalue(tmp[11])
if mixingmodel:
if nbins == 1:
@@ -783,12 +783,12 @@ def process_geometrycmds(geometry, G):
# Only process rough surfaces for this fractal volume
if tmp[12] == volume.ID:
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
if xs < 0 or xs > G.nx:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs * G.dx))
@@ -817,7 +817,7 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
if xs != volume.xs and xs != volume.xf:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
fractalrange = (rvalue(float(tmp[10]) / G.dx), rvalue(float(tmp[11]) / G.dx))
fractalrange = (roundvalue(float(tmp[10]) / G.dx), roundvalue(float(tmp[11]) / G.dx))
# xminus surface
if xs == volume.xs:
if fractalrange[0] < 0:
@@ -834,7 +834,7 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
if ys != volume.ys and ys != volume.yf:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
fractalrange = (rvalue(float(tmp[10]) / G.dy), rvalue(float(tmp[11]) / G.dy))
fractalrange = (roundvalue(float(tmp[10]) / G.dy), roundvalue(float(tmp[11]) / G.dy))
# yminus surface
if ys == volume.ys:
if fractalrange[0] < 0:
@@ -851,7 +851,7 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
if zs != volume.zs and zs != volume.zf:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' can only be used on the external surfaces of a fractal box')
fractalrange = (rvalue(float(tmp[10]) / G.dz), rvalue(float(tmp[11]) / G.dz))
fractalrange = (roundvalue(float(tmp[10]) / G.dz), roundvalue(float(tmp[11]) / G.dz))
# zminus surface
if zs == volume.zs:
if fractalrange[0] < 0:
@@ -890,12 +890,12 @@ def process_geometrycmds(geometry, G):
# Only process surfaces for this fractal volume
if tmp[8] == volume.ID:
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
depth = float(tmp[7])
if xs < 0 or xs > G.nx:
@@ -927,7 +927,7 @@ def process_geometrycmds(geometry, G):
# xplus surface
elif xf == volume.xf:
requestedsurface = 'xplus'
filldepthcells = rvalue(depth / G.dx)
filldepthcells = roundvalue(depth / G.dx)
filldepth = filldepthcells * G.dx
elif ys == yf:
@@ -941,7 +941,7 @@ def process_geometrycmds(geometry, G):
# yplus surface
elif yf == volume.yf:
requestedsurface = 'yplus'
filldepthcells = rvalue(depth / G.dy)
filldepthcells = roundvalue(depth / G.dy)
filldepth = filldepthcells * G.dy
elif zs == zf:
@@ -955,7 +955,7 @@ def process_geometrycmds(geometry, G):
# zplus surface
elif zf == volume.zf:
requestedsurface = 'zplus'
filldepthcells = rvalue(depth / G.dz)
filldepthcells = roundvalue(depth / G.dz)
filldepth = filldepthcells * G.dz
else:
@@ -1004,12 +1004,12 @@ def process_geometrycmds(geometry, G):
# Only process grass for this fractal volume
if tmp[11] == volume.ID:
xs = rvalue(float(tmp[1])/G.dx)
xf = rvalue(float(tmp[4])/G.dx)
ys = rvalue(float(tmp[2])/G.dy)
yf = rvalue(float(tmp[5])/G.dy)
zs = rvalue(float(tmp[3])/G.dz)
zf = rvalue(float(tmp[6])/G.dz)
xs = roundvalue(float(tmp[1])/G.dx)
xf = roundvalue(float(tmp[4])/G.dx)
ys = roundvalue(float(tmp[2])/G.dy)
yf = roundvalue(float(tmp[5])/G.dy)
zs = roundvalue(float(tmp[3])/G.dz)
zf = roundvalue(float(tmp[6])/G.dz)
numblades = int(tmp[10])
if xs < 0 or xs > G.nx:
@@ -1037,7 +1037,7 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
if xs != volume.xs and xs != volume.xf:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must specify external surfaces on a fractal box')
fractalrange = (rvalue(float(tmp[8]) / G.dx), rvalue(float(tmp[9]) / G.dx))
fractalrange = (roundvalue(float(tmp[8]) / G.dx), roundvalue(float(tmp[9]) / G.dx))
# xminus surface
if xs == volume.xs:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' grass can only be specified on surfaces in the positive axis direction')
@@ -1052,7 +1052,7 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
if ys != volume.ys and ys != volume.yf:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must specify external surfaces on a fractal box')
fractalrange = (rvalue(float(tmp[8]) / G.dy), rvalue(float(tmp[9]) / G.dy))
fractalrange = (roundvalue(float(tmp[8]) / G.dy), roundvalue(float(tmp[9]) / G.dy))
# yminus surface
if ys == volume.ys:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' grass can only be specified on surfaces in the positive axis direction')
@@ -1067,7 +1067,7 @@ def process_geometrycmds(geometry, G):
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' dimensions are not specified correctly')
if zs != volume.zs and zs != volume.zf:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' must specify external surfaces on a fractal box')
fractalrange = (rvalue(float(tmp[8]) / G.dz), rvalue(float(tmp[9]) / G.dz))
fractalrange = (roundvalue(float(tmp[8]) / G.dz), roundvalue(float(tmp[9]) / G.dz))
# zminus surface
if zs == volume.zs:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' grass can only be specified on surfaces in the positive axis direction')

88
gprMax/input_cmds_multiuse.py
查看文件

@@ -23,7 +23,7 @@ from gprMax.pml import CFSParameter, CFS
from gprMax.receivers import Rx
from gprMax.snapshots import Snapshot
from gprMax.sources import VoltageSource, HertzianDipole, MagneticDipole, TransmissionLine
from gprMax.utilities import rvalue
from gprMax.utilities import roundvalue
from gprMax.waveforms import Waveform
@@ -72,9 +72,9 @@ def process_multicmds(multicmds, G):
# Check polarity & position parameters
if tmp[0].lower() not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
positionx = rvalue(float(tmp[1])/G.dx)
positiony = rvalue(float(tmp[2])/G.dy)
positionz = rvalue(float(tmp[3])/G.dz)
positionx = roundvalue(float(tmp[1])/G.dx)
positiony = roundvalue(float(tmp[2])/G.dy)
positionz = roundvalue(float(tmp[3])/G.dz)
resistance = float(tmp[4])
if positionx < 0 or positionx >= G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' x-coordinate is not within the model domain')
@@ -135,9 +135,9 @@ def process_multicmds(multicmds, G):
# Check polarity & position parameters
if tmp[0].lower() not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
positionx = rvalue(float(tmp[1])/G.dx)
positiony = rvalue(float(tmp[2])/G.dy)
positionz = rvalue(float(tmp[3])/G.dz)
positionx = roundvalue(float(tmp[1])/G.dx)
positiony = roundvalue(float(tmp[2])/G.dy)
positionz = roundvalue(float(tmp[3])/G.dz)
if positionx < 0 or positionx >= G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' x-coordinate is not within the model domain')
if positiony < 0 or positiony >= G.ny:
@@ -194,9 +194,9 @@ def process_multicmds(multicmds, G):
# Check polarity & position parameters
if tmp[0].lower() not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
positionx = rvalue(float(tmp[1])/G.dx)
positiony = rvalue(float(tmp[2])/G.dy)
positionz = rvalue(float(tmp[3])/G.dz)
positionx = roundvalue(float(tmp[1])/G.dx)
positiony = roundvalue(float(tmp[2])/G.dy)
positionz = roundvalue(float(tmp[3])/G.dz)
if positionx < 0 or positionx >= G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' x-coordinate is not within the model domain')
if positiony < 0 or positiony >= G.ny:
@@ -253,9 +253,9 @@ def process_multicmds(multicmds, G):
# Check polarity & position parameters
if tmp[0].lower() not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
positionx = rvalue(float(tmp[1])/G.dx)
positiony = rvalue(float(tmp[2])/G.dy)
positionz = rvalue(float(tmp[3])/G.dz)
positionx = roundvalue(float(tmp[1])/G.dx)
positiony = roundvalue(float(tmp[2])/G.dy)
positionz = roundvalue(float(tmp[3])/G.dz)
resistance = float(tmp[4])
if positionx < 0 or positionx >= G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' x-coordinate is not within the model domain')
@@ -314,9 +314,9 @@ def process_multicmds(multicmds, G):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' has an incorrect number of parameters')
# Check position parameters
positionx = rvalue(float(tmp[0])/G.dx)
positiony = rvalue(float(tmp[1])/G.dy)
positionz = rvalue(float(tmp[2])/G.dz)
positionx = roundvalue(float(tmp[0])/G.dx)
positiony = roundvalue(float(tmp[1])/G.dy)
positionz = roundvalue(float(tmp[2])/G.dz)
if positionx < 0 or positionx >= G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' x-coordinate is not within the model domain')
if positiony < 0 or positiony >= G.ny:
@@ -354,15 +354,15 @@ def process_multicmds(multicmds, G):
if len(tmp) != 9:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly nine parameters')
xs = rvalue(float(tmp[0])/G.dx)
xf = rvalue(float(tmp[3])/G.dx)
ys = rvalue(float(tmp[1])/G.dy)
yf = rvalue(float(tmp[4])/G.dy)
zs = rvalue(float(tmp[2])/G.dz)
zf = rvalue(float(tmp[5])/G.dz)
dx = rvalue(float(tmp[6])/G.dx)
dy = rvalue(float(tmp[7])/G.dy)
dz = rvalue(float(tmp[8])/G.dz)
xs = roundvalue(float(tmp[0])/G.dx)
xf = roundvalue(float(tmp[3])/G.dx)
ys = roundvalue(float(tmp[1])/G.dy)
yf = roundvalue(float(tmp[4])/G.dy)
zs = roundvalue(float(tmp[2])/G.dz)
zf = roundvalue(float(tmp[5])/G.dz)
dx = roundvalue(float(tmp[6])/G.dx)
dy = roundvalue(float(tmp[7])/G.dy)
dz = roundvalue(float(tmp[8])/G.dz)
if xs < 0 or xs >= G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs))
@@ -403,20 +403,20 @@ def process_multicmds(multicmds, G):
if len(tmp) != 11:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly eleven parameters')
xs = rvalue(float(tmp[0])/G.dx)
xf = rvalue(float(tmp[3])/G.dx)
ys = rvalue(float(tmp[1])/G.dy)
yf = rvalue(float(tmp[4])/G.dy)
zs = rvalue(float(tmp[2])/G.dz)
zf = rvalue(float(tmp[5])/G.dz)
dx = rvalue(float(tmp[6])/G.dx)
dy = rvalue(float(tmp[7])/G.dy)
dz = rvalue(float(tmp[8])/G.dz)
xs = roundvalue(float(tmp[0])/G.dx)
xf = roundvalue(float(tmp[3])/G.dx)
ys = roundvalue(float(tmp[1])/G.dy)
yf = roundvalue(float(tmp[4])/G.dy)
zs = roundvalue(float(tmp[2])/G.dz)
zf = roundvalue(float(tmp[5])/G.dz)
dx = roundvalue(float(tmp[6])/G.dx)
dy = roundvalue(float(tmp[7])/G.dy)
dz = roundvalue(float(tmp[8])/G.dz)
# If real floating point value given
if '.' in tmp[9] or 'e' in tmp[9]:
if float(tmp[9]) > 0:
time = rvalue((float(tmp[9]) / G.dt)) + 1
time = roundvalue((float(tmp[9]) / G.dt)) + 1
else:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time value must be greater than zero')
# If number of iterations given
@@ -640,15 +640,15 @@ def process_multicmds(multicmds, G):
if len(tmp) != 11:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly eleven parameters')
xs = rvalue(float(tmp[0])/G.dx)
xf = rvalue(float(tmp[3])/G.dx)
ys = rvalue(float(tmp[1])/G.dy)
yf = rvalue(float(tmp[4])/G.dy)
zs = rvalue(float(tmp[2])/G.dz)
zf = rvalue(float(tmp[5])/G.dz)
dx = rvalue(float(tmp[6])/G.dx)
dy = rvalue(float(tmp[7])/G.dy)
dz = rvalue(float(tmp[8])/G.dz)
xs = roundvalue(float(tmp[0])/G.dx)
xf = roundvalue(float(tmp[3])/G.dx)
ys = roundvalue(float(tmp[1])/G.dy)
yf = roundvalue(float(tmp[4])/G.dy)
zs = roundvalue(float(tmp[2])/G.dz)
zf = roundvalue(float(tmp[5])/G.dz)
dx = roundvalue(float(tmp[6])/G.dx)
dy = roundvalue(float(tmp[7])/G.dy)
dz = roundvalue(float(tmp[8])/G.dz)
if xs < 0 or xs > G.nx:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower x-coordinate {} is not within the model domain'.format(xs * G.dx))

14
gprMax/snapshots.py
查看文件

@@ -21,7 +21,7 @@ import numpy as np
from struct import pack
from gprMax.constants import floattype
from gprMax.utilities import rvalue
from gprMax.utilities import roundvalue
class Snapshot:
@@ -82,12 +82,12 @@ class Snapshot:
self.filename = G.inputdirectory + self.filename + '_' + str(modelrun) + '.vti'
# Calculate number of cells according to requested sampling
self.vtk_xscells = rvalue(self.xs / self.dx)
self.vtk_xfcells = rvalue(self.xf / self.dx)
self.vtk_yscells = rvalue(self.ys / self.dy)
self.vtk_yfcells = rvalue(self.yf / self.dz)
self.vtk_zscells = rvalue(self.zs / self.dz)
self.vtk_zfcells = rvalue(self.zf / self.dz)
self.vtk_xscells = roundvalue(self.xs / self.dx)
self.vtk_xfcells = roundvalue(self.xf / self.dx)
self.vtk_yscells = roundvalue(self.ys / self.dy)
self.vtk_yfcells = roundvalue(self.yf / self.dz)
self.vtk_zscells = roundvalue(self.zs / self.dz)
self.vtk_zfcells = roundvalue(self.zf / self.dz)
vtk_hfield_offset = 3 * np.dtype(floattype).itemsize * (self.vtk_xfcells - self.vtk_xscells) * (self.vtk_yfcells - self.vtk_yscells) * (self.vtk_zfcells - self.vtk_zscells) + np.dtype(np.uint32).itemsize
# vtk_current_offset = 2 * vtk_hfield_offset