你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-08 07:24:19 +08:00
Create new geometry views/objects using grid views
这个提交包含在：
nmannall
@@ -0,0 +1,181 @@
|
||||
from io import TextIOWrapper
|
||||
from pathlib import Path
|
||||
from typing import Generic
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
from mpi4py import MPI
|
||||
from tqdm import tqdm
|
||||
|
||||
from gprMax import config
|
||||
from gprMax._version import __version__
|
||||
from gprMax.grid.fdtd_grid import FDTDGrid
|
||||
from gprMax.materials import Material
|
||||
from gprMax.output_controllers.grid_view import GridViewType, MPIGridView
|
||||
|
||||
|
||||
class GeometryObjects(Generic[GridViewType]):
|
||||
"""Geometry objects to be written to file."""
|
||||
|
||||
def __init__(self, grid_view: GridViewType, filename: str):
|
||||
"""
|
||||
Args:
|
||||
xs, xf, ys, yf, zs, zf: ints for extent of the volume in cells.
|
||||
filename: string for filename.
|
||||
"""
|
||||
self.grid_view = grid_view
|
||||
|
||||
# Set filenames
|
||||
parts = config.sim_config.input_file_path.with_suffix("").parts
|
||||
self.filename = Path(*parts[:-1], filename)
|
||||
self.filename_hdf5 = self.filename.with_suffix(".h5")
|
||||
self.filename_materials = Path(self.filename, f"{self.filename}_materials")
|
||||
self.filename_materials = self.filename_materials.with_suffix(".txt")
|
||||
|
||||
# Sizes of arrays to write necessary to update progress bar
|
||||
self.solidsize = np.prod(self.grid_view.size + 1) * np.dtype(np.uint32).itemsize
|
||||
self.rigidsize = 18 * np.prod(self.grid_view.size + 1) * np.dtype(np.int8).itemsize
|
||||
self.IDsize = 6 * np.prod(self.grid_view.size + 1) * np.dtype(np.uint32).itemsize
|
||||
self.datawritesize = self.solidsize + self.rigidsize + self.IDsize
|
||||
|
||||
@property
|
||||
def grid(self) -> FDTDGrid:
|
||||
return self.grid_view.grid
|
||||
|
||||
def write_metadata(self, file_handler: h5py.File, title: str):
|
||||
file_handler.attrs["gprMax"] = __version__
|
||||
file_handler.attrs["Title"] = title
|
||||
file_handler.attrs["dx_dy_dz"] = (self.grid.dx, self.grid.dy, self.grid.dz)
|
||||
|
||||
def output_material(self, material: Material, file: TextIOWrapper):
|
||||
file.write(
|
||||
f"#material: {material.er:g} {material.se:g} "
|
||||
f"{material.mr:g} {material.sm:g} {material.ID}\n"
|
||||
)
|
||||
if hasattr(material, "poles"):
|
||||
if "debye" in material.type:
|
||||
dispersionstr = "#add_dispersion_debye: " f"{material.poles:g} "
|
||||
for pole in range(material.poles):
|
||||
dispersionstr += f"{material.deltaer[pole]:g} " f"{material.tau[pole]:g} "
|
||||
elif "lorenz" in material.type:
|
||||
dispersionstr = f"#add_dispersion_lorenz: " f"{material.poles:g} "
|
||||
for pole in range(material.poles):
|
||||
dispersionstr += (
|
||||
f"{material.deltaer[pole]:g} "
|
||||
f"{material.tau[pole]:g} "
|
||||
f"{material.alpha[pole]:g} "
|
||||
)
|
||||
elif "drude" in material.type:
|
||||
dispersionstr = f"#add_dispersion_drude: " f"{material.poles:g} "
|
||||
for pole in range(material.poles):
|
||||
dispersionstr += f"{material.tau[pole]:g} " f"{material.alpha[pole]:g} "
|
||||
dispersionstr += material.ID
|
||||
file.write(dispersionstr + "\n")
|
||||
|
||||
def write_hdf5(self, title: str, G: FDTDGrid, pbar: tqdm):
|
||||
"""Writes a geometry objects file in HDF5 format.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
pbar: Progress bar class instance.
|
||||
"""
|
||||
|
||||
self.grid_view.initialise_materials()
|
||||
|
||||
ID = self.grid_view.get_ID()
|
||||
data = self.grid_view.get_solid().astype(np.int16)
|
||||
rigidE = self.grid_view.get_rigidE()
|
||||
rigidH = self.grid_view.get_rigidH()
|
||||
|
||||
ID = self.grid_view.map_to_view_materials(ID)
|
||||
data = self.grid_view.map_to_view_materials(data)
|
||||
|
||||
with h5py.File(self.filename_hdf5, "w") as fdata:
|
||||
self.write_metadata(fdata, title)
|
||||
|
||||
fdata["/data"] = data
|
||||
pbar.update(self.solidsize)
|
||||
|
||||
fdata["/rigidE"] = rigidE
|
||||
fdata["/rigidH"] = rigidH
|
||||
pbar.update(self.rigidsize)
|
||||
|
||||
fdata["/ID"] = ID
|
||||
pbar.update(self.IDsize)
|
||||
|
||||
# Write materials list to a text file
|
||||
with open(self.filename_materials, "w") as fmaterials:
|
||||
for material in self.grid_view.materials:
|
||||
self.output_material(material, fmaterials)
|
||||
|
||||
|
||||
class MPIGeometryObjects(GeometryObjects[MPIGridView]):
|
||||
def __init__(
|
||||
self,
|
||||
grid_view: MPIGridView,
|
||||
filename: str,
|
||||
comm: MPI.Comm = None,
|
||||
):
|
||||
super().__init__(grid_view, filename)
|
||||
|
||||
self.comm = comm
|
||||
|
||||
def write_hdf5(self, title: str, pbar: tqdm):
|
||||
"""Writes a geometry objects file in HDF5 format.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
pbar: Progress bar class instance.
|
||||
"""
|
||||
|
||||
self.grid_view.initialise_materials(self.comm)
|
||||
|
||||
ID = self.grid_view.get_ID()
|
||||
data = self.grid_view.get_solid().astype(np.int16)
|
||||
rigidE = self.grid_view.get_rigidE()
|
||||
rigidH = self.grid_view.get_rigidH()
|
||||
|
||||
ID = self.grid_view.map_to_view_materials(ID)
|
||||
data = self.grid_view.map_to_view_materials(data)
|
||||
|
||||
with h5py.File(self.filename_hdf5, "w", driver="mpio", comm=self.comm) as fdata:
|
||||
self.write_metadata(fdata, title)
|
||||
|
||||
dset_slice = (
|
||||
self.grid_view.get_slice(0),
|
||||
self.grid_view.get_slice(1),
|
||||
self.grid_view.get_slice(2),
|
||||
)
|
||||
|
||||
dset = fdata.create_dataset("/data", self.grid_view.global_size, dtype=data.dtype)
|
||||
dset[dset_slice] = data
|
||||
pbar.update(self.solidsize)
|
||||
|
||||
rigid_E_dataset = fdata.create_dataset(
|
||||
"/rigidE", (12, *self.grid_view.global_size), dtype=rigidE.dtype
|
||||
)
|
||||
rigid_E_dataset[:, dset_slice[0], dset_slice[1], dset_slice[2]] = rigidE
|
||||
|
||||
rigid_H_dataset = fdata.create_dataset(
|
||||
"/rigidH", (6, *self.grid_view.global_size), dtype=rigidH.dtype
|
||||
)
|
||||
rigid_H_dataset[:, dset_slice[0], dset_slice[1], dset_slice[2]] = rigidH
|
||||
pbar.update(self.rigidsize)
|
||||
|
||||
dset_slice = (
|
||||
self.grid_view.get_slice(0, upper_bound_exclusive=False),
|
||||
self.grid_view.get_slice(1, upper_bound_exclusive=False),
|
||||
self.grid_view.get_slice(2, upper_bound_exclusive=False),
|
||||
)
|
||||
|
||||
dset = fdata.create_dataset(
|
||||
"/ID", (6, *(self.grid_view.global_size + 1)), dtype=ID.dtype
|
||||
)
|
||||
dset[:, dset_slice[0], dset_slice[1], dset_slice[2]] = ID
|
||||
pbar.update(self.IDsize)
|
||||
|
||||
# Write materials list to a text file
|
||||
if self.grid_view.materials is not None:
|
||||
with open(self.filename_materials, "w") as materials_file:
|
||||
for material in self.grid_view.materials:
|
||||
self.output_material(material, materials_file)
|
||||
@@ -0,0 +1,171 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# GNU General Public License for more details.
|
||||
#
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
import logging
|
||||
|
||||
import numpy as np
|
||||
from mpi4py import MPI
|
||||
|
||||
from gprMax._version import __version__
|
||||
from gprMax.cython.geometry_outputs import get_line_properties
|
||||
from gprMax.grid.mpi_grid import MPIGrid
|
||||
from gprMax.output_controllers.geometry_views import GeometryView, GridViewType, Metadata
|
||||
from gprMax.output_controllers.grid_view import MPIGridView
|
||||
from gprMax.subgrids.grid import SubGridBaseGrid
|
||||
from gprMax.vtkhdf_filehandlers.vtk_unstructured_grid import VtkUnstructuredGrid
|
||||
from gprMax.vtkhdf_filehandlers.vtkhdf import VtkCellType
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
class GeometryViewLines(GeometryView[GridViewType]):
|
||||
"""Unstructured grid for a per-cell-edge geometry view."""
|
||||
|
||||
def prep_vtk(self):
|
||||
"""Prepares data for writing to VTKHDF file."""
|
||||
|
||||
ID = self.grid_view.get_ID()
|
||||
|
||||
x = np.arange(self.grid_view.nx + 1, dtype=np.float64)
|
||||
y = np.arange(self.grid_view.ny + 1, dtype=np.float64)
|
||||
z = np.arange(self.grid_view.nz + 1, dtype=np.float64)
|
||||
coords = np.meshgrid(x, y, z, indexing="ij")
|
||||
self.points = np.vstack(list(map(np.ravel, coords))).T
|
||||
self.points += self.grid_view.start
|
||||
self.points *= self.grid_view.step * self.grid.dl
|
||||
|
||||
# Add offset to subgrid geometry to correctly locate within main grid
|
||||
if isinstance(self.grid, SubGridBaseGrid):
|
||||
offset = [self.grid.i0, self.grid.j0, self.grid.k0]
|
||||
self.points += offset * self.grid.dl * self.grid.ratio
|
||||
|
||||
# Each point is the 'source' for 3 lines.
|
||||
# NB: Excluding points at the far edge of the geometry as those
|
||||
# are the 'source' for no lines
|
||||
n_lines = 3 * np.prod(self.grid_view.size)
|
||||
|
||||
self.cell_types = np.full(n_lines, VtkCellType.LINE)
|
||||
self.cell_offsets = np.arange(0, 2 * n_lines + 2, 2, dtype=np.intc)
|
||||
|
||||
self.connectivity, self.material_data = get_line_properties(
|
||||
n_lines, *self.grid_view.size, ID
|
||||
)
|
||||
|
||||
assert isinstance(self.connectivity, np.ndarray)
|
||||
assert isinstance(self.material_data, np.ndarray)
|
||||
|
||||
# Write information about any PMLs, sources, receivers
|
||||
self.metadata = Metadata(self.grid_view, averaged_materials=True, materials_only=True)
|
||||
|
||||
# Number of bytes of data to be written to file
|
||||
self.nbytes = (
|
||||
self.points.nbytes
|
||||
+ self.cell_types.nbytes
|
||||
+ self.connectivity.nbytes
|
||||
+ self.cell_offsets.nbytes
|
||||
+ self.material_data.nbytes
|
||||
)
|
||||
|
||||
def write_vtk(self):
|
||||
"""Writes geometry information to a VTKHDF file."""
|
||||
|
||||
# Write the VTK file
|
||||
with VtkUnstructuredGrid(
|
||||
self.filename,
|
||||
self.points,
|
||||
self.cell_types,
|
||||
self.connectivity,
|
||||
self.cell_offsets,
|
||||
) as f:
|
||||
f.add_cell_data("Material", self.material_data)
|
||||
self.metadata.write_to_vtkhdf(f)
|
||||
|
||||
|
||||
class MPIGeometryViewLines(GeometryViewLines[MPIGridView]):
|
||||
"""Image data for a per-cell geometry view."""
|
||||
|
||||
def __init__(self, grid_view: MPIGridView, filename: str, comm: MPI.Comm):
|
||||
super().__init__(grid_view, filename)
|
||||
|
||||
self.comm = comm
|
||||
|
||||
@property
|
||||
def grid(self) -> MPIGrid:
|
||||
return self.grid_view.grid
|
||||
|
||||
def prep_vtk(self):
|
||||
"""Prepares data for writing to VTKHDF file."""
|
||||
|
||||
ID = self.grid_view.get_ID()
|
||||
|
||||
x = np.arange(self.grid_view.gx + 1, dtype=np.float64)
|
||||
y = np.arange(self.grid_view.gy + 1, dtype=np.float64)
|
||||
z = np.arange(self.grid_view.gz + 1, dtype=np.float64)
|
||||
coords = np.meshgrid(x, y, z, indexing="ij")
|
||||
self.points = np.vstack(list(map(np.ravel, coords))).T
|
||||
self.points += self.grid_view.global_start
|
||||
self.points *= self.grid_view.step * self.grid.dl
|
||||
|
||||
# Add offset to subgrid geometry to correctly locate within main grid
|
||||
if isinstance(self.grid, SubGridBaseGrid):
|
||||
offset = [self.grid.i0, self.grid.j0, self.grid.k0]
|
||||
self.points += offset * self.grid.dl * self.grid.ratio
|
||||
|
||||
# Each point is the 'source' for 3 lines.
|
||||
# NB: Excluding points at the far edge of the geometry as those
|
||||
# are the 'source' for no lines
|
||||
n_lines = 3 * np.prod(self.grid_view.global_size)
|
||||
|
||||
self.cell_types = np.full(n_lines, VtkCellType.LINE)
|
||||
self.cell_offsets = np.arange(0, 2 * n_lines + 2, 2, dtype=np.intc)
|
||||
|
||||
self.connectivity, self.material_data = get_line_properties(
|
||||
n_lines, *self.grid_view.size, ID
|
||||
)
|
||||
|
||||
assert isinstance(self.connectivity, np.ndarray)
|
||||
assert isinstance(self.material_data, np.ndarray)
|
||||
|
||||
# Write information about any PMLs, sources, receivers
|
||||
self.metadata = Metadata(self.grid_view, averaged_materials=True, materials_only=True)
|
||||
|
||||
# Number of bytes of data to be written to file
|
||||
self.nbytes = (
|
||||
self.points.nbytes
|
||||
+ self.cell_types.nbytes
|
||||
+ self.connectivity.nbytes
|
||||
+ self.cell_offsets.nbytes
|
||||
+ self.material_data.nbytes
|
||||
)
|
||||
|
||||
# Use global material IDs rather than local IDs
|
||||
self.material_data = self.grid.local_to_global_material_id_map(self.material_data)
|
||||
|
||||
def write_vtk(self):
|
||||
"""Writes geometry information to a VTKHDF file."""
|
||||
|
||||
with VtkUnstructuredGrid(
|
||||
self.filename,
|
||||
self.points,
|
||||
self.cell_types,
|
||||
self.connectivity,
|
||||
self.cell_offsets,
|
||||
comm=self.comm,
|
||||
) as f:
|
||||
f.add_cell_data("Material", self.material_data, self.grid_view.offset)
|
||||
self.metadata.write_to_vtkhdf(f)
|
||||
@@ -0,0 +1,95 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# GNU General Public License for more details.
|
||||
#
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
import logging
|
||||
|
||||
import numpy as np
|
||||
from mpi4py import MPI
|
||||
|
||||
from gprMax._version import __version__
|
||||
from gprMax.grid.mpi_grid import MPIGrid
|
||||
from gprMax.output_controllers.geometry_views import GeometryView, GridViewType, Metadata
|
||||
from gprMax.output_controllers.grid_view import MPIGridView
|
||||
from gprMax.subgrids.grid import SubGridBaseGrid
|
||||
from gprMax.vtkhdf_filehandlers.vtk_image_data import VtkImageData
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
class GeometryViewVoxels(GeometryView[GridViewType]):
|
||||
"""Image data for a per-cell geometry view."""
|
||||
|
||||
def prep_vtk(self):
|
||||
"""Prepares data for writing to VTKHDF file."""
|
||||
|
||||
self.material_data = self.grid_view.get_solid()
|
||||
|
||||
if isinstance(self.grid, SubGridBaseGrid):
|
||||
self.origin = np.array(
|
||||
[
|
||||
(self.grid.i0 * self.grid.dx * self.grid.ratio),
|
||||
(self.grid.j0 * self.grid.dy * self.grid.ratio),
|
||||
(self.grid.k0 * self.grid.dz * self.grid.ratio),
|
||||
]
|
||||
)
|
||||
else:
|
||||
self.origin = self.grid_view.start * self.grid.dl
|
||||
|
||||
self.spacing = self.grid_view.step * self.grid.dl
|
||||
|
||||
# Write information about any PMLs, sources, receivers
|
||||
self.metadata = Metadata(self.grid_view)
|
||||
|
||||
self.nbytes = self.material_data.nbytes
|
||||
|
||||
def write_vtk(self):
|
||||
"""Writes geometry information to a VTKHDF file."""
|
||||
|
||||
with VtkImageData(self.filename, self.grid_view.size, self.origin, self.spacing) as f:
|
||||
f.add_cell_data("Material", self.material_data)
|
||||
self.metadata.write_to_vtkhdf(f)
|
||||
|
||||
|
||||
class MPIGeometryViewVoxels(GeometryViewVoxels[MPIGridView]):
|
||||
"""Image data for a per-cell geometry view."""
|
||||
|
||||
def __init__(self, grid_view: MPIGridView, filename: str, comm: MPI.Comm):
|
||||
super().__init__(grid_view, filename)
|
||||
|
||||
self.comm = comm
|
||||
|
||||
@property
|
||||
def grid(self) -> MPIGrid:
|
||||
return self.grid_view.grid
|
||||
|
||||
def prep_vtk(self):
|
||||
"""Prepares data for writing to VTKHDF file."""
|
||||
|
||||
super().prep_vtk()
|
||||
|
||||
# Use global material IDs rather than local IDs
|
||||
self.material_data = self.grid.local_to_global_material_id_map(self.material_data)
|
||||
|
||||
def write_vtk(self):
|
||||
"""Writes geometry information to a VTKHDF file."""
|
||||
|
||||
with VtkImageData(
|
||||
self.filename, self.grid_view.global_size, self.origin, self.spacing, comm=self.comm
|
||||
) as f:
|
||||
f.add_cell_data("Material", self.material_data, self.grid_view.offset)
|
||||
self.metadata.write_to_vtkhdf(f)
|
||||
@@ -0,0 +1,239 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
#
|
||||
# gprMax is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# gprMax is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# GNU General Public License for more details.
|
||||
#
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
import logging
|
||||
import sys
|
||||
from abc import abstractmethod
|
||||
from pathlib import Path
|
||||
from typing import Dict, Generic, List, Sequence, Tuple, Union
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
from tqdm import tqdm
|
||||
|
||||
import gprMax.config as config
|
||||
from gprMax._version import __version__
|
||||
from gprMax.grid.fdtd_grid import FDTDGrid
|
||||
from gprMax.output_controllers.grid_view import GridType, GridView, GridViewType
|
||||
from gprMax.receivers import Rx
|
||||
from gprMax.sources import Source
|
||||
from gprMax.utilities.utilities import get_terminal_width
|
||||
from gprMax.vtkhdf_filehandlers.vtkhdf import VtkHdfFile
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
def save_geometry_views(gvs: "List[GeometryView]"):
|
||||
"""Creates and saves geometryviews.
|
||||
|
||||
Args:
|
||||
gvs: list of all GeometryViews.
|
||||
"""
|
||||
|
||||
logger.info("")
|
||||
for i, gv in enumerate(gvs):
|
||||
gv.set_filename()
|
||||
gv.prep_vtk()
|
||||
pbar = tqdm(
|
||||
total=gv.nbytes,
|
||||
unit="byte",
|
||||
unit_scale=True,
|
||||
desc=f"Writing geometry view file {i + 1}/{len(gvs)}, {gv.filename.name}",
|
||||
ncols=get_terminal_width() - 1,
|
||||
file=sys.stdout,
|
||||
disable=not config.sim_config.general["progressbars"],
|
||||
)
|
||||
gv.write_vtk()
|
||||
pbar.update(gv.nbytes)
|
||||
pbar.close()
|
||||
|
||||
logger.info("")
|
||||
|
||||
|
||||
class GeometryView(Generic[GridViewType]):
|
||||
"""Base class for Geometry Views."""
|
||||
|
||||
FILE_EXTENSION = ".vtkhdf"
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
grid_view: GridViewType,
|
||||
filename: str,
|
||||
):
|
||||
"""
|
||||
Args:
|
||||
xs, xf, ys, yf, zs, zf: ints for extent of geometry view in cells.
|
||||
dx, dy, dz: ints for spatial discretisation of geometry view in cells.
|
||||
filename: string for filename.
|
||||
grid: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
self.grid_view = grid_view
|
||||
|
||||
self.filenamebase = filename
|
||||
self.nbytes = None
|
||||
|
||||
self.material_data = None
|
||||
self.materials = None
|
||||
|
||||
@property
|
||||
def grid(self) -> FDTDGrid:
|
||||
return self.grid_view.grid
|
||||
|
||||
def set_filename(self):
|
||||
"""Construct filename from user-supplied name and model number."""
|
||||
parts = config.get_model_config().output_file_path.parts
|
||||
self.filename = Path(
|
||||
*parts[:-1], self.filenamebase + config.get_model_config().appendmodelnumber
|
||||
).with_suffix(self.FILE_EXTENSION)
|
||||
|
||||
@abstractmethod
|
||||
def prep_vtk(self):
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def write_vtk(self):
|
||||
pass
|
||||
|
||||
|
||||
class Metadata(Generic[GridType]):
|
||||
"""Comments can be strings included in the header of XML VTK file, and are
|
||||
used to hold extra (gprMax) information about the VTK data.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
gv: GridView[GridType],
|
||||
averaged_materials: bool = False,
|
||||
materials_only: bool = False,
|
||||
):
|
||||
self.gv = gv
|
||||
self.averaged_materials = averaged_materials
|
||||
self.materials_only = materials_only
|
||||
|
||||
self.gprmax_version = __version__
|
||||
self.dx_dy_dz = self.grid.dl
|
||||
self.nx_ny_nz = np.array([self.grid.nx, self.grid.ny, self.grid.nz], dtype=np.intc)
|
||||
|
||||
self.materials = self.materials_comment()
|
||||
|
||||
# Write information on PMLs, sources, and receivers
|
||||
if not self.materials_only:
|
||||
# Information on PML thickness
|
||||
if self.grid.pmls["slabs"]:
|
||||
self.pml_thickness = self.pml_gv_comment()
|
||||
else:
|
||||
self.pml_thickness = None
|
||||
srcs = (
|
||||
self.grid.hertziandipoles
|
||||
+ self.grid.magneticdipoles
|
||||
+ self.grid.voltagesources
|
||||
+ self.grid.transmissionlines
|
||||
)
|
||||
if srcs:
|
||||
self.source_ids, self.source_positions = self.srcs_rx_gv_comment(srcs)
|
||||
else:
|
||||
self.source_ids = None
|
||||
self.source_positions = None
|
||||
if self.grid.rxs:
|
||||
self.receiver_ids, self.receiver_positions = self.srcs_rx_gv_comment(self.grid.rxs)
|
||||
else:
|
||||
self.receiver_ids = None
|
||||
self.receiver_positions = None
|
||||
|
||||
@property
|
||||
def grid(self) -> GridType:
|
||||
return self.gv.grid
|
||||
|
||||
def write_to_vtkhdf(self, file_handler: VtkHdfFile):
|
||||
file_handler.add_field_data(
|
||||
"gprMax_version",
|
||||
self.gprmax_version,
|
||||
dtype=h5py.string_dtype(),
|
||||
)
|
||||
file_handler.add_field_data("dx_dy_dz", self.dx_dy_dz)
|
||||
file_handler.add_field_data("nx_ny_nz", self.nx_ny_nz)
|
||||
|
||||
file_handler.add_field_data(
|
||||
"material_ids",
|
||||
self.materials,
|
||||
dtype=h5py.string_dtype(),
|
||||
)
|
||||
|
||||
if not self.materials_only:
|
||||
if self.pml_thickness is not None:
|
||||
file_handler.add_field_data("pml_thickness", self.pml_thickness)
|
||||
|
||||
if self.source_ids is not None and self.source_positions is not None:
|
||||
file_handler.add_field_data(
|
||||
"source_ids", self.source_ids, dtype=h5py.string_dtype()
|
||||
)
|
||||
file_handler.add_field_data("sources", self.source_positions)
|
||||
|
||||
if self.receiver_ids is not None and self.receiver_positions is not None:
|
||||
file_handler.add_field_data(
|
||||
"receiver_ids", self.receiver_ids, dtype=h5py.string_dtype()
|
||||
)
|
||||
file_handler.add_field_data("receivers", self.receiver_positions)
|
||||
|
||||
def pml_gv_comment(self) -> List[int]:
|
||||
grid = self.grid
|
||||
|
||||
# Only render PMLs if they are in the geometry view
|
||||
thickness: Dict[str, int] = grid.pmls["thickness"]
|
||||
gv_pml_depth = dict.fromkeys(thickness, 0)
|
||||
|
||||
if self.gv.xs < thickness["x0"]:
|
||||
gv_pml_depth["x0"] = thickness["x0"] - self.gv.xs
|
||||
if self.gv.ys < thickness["y0"]:
|
||||
gv_pml_depth["y0"] = thickness["y0"] - self.gv.ys
|
||||
if thickness["z0"] - self.gv.zs > 0:
|
||||
gv_pml_depth["z0"] = thickness["z0"] - self.gv.zs
|
||||
if self.gv.xf > grid.nx - thickness["xmax"]:
|
||||
gv_pml_depth["xmax"] = self.gv.xf - (grid.nx - thickness["xmax"])
|
||||
if self.gv.yf > grid.ny - thickness["ymax"]:
|
||||
gv_pml_depth["ymax"] = self.gv.yf - (grid.ny - thickness["ymax"])
|
||||
if self.gv.zf > grid.nz - thickness["zmax"]:
|
||||
gv_pml_depth["zmax"] = self.gv.zf - (grid.nz - thickness["zmax"])
|
||||
|
||||
return list(gv_pml_depth.values())
|
||||
|
||||
def srcs_rx_gv_comment(
|
||||
self, srcs: Union[Sequence[Source], List[Rx]]
|
||||
) -> Tuple[List[str], npt.NDArray[np.float32]]:
|
||||
"""Used to name sources and/or receivers."""
|
||||
names: List[str] = []
|
||||
positions: npt.NDArray[np.float32] = np.empty((len(srcs), 3))
|
||||
for index, src in enumerate(srcs):
|
||||
position = src.coord * self.grid.dl
|
||||
names.append(src.ID)
|
||||
positions[index] = position
|
||||
|
||||
return names, positions
|
||||
|
||||
def dx_dy_dz_comment(self) -> npt.NDArray[np.float64]:
|
||||
return self.grid.dl
|
||||
|
||||
def nx_ny_nz_comment(self) -> npt.NDArray[np.intc]:
|
||||
return np.array([self.grid.nx, self.grid.ny, self.grid.nz], dtype=np.intc)
|
||||
|
||||
def materials_comment(self) -> List[str]:
|
||||
if not self.averaged_materials:
|
||||
return [m.ID for m in self.grid.materials if m.type != "dielectric-smoothed"]
|
||||
else:
|
||||
return [m.ID for m in self.grid.materials]
|
||||
@@ -0,0 +1,253 @@
|
||||
from typing import Generic
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
from mpi4py import MPI
|
||||
from typing_extensions import TypeVar
|
||||
|
||||
from gprMax.grid.fdtd_grid import FDTDGrid
|
||||
from gprMax.grid.mpi_grid import MPIGrid
|
||||
from gprMax.materials import Material
|
||||
|
||||
GridType = TypeVar("GridType", bound=FDTDGrid, default="FDTDGrid")
|
||||
GridViewType = TypeVar("GridViewType", bound="GridView[GridType]", default="GridView")
|
||||
|
||||
|
||||
class GridView(Generic[GridType]):
|
||||
def __init__(
|
||||
self,
|
||||
grid: GridType,
|
||||
xs: int,
|
||||
ys: int,
|
||||
zs: int,
|
||||
xf: int,
|
||||
yf: int,
|
||||
zf: int,
|
||||
dx: int = 1,
|
||||
dy: int = 1,
|
||||
dz: int = 1,
|
||||
):
|
||||
"""
|
||||
Args:
|
||||
xs, xf, ys, yf, zs, zf: ints for extent of the volume in cells.
|
||||
filename: string for filename.
|
||||
"""
|
||||
|
||||
self.start = np.array([xs, ys, zs], dtype=np.int32)
|
||||
self.stop = np.array([xf, yf, zf], dtype=np.int32)
|
||||
self.step = np.array([dx, dy, dz], dtype=np.int32)
|
||||
self.size = np.ceil((self.stop - self.start) / self.step).astype(np.int32)
|
||||
|
||||
self.grid = grid
|
||||
|
||||
self._ID = None
|
||||
|
||||
# Properties for backwards compatibility
|
||||
@property
|
||||
def xs(self) -> int:
|
||||
return self.start[0]
|
||||
|
||||
@property
|
||||
def ys(self) -> int:
|
||||
return self.start[1]
|
||||
|
||||
@property
|
||||
def zs(self) -> int:
|
||||
return self.start[2]
|
||||
|
||||
@property
|
||||
def xf(self) -> int:
|
||||
return self.stop[0]
|
||||
|
||||
@property
|
||||
def yf(self) -> int:
|
||||
return self.stop[1]
|
||||
|
||||
@property
|
||||
def zf(self) -> int:
|
||||
return self.stop[2]
|
||||
|
||||
@property
|
||||
def dx(self) -> int:
|
||||
return self.step[0]
|
||||
|
||||
@property
|
||||
def dy(self) -> int:
|
||||
return self.step[1]
|
||||
|
||||
@property
|
||||
def dz(self) -> int:
|
||||
return self.step[2]
|
||||
|
||||
@property
|
||||
def nx(self) -> int:
|
||||
return self.size[0]
|
||||
|
||||
@property
|
||||
def ny(self) -> int:
|
||||
return self.size[1]
|
||||
|
||||
@property
|
||||
def nz(self) -> int:
|
||||
return self.size[2]
|
||||
|
||||
def get_slice(self, index: int, upper_bound_exclusive: bool = True):
|
||||
if upper_bound_exclusive:
|
||||
stop = self.stop[index]
|
||||
else:
|
||||
stop = self.stop[index] + self.step[index]
|
||||
|
||||
return slice(self.start[index], stop, self.step[index])
|
||||
|
||||
def slice_array(self, array: npt.NDArray, upper_bound_exclusive: bool = True):
|
||||
return np.ascontiguousarray(
|
||||
array[
|
||||
...,
|
||||
self.get_slice(0, upper_bound_exclusive),
|
||||
self.get_slice(1, upper_bound_exclusive),
|
||||
self.get_slice(2, upper_bound_exclusive),
|
||||
]
|
||||
)
|
||||
|
||||
def initialise_materials(self):
|
||||
ID = self.get_ID(force_refresh=True)
|
||||
|
||||
materials_in_grid_view = np.unique(ID)
|
||||
|
||||
self.materials = np.array(self.grid.materials, dtype=Material)[materials_in_grid_view]
|
||||
self.materials.sort()
|
||||
|
||||
# Create map from material ID to 0 - number of materials
|
||||
materials_map = {material.numID: index for index, material in enumerate(self.materials)}
|
||||
self.map_materials_func = np.vectorize(lambda id: materials_map[id])
|
||||
|
||||
def map_to_view_materials(self, array: npt.NDArray):
|
||||
return self.map_materials_func(array)
|
||||
|
||||
def get_ID(self, force_refresh=False):
|
||||
if self._ID is None or force_refresh:
|
||||
self._ID = self.slice_array(self.grid.ID, upper_bound_exclusive=False)
|
||||
return self._ID
|
||||
|
||||
def get_solid(self):
|
||||
return self.slice_array(self.grid.solid)
|
||||
|
||||
def get_rigidE(self):
|
||||
return self.slice_array(self.grid.rigidE)
|
||||
|
||||
def get_rigidH(self):
|
||||
return self.slice_array(self.grid.rigidH)
|
||||
|
||||
def get_Ex(self):
|
||||
return self.slice_array(self.grid.Ex, upper_bound_exclusive=False)
|
||||
|
||||
def get_Ey(self):
|
||||
return self.slice_array(self.grid.Ey, upper_bound_exclusive=False)
|
||||
|
||||
def get_Ez(self):
|
||||
return self.slice_array(self.grid.Ez, upper_bound_exclusive=False)
|
||||
|
||||
def get_Hx(self):
|
||||
return self.slice_array(self.grid.Hx, upper_bound_exclusive=False)
|
||||
|
||||
def get_Hy(self):
|
||||
return self.slice_array(self.grid.Hy, upper_bound_exclusive=False)
|
||||
|
||||
def get_Hz(self):
|
||||
return self.slice_array(self.grid.Hz, upper_bound_exclusive=False)
|
||||
|
||||
|
||||
class MPIGridView(GridView[MPIGrid]):
|
||||
def __init__(
|
||||
self,
|
||||
grid: MPIGrid,
|
||||
xs: int,
|
||||
ys: int,
|
||||
zs: int,
|
||||
xf: int,
|
||||
yf: int,
|
||||
zf: int,
|
||||
dx: int = 1,
|
||||
dy: int = 1,
|
||||
dz: int = 1,
|
||||
):
|
||||
super().__init__(grid, xs, ys, zs, xf, yf, zf, dx, dy, dz)
|
||||
|
||||
self.global_size = self.size
|
||||
|
||||
# Calculate start for the local grid
|
||||
self.global_start = self.start
|
||||
self.start = self.grid.global_to_local_coordinate(self.start)
|
||||
|
||||
# Bring start into the local grid (and not in the negative halo)
|
||||
# local_start must still be aligned with the provided step.
|
||||
self.start = np.where(
|
||||
self.start < self.grid.negative_halo_offset,
|
||||
self.grid.negative_halo_offset
|
||||
+ ((self.start - self.grid.negative_halo_offset) % self.step),
|
||||
self.start,
|
||||
)
|
||||
|
||||
# Calculate stop for the local grid
|
||||
self.stop = self.grid.global_to_local_coordinate(self.stop)
|
||||
|
||||
self.has_positive_neighbour = self.stop > self.grid.size
|
||||
|
||||
# Limit stop such that it is at most one step beyond the max
|
||||
# index of the grid. As stop is the upper bound, it is
|
||||
# exclusive, meaning when used to slice an array (with the
|
||||
# provided step), the last element accessed will one step below
|
||||
# stop.
|
||||
# Note: using self.grid.size as an index in any dimension would
|
||||
# fall in the positive halo (this counts as outside the local
|
||||
# grid).
|
||||
self.stop = np.where(
|
||||
self.has_positive_neighbour,
|
||||
self.grid.size + ((self.stop - self.grid.size) % self.step),
|
||||
self.stop,
|
||||
)
|
||||
|
||||
# Calculate offset for the local grid view
|
||||
self.offset = self.grid.local_to_global_coordinate(self.start) - self.global_start
|
||||
|
||||
# Update local size
|
||||
self.size = self.stop - self.start
|
||||
|
||||
@property
|
||||
def gx(self) -> int:
|
||||
return self.global_size[0]
|
||||
|
||||
@property
|
||||
def gy(self) -> int:
|
||||
return self.global_size[1]
|
||||
|
||||
@property
|
||||
def gz(self) -> int:
|
||||
return self.global_size[2]
|
||||
|
||||
def get_slice(self, index: int, upper_bound_exclusive: bool = True):
|
||||
if upper_bound_exclusive:
|
||||
stop = self.stop[index]
|
||||
else:
|
||||
# Make slice of array one step larger if this rank does not
|
||||
# have a positive neighbour
|
||||
stop = np.where(
|
||||
self.has_positive_neighbour,
|
||||
self.stop,
|
||||
self.stop + self.step,
|
||||
)
|
||||
|
||||
return slice(self.start[index], stop, self.step[index])
|
||||
|
||||
def initialise_materials(self, comm: MPI.Comm):
|
||||
ID = self.get_ID(force_refresh=True)
|
||||
|
||||
local_material_ids = np.unique(ID)
|
||||
local_materials = np.array(self.grid.materials, dtype=Material)[local_material_ids]
|
||||
|
||||
self.materials, material_id_map = self.grid.remap_material_ids(
|
||||
local_materials.tolist(), comm
|
||||
)
|
||||
|
||||
# Create map from material ID to 0 - number of materials
|
||||
self.map_materials_func = np.vectorize(lambda id: material_id_map(id))
|
||||
在新工单中引用
屏蔽一个用户