你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 04:56:51 +08:00
Merge remote-tracking branch 'origin/master'
这个提交包含在：
Craig Warren
@@ -11,6 +11,7 @@ dependencies:
|
||||
- mpi4py
|
||||
- numpy
|
||||
- scipy
|
||||
- lxml
|
||||
|
||||
- pip:
|
||||
- psutil
|
||||
|
||||
@@ -27,8 +27,6 @@ AustinMan and AustinWoman (http://bit.ly/AustinMan) are open source electromagne
|
||||
|
||||
FDTD geometry mesh showing the head of the AustinMan model (2x2x2mm :math:`^3`).
|
||||
|
||||
The AustinMan and AustinWoman models are not currently included in the user libraries sub-package, however they can be downloaded from http://bit.ly/AustinMan.
|
||||
|
||||
The following whole body models are available.
|
||||
|
||||
=========== ========================== ==================
|
||||
@@ -44,16 +42,25 @@ AustinWoman 2x2x2 342 x 187 x 865
|
||||
AustinWoman 1x1x1 683 x 374 x 1730
|
||||
=========== ========================== ==================
|
||||
|
||||
Package overview
|
||||
================
|
||||
|
||||
.. code-block:: none
|
||||
|
||||
AustinManWoman_gprMax_materials.txt
|
||||
AustinManWoman_gprMax_materials_dispersive.txt
|
||||
head_only_hdf5.py
|
||||
|
||||
* ``AustinManWoman_gprMax_materials.txt`` is a text file containing non-dispersive material properties at 900 MHz (http://niremf.ifac.cnr.it/tissprop/).
|
||||
* ``AustinManWoman_gprMax_materials_dispersive.txt`` is a text file containing dispersive material properties using a 3-pole Debye model (http://dx.doi.org/10.1109/LMWC.2011.2180371). Note the main body tissues are described using a 3-pole Debye model, nbut ot all materials have a dispersive description.
|
||||
* ``head_only_hdf5.py`` is a script to assist with creating a model of only the head from a full body AustinMan/Woman model.
|
||||
|
||||
How to use the models
|
||||
=====================
|
||||
|
||||
From http://bit.ly/AustinMan:
|
||||
The AustinMan and AustinWoman models themselves are not included in the user libraries sub-package.
|
||||
|
||||
* Download a HDF5 file (.h5) of AustinMan or AustinWoman at the resolution you wish to use.
|
||||
* Download a text file of material descriptions for gprMax, either
|
||||
|
||||
* ``AustinManWoman_gprMax_materials.txt`` for non-dispersive material properties at 900 MHz (http://niremf.ifac.cnr.it/tissprop/).
|
||||
* ``AustinManWoman_gprMax_materials_dispersive.txt`` for dispersive material properties that feature a 3-pole Debye model (http://dx.doi.org/10.1109/LMWC.2011.2180371). Not all materials have a dispersive description.
|
||||
* Download a HDF5 file (.h5) of AustinMan or AustinWoman at the resolution you wish to use from http://bit.ly/AustinMan
|
||||
|
||||
To insert either AustinMan or AustinWoman models into a simulation use the ``#geometry_objects_file``.
|
||||
|
||||
@@ -66,7 +73,7 @@ To insert a 2x2x2mm :math:`^3` AustinMan with the lower left corner 40mm from th
|
||||
|
||||
#geometry_objects_file: 0.04 0.04 0.04 ../user_libs/AustinManWoman/AustinMan_v2.3_2x2x2.h5 ../user_libs/AustinManWoman/AustinManWoman_gprMax_materials_dispersive.txt
|
||||
|
||||
For further information on the `#geometry_objects_file` see the section on object contruction commands in the :ref:`Input commands section <commands>`.
|
||||
For further information on the ``#geometry_objects_file`` see the section on object contruction commands in the :ref:`Input commands section <commands>`.
|
||||
|
||||
.. figure:: images/user_libs/AustinMan.png
|
||||
:width: 300 px
|
||||
|
||||
@@ -1 +1,4 @@
|
||||
from gprMax._version import __version__
|
||||
from ._version import __version__
|
||||
from .gprMax import api as run
|
||||
|
||||
__name__ = 'gprMax'
|
||||
|
||||
@@ -19,13 +19,14 @@
|
||||
import os, sys
|
||||
import numpy as np
|
||||
from struct import pack
|
||||
from gprMax.xdmf import write_output_file
|
||||
|
||||
from gprMax.utilities import round_value
|
||||
|
||||
|
||||
class GeometryView(object):
|
||||
"""Views of the geometry of the model."""
|
||||
|
||||
|
||||
if sys.byteorder == 'little':
|
||||
byteorder = 'LittleEndian'
|
||||
else:
|
||||
@@ -39,7 +40,7 @@ class GeometryView(object):
|
||||
filename (str): Filename to save to.
|
||||
type (str): Either 'n' for a per cell geometry view, or 'f' for a per cell edge geometry view.
|
||||
"""
|
||||
|
||||
|
||||
self.xs = xs
|
||||
self.ys = ys
|
||||
self.zs = zs
|
||||
@@ -55,26 +56,30 @@ class GeometryView(object):
|
||||
self.basefilename = filename
|
||||
self.type = type
|
||||
|
||||
def write_xdmf(self, modelrun, numbermodelruns, G):
|
||||
filename = self.filename[:-4]
|
||||
write_output_file(filename, G)
|
||||
|
||||
def write_vtk(self, modelrun, numbermodelruns, G):
|
||||
"""Writes the geometry information to a VTK file. Either ImageData (.vti) for a per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view.
|
||||
|
||||
|
||||
N.B. No Python 3 support for VTK at time of writing (03/2015)
|
||||
|
||||
|
||||
Args:
|
||||
modelrun (int): Current model run number.
|
||||
numbermodelruns (int): Total number of model runs.
|
||||
G (class): Grid class instance - holds essential parameters describing the model.
|
||||
"""
|
||||
|
||||
|
||||
# Construct filename from user-supplied name and model run number
|
||||
if numbermodelruns == 1:
|
||||
self.filename = os.path.abspath(os.path.join(G.inputdirectory, self.basefilename))
|
||||
else:
|
||||
self.filename = os.path.abspath(os.path.join(G.inputdirectory, self.basefilename + str(modelrun)))
|
||||
|
||||
|
||||
if self.type == 'n':
|
||||
self.filename += '.vti'
|
||||
|
||||
|
||||
# Calculate number of cells according to requested sampling for geometry view
|
||||
self.vtk_xscells = round_value(self.xs / self.dx)
|
||||
self.vtk_xfcells = round_value(self.xf / self.dx)
|
||||
@@ -85,7 +90,7 @@ class GeometryView(object):
|
||||
self.vtk_nxcells = self.vtk_xfcells - self.vtk_xscells
|
||||
self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells
|
||||
self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells
|
||||
|
||||
|
||||
# Create arrays and add numeric IDs for PML, sources and receivers (0 is not set, 1 is PML, srcs and rxs numbered thereafter)
|
||||
self.srcs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
|
||||
self.rxs = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
|
||||
@@ -95,10 +100,10 @@ class GeometryView(object):
|
||||
self.srcs_pml[src.xcoord, src.ycoord, src.zcoord] = index + 2
|
||||
for index, rx in enumerate(G.rxs):
|
||||
self.rxs[rx.xcoord, rx.ycoord, rx.zcoord] = index + 1
|
||||
|
||||
|
||||
vtk_srcs_pml_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
|
||||
vtk_rxs_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize + (np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
|
||||
|
||||
|
||||
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'))
|
||||
@@ -110,7 +115,7 @@ class GeometryView(object):
|
||||
f.write('<DataArray type="Int8" Name="Receivers" format="appended" offset="{}" />\n'.format(vtk_rxs_offset).encode('utf-8'))
|
||||
f.write('</CellData>\n'.encode('utf-8'))
|
||||
f.write('</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
||||
|
||||
|
||||
# Write material IDs
|
||||
datasize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||
# Write number of bytes of appended data as UInt32
|
||||
@@ -119,7 +124,7 @@ class GeometryView(object):
|
||||
for j in range(self.ys, self.yf, self.dy):
|
||||
for i in range(self.xs, self.xf, self.dx):
|
||||
f.write(pack('I', G.solid[i, j, k]))
|
||||
|
||||
|
||||
# Write source/PML IDs
|
||||
datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||
f.write(pack('I', datasize))
|
||||
@@ -127,7 +132,7 @@ class GeometryView(object):
|
||||
for j in range(self.ys, self.yf, self.dy):
|
||||
for i in range(self.xs, self.xf, self.dx):
|
||||
f.write(pack('b', self.srcs_pml[i, j, k]))
|
||||
|
||||
|
||||
# Write receiver IDs
|
||||
datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||
f.write(pack('I', datasize))
|
||||
@@ -137,34 +142,34 @@ class GeometryView(object):
|
||||
f.write(pack('b', self.rxs[i, j, k]))
|
||||
|
||||
f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
|
||||
|
||||
|
||||
self.write_gprmax_info(f, G)
|
||||
|
||||
elif self.type == 'f':
|
||||
self.filename += '.vtp'
|
||||
|
||||
|
||||
vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
|
||||
vtk_numpoint_components = 3
|
||||
vtk_numlines = 2 * self.nx * self.ny + 2 * self.ny * self.nz + 2 * self.nx * self.nz + 3 * self.nx * self.ny * self.nz + self.nx + self.ny + self.nz
|
||||
vtk_numline_components = 2;
|
||||
|
||||
|
||||
vtk_connectivity_offset = round_value((vtk_numpoints * vtk_numpoint_components * np.dtype(np.float32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||
vtk_offsets_offset = round_value(vtk_connectivity_offset + (vtk_numlines * vtk_numline_components * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||
vtk_materials_offset = round_value(vtk_offsets_offset + (vtk_numlines * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||
|
||||
|
||||
with open(self.filename, 'wb') as f:
|
||||
f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
|
||||
f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
|
||||
f.write('<PolyData>\n<Piece NumberOfPoints="{}" NumberOfVerts="0" NumberOfLines="{}" NumberOfStrips="0" NumberOfPolys="0">\n'.format(vtk_numpoints, vtk_numlines).encode('utf-8'))
|
||||
|
||||
|
||||
f.write('<Points>\n<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="0" />\n</Points>\n'.encode('utf-8'))
|
||||
f.write('<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{}" />\n'.format(vtk_connectivity_offset).encode('utf-8'))
|
||||
f.write('<DataArray type="UInt32" Name="offsets" format="appended" offset="{}" />\n</Lines>\n'.format(vtk_offsets_offset).encode('utf-8'))
|
||||
|
||||
|
||||
f.write('<CellData Scalars="Material">\n'.encode('utf-8'))
|
||||
f.write('<DataArray type="UInt32" Name="Material" format="appended" offset="{}" />\n'.format(vtk_materials_offset).encode('utf-8'))
|
||||
f.write('</CellData>\n'.encode('utf-8'))
|
||||
|
||||
|
||||
f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
||||
|
||||
# Write points
|
||||
@@ -222,7 +227,7 @@ class GeometryView(object):
|
||||
for j in range(self.ys, self.yf + 1):
|
||||
for k in range(self.zs, self.zf + 1):
|
||||
f.write(pack('I', G.ID[0, i, j, k]))
|
||||
|
||||
|
||||
for i in range(self.xs, self.xf + 1):
|
||||
for j in range(self.ys, self.yf):
|
||||
for k in range(self.zs, self.zf + 1):
|
||||
@@ -240,7 +245,7 @@ class GeometryView(object):
|
||||
|
||||
def write_gprmax_info(self, f, G, materialsonly=False):
|
||||
"""Writes gprMax specific information relating material, source, and receiver names to numeric identifiers.
|
||||
|
||||
|
||||
Args:
|
||||
f (filehandle): VTK file.
|
||||
G (class): Grid class instance - holds essential parameters describing the model.
|
||||
@@ -257,5 +262,3 @@ class GeometryView(object):
|
||||
for index, rx in enumerate(G.rxs):
|
||||
f.write('<Receivers name="{}">{}</Receivers>\n'.format(rx.ID, index + 1).encode('utf-8'))
|
||||
f.write('</gprMax>\n'.encode('utf-8'))
|
||||
|
||||
|
||||
|
||||
143
gprMax/gprMax.py
143
gprMax/gprMax.py
@@ -18,33 +18,38 @@
|
||||
|
||||
"""gprMax.gprMax: provides entry point main()."""
|
||||
|
||||
import argparse, datetime, itertools, os, psutil, sys
|
||||
from time import perf_counter
|
||||
from enum import Enum
|
||||
import argparse
|
||||
import datetime
|
||||
import itertools
|
||||
import os
|
||||
import psutil
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
import gprMax
|
||||
from gprMax.constants import c, e0, m0, z0, floattype
|
||||
from gprMax.exceptions import GeneralError
|
||||
from gprMax.fields_update import update_electric, update_magnetic, update_electric_dispersive_multipole_A, update_electric_dispersive_multipole_B, update_electric_dispersive_1pole_A, update_electric_dispersive_1pole_B
|
||||
from gprMax.grid import FDTDGrid, dispersion_check
|
||||
from gprMax.input_cmds_geometry import process_geometrycmds
|
||||
from gprMax.input_cmds_file import process_python_include_code, write_processed_file, check_cmd_names
|
||||
from gprMax.input_cmds_multiuse import process_multicmds
|
||||
from gprMax.input_cmds_singleuse import process_singlecmds
|
||||
from gprMax.materials import Material
|
||||
from gprMax.writer_hdf5 import prepare_hdf5, write_hdf5
|
||||
from gprMax.pml import build_pmls, update_electric_pml, update_magnetic_pml
|
||||
from gprMax.utilities import update_progress, logo, human_size
|
||||
from gprMax.yee_cell_build import build_electric_components, build_magnetic_components
|
||||
from time import perf_counter
|
||||
from enum import Enum
|
||||
from ._version import __version__
|
||||
from .constants import c, e0, m0, z0
|
||||
from .exceptions import GeneralError
|
||||
from .fields_update import update_electric, update_magnetic, update_electric_dispersive_multipole_A, update_electric_dispersive_multipole_B, update_electric_dispersive_1pole_A, update_electric_dispersive_1pole_B
|
||||
from .grid import FDTDGrid, dispersion_check
|
||||
from .input_cmds_geometry import process_geometrycmds
|
||||
from .input_cmds_file import process_python_include_code, write_processed_file, check_cmd_names
|
||||
from .input_cmds_multiuse import process_multicmds
|
||||
from .input_cmds_singleuse import process_singlecmds
|
||||
from .materials import Material
|
||||
from .writer_hdf5 import prepare_hdf5, write_hdf5
|
||||
from .pml import build_pmls, update_electric_pml, update_magnetic_pml
|
||||
from .utilities import update_progress, logo, human_size
|
||||
from .yee_cell_build import build_electric_components, build_magnetic_components
|
||||
|
||||
|
||||
def main():
|
||||
"""This is the main function for gprMax."""
|
||||
|
||||
|
||||
# Print gprMax logo, version, and licencing/copyright information
|
||||
logo(gprMax.__version__ + ' (Bowmore)')
|
||||
logo(__version__ + ' (Bowmore)')
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(prog='gprMax', description='Electromagnetic modelling software based on the Finite-Difference Time-Domain (FDTD) method')
|
||||
@@ -57,10 +62,35 @@ def main():
|
||||
parser.add_argument('--write-processed', action='store_true', default=False, help='write an input file after any Python code and include commands in the original input file have been processed')
|
||||
parser.add_argument('--opt-taguchi', action='store_true', default=False, help='optimise parameters using the Taguchi optimisation method')
|
||||
args = parser.parse_args()
|
||||
|
||||
run_main(args)
|
||||
|
||||
|
||||
def api(inputfile, n=1, mpi=False, benchmark=False, geometry_only=False, geometry_fixed=False, write_processed=False, opt_taguchi=False):
|
||||
"""If you have installed gprMax as a module this is the entry point"""
|
||||
class ImportArguments:
|
||||
pass
|
||||
|
||||
args = ImportArguments()
|
||||
|
||||
args.inputfile = inputfile
|
||||
args.n = n
|
||||
args.mpi = mpi
|
||||
args.benchmark = benchmark
|
||||
args.geometry_only = geometry_only
|
||||
args.geometry_fixed = geometry_fixed
|
||||
args.write_processed = write_processed
|
||||
args.opt_taguchi = opt_taguchi
|
||||
|
||||
run_main(args)
|
||||
|
||||
|
||||
def run_main(args):
|
||||
|
||||
numbermodelruns = args.n
|
||||
inputdirectory = os.path.dirname(os.path.abspath(args.inputfile))
|
||||
inputfile = os.path.abspath(os.path.join(inputdirectory, os.path.basename(args.inputfile)))
|
||||
|
||||
|
||||
# Create a separate namespace that users can access in any Python code blocks in the input file
|
||||
usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'input_directory': inputdirectory}
|
||||
|
||||
@@ -93,7 +123,7 @@ def main():
|
||||
|
||||
def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None):
|
||||
"""Run standard simulation - models are run one after another and each model is parallelised with OpenMP
|
||||
|
||||
|
||||
Args:
|
||||
args (dict): Namespace with command line arguments
|
||||
numbermodelruns (int): Total number of model runs.
|
||||
@@ -118,20 +148,20 @@ def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
||||
|
||||
def run_benchmark_sim(args, inputfile, usernamespace):
|
||||
"""Run standard simulation in benchmarking mode - models are run one after another and each model is parallelised with OpenMP
|
||||
|
||||
|
||||
Args:
|
||||
args (dict): Namespace with command line arguments
|
||||
inputfile (str): Name of the input file to open.
|
||||
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
||||
"""
|
||||
|
||||
|
||||
# Number of threads to test - start from max physical CPU cores and divide in half until 1
|
||||
thread = psutil.cpu_count(logical=False)
|
||||
threads = [thread]
|
||||
while not thread%2:
|
||||
thread /= 2
|
||||
threads.append(int(thread))
|
||||
|
||||
|
||||
benchtimes = np.zeros(len(threads))
|
||||
|
||||
numbermodelruns = len(threads)
|
||||
@@ -151,7 +181,7 @@ def run_benchmark_sim(args, inputfile, usernamespace):
|
||||
|
||||
def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None):
|
||||
"""Run mixed mode MPI/OpenMP simulation - MPI task farm for models with each model parallelised with OpenMP
|
||||
|
||||
|
||||
Args:
|
||||
args (dict): Namespace with command line arguments
|
||||
numbermodelruns (int): Total number of model runs.
|
||||
@@ -181,7 +211,7 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
||||
data = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=status)
|
||||
source = status.Get_source()
|
||||
tag = status.Get_tag()
|
||||
|
||||
|
||||
if tag == tags.READY.value: # Worker is ready, so send it a task
|
||||
if modelrun < numbermodelruns + 1:
|
||||
comm.send(modelrun, dest=source, tag=tags.START.value)
|
||||
@@ -189,10 +219,10 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
||||
modelrun += 1
|
||||
else:
|
||||
comm.send(None, dest=source, tag=tags.EXIT.value)
|
||||
|
||||
|
||||
elif tag == tags.DONE.value:
|
||||
print('Worker {}: completed.'.format(source))
|
||||
|
||||
|
||||
elif tag == tags.EXIT.value:
|
||||
print('Worker {}: exited.'.format(source))
|
||||
closedworkers += 1
|
||||
@@ -203,7 +233,7 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
||||
comm.send(None, dest=0, tag=tags.READY.value)
|
||||
modelrun = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) # Receive a model number to run from the master
|
||||
tag = status.Get_tag()
|
||||
|
||||
|
||||
# Run a model
|
||||
if tag == tags.START.value:
|
||||
if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
|
||||
@@ -213,10 +243,10 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
||||
modelusernamespace.update({'optparams': tmp})
|
||||
else:
|
||||
modelusernamespace = usernamespace
|
||||
|
||||
|
||||
run_model(args, modelrun, numbermodelruns, inputfile, modelusernamespace)
|
||||
comm.send(None, dest=0, tag=tags.DONE.value)
|
||||
|
||||
|
||||
elif tag == tags.EXIT.value:
|
||||
break
|
||||
|
||||
@@ -225,39 +255,39 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
||||
|
||||
def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
"""Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop.
|
||||
|
||||
|
||||
Args:
|
||||
args (dict): Namespace with command line arguments
|
||||
modelrun (int): Current model run number.
|
||||
numbermodelruns (int): Total number of model runs.
|
||||
inputfile (str): Name of the input file to open.
|
||||
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
||||
|
||||
|
||||
Returns:
|
||||
tsolve (int): Length of time (seconds) of main FDTD calculations
|
||||
"""
|
||||
|
||||
|
||||
# Monitor memory usage
|
||||
p = psutil.Process()
|
||||
|
||||
|
||||
# Declare variable to hold FDTDGrid class
|
||||
global G
|
||||
|
||||
|
||||
# Normal model reading/building process; bypassed if geometry information to be reused
|
||||
if not 'G' in globals():
|
||||
print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile))
|
||||
|
||||
|
||||
# Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
|
||||
usernamespace['current_model_run'] = modelrun
|
||||
print('Constants/variables available for Python scripting: {}\n'.format(usernamespace))
|
||||
|
||||
|
||||
# Read input file and process any Python or include commands
|
||||
processedlines = process_python_include_code(inputfile, usernamespace)
|
||||
|
||||
|
||||
# Write a file containing the input commands after Python or include commands have been processed
|
||||
if args.write_processed:
|
||||
write_processed_file(inputfile, modelrun, numbermodelruns, processedlines)
|
||||
|
||||
|
||||
# Check validity of command names and that essential commands are present
|
||||
singlecmds, multicmds, geometry = check_cmd_names(processedlines)
|
||||
|
||||
@@ -294,7 +324,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
|
||||
# Build the PML and calculate initial coefficients
|
||||
build_pmls(G)
|
||||
|
||||
|
||||
# Build the model, i.e. set the material properties (ID) for every edge of every Yee cell
|
||||
tbuildstart = perf_counter()
|
||||
build_electric_components(G.solid, G.rigidE, G.ID, G)
|
||||
@@ -305,7 +335,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
# Process any voltage sources (that have resistance) to create a new material at the source location
|
||||
for voltagesource in G.voltagesources:
|
||||
voltagesource.create_material(G)
|
||||
|
||||
|
||||
# Initialise arrays of update coefficients to pass to update functions
|
||||
G.initialise_std_update_coeff_arrays()
|
||||
|
||||
@@ -319,22 +349,22 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
print('ID\tName\t\tProperties')
|
||||
print('{}'.format('-'*50))
|
||||
for material in G.materials:
|
||||
|
||||
|
||||
# Calculate update coefficients for material
|
||||
material.calculate_update_coeffsE(G)
|
||||
material.calculate_update_coeffsH(G)
|
||||
|
||||
|
||||
# Store all update coefficients together
|
||||
G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
|
||||
G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
|
||||
|
||||
|
||||
# Store coefficients for any dispersive materials
|
||||
if Material.maxpoles != 0:
|
||||
z = 0
|
||||
for pole in range(Material.maxpoles):
|
||||
G.updatecoeffsdispersive[material.numID, z:z+3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole]
|
||||
z += 3
|
||||
|
||||
|
||||
if G.messages:
|
||||
if material.deltaer and material.tau:
|
||||
tmp = 'delta_epsr={}, tau={} secs; '.format(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer), ', '.join('{:g}'.format(tau) for tau in material.tau))
|
||||
@@ -352,10 +382,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution))
|
||||
|
||||
# If geometry information to be reused between model runs
|
||||
else:
|
||||
else:
|
||||
# Clear arrays for field components
|
||||
G.initialise_field_arrays()
|
||||
|
||||
|
||||
# Clear arrays for fields in PML
|
||||
for pml in G.pmls:
|
||||
pml.initialise_field_arrays()
|
||||
@@ -377,7 +407,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
receiver.xcoord = receiver.xcoordbase + (modelrun - 1) * G.rxstepx
|
||||
receiver.ycoord = receiver.ycoordbase + (modelrun - 1) * G.rxstepy
|
||||
receiver.zcoord = receiver.zcoordbase + (modelrun - 1) * G.rxstepz
|
||||
|
||||
|
||||
# Write files for any geometry views
|
||||
if not G.geometryviews and args.geometry_only:
|
||||
raise GeneralError('No geometry views found.')
|
||||
@@ -385,16 +415,17 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
tgeostart = perf_counter()
|
||||
for geometryview in G.geometryviews:
|
||||
geometryview.write_vtk(modelrun, numbermodelruns, G)
|
||||
geometryview.write_xdmf(modelrun, numbermodelruns, G)
|
||||
tgeoend = perf_counter()
|
||||
print('\nGeometry file(s) written in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tgeoend - tgeostart))))
|
||||
|
||||
# Run simulation (if not doing geometry only)
|
||||
if not args.geometry_only:
|
||||
|
||||
|
||||
# Prepare any snapshot files
|
||||
for snapshot in G.snapshots:
|
||||
snapshot.prepare_vtk_imagedata(modelrun, numbermodelruns, G)
|
||||
|
||||
|
||||
# Prepare output file
|
||||
inputfileparts = os.path.splitext(inputfile)
|
||||
if numbermodelruns == 1:
|
||||
@@ -415,10 +446,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
for timestep in range(G.iterations):
|
||||
if timestep == 0:
|
||||
tstepstart = perf_counter()
|
||||
|
||||
|
||||
# Write field outputs to file
|
||||
write_hdf5(f, timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
|
||||
|
||||
|
||||
# Write any snapshots to file
|
||||
for snapshot in G.snapshots:
|
||||
if snapshot.time == timestep + 1:
|
||||
@@ -451,7 +482,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
|
||||
# Increment absolute time value
|
||||
abstime += 0.5 * G.dt
|
||||
|
||||
|
||||
# Update magnetic field components
|
||||
update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
|
||||
|
||||
@@ -466,7 +497,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
|
||||
# Increment absolute time value
|
||||
abstime += 0.5 * G.dt
|
||||
|
||||
|
||||
# Calculate time for two iterations, used to estimate overall runtime
|
||||
if timestep == 1:
|
||||
tstepend = perf_counter()
|
||||
@@ -476,7 +507,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||
sys.stdout.flush()
|
||||
elif timestep > 1:
|
||||
update_progress((timestep + 1) / G.iterations)
|
||||
|
||||
|
||||
# Close output file
|
||||
f.close()
|
||||
|
||||
|
||||
@@ -23,9 +23,34 @@ from gprMax.constants import c, floattype, complextype
|
||||
from gprMax.materials import Material
|
||||
|
||||
|
||||
class FDTDGrid(object):
|
||||
class Grid():
|
||||
|
||||
def __init__(self, grid):
|
||||
self.nx = grid.shape[0]
|
||||
self.ny = grid.shape[1]
|
||||
self.nz = grid.shape[2]
|
||||
self.grid = grid
|
||||
|
||||
def n_edges(self):
|
||||
l = self.nx
|
||||
m = self.ny
|
||||
n = self.nz
|
||||
e = (l * m * (n - 1)) + (m * n * (l - 1)) + (l * n * (m - 1))
|
||||
return e
|
||||
|
||||
def n_nodes(self):
|
||||
return self.nx * self.ny * self.nz
|
||||
|
||||
def n_cells(self):
|
||||
return (self.nx - 1) * (self.ny - 1) * (self.nz - 1)
|
||||
|
||||
def get(self, i, j, k):
|
||||
return self.grid[i, j, k]
|
||||
|
||||
|
||||
class FDTDGrid(Grid):
|
||||
"""Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters."""
|
||||
|
||||
|
||||
def __init__(self):
|
||||
self.inputfilename = ''
|
||||
self.inputdirectory = ''
|
||||
@@ -62,7 +87,7 @@ class FDTDGrid(object):
|
||||
self.rxstepy = 0
|
||||
self.rxstepz = 0
|
||||
self.snapshots = []
|
||||
|
||||
|
||||
def initialise_geometry_arrays(self):
|
||||
"""Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid);
|
||||
and an array for cell edge IDs (ID). Solid and ID arrays are initialised to free_space (one); rigid arrays to allow dielectric smoothing (zero).
|
||||
@@ -72,7 +97,7 @@ class FDTDGrid(object):
|
||||
self.rigidH = np.zeros((6, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.int8)
|
||||
self.IDlookup = {'Ex': 0, 'Ey': 1, 'Ez': 2, 'Hx': 3, 'Hy': 4, 'Hz': 5}
|
||||
self.ID = np.ones((6, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.uint32)
|
||||
|
||||
|
||||
def initialise_field_arrays(self):
|
||||
"""Initialise arrays for the electric and magnetic field components."""
|
||||
self.Ex = np.zeros((self.nx, self.ny + 1, self.nz + 1), dtype=floattype)
|
||||
@@ -81,7 +106,7 @@ class FDTDGrid(object):
|
||||
self.Hx = np.zeros((self.nx + 1, self.ny, self.nz), dtype=floattype)
|
||||
self.Hy = np.zeros((self.nx, self.ny + 1, self.nz), dtype=floattype)
|
||||
self.Hz = np.zeros((self.nx, self.ny, self.nz + 1), dtype=floattype)
|
||||
|
||||
|
||||
def initialise_std_update_coeff_arrays(self):
|
||||
"""Initialise arrays for storing update coefficients."""
|
||||
self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=floattype)
|
||||
@@ -97,32 +122,32 @@ class FDTDGrid(object):
|
||||
|
||||
def dispersion_check(G):
|
||||
"""Check for potential numerical dispersion. Is the smallest wavelength present in the simulation discretised by at least a factor of 10
|
||||
|
||||
|
||||
Args:
|
||||
G (class): Grid class instance - holds essential parameters describing the model.
|
||||
|
||||
|
||||
Returns:
|
||||
resolution (float): Potential numerical dispersion
|
||||
"""
|
||||
|
||||
|
||||
# Minimum number of spatial steps to resolve smallest wavelength
|
||||
resolvedsteps = 10
|
||||
|
||||
|
||||
# Find maximum frequency
|
||||
maxfreqs = []
|
||||
for waveform in G.waveforms:
|
||||
|
||||
if waveform.type == 'sine' or waveform.type == 'contsine':
|
||||
maxfreqs.append(4 * waveform.freq)
|
||||
|
||||
|
||||
elif waveform.type =='impulse':
|
||||
pass
|
||||
|
||||
|
||||
else:
|
||||
# User-defined waveform
|
||||
if waveform.type == 'user':
|
||||
waveformvalues = waveform.uservalues
|
||||
|
||||
|
||||
# Built-in waveform
|
||||
else:
|
||||
time = np.linspace(0, 1, G.iterations)
|
||||
@@ -133,31 +158,31 @@ def dispersion_check(G):
|
||||
while not timeiter.finished:
|
||||
waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt)
|
||||
timeiter.iternext()
|
||||
|
||||
|
||||
# Calculate magnitude of frequency spectra of waveform
|
||||
power = 20 * np.log10(np.abs(np.fft.fft(waveformvalues))**2)
|
||||
freqs = np.fft.fftfreq(power.size, d=G.dt)
|
||||
|
||||
# Shift powers so that frequency with maximum power is at zero decibels
|
||||
power -= np.amax(power)
|
||||
|
||||
|
||||
# Set maximum frequency to -60dB from maximum power
|
||||
freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
|
||||
maxfreqs.append(freqs[freq])
|
||||
|
||||
if maxfreqs:
|
||||
maxfreq = max(maxfreqs)
|
||||
|
||||
|
||||
# Find minimum wavelength
|
||||
ers = [material.er for material in G.materials]
|
||||
maxer = max(ers)
|
||||
|
||||
# Minimum velocity
|
||||
minvelocity = c / np.sqrt(maxer)
|
||||
|
||||
|
||||
# Minimum wavelength
|
||||
minwavelength = minvelocity / maxfreq
|
||||
|
||||
|
||||
# Resolution of minimum wavelength
|
||||
resolution = minwavelength / resolvedsteps
|
||||
|
||||
@@ -169,10 +194,10 @@ def dispersion_check(G):
|
||||
|
||||
def get_other_directions(direction):
|
||||
"""Return the two other directions from x, y, z given a single direction
|
||||
|
||||
|
||||
Args:
|
||||
direction (str): Component x, y or z
|
||||
|
||||
|
||||
Returns:
|
||||
(tuple): Two directions from x, y, z
|
||||
"""
|
||||
@@ -184,7 +209,7 @@ def get_other_directions(direction):
|
||||
|
||||
def Ix(x, y, z, Hy, Hz, G):
|
||||
"""Calculates the x-component of current at a grid position.
|
||||
|
||||
|
||||
Args:
|
||||
x, y, z (float): Coordinates of position in grid.
|
||||
Hy, Hz (memory view): numpy array of magnetic field values.
|
||||
@@ -194,14 +219,14 @@ def Ix(x, y, z, Hy, Hz, G):
|
||||
if y == 0 or z == 0:
|
||||
Ix = 0
|
||||
return Ix
|
||||
|
||||
|
||||
else:
|
||||
Ix = G.dy * (Hy[x, y, z - 1] - Hy[x, y, z]) + G.dz * (Hz[x, y, z] - Hz[x, y - 1, z])
|
||||
return Ix
|
||||
|
||||
def Iy(x, y, z, Hx, Hz, G):
|
||||
"""Calculates the y-component of current at a grid position.
|
||||
|
||||
|
||||
Args:
|
||||
x, y, z (float): Coordinates of position in grid.
|
||||
Hx, Hz (memory view): numpy array of magnetic field values.
|
||||
@@ -211,14 +236,14 @@ def Iy(x, y, z, Hx, Hz, G):
|
||||
if x == 0 or z == 0:
|
||||
Iy = 0
|
||||
return Iy
|
||||
|
||||
|
||||
else:
|
||||
Iy = G.dx * (Hx[x, y, z] - Hx[x, y, z - 1]) + G.dz * (Hz[x - 1, y, z] - Hz[x, y, z])
|
||||
return Iy
|
||||
|
||||
def Iz(x, y, z, Hx, Hy, G):
|
||||
"""Calculates the z-component of current at a grid position.
|
||||
|
||||
|
||||
Args:
|
||||
x, y, z (float): Coordinates of position in grid.
|
||||
Hx, Hy (memory view): numpy array of magnetic field values.
|
||||
@@ -228,7 +253,7 @@ def Iz(x, y, z, Hx, Hy, G):
|
||||
if x == 0 or y == 0:
|
||||
Iz = 0
|
||||
return Iz
|
||||
|
||||
|
||||
else:
|
||||
Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z])
|
||||
return Iz
|
||||
|
||||
@@ -17,26 +17,25 @@
|
||||
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
import gprMax
|
||||
from ._version import __version__
|
||||
from gprMax.constants import floattype
|
||||
from gprMax.grid import Ix, Iy, Iz
|
||||
|
||||
|
||||
def prepare_hdf5(outputfile, G):
|
||||
"""Prepares an output file in HDF5 format for writing.
|
||||
|
||||
|
||||
Args:
|
||||
outputfile (str): Name of the output file.
|
||||
G (class): Grid class instance - holds essential parameters describing the model.
|
||||
|
||||
|
||||
Returns:
|
||||
f (file object): File object for the file to be written to.
|
||||
f (file object): File object for the file to be written to.
|
||||
"""
|
||||
|
||||
f = h5py.File(outputfile, 'w')
|
||||
f.attrs['gprMax'] = gprMax.__version__
|
||||
f.attrs['gprMax'] = __version__
|
||||
f.attrs['Title'] = G.title
|
||||
f.attrs['Iterations'] = G.iterations
|
||||
f.attrs['nx, ny, nz'] = (G.nx, G.ny, G.nz)
|
||||
@@ -47,7 +46,7 @@ def prepare_hdf5(outputfile, G):
|
||||
f.attrs['nrx'] = len(G.rxs)
|
||||
f.attrs['srcsteps'] = (G.srcstepx, G.srcstepy, G.srcstepz)
|
||||
f.attrs['rxsteps'] = (G.rxstepx, G.rxstepy, G.rxstepz)
|
||||
|
||||
|
||||
|
||||
# Create group for sources (except transmission lines); add type and positional data attributes
|
||||
srclist = G.voltagesources + G.hertziandipoles + G.magneticdipoles
|
||||
@@ -55,7 +54,7 @@ def prepare_hdf5(outputfile, G):
|
||||
grp = f.create_group('/srcs/src' + str(srcindex + 1))
|
||||
grp.attrs['Type'] = type(src).__name__
|
||||
grp.attrs['Position'] = (src.xcoord * G.dx, src.ycoord * G.dy, src.zcoord * G.dz)
|
||||
|
||||
|
||||
# Create group for transmission lines; add positional data, line resistance and line discretisation attributes; initialise arrays for line voltages and currents
|
||||
if G.transmissionlines:
|
||||
for tlindex, tl in enumerate(G.transmissionlines):
|
||||
@@ -68,7 +67,7 @@ def prepare_hdf5(outputfile, G):
|
||||
grp['Iinc'] = tl.Iinc
|
||||
grp.create_dataset('Vtotal', (G.iterations, ), dtype=floattype)
|
||||
grp.create_dataset('Itotal', (G.iterations, ), dtype=floattype)
|
||||
|
||||
|
||||
# Create group and add positional data and initialise field component arrays for receivers
|
||||
for rxindex, rx in enumerate(G.rxs):
|
||||
grp = f.create_group('/rxs/rx' + str(rxindex + 1))
|
||||
@@ -83,7 +82,7 @@ def prepare_hdf5(outputfile, G):
|
||||
|
||||
def write_hdf5(f, timestep, Ex, Ey, Ez, Hx, Hy, Hz, G):
|
||||
"""Writes field component values to an output file in HDF5 format.
|
||||
|
||||
|
||||
Args:
|
||||
f (file object): File object for the file to be written to.
|
||||
timestep (int): Current iteration number.
|
||||
@@ -116,4 +115,3 @@ def write_hdf5(f, timestep, Ex, Ey, Ez, Hx, Hy, Hz, G):
|
||||
for tlindex, tl in enumerate(G.transmissionlines):
|
||||
f['/tls/tl' + str(tlindex + 1) + '/Vtotal'][timestep] = tl.voltage[tl.antpos]
|
||||
f['/tls/tl' + str(tlindex + 1) + '/Itotal'][timestep] = tl.current[tl.antpos]
|
||||
|
||||
|
||||
327
gprMax/xdmf.py
普通文件
327
gprMax/xdmf.py
普通文件
@@ -0,0 +1,327 @@
|
||||
import numpy as np
|
||||
import h5py
|
||||
from lxml import etree
|
||||
from gprMax.grid import Grid
|
||||
import copy
|
||||
|
||||
|
||||
class Edges:
|
||||
|
||||
def __init__(self, grid):
|
||||
|
||||
"""
|
||||
Class to define some connectivity of for an n x l x m
|
||||
grid
|
||||
"""
|
||||
self.total_edges = grid.n_edges()
|
||||
self.edges = np.zeros((self.total_edges, 2), np.float32)
|
||||
self.edge_count = 0
|
||||
self.grid = grid
|
||||
"""
|
||||
Adds the the edge specified by in_node and the i,j,k position of the outnode
|
||||
"""
|
||||
def add_edge(self, in_label, i, j, k):
|
||||
|
||||
out_label = self.grid.get(i, j, k)
|
||||
self.edges[self.edge_count] = np.array([in_label, out_label])
|
||||
self.edge_count += 1
|
||||
|
||||
|
||||
class Coordinates:
|
||||
|
||||
def __init__(self, grid):
|
||||
self.total_coordinates = grid.nx * grid.ny * grid.nz
|
||||
self.coordinate_count = 0
|
||||
self.coordinates = np.zeros((self.total_coordinates, 3), np.float32)
|
||||
|
||||
def add_coordinate(self, x, y, z):
|
||||
self.coordinates[self.coordinate_count] = np.array([x, y, z])
|
||||
self.coordinate_count += 1
|
||||
|
||||
|
||||
def hexCellPicker(grid, i, j, k):
|
||||
"""
|
||||
This is the ordering of nodes in the hexahedron cell.
|
||||
|
||||
7 --------- 6
|
||||
/ /|
|
||||
4 --------- 5 2
|
||||
| 3 | /
|
||||
| / |/
|
||||
0 --------- 1
|
||||
|
||||
0 1 2 3 4 5 6 7
|
||||
"""
|
||||
|
||||
cell = [
|
||||
grid[i][j][k],
|
||||
# 1
|
||||
grid[i + 1][j][k],
|
||||
# 2
|
||||
grid[i + 1][j + 1][k],
|
||||
# 3
|
||||
grid[i][j + 1][k],
|
||||
# 4
|
||||
grid[i][j][k + 1],
|
||||
# 5
|
||||
grid[i + 1][j][k + 1],
|
||||
# 6
|
||||
grid[i + 1][j + 1][k + 1],
|
||||
# 7
|
||||
grid[i][j + 1][k + 1]
|
||||
]
|
||||
|
||||
return cell
|
||||
|
||||
|
||||
class Solids:
|
||||
|
||||
def __init__(self, fdtd_grid):
|
||||
self.count = 0
|
||||
self.fdtd_grid = fdtd_grid
|
||||
self.total_solids = fdtd_grid.n_cells()
|
||||
self.solids = np.zeros((self.total_solids), np.float32)
|
||||
|
||||
def add_solid(self, i, j, k):
|
||||
|
||||
self.solids[self.count] = self.fdtd_grid.solid[i][j][k]
|
||||
self.count += 1
|
||||
|
||||
|
||||
class SolidLabels():
|
||||
|
||||
def __init__(self, label_grid):
|
||||
self.count = 0
|
||||
self.label_grid = label_grid
|
||||
self.total_solids = label_grid.n_cells()
|
||||
self.solid_labels = np.zeros((self.total_solids, 8), np.float32)
|
||||
|
||||
def add(self, i, j, k):
|
||||
|
||||
solid_labels = hexCellPicker(self.label_grid.grid, i, j, k)
|
||||
self.solid_labels[self.count] = solid_labels
|
||||
self.count += 1
|
||||
|
||||
|
||||
class Materials:
|
||||
|
||||
def __init__(self, fdtd_grid):
|
||||
self.fdtd_grid = fdtd_grid
|
||||
self.n_edges = fdtd_grid.n_edges()
|
||||
self.materials = np.zeros((self.n_edges), np.float32)
|
||||
self.material_count = 0
|
||||
|
||||
# direction x->0 y->1 z->2
|
||||
def add_material(self, i, j, k, direction):
|
||||
|
||||
material = self.fdtd_grid.ID[direction, i, j, k]
|
||||
self.materials[self.material_count] = material
|
||||
|
||||
self.material_count += 1
|
||||
|
||||
|
||||
def process_grid(fdtd_grid):
|
||||
|
||||
# Dimensions of the problem domain.
|
||||
nx = fdtd_grid.nx
|
||||
ny = fdtd_grid.ny
|
||||
nz = fdtd_grid.nz
|
||||
|
||||
# label each node in the space
|
||||
labels = np.arange(nx * ny * nz).reshape(nx, ny, nz)
|
||||
|
||||
label_grid = Grid(labels)
|
||||
|
||||
# Edges define the connectivity of the grid.
|
||||
edges = Edges(label_grid)
|
||||
|
||||
# Material for each edge
|
||||
edge_materials = Materials(fdtd_grid)
|
||||
|
||||
# Define coordinates for each node
|
||||
coordinates = Coordinates(fdtd_grid)
|
||||
|
||||
# Material for each solid
|
||||
solids = Solids(fdtd_grid)
|
||||
|
||||
# Connectivity for hexhahedron grid
|
||||
solid_labels = SolidLabels(label_grid)
|
||||
|
||||
i_max = nx - 1
|
||||
j_max = ny - 1
|
||||
k_max = nz - 1
|
||||
|
||||
for i, ix in enumerate(labels):
|
||||
for j, jx in enumerate(ix):
|
||||
for k, kx in enumerate(jx):
|
||||
|
||||
label = labels[i][j][k]
|
||||
|
||||
# Each vertex can have varying numbers of edges
|
||||
|
||||
# Type 1 vertex
|
||||
if i < i_max and j < j_max and k < k_max:
|
||||
edges.add_edge(label, i + 1, j, k)
|
||||
edges.add_edge(label, i, j + 1, k)
|
||||
edges.add_edge(label, i, j, k + 1)
|
||||
|
||||
edge_materials.add_material(i, j, k, 0)
|
||||
edge_materials.add_material(i, j, k, 1)
|
||||
edge_materials.add_material(i, j, k, 2)
|
||||
|
||||
# Only this node can support a cell
|
||||
solids.add_solid(i, j, k)
|
||||
solid_labels.add(i, j, k)
|
||||
|
||||
# Type 2 vertex
|
||||
elif i < i_max and j == j_max and k == k_max:
|
||||
edges.add_edge(label, i + 1, j, k)
|
||||
edge_materials.add_material(i, j, k, 0)
|
||||
|
||||
# Type 7 vertex
|
||||
elif i < i_max and j == j_max and k < k_max:
|
||||
edges.add_edge(label, i + 1, j, k)
|
||||
edges.add_edge(label, i, j, k + 1)
|
||||
edge_materials.add_material(i, j, k, 0)
|
||||
edge_materials.add_material(i, j, k, 2)
|
||||
|
||||
# Type 6 vertex
|
||||
elif i == i_max and j == j_max and k < k_max:
|
||||
edges.add_edge(label, i, j, k + 1)
|
||||
edge_materials.add_material(i, j, k, 2)
|
||||
|
||||
# Type 5 vertex
|
||||
elif i == i_max and j < j_max and k < k_max:
|
||||
edges.add_edge(label, i, j, k + 1)
|
||||
edges.add_edge(label, i, j + 1, k)
|
||||
edge_materials.add_material(i, j, k, 2)
|
||||
edge_materials.add_material(i, j, k, 1)
|
||||
|
||||
# Type 4 vertex
|
||||
elif i == i_max and j < j_max and k == k_max:
|
||||
edges.add_edge(label, i, j + 1, k)
|
||||
edge_materials.add_material(i, j, k, 1)
|
||||
|
||||
# Type 8 vertex
|
||||
elif i < i_max and j < j_max and k == k_max:
|
||||
edges.add_edge(label, i, j + 1, k)
|
||||
edges.add_edge(label, i + 1, j, k)
|
||||
edge_materials.add_material(i, j, k, 1)
|
||||
edge_materials.add_material(i, j, k, 0)
|
||||
|
||||
# Type 3 vertex
|
||||
# Has no new connectivity
|
||||
elif i == i_max and j == j_max and k == k_max:
|
||||
pass
|
||||
else:
|
||||
print('oh no')
|
||||
|
||||
# Add the coordinates
|
||||
coordinates.add_coordinate(i, j, k)
|
||||
|
||||
return {
|
||||
'coordinates': coordinates,
|
||||
'solids': solids,
|
||||
'solid_labels': solid_labels,
|
||||
'edges': edges,
|
||||
'edge_materials': edge_materials,
|
||||
}
|
||||
|
||||
|
||||
def write_output_file(filename, grid):
|
||||
|
||||
data = process_grid(grid)
|
||||
data['filename'] = filename
|
||||
data['xml_doc'] = create_xdmf_markup(data)
|
||||
|
||||
write_H5file(data)
|
||||
write_xml_doc(data)
|
||||
|
||||
|
||||
def write_xml_doc(options):
|
||||
#write xml to file
|
||||
with open(options['filename'] + '.xdmf', 'wb') as xdmf_f:
|
||||
xdmf_f.write(options['xml_doc'])
|
||||
|
||||
|
||||
def write_H5file(options):
|
||||
|
||||
f = h5py.File(options['filename'] + '.h5', "w")
|
||||
coords = f.create_group("mesh")
|
||||
coords.create_dataset('coordinates', data=options['coordinates'].coordinates)
|
||||
coords.create_dataset('connectivity', data=options['edges'].edges)
|
||||
coords.create_dataset('solid_connectivity', data=options['solid_labels'].solid_labels)
|
||||
data = f.create_group("data")
|
||||
data.create_dataset('materials', data=options['edge_materials'].materials)
|
||||
data.create_dataset('solids', data=options['solids'].solids)
|
||||
|
||||
|
||||
def create_xdmf_markup(options):
|
||||
|
||||
# Write the XDMF markup for edge style grid
|
||||
xdmf_el = etree.Element("Xdmf", Version="2.0")
|
||||
|
||||
domain_el = etree.Element("Domain")
|
||||
xdmf_el.append(domain_el)
|
||||
|
||||
grid_el = etree.Element("Grid", Name="Edges", GridType="Uniform")
|
||||
domain_el.append(grid_el)
|
||||
|
||||
# Create the grid node
|
||||
topology_el = etree.Element("Topology", TopologyType="Polyline", NumberOfElements=str(options['edges'].total_edges))
|
||||
grid_el.append(topology_el)
|
||||
|
||||
topology_dimensions = "{} 2".format(options['edges'].total_edges)
|
||||
top_data_el = etree.Element("DataItem", Dimensions=topology_dimensions, NumberType="Float", Precision="8", Format="HDF")
|
||||
top_data_el.text = "{}:/mesh/connectivity".format(options['filename'] + '.h5')
|
||||
topology_el.append(top_data_el)
|
||||
|
||||
# Create the Geometry node
|
||||
geometry_el = etree.Element("Geometry", GeometryType="XYZ")
|
||||
grid_el.append(geometry_el)
|
||||
|
||||
# Create the origin coordinates
|
||||
coordinates_dimensions = "{} 3".format(options['coordinates'].total_coordinates)
|
||||
origin_el = etree.Element("DataItem", Dimensions=coordinates_dimensions, NumberType="Float", Precision="8", Format="HDF")
|
||||
origin_el.text = "{}:/mesh/coordinates".format(options['filename'] + '.h5')
|
||||
geometry_el.append(origin_el)
|
||||
|
||||
# Create the materials attribute
|
||||
attr_el = etree.Element("Attribute", Center="Cell", Name="Edge_Materials")
|
||||
grid_el.append(attr_el)
|
||||
|
||||
materials_dimensions = "{} 1".format(options['edge_materials'].material_count)
|
||||
materials_el = etree.Element("DataItem", Dimensions=materials_dimensions, NumberType="Float", Precision="8", Format="HDF")
|
||||
materials_el.text = "{}:/data/materials".format(options['filename'] + '.h5')
|
||||
attr_el.append(materials_el)
|
||||
|
||||
v_grid_el = etree.Element("Grid", Name="Voxel", GridType="Uniform")
|
||||
domain_el.append(v_grid_el)
|
||||
|
||||
n_solids = str(options['solids'].solids.size)
|
||||
v_topology_el = etree.Element("Topology", TopologyType="Hexahedron", NumberOfElements=str(options['solids'].solids.size))
|
||||
v_grid_el.append(v_topology_el)
|
||||
|
||||
solid_label_d = "{} {}".format(n_solids, 8)
|
||||
solid_labels_el = etree.Element("DataItem", Dimensions=solid_label_d, Format="HDF")
|
||||
solid_labels_el.text = "{}:/mesh/solid_connectivity".format(options['filename'] + '.h5')
|
||||
v_topology_el.append(solid_labels_el)
|
||||
|
||||
# Same geometry as edges
|
||||
v_grid_el.append(copy.deepcopy(geometry_el))
|
||||
|
||||
v_attr = etree.Element("Attribute", Name="Voxel_Materials", Center="Cell")
|
||||
v_grid_el.append(v_attr)
|
||||
|
||||
d4 = etree.Element("DataItem", Format="HDF", NumberType="Float", Precision="4", Dimensions=str(options['solids'].solids.size))
|
||||
d4.text = "{}:/data/solids".format(options['filename'] + '.h5')
|
||||
v_attr.append(d4)
|
||||
|
||||
# Define a doctype - useful for parsers
|
||||
doc_type = '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>'
|
||||
|
||||
# Serialize elements
|
||||
xml_doc = etree.tostring(xdmf_el, xml_declaration=True,
|
||||
encoding="utf-8", doctype=doc_type, pretty_print=True)
|
||||
|
||||
return xml_doc
|
||||
37
setup.py
37
setup.py
@@ -27,28 +27,30 @@ try:
|
||||
except ImportError:
|
||||
raise ImportError('gprMax requires the NumPy package.')
|
||||
|
||||
import glob, os, shutil, sys
|
||||
import glob
|
||||
import os
|
||||
import shutil
|
||||
import sys
|
||||
import re
|
||||
|
||||
# Importing _version__.py before building can cause issues.
|
||||
with open('gprMax/_version.py', 'r') as fd:
|
||||
version = re.search(r'^__version__\s*=\s*[\'"]([^\'"]*)[\'"]',
|
||||
fd.read(), re.MULTILINE).group(1)
|
||||
|
||||
"""
|
||||
Parse package name from init file. Importing __init__.py / gprMax will
|
||||
break as gprMax depends on compiled .pyx files.
|
||||
"""
|
||||
with open('gprMax/__init__.py', 'r') as fd:
|
||||
packagename = re.search(r'^__name__\s*=\s*[\'"]([^\'"]*)[\'"]',
|
||||
fd.read(), re.MULTILINE).group(1)
|
||||
|
||||
# Python version
|
||||
if sys.version_info[:2] < (3, 4):
|
||||
print('gprMax requires Python 3.4 or newer')
|
||||
sys.exit(-1)
|
||||
|
||||
import gprMax
|
||||
|
||||
# Package name
|
||||
packagename = gprMax.__name__
|
||||
|
||||
# Package version
|
||||
version = gprMax.__version__
|
||||
|
||||
# Process 'install' command line argument
|
||||
if 'install' in sys.argv:
|
||||
print("'install' is not required for this package, running 'build_ext --inplace' instead.")
|
||||
sys.argv.remove('install')
|
||||
sys.argv.append('build_ext')
|
||||
sys.argv.append('--inplace')
|
||||
|
||||
# Process 'build' command line argument
|
||||
if 'build' in sys.argv:
|
||||
print("'build' is not required for this package, running 'build_ext --inplace' instead.")
|
||||
@@ -118,7 +120,7 @@ elif sys.platform == 'darwin':
|
||||
compile_args = ['-O3', '-w', '-fstrict-aliasing', '-fno-common', '-fopenmp']
|
||||
linker_args = ['-fopenmp']
|
||||
extra_objects = []
|
||||
|
||||
|
||||
# If user is cwarren do static linking as this is for binaries uploaded to GitHub
|
||||
# if os.getlogin() == 'cwarren':
|
||||
# linker_args = ['-static-libgcc']
|
||||
@@ -176,4 +178,5 @@ setup(name=packagename,
|
||||
'Topic :: Scientific/Engineering'
|
||||
],
|
||||
ext_modules=extensions,
|
||||
packages=['gprMax'],
|
||||
include_dirs=[np.get_include()])
|
||||
|
||||
@@ -0,0 +1,58 @@
|
||||
## Tissue properties
|
||||
## Non-dispersive properties (900MHz) from http://niremf.ifac.cnr.it/tissprop/htmlclie/htmlclie.php
|
||||
#material: 1 0 1 0 Air
|
||||
#material: 44.77 0.6961 1 0 Aorta
|
||||
#material: 18.94 0.3831 1 0 Bladder
|
||||
#material: 61.36 1.538 1 0 Blood
|
||||
#material: 44.77 0.6961 1 0 BloodVessel
|
||||
#material: 68.9 1.636 1 0 BodyFluid
|
||||
#material: 20.79 0.34 1 0 BoneCancellous
|
||||
#material: 12.45 0.1433 1 0 BoneCortical
|
||||
#material: 5.504 0.04021 1 0 BoneMarrow
|
||||
#material: 52.73 0.9423 1 0 BrainGreyMatter
|
||||
#material: 38.89 0.5908 1 0 BrainWhiteMatter
|
||||
#material: 5.424 0.04896 1 0 BreastFat
|
||||
#material: 42.65 0.7824 1 0 Cartilage
|
||||
#material: 49.44 1.263 1 0 Cerebellum
|
||||
#material: 68.64 2.413 1 0 CerebroSpinalFluid
|
||||
#material: 49.78 0.9591 1 0 Cervix
|
||||
#material: 57.94 1.08 1 0 Colon
|
||||
#material: 55.23 1.394 1 0 Cornea
|
||||
#material: 65.06 1.187 1 0 Duodenum
|
||||
#material: 44.43 0.9611 1 0 Dura
|
||||
#material: 55.27 1.167 1 0 EyeSclera
|
||||
#material: 5.462 0.05104 1 0 Fat
|
||||
#material: 59.14 1.257 1 0 GallBladder
|
||||
#material: 70.19 1.838 1 0 GallBladderBile
|
||||
#material: 59.68 1.038 1 0 Gland
|
||||
#material: 59.89 1.23 1 0 Heart
|
||||
#material: 58.67 1.392 1 0 Kidney
|
||||
#material: 46.57 0.7934 1 0 Lens
|
||||
#material: 46.83 0.855 1 0 Liver
|
||||
#material: 51.42 0.858 1 0 LungDeflated
|
||||
#material: 22 0.4567 1 0 LungInflated
|
||||
#material: 59.68 1.038 1 0 Lymph
|
||||
#material: 46.08 0.8447 1 0 MucousMembrane
|
||||
#material: 55.03 0.9429 1 0 Muscle
|
||||
#material: 12.45 0.1433 1 0 Nail
|
||||
#material: 32.53 0.5737 1 0 Nerve
|
||||
#material: 65.06 1.187 1 0 Oesophagus
|
||||
#material: 50.47 1.29 1 0 Ovary
|
||||
#material: 59.68 1.038 1 0 Pancreas
|
||||
#material: 60.55 1.21 1 0 Prostate
|
||||
#material: 55.27 1.167 1 0 Retina
|
||||
#material: 41.41 0.8667 1 0 SkinDry
|
||||
#material: 46.08 0.8447 1 0 SkinWet
|
||||
#material: 59.49 2.165 1 0 SmallIntestine
|
||||
#material: 32.53 0.5737 1 0 SpinalCord
|
||||
#material: 57.18 1.273 1 0 Spleen
|
||||
#material: 65.06 1.187 1 0 Stomach
|
||||
#material: 45.83 0.7184 1 0 Tendon
|
||||
#material: 60.55 1.21 1 0 Testis
|
||||
#material: 59.68 1.038 1 0 Thymus
|
||||
#material: 59.68 1.038 1 0 Thyroid
|
||||
#material: 55.27 0.9363 1 0 Tongue
|
||||
#material: 12.45 0.1433 1 0 Tooth
|
||||
#material: 42.01 0.7711 1 0 Trachea
|
||||
#material: 61.12 1.27 1 0 Uterus
|
||||
#material: 68.9 1.636 1 0 VitreousHumor
|
||||
@@ -0,0 +1,124 @@
|
||||
## Tissue properties
|
||||
## 3-Pole Debye model (1MHz-100GHz) from http://dx.doi.org/10.1109/LMWC.2011.2180371
|
||||
## Non-dispersive properties (900MHz) from http://niremf.ifac.cnr.it/tissprop/htmlclie/htmlclie.php
|
||||
|
||||
#material: 1 0 1 0 Air
|
||||
#material: 44.77 0.6961 1 0 Aorta
|
||||
|
||||
#material: 3.4969 0.2261 1 0 Bladder
|
||||
#add_dispersion_debye: 3 370.48 9.1145e-8 21.866 4.6778e-10 15.38 0.94981e-11 Bladder
|
||||
|
||||
#material: 7.5425 0.73265 1 0 Blood
|
||||
#add_dispersion_debye: 3 3529.4 8.8589e-8 110.21 5.4229e-9 53.637 0.89892 Blood
|
||||
|
||||
#material: 44.77 0.6961 1 0 BloodVessel
|
||||
#material: 68.9 1.636 1 0 BodyFluid
|
||||
|
||||
## Avg. Bone used for BoneCancellous & BoneMarrow
|
||||
#material: 3.8518 0.095687 1 0 BoneCancellous
|
||||
#add_dispersion_debye: 3 125.82 3.2802e-8 18.125 2.4901e-9 11.748 1.379e-11 BoneCancellous
|
||||
|
||||
#material: 3.4498 0.022411 1 0 BoneCortical
|
||||
#add_dispersion_debye: 3 118.19 4.6747e-8 9.5114 2.2822e-9 8.807 1.3731e-11 BoneCortical
|
||||
|
||||
## Avg. Bone used for BoneCancellous & BoneMarrow
|
||||
#material: 3.8518 0.095687 1 0 BoneMarrow
|
||||
#add_dispersion_debye: 3 125.82 3.2802e-8 18.125 2.4901e-9 11.748 1.379e-11 BoneMarrow
|
||||
|
||||
## Avg. Brain used for BrainGreyMatter & BrainWhiteMatter
|
||||
#material: 6.9992 0.17911 1 0 BrainGreyMatter
|
||||
#add_dispersion_debye: 3 787.41 2.2589e-8 142.26 2.5012e-9 40.264 0.91939e-11 BrainGreyMatter
|
||||
|
||||
#material: 6.9992 0.17911 1 0 BrainWhiteMatter
|
||||
#add_dispersion_debye: 3 787.41 2.2589e-8 142.26 2.5012e-9 40.264 0.91939e-11 BrainWhiteMatter
|
||||
|
||||
#material: 5.424 0.04896 1 0 BreastFat
|
||||
#material: 42.65 0.7824 1 0 Cartilage
|
||||
|
||||
#material: 6.9992 0.17911 1 0 Cerebellum
|
||||
#add_dispersion_debye: 3 787.41 2.2589e-8 142.26 2.5012e-9 40.264 0.91939e-11 Cerebellum
|
||||
|
||||
#material: 5.7772 2.0002 1 0 CerebroSpinalFluid
|
||||
#add_dispersion_debye: 3 40.84 0.15293e-8 27.951 0.014519e-9 34.334 0.50631e-11 CerebroSpinalFluid
|
||||
|
||||
#material: 49.78 0.9591 1 0 Cervix
|
||||
#material: 57.94 1.08 1 0 Colon
|
||||
#material: 55.23 1.394 1 0 Cornea
|
||||
|
||||
#material: 7.8582 0.55004 1 0 Duodenum
|
||||
#add_dispersion_debye: 3 1674.5 6.273e-8 74.699 4.4123e-9 56.797 0.85287e-11 Duodenum
|
||||
|
||||
#material: 44.43 0.9611 1 0 Dura
|
||||
#material: 55.27 1.167 1 0 EyeSclera
|
||||
|
||||
#material: 3.261 0.043999 1 0 Fat
|
||||
#add_dispersion_debye: 3 34.522 1.3971e-8 2.2452 0.04806e-9 6.2928 0.60675e-11 Fat
|
||||
|
||||
#material: 4.7563 0.90026 1 0 GallBladder
|
||||
#add_dispersion_debye: 3 40.681 0.1583e-8 20.066 0.01249e-9 33.876 0.57641e-11 GallBladder
|
||||
|
||||
#material: 70.19 1.838 1 0 GallBladderBile
|
||||
#material: 59.68 1.038 1 0 Gland
|
||||
|
||||
#material: 7.8989 0.27449 1 0 Heart
|
||||
#add_dispersion_debye: 3 1935.3 7.9786e-8 138.97 3.2768e-9 50.762 0.93944e-11 Heart
|
||||
|
||||
#material: 7.8691 0.23004 1 0 Kidney
|
||||
#add_dispersion_debye: 3 2025 6.5466e-8 167.3 3.0341e-9 48.945 0.96655e-11 Kidney
|
||||
|
||||
## Avg. Lens
|
||||
#material: 5.8421 0.27596 1 0 Lens
|
||||
#add_dispersion_debye: 3 630.99 5.9556e-8 62.733 3.6647e-9 30.811 0.89718e-11 Lens
|
||||
|
||||
#material: 6.809 0.14243 1 0 Liver
|
||||
#add_dispersion_debye: 3 1529.8 8.4314e-8 112.04 3.7562e-9 39.593 1.0293e-11 Liver
|
||||
|
||||
## Avg. Lung used for LungDeflated & LungInflated
|
||||
#material: 5.5014 0.21296 1 0 LungDeflated
|
||||
#add_dispersion_debye: 3 899.49 6.9396e-8 62.733 3.6647e-9 30.811 0.89718e-11 LungDeflated
|
||||
|
||||
#material: 5.5014 0.21296 1 0 LungInflated
|
||||
#add_dispersion_debye: 3 899.49 6.9396e-8 62.733 3.6647e-9 30.811 0.89718e-11 LungInflated
|
||||
|
||||
#material: 59.68 1.038 1 0 Lymph
|
||||
#material: 46.08 0.8447 1 0 MucousMembrane
|
||||
|
||||
## Avg. Muscle
|
||||
#material: 7.5933 0.42066 1 0 Muscle
|
||||
#add_dispersion_debye: 3 2646.6 13.044e-8 87.588 6.1808e-9 47.546 0.79767e-11 Muscle
|
||||
|
||||
#material: 12.45 0.1433 1 0 Nail
|
||||
#material: 32.53 0.5737 1 0 Nerve
|
||||
#material: 65.06 1.187 1 0 Oesophagus
|
||||
#material: 50.47 1.29 1 0 Ovary
|
||||
#material: 59.68 1.038 1 0 Pancreas
|
||||
#material: 60.55 1.21 1 0 Prostate
|
||||
#material: 55.27 1.167 1 0 Retina
|
||||
|
||||
## Avg. Skin used from SkinDry & SkinWet
|
||||
#material: 5.8661 0.092415 1 0 SkinDry
|
||||
#add_dispersion_debye: 3 1123 5.4817e-8 144.69 4.0078e-9 37.434 0.88726e-11 SkinDry
|
||||
|
||||
#material: 5.8661 0.092415 1 0 SkinWet
|
||||
#add_dispersion_debye: 3 1123 5.4817e-8 144.69 4.0078e-9 37.434 0.88726e-11 SkinWet
|
||||
|
||||
#material: 7.7811 0.6884 1 0 SmallIntestine
|
||||
#add_dispersion_debye: 3 6590.7 9.123e-8 211.31 4.2003e-9 50.6 0.92249e-11 SmallIntestine
|
||||
|
||||
#material: 32.53 0.5737 1 0 SpinalCord
|
||||
#material: 57.18 1.273 1 0 Spleen
|
||||
|
||||
#material: 7.8582 0.55004 1 0 Stomach
|
||||
#add_dispersion_debye: 3 1674.5 6.273e-8 74.699 4.4123e-9 56.797 0.85287e-11 Stomach
|
||||
|
||||
#material: 45.83 0.7184 1 0 Tendon
|
||||
#material: 60.55 1.21 1 0 Testis
|
||||
#material: 59.68 1.038 1 0 Thymus
|
||||
#material: 59.68 1.038 1 0 Thyroid
|
||||
#material: 55.27 0.9363 1 0 Tongue
|
||||
#material: 12.45 0.1433 1 0 Tooth
|
||||
#material: 42.01 0.7711 1 0 Trachea
|
||||
#material: 61.12 1.27 1 0 Uterus
|
||||
|
||||
#material: 4.0041 1.5001 1 0 VitreousHumor
|
||||
#add_dispersion_debye: 3 21.766 12.46e-8 1.364 15.742e-9 65.024 0.72371e-11 VitreousHumor
|
||||
@@ -0,0 +1,23 @@
|
||||
import argparse, os
|
||||
import h5py
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(description='Writes a HDF5 file of AustinMan or AustinWoman head only.', usage='python head_only_hdf5 filename')
|
||||
parser.add_argument('filename', help='name and path to (HDF5) file containing AustinMan or AustinWoman model')
|
||||
args = parser.parse_args()
|
||||
|
||||
# Read full body HDF5 file
|
||||
f = h5py.File(args.filename, 'r')
|
||||
dx_dy_dz = f.attrs['dx, dy, dz']
|
||||
data = f['/data'][:,:,:]
|
||||
|
||||
# Define head as last 1/8 of total body height
|
||||
nzhead = 7 * int(data.shape[2] / 8)
|
||||
|
||||
print('Dimensions of head model: {:g} x {:g} x {:g} cells'.format(data.shape[0], data.shape[1], data.shape[2] - nzhead))
|
||||
|
||||
# Write HDF5 file
|
||||
headfile = os.path.splitext(args.filename)[0] + '_head.h5'
|
||||
f = h5py.File(headfile, 'w')
|
||||
f.attrs['dx, dy, dz'] = dx_dy_dz
|
||||
f['/data'] = data[:,:,nzhead:data.shape[2]]
|
||||
在新工单中引用
屏蔽一个用户